Spatial transcriptomics reveal heterogeneous cell‒cell interactions among brain regions in cuprizone model consistent with multiple sclerosis lesions.

1. #' Scissor: Single-Cell Identification of Subpopulations with bulk Sample phenOtype coRrelation
2. #'
3. #' \code{Scissor} is a novel approach that utilizes the phenotypes, such as disease stage, tumor metastasis, treatment response, and survival
4. #' outcomes, collected from bulk assays to identify the most highly phenotype-associated cell subpopulations from single-cell data.
5. #'
6. #' Scissor is a novel algorithm to identify cell subpopulations from single-cell data that are most highly associated with the given phenotypes.
7. #' The three data sources for Scissor inputs are a single-cell expression matrix, a bulk expression matrix, and a phenotype of interest.
8. #' The phenotype annotation of each bulk sample can be a continuous dependent variable, binary group indicator vector, or clinical survival data.
9. #' The key step of Scissor is to quantify the similarity between the single-cell and bulk samples by Pearson correlation for each pair of cells and bulk samples.
10. #' After this, Scissor optimizes a regression model on the correlation matrix with the sample phenotype.
11. #' The selection of the regression model depends on the type of the input phenotype, i.e., linear regression for continuous variables,
12. #' logistic regression for dichotomous variables, and Cox regression for clinical survival data.
13. #' Based on the signs of the estimated regression coefficients, the cells with non-zero coefficients can be indicated as
14. #' Scissor positive (Scissor+) cells and Scissor negative (Scissor-) cells, which are positively and negatively associated
15. #' with the phenotype of interest, respectively.
16. #'
17. #' @param bulk_dataset Bulk expression matrix of related disease. Each row represents a gene and each column represents a sample.
18. #' @param sc_dataset Single-cell RNA-seq expression matrix of related disease. Each row represents a gene and each column represents a sample.
19. #' A Seurat object that contains the preprocessed data and constructed network is preferred. Otherwise, a cell-cell similarity network is
20. #' constructed based on the input matrix.
21. #' @param phenotype Phenotype annotation of each bulk sample. It can be a continuous dependent variable,
22. #' binary group indicator vector, or clinical survival data:
23. #' \itemize{
24. #' \item Continuous dependent variable. Should be a quantitative vector for \code{family = gaussian}.
25. #' \item Binary group indicator vector. Should be either a 0-1 encoded vector or a factor with two levels for \code{family = binomial}.
26. #' \item Clinical survival data. Should be a two-column matrix with columns named 'time' and 'status'. The latter is a binary variable,
27. #' with '1' indicating event (e.g.recurrence of cancer or death), and '0' indicating right censored.
28. #' The function \code{Surv()} in package survival produces such a matrix.
29. #' }
30. #' @param tag Names for each phenotypic group. Used for linear and logistic regressions only.
31. #' @param alpha Parameter used to balance the effect of the l1 norm and the network-based penalties. It can be a number or a searching vector.
32. #' If \code{alpha = NULL}, a default searching vector is used. The range of alpha is in \code{[0,1]}. A larger alpha lays more emphasis on the l1 norm.
33. #' @param cutoff Cutoff for the percentage of the Scissor selected cells in total cells. This parameter is used to restrict the number of the
34. #' Scissor selected cells. A cutoff less than \code{50\%} (default \code{20\%}) is recommended depending on the input data.
35. #' @param family Response type for the regression model. It depends on the type of the given phenotype and
36. #' can be \code{family = gaussian} for linear regression, \code{family = binomial} for classification, or \code{family = cox} for Cox regression.
37. #' @param Save_file File name for saving the preprocessed regression inputs into a RData.
38. #' @param Load_file File name for loading the preprocessed regression inputs. It can help to tune the model parameter \code{alpha}.
39. #' Please see Scissor Tutorial for more details.
40. #'
41. #' @return This function returns a list with the following components:
42. #' \item{para}{A list contains the final model parameters.}
43. #' \item{Coefs}{The regression coefficient for each cell.}
44. #' \item{Scissor_pos}{The cell IDs of Scissor+ cells.}
45. #' \item{Scissor_neg}{The cell IDs of Scissor- cells.}
46. #'
47. #' @references Duanchen Sun and Zheng Xia (2021): Phenotype-guided subpopulation identification from single-cell sequencing data. Nature Biotechnology.
48. #' @import Seurat Matrix preprocessCore
49. #' @export
50. Scissor <- function(bulk_dataset, sc_dataset, phenotype, tag = NULL,
51. alpha = NULL, cutoff = 0.2, family = c("gaussian","binomial","cox"),
52. Save_file = "Scissor_inputs.RData", Load_file = NULL){
53. library(Seurat)
54. library(Matrix)
55. library(preprocessCore)
58. if (is.null(Load_file)){
59. common <- intersect(rownames(bulk_dataset), rownames(sc_dataset))
60. if (length(common) == 0) {
61. stop("There is no common genes between the given single-cell and bulk samples.")
62. }
63. if (class(sc_dataset) == "Seurat"){
64. sc_exprs <- as.matrix(sc_dataset@assays$RNA@data)
65. network <- as.matrix(sc_dataset@graphs$RNA_snn)
66. }else{
67. sc_exprs <- as.matrix(sc_dataset)
68. Seurat_tmp <- CreateSeuratObject(sc_dataset)
69. Seurat_tmp <- FindVariableFeatures(Seurat_tmp, selection.method = "vst", verbose = F)
70. Seurat_tmp <- ScaleData(Seurat_tmp, verbose = F)
71. Seurat_tmp <- RunPCA(Seurat_tmp, features = VariableFeatures(Seurat_tmp), verbose = F)
72. Seurat_tmp <- FindNeighbors(Seurat_tmp, dims = 1:10, verbose = F)
73. network <- as.matrix(Seurat_tmp@graphs$RNA_snn)
74. }
75. diag(network) <- 0
76. network[which(network != 0)] <- 1
78. dataset0 <- cbind(bulk_dataset[common,], sc_exprs[common,]) # Dataset before quantile normalization.
79. dataset1 <- normalize.quantiles(dataset0) # Dataset after quantile normalization.
80. rownames(dataset1) <- rownames(dataset0)
81. colnames(dataset1) <- colnames(dataset0)
83. Expression_bulk <- dataset1[,1:ncol(bulk_dataset)]
84. Expression_cell <- dataset1[,(ncol(bulk_dataset) + 1):ncol(dataset1)]
85. X <- cor(Expression_bulk, Expression_cell)
87. quality_check <- quantile(X)
88. print("|**************************************************|")
89. print("Performing quality-check for the correlations")
90. print("The five-number summary of correlations:")
91. print(quality_check)
92. print("|**************************************************|")
93. if (quality_check[3] < 0.01){
94. warning("The median correlation between the single-cell and bulk samples is relatively low.")
95. }
96. if (family == "binomial"){
97. Y <- as.numeric(phenotype)
98. z <- table(Y)
99. if (length(z) != length(tag)){
100. stop("The length differs between tags and phenotypes. Please check Scissor inputs and selected regression type.")
101. }else{
102. print(sprintf("Current phenotype contains %d %s and %d %s samples.", z[1], tag[1], z[2], tag[2]))
103. print("Perform logistic regression on the given phenotypes:")
104. }
105. }
106. if (family == "gaussian"){
107. Y <- as.numeric(phenotype)
108. z <- table(Y)
109. if (length(z) != length(tag)){
110. stop("The length differs between tags and phenotypes. Please check Scissor inputs and selected regression type.")
111. }else{
112. tmp <- paste(z, tag)
113. print(paste0("Current phenotype contains ", paste(tmp[1:(length(z)-1)], collapse = ", "), ", and ", tmp[length(z)], " samples."))
114. print("Perform linear regression on the given phenotypes:")
115. }
116. }
117. if (family == "cox"){
118. Y <- as.matrix(phenotype)
119. if (ncol(Y) != 2){
120. stop("The size of survival data is wrong. Please check Scissor inputs and selected regression type.")
121. }else{
122. print("Perform cox regression on the given clinical outcomes:")
123. }
124. }
125. save(X, Y, network, Expression_bulk, Expression_cell, file = Save_file)
126. }else{
127. load(Load_file)
128. }
130. if (is.null(alpha)){
131. alpha <- c(0.005, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
132. }
133. for (i in 1:length(alpha)){
134. set.seed(123)
135. fit0 <- APML1(X, Y, family = family, penalty = "Net", alpha = alpha[i], Omega = network, nlambda = 100, nfolds = min(10,nrow(X)))
136. fit1 <- APML1(X, Y, family = family, penalty = "Net", alpha = alpha[i], Omega = network, lambda = fit0$lambda.min)
137. if (family == "binomial"){
138. Coefs <- as.numeric(fit1$Beta[2:(ncol(X)+1)])
139. }else{
140. Coefs <- as.numeric(fit1$Beta)
141. }
142. Cell1 <- colnames(X)[which(Coefs > 0)]
143. Cell2 <- colnames(X)[which(Coefs < 0)]
144. percentage <- (length(Cell1) + length(Cell2)) / ncol(X)
145. print(sprintf("alpha = %s", alpha[i]))
146. print(sprintf("Scissor identified %d Scissor+ cells and %d Scissor- cells.", length(Cell1), length(Cell2)))
147. print(sprintf("The percentage of selected cell is: %s%%", formatC(percentage*100, format = 'f', digits = 3)))
149. if (percentage < cutoff){
150. break
151. }
152. cat("\n")
153. }
154. print("|**************************************************|")
156. return(list(para = list(alpha = alpha[i], lambda = fit0$lambda.min, family = family),
157. Coefs = Coefs,
158. Scissor_pos = Cell1,
159. Scissor_neg = Cell2))
160. }

Scissor.R at commit 311560a, under GPL-3.0 · at the source