Browse code

version 5

miccec authored on 02/03/2020 18:19:17
Showing19 changed files

... ...
@@ -1,6 +1,6 @@
1 1
 Package: scTHI
2 2
 Title:  Indentification of  significantly activated ligand-receptor interactions across clusters of cells from single-cell RNA sequencing data. 
3
-Version:  0.99.4
3
+Version:  0.99.5:
4 4
 Authors@R: c(person("Francesca Pia", "Caruso", email = "francescapia.caruso@gmail.com", role = c("aut")),
5 5
 	person("Michele", "Ceccarelli", email = "m.ceccarelli@gmail.com", role = c("aut", "cre")))
6 6
 Description: scTHI is an R package to identify active pairs of ligand-receptors from single cells in order to study,among others, tumor-host interactions. scTHI contains a set of signatures to classify cells from the tumor microenvironment.
... ...
@@ -1,3 +1,3 @@
1
-Changes in version 0.99.4 (2020-03-02)
1
+Changes in version 0.99.5 (2020-03-02)
2 2
 + Submitted to Bioconductor
3 3
 
4 4
new file mode 100644
... ...
@@ -0,0 +1,152 @@
1
+#' TME_classification
2
+#'
3
+#' The function allows the user to classify non-tumor cells in tumor
4
+#' microenvironment. It implements the Mann-Whitney-Wilcoxon Gene
5
+#' Set Test
6
+#' (MWW-GST) algorithm and tests for each cell the enrichment of
7
+#' a collection
8
+#' of signatures of different cell types.
9
+#' @param expMat Gene expression matrix where rows are genes
10
+#'   presented with
11
+#'   Hugo Symbols and columns are cells. Gene expression values
12
+#'   should be normalized counts.
13
+#' @param minLenGeneSet Minimum gene set length
14
+#' @param pvalFilter Logical, if TRUE results will be filtered
15
+#' for p-Value.
16
+#'   Defoult is FALSE.
17
+#' @param alternative a character string specifying the alternative
18
+#' hypothesis of wilcoxon test, must be one of "two.sided" (default),
19
+#' "greater" or "less".
20
+#' @param fdrFilter Logical, if TRUE results will be filtered for FDR.
21
+#' @param pvalCutoff Numeric p-Value (or FDR) threshold. Gene set with
22
+#'   p-Value (or FDR) greater than pvalCutoff will be discarded
23
+#'   (default is 0.01).
24
+#' @param nesCutoff Numeric threshold. Gene set with NES greater than
25
+#'   nesCutoff will be discarded (default is 0.58)
26
+#' @param nNES Default is 0.58, so each cell is classified with
27
+#'  a specific
28
+#'   phenotype based on the first significant enriched gene set.
29
+#' @examples
30
+#' data(scExample)
31
+#' Class <- TME_classification(scExample)
32
+#' @return A list with two items: Class (character) and ClassLegend
33
+#'   (character)
34
+#' @export
35
+#'
36
+#' TME_classification
37
+
38
+TME_classification <- function(expMat,
39
+                               minLenGeneSet = 10,
40
+                               alternative = "two.sided",
41
+                               pvalFilter = FALSE,
42
+                               fdrFilter = TRUE,
43
+                               pvalCutoff = 0.01,
44
+                               nesCutoff = 0.58,
45
+                               nNES = 1) {
46
+  
47
+  zerogenes <- rowSums(expMat == 0)
48
+  expMat <- expMat[zerogenes != ncol(expMat), ]
49
+  means <- rowMeans(expMat)
50
+  sds <- apply(expMat, 1, sd)
51
+  E <- apply(expMat, 2, function(x) {
52
+    (x - means) / sds
53
+  })
54
+  
55
+  E_ <- apply(E, 2, rank)
56
+  tmp <- lapply(signatures, function(x) {
57
+    sum(x %in% rownames(E))
58
+  })
59
+  signatures_ <- signatures[tmp > minLenGeneSet]
60
+  
61
+  ans <- vector("list", length(signatures_))
62
+  names(ans) <- names(signatures_)
63
+  ans <- lapply(signatures_, function(S) {
64
+    geneSet <- S
65
+    gs <- geneSet[which(geneSet %in% rownames(E))]
66
+    outside_gs <- setdiff(rownames(E), gs)
67
+    nx <- length(gs)
68
+    ny <- length(outside_gs)
69
+    Tstat <- colSums(E_[outside_gs, ])
70
+    Ustat <- nx * ny + ny * (ny + 1) / 2 - Tstat
71
+    mu <- nx * ny / 2
72
+    sigma <- sqrt(mu * (nx + ny + 1) / 6)
73
+    zValue <- Ustat - mu
74
+    correction <- vapply(zValue, function(x) {
75
+      switch(
76
+        alternative,
77
+        two.sided = sign(x) * 0.5,
78
+        greater = 0.5,
79
+        less = -0.5
80
+      )
81
+    },numeric(1))
82
+    zValue <- (zValue - correction) / sigma
83
+    pValue <- vapply(zValue, function(x) {
84
+      switch(
85
+        alternative,
86
+        less = 1 - pnorm(x),
87
+        greater = pnorm(-x),
88
+        two.sided = 2 * pnorm(-abs(x))
89
+      )
90
+    },numeric(1))
91
+    nes <- Ustat / nx / ny
92
+    pu <- nes / (1 - nes)
93
+    log.pu <- log2(pu)
94
+    names(Ustat) <- colnames(E_)
95
+    names(pValue) <- colnames(E_)
96
+    names(nes) <- colnames(E_)
97
+    names(pu) <- colnames(E_)
98
+    names(log.pu) <- colnames(E_)
99
+    #ans[[i]] <- 
100
+    list(
101
+      statistic = Ustat,
102
+      p.value = pValue,
103
+      nes = nes,
104
+      pu = pu,
105
+      log.pu = log.pu
106
+    )})
107
+  
108
+  
109
+  NES <- vapply(ans, function(x) {
110
+    cbind(x$log.pu)
111
+  },numeric(ncol(expMat)))
112
+  NES <- t(NES)
113
+  colnames(NES) <- colnames(expMat)
114
+  pValue <- vapply(ans, function(x) {
115
+    cbind(x$p.value)
116
+  },numeric(ncol(expMat)))
117
+  pValue <- t(pValue)
118
+  colnames(pValue) <- colnames(expMat)
119
+  
120
+  FDR <- apply(pValue, 2, function(x) {
121
+    p.adjust(x, method = "fdr")
122
+  })
123
+  if (pvalFilter == TRUE) {
124
+    NES[pValue >= pvalCutoff] <- 0
125
+  }
126
+  if (fdrFilter == TRUE) {
127
+    NES[FDR >= pvalCutoff] <- 0
128
+  }
129
+  NES[NES < nesCutoff] <- 0
130
+  if (nNES == 1) {
131
+    Class <- apply(NES, 2, function(x) {
132
+      names(which.max(x))
133
+    })
134
+  }
135
+  if (nNES != 1) {
136
+    Class <- apply(NES, 2, function(x) {
137
+      names(sort(x, decreasing = TRUE))[nNES]
138
+    })
139
+  }
140
+  Class[colSums(NES != 0) < nNES] <- "nc"
141
+  phenotype <- signaturesColors[Class, ]
142
+  rownames(phenotype) <- names(Class)
143
+  Class <- phenotype$ALLPhenotypeFinal
144
+  names(Class) <- rownames(phenotype)
145
+  ClassLegend <- phenotype$Color
146
+  names(ClassLegend) <- phenotype$ALLPhenotypeFinal
147
+  ClassLegend <- ClassLegend[!duplicated(ClassLegend)]
148
+  #print(sort(table(Class), decreasing = TRUE))
149
+  Classification <- list(Class, ClassLegend)
150
+  names(Classification) <- c("Class", "ClassLegend")
151
+  return(Classification)
152
+}
0 153
\ No newline at end of file
1 154
new file mode 100644
... ...
@@ -0,0 +1,54 @@
1
+#' TME_plot
2
+#'
3
+#' Generates a plot on the t-SNE coordinates, labeling cells by TME
4
+#' classification.
5
+#'
6
+#' @param tsneData X and y coordinates of points in the plot.
7
+#' @param Class Object returned by TME_classification function.
8
+#' @param cexPoint Set the point size.
9
+#' @examples
10
+#' data(scExample)
11
+#' result <-  scTHI_score(scExample,
12
+#'            cellCusterA = colnames(scExample)[1:30],
13
+#'            cellCusterB = colnames(scExample)[31:100],
14
+#'            cellCusterAName = "ClusterA",
15
+#'            cellCusterBName = "ClusterB", filterCutoff = 0,
16
+#'            pvalueCutoff = 1, nPermu = 100, ncore = 8)
17
+#' result <- scTHI_runTsne(result)
18
+#' Class <- TME_classification(scExample)
19
+#' TME_plot(tsneData = result$tsneData, Class)
20
+#' @return  None
21
+#' @export
22
+#'
23
+#' TME_plot
24
+
25
+TME_plot <- function(tsneData, Class, cexPoint = 0.8) {
26
+  
27
+  ClassColor <- Class$ClassLegend[Class$Class]
28
+  names(ClassColor) <- names(Class$Class)
29
+  
30
+  TmeColors <- rep("gray80", nrow(tsneData))
31
+  names(TmeColors) <- rownames(tsneData)
32
+  TmeColors[names(ClassColor)] <- ClassColor
33
+  par(oma = c(0, 0, 0, 6))
34
+  plot(
35
+    tsneData[, 1:2],
36
+    pch = 16,
37
+    cex = cexPoint,
38
+    col = TmeColors[rownames(tsneData)],
39
+    main = "",
40
+    xlab = "tSNE 1",
41
+    ylab = "tSNE 2"
42
+  )
43
+  legend(
44
+    par("usr")[2],
45
+    par("usr")[4],
46
+    bty = "n",
47
+    xpd = NA,
48
+    legend = names(Class$ClassLegend),
49
+    col = Class$ClassLegend,
50
+    cex = 0.9,
51
+    pch = 15,
52
+    box.lty = 0
53
+  )
54
+}
0 55
new file mode 100644
... ...
@@ -0,0 +1,5 @@
1
+#' @importFrom grDevices colorRampPalette hcl
2
+#' @importFrom graphics barplot legend par plot text
3
+#' @importFrom stats density p.adjust quantile sd pnorm
4
+#' @importFrom BiocParallel SnowParam bplapply
5
+NULL
0 6
\ No newline at end of file
1 7
new file mode 100644
... ...
@@ -0,0 +1,66 @@
1
+#' scTHI_plotCluster
2
+#'
3
+#' Graphs the output of scTHI_runTsne, labeling cells by clusters.
4
+#' @param scTHIresult scTHI object.
5
+#' @param cexPoint Set the point size.
6
+#' @param legendPos Character string to custom the legend position.
7
+#' @examples
8
+#' data(scExample)
9
+#' result <-  scTHI_score(scExample,
10
+#'                        cellCusterA = colnames(scExample)[1:30],
11
+#'                        cellCusterB = colnames(scExample)[31:100],
12
+#'                        cellCusterAName = "ClusterA",
13
+#'                        cellCusterBName = "ClusterB", filterCutoff = 0,
14
+#'                        pvalueCutoff = 1, nPermu = 100, ncore = 8)
15
+#' result <- scTHI_runTsne(result)
16
+#' scTHI_plotCluster(result)
17
+#' @return None
18
+#' @export
19
+#'
20
+#' scTHI_plotCluster
21
+
22
+scTHI_plotCluster <- function(scTHIresult,
23
+                              cexPoint = 0.8,
24
+                              legendPos = c(
25
+                                "topleft", "topright",
26
+                                "bottomright", "bottomleft"
27
+                              )) {
28
+  
29
+  tsneData <- scTHIresult$tsneData
30
+  legendPos <- match.arg(legendPos)
31
+  
32
+  nCluster <- length(scTHIresult) - 3
33
+  ggplot_colors <- function(n) {
34
+    hues <- seq(15, 375, length = n + 1)
35
+    hcl(h = hues, l = 65, c = 100)[1:n]
36
+  }
37
+  
38
+  Colors <- ggplot_colors(n = nCluster)
39
+  names(Colors) <- names(scTHIresult)[4:length(scTHIresult)]
40
+  
41
+  ClusterColors <- rep("gray80", nrow(tsneData))
42
+  names(ClusterColors) <- rownames(tsneData)
43
+  
44
+  ClusterColors[scTHIresult[[names(Colors)[1]]]] <- Colors[1]
45
+  ClusterColors[scTHIresult[[names(Colors)[2]]]] <- Colors[2]
46
+  
47
+  plot(
48
+    tsneData[, 1:2],
49
+    pch = 16,
50
+    cex = cexPoint,
51
+    col = ClusterColors[rownames(tsneData)],
52
+    main = "",
53
+    xlab = "tSNE 1",
54
+    ylab = "tSNE 2"
55
+  )
56
+  legend(
57
+    legendPos,
58
+    inset = .02,
59
+    legend = names(Colors),
60
+    col = Colors,
61
+    cex = 1,
62
+    pch = 15,
63
+    box.lty = 0,
64
+    bg = "transparent"
65
+  )
66
+}
0 67
\ No newline at end of file
1 68
new file mode 100644
... ...
@@ -0,0 +1,199 @@
1
+colByValue_ <-
2
+  function(x,
3
+           col,
4
+           range = NA,
5
+           breaks = NA,
6
+           cex.axis = 2,
7
+           las = 1,
8
+           ...) {
9
+    if (is.vector(x)) {
10
+      x <- as.matrix(x)
11
+    }
12
+    if (is.na(range[1])) {
13
+      
14
+    }
15
+    else {
16
+      x[x < range[1]] <- range[1]
17
+      x[x > range[2]] <- range[2]
18
+    }
19
+    if (is.na(breaks[1])) {
20
+      ff <- seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE),
21
+                length = length(col)
22
+      )
23
+      bg2 <- apply(
24
+        as.matrix(as.numeric(unlist(x))), 1,
25
+        function(x) {
26
+          rank(c(ff, x), ties.method = "min")[length(col) + 1]
27
+        }
28
+      )
29
+      dens <- matrix(bg2, nrow(x), ncol(x))
30
+      result <- matrix(col[dens], nrow = nrow(x), ncol = ncol(x))
31
+      row.names(result) <- row.names(x)
32
+      # image(x = 1:2, y = as.matrix(ff), z = t(ff), col = col,
33
+      #       xaxt = "n", ylab = "", las = las, xlab = "",
34
+      #       xlim = c(1, 4), bty = "n", ...)
35
+      return(result)
36
+    }
37
+    else {
38
+      temp <- cut(as.numeric(unlist(x)),
39
+                  breaks = breaks,
40
+                  include.lowest = TRUE
41
+      )
42
+      if (length(col) != length(levels(temp))) {
43
+        stop("length:col != length: cut result")
44
+      }
45
+      result <- matrix(col[as.numeric(temp)],
46
+                       nrow = nrow(x),
47
+                       ncol = ncol(x)
48
+      )
49
+      row.names(result) <- row.names(x)
50
+      return(result)
51
+    }
52
+  }
53
+
54
+
55
+
56
+getcolors <- function(genesToplot, expMat, 
57
+                      tsneData){
58
+  
59
+  ans <- lapply(genesToplot, function(x){
60
+    tmp_exp <-
61
+      t(expMat[x, colnames(expMat), drop = FALSE])
62
+    mean_exp_samples <- mean(tmp_exp)
63
+    
64
+    ##### range gene x gene
65
+    range_exp <- range(tmp_exp)
66
+    bre <- seq(
67
+      min(tmp_exp[, x]) - 0.01,
68
+      max(tmp_exp[, x]) + 0.01, 0.01
69
+    )
70
+    colori <-
71
+      colByValue_(
72
+        tmp_exp,
73
+        col = colorRampPalette(c(
74
+          "gray80", "coral",
75
+          "red"
76
+        ))(length(bre) - 1),
77
+        breaks = bre,
78
+        cex.axis = 0.8
79
+      )
80
+    colori <- colori[rownames(tsneData), , drop = FALSE]
81
+    tmp <- list(colori, range_exp, mean_exp_samples)
82
+    names(tmp) <- c("colori", "range_exp", "mean_exp_samples")
83
+    tmp})
84
+  
85
+  return(ans)
86
+}
87
+
88
+
89
+#' scTHI_plotPairs
90
+#'
91
+#' Generates a plot on the t-SNE coordinates to
92
+#' show the expression levels
93
+#' of an interaction pair of interest. Each cell is
94
+#' colored according to the
95
+#' corresponding
96
+#' gene expression value.
97
+#' @param scTHIresult scTHI object.
98
+#' @param cexPoint Set the point size.
99
+#' @param interactionToplot Interaction pair to plot.
100
+#' @examples
101
+#' data(scExample)
102
+#' result <-  scTHI_score(scExample,
103
+#'                  cellCusterA = colnames(scExample)[1:30],
104
+#'                  cellCusterB = colnames(scExample)[31:100],
105
+#'                  cellCusterAName = "ClusterA",
106
+#'                  cellCusterBName = "ClusterB", filterCutoff = 0,
107
+#'                  pvalueCutoff = 1, nPermu = 100, ncore = 8)
108
+#' result <- scTHI_runTsne(result)
109
+#' scTHI_plotPairs(result,interactionToplot = "CXCL12_CD4")
110
+
111
+#' @return None
112
+#' @export
113
+#'
114
+#' scTHI_plotPairs
115
+scTHI_plotPairs <-
116
+  function(scTHIresult,
117
+           cexPoint = 0.8,
118
+           interactionToplot) {
119
+   
120
+    tsneData <- scTHIresult$tsneData
121
+    result <- scTHIresult$result
122
+    expMat <- scTHIresult$expMat
123
+    
124
+    genesToplotA <- result[interactionToplot, c("partnerA")]
125
+    genesToplotA <- unlist(strsplit(genesToplotA, ":"))
126
+    genesToplotB <- result[interactionToplot, c("partnerB")]
127
+    genesToplotB <- unlist(strsplit(genesToplotB, ":"))
128
+    
129
+    #################  create list colori
130
+    list_colori_A <- getcolors(genesToplotA, expMat, tsneData)
131
+    names(list_colori_A) <- genesToplotA
132
+    
133
+    list_colori_B <- getcolors(genesToplotB, expMat, tsneData)
134
+    names(list_colori_B) <- genesToplotB
135
+    
136
+    ################# plot
137
+    par(
138
+      mfrow = c(1, length(c(
139
+        genesToplotA, genesToplotB
140
+      ))),
141
+      mar = c(6, 4, 3, 2),
142
+      cex.main = 1.5,
143
+      cex.sub = 1.2,
144
+      col.main = "black",
145
+      col.sub = "gray30",
146
+      font.main = 2,
147
+      font.sub = 3
148
+    )
149
+    
150
+    lapply(names(list_colori_A), function(x) {
151
+      colori <- list_colori_A[[x]]$colori
152
+      range_exp <- round(list_colori_A[[x]]$range_exp, digits = 2)
153
+      mean_exp_samples <-
154
+        round(list_colori_A[[x]]$mean_exp_samples, digits = 2)
155
+      plot(
156
+        tsneData[, 1:2],
157
+        pch = 16,
158
+        cex = cexPoint,
159
+        col = colori,
160
+        main = x,
161
+        xlab = "",
162
+        ylab = "",
163
+        sub = paste0(
164
+          "Range: ",
165
+          range_exp[1],
166
+          " to ",
167
+          range_exp[2],
168
+          ", MeanExp: ",
169
+          mean_exp_samples
170
+        )
171
+      )
172
+    })
173
+    
174
+    lapply(names(list_colori_B), function(x) {
175
+      colori <- list_colori_B[[x]]$colori
176
+      range_exp <- round(list_colori_B[[x]]$range_exp, digits = 2)
177
+      mean_exp_samples <-
178
+        round(list_colori_B[[x]]$mean_exp_samples, digits = 2)
179
+      plot(
180
+        tsneData[, 1:2],
181
+        pch = 16,
182
+        cex = cexPoint,
183
+        col = colori,
184
+        main = x,
185
+        xlab = "",
186
+        ylab = "",
187
+        sub = paste0(
188
+          "Range: ",
189
+          range_exp[1],
190
+          " to ",
191
+          range_exp[2],
192
+          ", MeanExp: ",
193
+          mean_exp_samples
194
+        )
195
+      )
196
+    })
197
+  
198
+    par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)
199
+  }
0 200
\ No newline at end of file
1 201
new file mode 100644
... ...
@@ -0,0 +1,104 @@
1
+#' scTHI_plotResult
2
+#'
3
+#' Creates barplots of scTHI_score results.
4
+#' @param scTHIresult scTHI object.
5
+#' @param cexNames Size of names in barplot.
6
+#' @param plotType Type of plot to be generated. Default is
7
+#' "score", can be  also "pair". The "score" option will generate
8
+#' a barplot for each resulted interaction pair, representing
9
+#' the calculated interaction score
10
+#'   and the related p-Value.The "pair" option will generate
11
+#'   two barplot for each resulted interaction pair, representing
12
+#'   the percentage of cells of each cluster expressing partnerA
13
+#'   and partnerB gene, respectively.
14
+#' @param nRes Number of pairs to plot (all if NULL).
15
+#' @examples
16
+#' data(scExample)
17
+#' result <-  scTHI_score(scExample,
18
+#'                        cellCusterA = colnames(scExample)[1:30],
19
+#'                        cellCusterB = colnames(scExample)[31:100],
20
+#'                        cellCusterAName = "ClusterA",
21
+#'                        cellCusterBName = "ClusterB", filterCutoff = 0,
22
+#'                        pvalueCutoff = 1, nPermu = 100, ncore = 8)
23
+#'
24
+#' scTHI_plotResult(result, plotType = "score")
25
+#' scTHI_plotResult(result, plotType = "pair")
26
+#' @return None
27
+#' @export
28
+#'
29
+#' scTHI_plotResult
30
+
31
+scTHI_plotResult <- function(scTHIresult,
32
+                             cexNames = 0.8,
33
+                             plotType = c("score", "pair"),
34
+                             nRes = NULL) {
35
+  
36
+  result <- scTHIresult$result
37
+  if (!is.null(nRes)) {
38
+    result <- result[1:nRes, ]
39
+  }
40
+  
41
+  
42
+  if (plotType == "score") {
43
+    tmp <- as.matrix(result[, c("SCORE", "pValue")])
44
+    pvalues <- apply(tmp, 1, function(x) {
45
+      if (x["pValue"] > 0.05) "ns" else if (x["pValue"] <= 0.0001){ 
46
+        "****"
47
+      }else if (x["pValue"] <= 0.001) {
48
+        "***"
49
+      }else if (x["pValue"] <= 0.01) {
50
+        "**"
51
+      }else if (x["pValue"] <= 0.05) {
52
+        "*"
53
+      }
54
+      
55
+    })
56
+    
57
+    par(mar = c(5, 10, 4, 2) + 0.1)
58
+    barplot(
59
+      t(tmp[, "SCORE"]),
60
+      beside = TRUE,
61
+      xlim = c(0, 1.1),
62
+      xlab = "scTHI Score",
63
+      col = c("lightseagreen"),
64
+      cex.names = cexNames,
65
+      horiz = TRUE,
66
+      las = 2
67
+    )
68
+    text(
69
+      x = 1.05,
70
+      y = seq(1.5, by = 2, length.out = nrow(tmp)),
71
+      labels = pvalues,
72
+      cex = 1.2
73
+    )
74
+    par(mar = c(5, 4, 4, 2) + 0.1)
75
+  }
76
+  
77
+  if (plotType == "pair") {
78
+    tmp <- as.matrix(result[, grep("rnk", colnames(result))])
79
+    par(mar = c(5, 10, 4, 2) + 0.1, oma = c(0, 0, 0, 6))
80
+    barplot(
81
+      t(tmp),
82
+      beside = TRUE,
83
+      xlim = c(0, 1.1),
84
+      xlab = "% of Cells",
85
+      col = c("#F8766D", "#00BFC4"),
86
+      cex.names = cexNames,
87
+      horiz = TRUE,
88
+      las = 2
89
+    )
90
+    legend(
91
+      par("usr")[2],
92
+      par("usr")[4],
93
+      bty = "n",
94
+      xpd = NA,
95
+      legend = c("ClusterA", "ClusterB"),
96
+      col = c("#F8766D", "#00BFC4"),
97
+      cex = 1,
98
+      pch = 15,
99
+      box.lty = 0,
100
+      bg = "transparent"
101
+    )
102
+    par(mar = c(5, 4, 4, 2) + 0.1, oma = c(0, 0, 0, 0))
103
+  }
104
+}
0 105
\ No newline at end of file
1 106
new file mode 100644
... ...
@@ -0,0 +1,55 @@
1
+#' scTHI_runTsne
2
+#'
3
+#' Runs t-SNE dimensionality reduction on selected features based on Rtsne
4
+#' package.
5
+#' @param scTHIresult scTHI object.
6
+#' @examples
7
+#' data(scExample)
8
+#' result <-  scTHI_score(scExample,
9
+#'                        cellCusterA = colnames(scExample)[1:30],
10
+#'                        cellCusterB = colnames(scExample)[31:100],
11
+#'                        cellCusterAName = "ClusterA",
12
+#'                        cellCusterBName = "ClusterB", filterCutoff = 0,
13
+#'                        pvalueCutoff = 1, nPermu = 100, ncore = 8)
14
+#' result <- scTHI_runTsne(result)
15
+#' @return The same object as scTHI_score with a fifth item tsneData
16
+#'   (data.frame)
17
+#' @export
18
+#'
19
+#' scTHI_runTsne
20
+
21
+scTHI_runTsne <- function(scTHIresult) {
22
+
23
+  expMat <- scTHIresult$expMat
24
+  
25
+  ## create tsne
26
+  eta <- .1
27
+  filter <- apply(expMat, 1, function(x) {
28
+    sum(quantile(x,
29
+                 probs = c(1 - eta, eta)
30
+    ) * c(1, -1))
31
+  })
32
+  # fivenum(filter)
33
+  # plot(density(log2(filter)))
34
+  cutoff <- density(log2(filter))
35
+  cutoff <- data.frame(x = cutoff$x, y = cutoff$y)
36
+  cutoff <- cutoff$x[which.max(cutoff$y)]
37
+  foldChange <- 2^cutoff
38
+  # sum(filter > foldChange) # 4567
39
+  variableGenes <- names(filter)[filter > foldChange]
40
+  expMat <- expMat[variableGenes, ]
41
+  
42
+  requireNamespace("Rtsne")
43
+  expMatT <- t(expMat)
44
+  tsne_out <- Rtsne::Rtsne(expMatT)
45
+  tsneData <- data.frame(
46
+    x = tsne_out$Y[, 1],
47
+    y = tsne_out$Y[, 2],
48
+    Sample = colnames(expMat),
49
+    stringsAsFactors = FALSE
50
+  )
51
+  rownames(tsneData) <- colnames(expMat)
52
+  scTHIresult <- c(list(tsneData), scTHIresult)
53
+  names(scTHIresult)[1] <- "tsneData"
54
+  return(scTHIresult)
55
+}
0 56
\ No newline at end of file
1 57
new file mode 100644
... ...
@@ -0,0 +1,576 @@
1
+checkInteractions <- function(expMat){
2
+  tmp_check <-
3
+    interaction_table[, c(
4
+      "partnerA1",
5
+      "partnerA2",
6
+      "partnerA3",
7
+      "partnerB1",
8
+      "partnerB2",
9
+      "partnerB3"
10
+    )]
11
+  tmp_check <- unique(tmp_check[!is.na(tmp_check)])
12
+  if (sum(tmp_check %in% rownames(expMat)) != length(tmp_check)) {
13
+    tmp_genes <- setdiff(tmp_check, rownames(expMat))
14
+    tmp <- which(
15
+      interaction_table$partnerA1 %in% tmp_genes |
16
+        interaction_table$partnerA2 %in% tmp_genes |
17
+        interaction_table$partnerA3 %in% tmp_genes |
18
+        interaction_table$partnerB1 %in% tmp_genes |
19
+        interaction_table$partnerB2 %in% tmp_genes |
20
+        interaction_table$partnerB3 %in% tmp_genes
21
+    )
22
+    interaction_table <- interaction_table[-tmp, ]
23
+    message("Warning: Not all interaction genes are in expMat")
24
+  }
25
+  if (nrow(interaction_table) == 0) {
26
+    stop("ERROR: No interaction genes to test")
27
+  }
28
+  return(interaction_table)
29
+}
30
+
31
+rankingValues <- function(x, topRank = 10){
32
+  n <- round(nrow(x) * topRank / 100)
33
+  ddata_ <- apply(x, 2, function(x) {
34
+    rank(-x)
35
+  })
36
+  ddata_ <- apply(ddata_, 2, function(x) {
37
+    x <= n
38
+  })
39
+  return(ddata_)
40
+}
41
+
42
+getScore <- function(expMat, interaction_table, cellCuster, 
43
+                     genes1columns, genes2columns,autocrineEffect){
44
+  ddata <- expMat[, cellCuster]
45
+  ddata_ <- rankingValues(ddata)
46
+  
47
+  rnk <- rep(0, nrow(interaction_table))
48
+  expValue <- rep(0, nrow(interaction_table))
49
+  names(rnk) <- rownames(interaction_table)
50
+  names(expValue) <- rownames(interaction_table)
51
+  
52
+  ans <-  apply(interaction_table, 1, function(x){
53
+    ggenes1 <- x[genes1columns]
54
+    ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
55
+    tmp_genes1 <- ddata_[ggenes1, , drop = FALSE]
56
+    tmp1 <- colSums(tmp_genes1 == TRUE)
57
+    
58
+    if (autocrineEffect == FALSE) {
59
+      ggenes2 <-
60
+        x[genes2columns]
61
+      ggenes2 <- as.vector(ggenes2[!is.na(ggenes2)])
62
+      tmp_genes2 <- ddata_[ggenes2, , drop = FALSE]
63
+      tmp2 <- colSums(tmp_genes2 == FALSE)
64
+      tmp <- tmp1 + tmp2
65
+      
66
+      rnk<- round(sum(tmp == (nrow(tmp_genes1) +
67
+                                nrow(tmp_genes2))) / ncol(tmp_genes1), digits = 2)
68
+      
69
+      expValue <-
70
+        paste(round(rowMeans(ddata[ggenes1, cellCuster,
71
+                                   drop = FALSE
72
+                                   ]), digits = 2), collapse = " # ")
73
+      
74
+    }
75
+    if (autocrineEffect == TRUE)  {
76
+      tmp <- tmp1
77
+      rnk <-
78
+        round(sum(tmp == nrow(tmp_genes1)) / ncol(tmp_genes1),
79
+              digits = 2
80
+        )
81
+      expValue <-
82
+        paste(round(rowMeans(ddata[ggenes1, cellCuster,
83
+                                   drop = FALSE
84
+                                   ]), digits = 2), collapse = " # ")
85
+    }
86
+    list(rnk,expValue)
87
+  })
88
+  rresult <- data.frame(
89
+    rnk = unlist(lapply(ans, function(x) x[[1]])),
90
+    expValue = unlist(lapply(ans, function(x) x[[2]]))
91
+  )
92
+  
93
+  return(rresult)
94
+}
95
+
96
+
97
+#' scTHI_score
98
+#'
99
+#' This function allows the user to compute a score for a set of
100
+#' ligand-receptor pairs, from a single cell gene expression matrix,
101
+#' and detect specific Tumor-Host interactions. You must specify at
102
+#' least two clusters of cells (for example tumor cells and immune
103
+#' cells).
104
+#' @param expMat ScRNA-seq gene expression matrix where rows are genes
105
+#'   presented with Hugo Symbols and columns are cells. Gene expression
106
+#'   values should be counts or normalized counts.
107
+#' @param cellCusterA Vector of columns of expMat that belong to the
108
+#' first cluster.
109
+#' @param cellCusterB Vector of columns of expMat that belong to the
110
+#' second cluster.
111
+#' @param cellCusterAName A character string labeling the clusterA.
112
+#' @param cellCusterBName A character string labeling the clusterB.
113
+#' @param topRank Filter threshold. Set to 10 (default) means that
114
+#' each gene of the interaction pair will be considered as expressed
115
+#' in a cell if it's in the top rank 10 percent.
116
+#' @param autocrineEffect if TRUE remove the paracrine filter
117
+#' @param fileNameBase Project name.
118
+#' @param filterCutoff Score threshold (default is 0.50). For each
119
+#'   interaction pair, if the score calculated (for the partnerA
120
+#'   or partnerB)
121
+#'   will be less than filterCutoff the interaction pair will be
122
+#'   discarded.
123
+#' @param PValue Logical, set to TRUE (default) compute statistical
124
+#'   iterations. If p.value < 0.05, the value will be returned.
125
+#' @param pvalueCutoff cutoff of the p-value
126
+#' @param nPermu Number of iterations to perform (default is 1000).
127
+#' @param ncore Number of processors to use.
128
+#' @keywords interaction
129
+#' @examples
130
+#'
131
+#' ####################### example of scTHI_score
132
+#' data(scExample)
133
+#' result <-  scTHI_score(scExample,
134
+#'       cellCusterA = colnames(scExample)[1:30],
135
+#'       cellCusterB = colnames(scExample)[31:100],
136
+#'       cellCusterAName = "ClusterA",
137
+#'       cellCusterBName = "ClusterB", filterCutoff = 0,
138
+#'      pvalueCutoff = 1, nPermu = 100, ncore = 8)
139
+#'
140
+#' @return A list of results, with four items: result (data.frame),
141
+#'   expMat (matrix), clusterA (character),  clusterA (character)
142
+#' @export
143
+#' scTHI_score
144
+
145
+scTHI_score <-
146
+  function(expMat,
147
+           cellCusterA,
148
+           cellCusterB,
149
+           cellCusterAName,
150
+           cellCusterBName,
151
+           topRank = 10,
152
+           autocrineEffect = TRUE,
153
+           fileNameBase = "scTHI",
154
+           filterCutoff = 0.5,
155
+           PValue = TRUE,
156
+           pvalueCutoff = 0.05,
157
+           nPermu = 1000,
158
+           ncore = 8) {
159
+    
160
+   
161
+    ############ check all interaction genes in expMat ############
162
+    interaction_table <- checkInteractions(expMat)
163
+    
164
+    #################### scTHI score ##############################
165
+    message(paste(
166
+      "Computing score for",
167
+      nrow(interaction_table),
168
+      "interaction pairs..."
169
+    ))
170
+    
171
+    ##### ddataA
172
+    resultA <- getScore(expMat, interaction_table = interaction_table,
173
+                        cellCuster = cellCusterA,
174
+                        genes1columns= c("partnerA1", "partnerA2", 
175
+                                         "partnerA3"),
176
+                        genes2columns=c("partnerB1", "partnerB2", 
177
+                                        "partnerB3"),
178
+                        autocrineEffect)
179
+    colnames(resultA) <-
180
+      paste0(c("rnkPartnerA", "expValueA"), "_", cellCusterAName)
181
+    print(paste("Computed ranked values for partner A"))
182
+    
183
+    ###### ddataB
184
+    resultB <- getScore(expMat, interaction_table = interaction_table,
185
+                        cellCuster = cellCusterB,
186
+                        genes1columns= c("partnerB1", "partnerB2", 
187
+                                         "partnerB3"),
188
+                        genes2columns=c("partnerA1", "partnerA2", 
189
+                                        "partnerA3"),
190
+                        autocrineEffect)
191
+    colnames(resultB) <-
192
+      paste0(c("rnkPartnerB", "expValueB"), "_", cellCusterBName)
193
+    print(paste("Computed ranked values for partner B"))
194
+    
195
+    
196
+    #################### merge results #################
197
+    SCORE <- rowMeans(cbind(resultA$rnkPartnerA, resultB$rnkPartnerB))
198
+    result <- data.frame(interaction_table, resultA, resultB, SCORE)
199
+    
200
+    ############### filter low-score int ###############
201
+    message("Removing low score interactions...")
202
+    interestColumn <- c(
203
+      paste0("rnkPartnerA_", cellCusterAName),
204
+      paste0("rnkPartnerB_", cellCusterBName)
205
+    )
206
+    tmp <- rowSums(result[, interestColumn] < filterCutoff)
207
+    result[tmp != 0, "SCORE"] <- NA
208
+    result <- result[!is.na(result$SCORE), ]
209
+    result <- result[order(result$SCORE, decreasing = TRUE), ]
210
+    result <- data.frame(
211
+      interationPair = rownames(result),
212
+      result,
213
+      stringsAsFactors = FALSE
214
+    )
215
+    columnToremove <- c(
216
+      "partnerA1",
217
+      "partnerA2",
218
+      "partnerA3",
219
+      "partnerB1",
220
+      "partnerB2",
221
+      "partnerB3"
222
+    )
223
+    result <- result[, setdiff(colnames(result), columnToremove)]
224
+    
225
+    if (nrow(result) == 0) {
226
+      stop("No interaction pair exceed the score filterCutoff")
227
+    }
228
+    
229
+    scTHIresult <- list(result, expMat, cellCusterA, cellCusterB)
230
+    names(scTHIresult) <-
231
+      c("result", "expMat", cellCusterAName, cellCusterBName)
232
+    
233
+    ###################### Permutation #########
234
+    if (PValue == TRUE) {
235
+      message("Computing permutation....")
236
+      interaction_table <- interaction_table[rownames(result), ]
237
+      columnToadd <- c(
238
+        "partnerA1",
239
+        "partnerA2",
240
+        "partnerA3",
241
+        "partnerB1",
242
+        "partnerB2",
243
+        "partnerB3"
244
+      )
245
+      result <- data.frame(interaction_table[, columnToadd], result,
246
+                           stringsAsFactors = FALSE
247
+      )
248
+      
249
+      if (length(grep("simple", result$interactionType)) != 0) {
250
+        sample_pairInteraction <-
251
+          rownames(result)[result$interactionType == "simple"][1]
252
+        
253
+        param <- SnowParam(workers = ncore, type = "SOCK")
254
+        if (autocrineEffect == FALSE) {
255
+          paracrineFun <- function(dummy, rankingValues) {
256
+            ddataA <- expMat[, cellCusterA]
257
+            ddataA <-
258
+              apply(ddataA, 2, function(x) {
259
+                sample(x, replace = FALSE)
260
+              })
261
+            rownames(ddataA) <- rownames(expMat)
262
+            ddataA_ <- rankingValues(ddataA)
263
+            
264
+            ggenes1 <-
265
+              result[sample_pairInteraction, c(
266
+                "partnerA1", "partnerA2",
267
+                "partnerA3"
268
+              )]
269
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
270
+            tmp_genes1 <- ddataA_[ggenes1, , drop = FALSE]
271
+            tmp1 <- colSums(tmp_genes1 == TRUE)
272
+            ggenes2 <-
273
+              result[sample_pairInteraction, c(
274
+                "partnerB1", "partnerB2",
275
+                "partnerB3"
276
+              )]
277
+            ggenes2 <- as.vector(ggenes2[!is.na(ggenes2)])
278
+            tmp_genes2 <- ddataA_[ggenes2, , drop = FALSE]
279
+            tmp2 <- colSums(tmp_genes2 == FALSE)
280
+            tmp <- tmp1 + tmp2
281
+            rnk_permuA <- round(sum(tmp == (
282
+              nrow(tmp_genes1) +
283
+                nrow(tmp_genes2)
284
+            )) / ncol(tmp_genes1), digits = 2)
285
+            
286
+            
287
+            ddataB <- expMat[, cellCusterB]
288
+            ddataB <-
289
+              apply(ddataB, 2, function(x) {
290
+                sample(x, replace = FALSE)
291
+              })
292
+            rownames(ddataB) <- rownames(expMat)
293
+            ddataB_ <- rankingValues(ddataB)
294
+            
295
+            ggenes1 <-
296
+              result[sample_pairInteraction, c(
297
+                "partnerB1", "partnerB2",
298
+                "partnerB3"
299
+              )]
300
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
301
+            tmp_genes1 <- ddataB_[ggenes1, , drop = FALSE]
302
+            tmp1 <- colSums(tmp_genes1 == TRUE)
303
+            
304
+            ggenes2 <-
305
+              result[sample_pairInteraction, c(
306
+                "partnerA1", "partnerA2",
307
+                "partnerA3"
308
+              )]
309
+            ggenes2 <- as.vector(ggenes2[!is.na(ggenes2)])
310
+            tmp_genes2 <- ddataB_[ggenes2, , drop = FALSE]
311
+            tmp2 <- colSums(tmp_genes2 == FALSE)
312
+            tmp <- tmp1 + tmp2
313
+            
314
+            rnk_permuB <- round(sum(tmp == (
315
+              nrow(tmp_genes1) +
316
+                nrow(tmp_genes2)
317
+            )) / ncol(tmp_genes1), digits = 2)
318
+            rnkInteraction <- mean(rnk_permuA, rnk_permuB)
319
+            return(rnkInteraction)
320
+          }
321
+          ans <- bplapply(1:nPermu,
322
+                          FUN = paracrineFun,
323
+                          rankingValues = rankingValues,
324
+                          BPPARAM = param
325
+                          
326
+          )
327
+          ans <- unlist(ans)
328
+          permutatedPvalue_simple <- ans
329
+        }
330
+        if (autocrineEffect == TRUE) {
331
+          autocrineFun <- function(dummy, rankingValues) {
332
+            ddataA <- expMat[, cellCusterA]
333
+            ddataA <-
334
+              apply(ddataA, 2, function(x) {
335
+                sample(x, replace = FALSE)
336
+              })
337
+            rownames(ddataA) <- rownames(expMat)
338
+            ddataA_ <- rankingValues(ddataA)
339
+            
340
+            ggenes1 <-
341
+              result[sample_pairInteraction, c(
342
+                "partnerA1", "partnerA2",
343
+                "partnerA3"
344
+              )]
345
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
346
+            tmp_genes1 <- ddataA_[ggenes1, , drop = FALSE]
347
+            tmp1 <- colSums(tmp_genes1 == TRUE)
348
+            tmp <- tmp1
349
+            rnk_permuA <-
350
+              round(sum(tmp == nrow(tmp_genes1)) / ncol(tmp_genes1),
351
+                    digits = 2
352
+              )
353
+            
354
+            
355
+            ddataB <- expMat[, cellCusterB]
356
+            ddataB <-
357
+              apply(ddataB, 2, function(x) {
358
+                sample(x, replace = FALSE)
359
+              })
360
+            rownames(ddataB) <- rownames(expMat)
361
+            ddataB_ <- rankingValues(ddataB)
362
+            
363
+            ggenes1 <-
364
+              result[sample_pairInteraction, c(
365
+                "partnerB1", "partnerB2",
366
+                "partnerB3"
367
+              )]
368
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
369
+            tmp_genes1 <- ddataB_[ggenes1, , drop = FALSE]
370
+            tmp1 <- colSums(tmp_genes1 == TRUE)
371
+            tmp <- tmp1
372
+            
373
+            rnk_permuB <-
374
+              round(sum(tmp == nrow(tmp_genes1)) / ncol(tmp_genes1),
375
+                    digits = 2
376
+              )
377
+            rnkInteraction <- mean(rnk_permuA, rnk_permuB)
378
+            return(rnkInteraction)
379
+          }
380
+          ans <- bplapply(1:nPermu,
381
+                          FUN = autocrineFun,
382
+                          rankingValues = rankingValues,
383
+                          BPPARAM = param
384
+                          
385
+          )
386
+          ans <- unlist(ans)
387
+          permutatedPvalue_simple <- ans
388
+        }
389
+      }
390
+      
391
+      if (length(grep("complex", result$interactionType)) != 0) {
392
+        sample_pairInteraction <-
393
+          rownames(result)[result$interactionType == "complex"][1]
394
+        
395
+        
396
+        param <- SnowParam(workers = ncore, type = "SOCK")
397
+        if (autocrineEffect == FALSE) {
398
+          paracrineFun2 <- function(dummy, rankingValues) {
399
+            ddataA <- expMat[, cellCusterA]
400
+            ddataA <-
401
+              apply(ddataA, 2, function(x) {
402
+                sample(x, replace = FALSE)
403
+              })
404
+            rownames(ddataA) <- rownames(expMat)
405
+            ddataA_ <- rankingValues(ddataA)
406
+            
407
+            ggenes1 <-
408
+              result[sample_pairInteraction, c(
409
+                "partnerA1", "partnerA2",
410
+                "partnerA3"
411
+              )]
412
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
413
+            tmp_genes1 <- ddataA_[ggenes1, , drop = FALSE]
414
+            tmp1 <- colSums(tmp_genes1 == TRUE)
415
+            
416
+            ggenes2 <-
417
+              result[sample_pairInteraction, c(
418
+                "partnerB1", "partnerB2",
419
+                "partnerB3"
420
+              )]
421
+            ggenes2 <- as.vector(ggenes2[!is.na(ggenes2)])
422
+            tmp_genes2 <- ddataA_[ggenes2, , drop = FALSE]
423
+            tmp2 <- colSums(tmp_genes2 == FALSE)
424
+            
425
+            tmp <- tmp1 + tmp2
426
+            rnk_permuA <- round(sum(tmp == (
427
+              nrow(tmp_genes1) +
428
+                nrow(tmp_genes2)
429
+            )) / ncol(tmp_genes1),
430
+            digits = 2
431
+            )
432
+            
433
+            
434
+            ddataB <- expMat[, cellCusterB]
435
+            ddataB <-
436
+              apply(ddataB, 2, function(x) {
437
+                sample(x, replace = FALSE)
438
+              })
439
+            rownames(ddataB) <- rownames(expMat)
440
+            ddataB_ <- rankingValues(ddataB)
441
+            
442
+            ggenes1 <-
443
+              result[sample_pairInteraction, c(
444
+                "partnerB1", "partnerB2",
445
+                "partnerB3"
446
+              )]
447
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
448
+            tmp_genes1 <- ddataB_[ggenes1, , drop = FALSE]
449
+            tmp1 <- colSums(tmp_genes1 == TRUE)
450
+            
451
+            ggenes2 <-
452
+              result[sample_pairInteraction, c(
453
+                "partnerA1", "partnerA2",
454
+                "partnerA3"
455
+              )]
456
+            ggenes2 <- as.vector(ggenes2[!is.na(ggenes2)])
457
+            tmp_genes2 <- ddataB_[ggenes2, , drop = FALSE]
458
+            tmp2 <- colSums(tmp_genes2 == FALSE)
459
+            tmp <- tmp1 + tmp2
460
+            
461
+            rnk_permuB <- round(sum(tmp == (
462
+              nrow(tmp_genes1) +
463
+                nrow(tmp_genes2)
464
+            )) / ncol(tmp_genes1), digits = 2)
465
+            rnkInteraction <- mean(rnk_permuA, rnk_permuB)
466
+            return(rnkInteraction)
467
+          }
468
+          ans <- bplapply(1:nPermu,
469
+                          FUN = paracrineFun2,
470
+                          rankingValues = rankingValues,
471
+                          BPPARAM = param
472
+                          
473
+          )
474
+          ans <- unlist(ans)
475
+          permutatedPvalue_complex <- ans
476
+        }
477
+        if (autocrineEffect == TRUE) {
478
+          autocrineFun2 <- function(dummy, rankingValues) {
479
+            ddataA <- expMat[, cellCusterA]
480
+            ddataA <-
481
+              apply(ddataA, 2, function(x) {
482
+                sample(x, replace = FALSE)
483
+              })
484
+            rownames(ddataA) <- rownames(expMat)
485
+            ddataA_ <- rankingValues(ddataA)
486
+            
487
+            ggenes1 <-
488
+              result[sample_pairInteraction, c(
489
+                "partnerA1", "partnerA2",
490
+                "partnerA3"
491
+              )]
492
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
493
+            tmp_genes1 <- ddataA_[ggenes1, , drop = FALSE]
494
+            tmp1 <- colSums(tmp_genes1 == TRUE)
495
+            tmp <- tmp1
496
+            rnk_permuA <-
497
+              round(sum(tmp == nrow(tmp_genes1)) / ncol(tmp_genes1),
498
+                    digits = 2
499
+              )
500
+            
501
+            
502
+            ddataB <- expMat[, cellCusterB]
503
+            ddataB <-
504
+              apply(ddataB, 2, function(x) {
505
+                sample(x, replace = FALSE)
506
+              })
507
+            rownames(ddataB) <- rownames(expMat)
508
+            ddataB_ <- rankingValues(ddataB)
509
+            
510
+            ggenes1 <-
511
+              result[sample_pairInteraction, c(
512
+                "partnerB1", "partnerB2",
513
+                "partnerB3"
514
+              )]
515
+            ggenes1 <- as.vector(ggenes1[!is.na(ggenes1)])
516
+            tmp_genes1 <- ddataB_[ggenes1, , drop = FALSE]
517
+            tmp1 <- colSums(tmp_genes1 == TRUE)
518
+            tmp <- tmp1
519
+            
520
+            rnk_permuB <-
521
+              round(sum(tmp == nrow(tmp_genes1)) / ncol(tmp_genes1),
522
+                    digits = 2
523
+              )
524
+            rnkInteraction <- mean(rnk_permuA, rnk_permuB)
525
+            return(rnkInteraction)
526
+          }
527
+          ans <- bplapply(1:nPermu,
528
+                          FUN = autocrineFun2,
529
+                          rankingValues = rankingValues,
530
+                          BPPARAM = param
531
+                          
532
+          )
533
+          ans <- unlist(ans)
534
+          permutatedPvalue_complex <- ans
535
+        }
536
+      }
537
+      
538
+      
539
+      
540
+      r <- apply(result, 1, function (x) {
541
+        if (x["interactionType"] == "simple") {
542
+          ans <-    sum(permutatedPvalue_simple > x["SCORE"]) / nPermu
543
+        }
544
+        if (x["interactionType"] == "complex") {
545
+          ans <- sum(permutatedPvalue_complex > x["SCORE"]) / nPermu
546
+        }
547
+        ans
548
+      })
549
+      SCOREpValue <- unlist(r)
550
+      FDR <- round(p.adjust(SCOREpValue, method = "fdr"), digits = 3)
551
+      result <- data.frame(
552
+        result,
553
+        pValue = SCOREpValue,
554
+        FDR = FDR,
555
+        stringsAsFactors = FALSE
556
+      )
557
+      
558
+      result <- result[result$pValue <= pvalueCutoff, ]
559
+      columnToremove <- c(
560
+        "partnerA1",
561
+        "partnerA2",
562
+        "partnerA3",
563
+        "partnerB1",
564
+        "partnerB2",
565
+        "partnerB3"
566
+      )
567
+      result <- result[, setdiff(colnames(result), columnToremove)]
568
+      scTHIresult$result <- result
569
+      # save(scTHIresult, file = paste0(fileNameBase,
570
+      # "_", cellCusterAName,
571
+      #  "&", cellCusterBName, ".RData"))
572
+    }
573
+    message("Interaction pairs detected:")
574
+    print(rownames(scTHIresult$result))
575
+    return(scTHIresult)
576
+  }
0 577
new file mode 100644
... ...
@@ -0,0 +1,78 @@
1
+---
2
+title: "External data"
3
+output: html_document
4
+---
5
+
6
+```{r setup, include=FALSE}
7
+knitr::opts_chunk$set(echo = TRUE)
8
+```
9
+
10
+# Interaction table
11
+
12
+The object `interaction_table` contains a list of 2,548 pairs of ligand and receptors. Data has been collected curating publicly available resources from:
13
+
14
+ * Ramilowski JA, Goldberg T, Harshbarger J, Kloppmann E, Lizio M, Satagopam VP, et al. A draft
15
+network of ligand-receptor-mediated multicellular signalling in human. Nat Commun. 2015;6:7866.
16
+doi:10.1038/ncomms8866.
17
+
18
+ * Vento-Tormo R, Efremova M, Botting RA, Turco MY, Vento-Tormo M, Meyer KB, et al. Single-cell
19
+reconstruction of the early maternal-fetal interface in humans. Nature. 2018;563:347–53.
20
+doi:10.1038/s41586-018-0698-6.
21
+
22
+The curated list is composed of known and novel literature-supported interactions and includes both heteromeric and monomeric ligands/receptors mainly related to chemokine, cytokine, growth factors, integrin, TGF and TNF family members, semaphorins, ephrins, Wnt and Notch signalings.. 
23
+
24
+
25
+
26
+# Signatures
27
+
28
+The object `signatures` contains a list of 295 cell-type specific signatures, including the Immune and Central Nervous Systems. Data has been collected integrating a cureted collection of marker genes and a set of signatures available from public databases and published studies, as:
29
+
30
+ * a compendium of 64 human cell types signatures including lymphoid, myeloid, stromal, tissue-specific and stem cells, collected from FANTOM5, ENCODE, Blueprint and Gene Expression Omnibus (GEO) data portals; 
31
+ 
32
+ * a set of markers for 30 immune cell types, including myeloid and lymphoid subpopulations identified from Peripheral Blood Mononuclear Cells (PBMCs) [1]; 
33
+ 
34
+ * a set of Central Nervous System cell signatures including astrocytes, neuron, oligodendrocytes, microglia, and endothelial cells [2];
35
+ 
36
+ * a set of 53 signatures corresponding to 26 different cell types [3-6]
37
+ 
38
+ * two gene expression programs related to microglia and bone marrow-derived macrophages in gliomas [7]. 
39
+
40
+
41
+The object `signaturesColors` is a data.frame including a list of colors for `signatures`.
42
+
43
+
44
+##References
45
+
46
+[1] Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data
47
+across different conditions, technologies, and species. Nat Biotechnol. 2018;36:411–20.
48
+doi:10.1038/nbt.4096.
49
+
50
+
51
+[2] Zhang Y, Chen K, Sloan SA, Bennett ML, Scholze AR, O’Keeffe S, et al. An RNA-sequencing
52
+transcriptome and splicing database of glia, neurons, and vascular cells of the cerebral cortex. J
53
+Neurosci. 2014;34:11929–47. doi:10.1523/JNEUROSCI.1860-14.2014.
54
+
55
+
56
+[3] Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal
57
+dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity.
58
+2013;39:782–95. doi:10.1016/j.immuni.2013.10.003.
59
+
60
+
61
+[4] Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer
62
+Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of
63
+Response to Checkpoint Blockade. Cell Rep. 2017;18:248–62. doi:10.1016/j.celrep.2016.12.019.
64
+
65
+
66
+[5] Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors
67
+associated with local immune cytolytic activity. Cell. 2015;160:48–61. doi:10.1016/j.cell.2014.12.033.
68
+
69
+
70
+[6] Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, et al. Dissecting the
71
+multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352:189–96.
72
+doi:10.1126/science.aad0501.
73
+
74
+
75
+[7] Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, et al. Decoupling genetics,
76
+lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. 2017;355.
77
+doi:10.1126/science.aai8478.
78
+
0 79
new file mode 100644
... ...
@@ -0,0 +1,55 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/TME_classification.R
3
+\name{TME_classification}
4
+\alias{TME_classification}
5
+\title{TME_classification}
6
+\usage{
7
+TME_classification(expMat, minLenGeneSet = 10,
8
+  alternative = "two.sided", pvalFilter = FALSE, fdrFilter = TRUE,
9
+  pvalCutoff = 0.01, nesCutoff = 0.58, nNES = 1)
10
+}
11
+\arguments{
12
+\item{expMat}{Gene expression matrix where rows are genes
13
+presented with
14
+Hugo Symbols and columns are cells. Gene expression values
15
+should be normalized counts.}
16
+
17
+\item{minLenGeneSet}{Minimum gene set length}
18
+
19
+\item{alternative}{a character string specifying the alternative
20
+hypothesis of wilcoxon test, must be one of "two.sided" (default),
21
+"greater" or "less".}
22
+
23
+\item{pvalFilter}{Logical, if TRUE results will be filtered
24
+for p-Value.
25
+  Defoult is FALSE.}
26
+
27
+\item{fdrFilter}{Logical, if TRUE results will be filtered for FDR.}
28
+
29
+\item{pvalCutoff}{Numeric p-Value (or FDR) threshold. Gene set with
30
+p-Value (or FDR) greater than pvalCutoff will be discarded
31
+(default is 0.01).}
32
+
33
+\item{nesCutoff}{Numeric threshold. Gene set with NES greater than
34
+nesCutoff will be discarded (default is 0.58)}
35
+
36
+\item{nNES}{Default is 0.58, so each cell is classified with
37
+a specific
38
+ phenotype based on the first significant enriched gene set.}
39
+}
40
+\value{
41
+A list with two items: Class (character) and ClassLegend
42
+  (character)
43
+}
44
+\description{
45
+The function allows the user to classify non-tumor cells in tumor
46
+microenvironment. It implements the Mann-Whitney-Wilcoxon Gene
47
+Set Test
48
+(MWW-GST) algorithm and tests for each cell the enrichment of
49
+a collection
50
+of signatures of different cell types.
51
+}
52
+\examples{
53
+data(scExample)
54
+Class <- TME_classification(scExample)
55
+}
0 56
new file mode 100644
... ...
@@ -0,0 +1,34 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/TME_plot.R
3
+\name{TME_plot}
4
+\alias{TME_plot}
5
+\title{TME_plot}
6
+\usage{
7
+TME_plot(tsneData, Class, cexPoint = 0.8)
8
+}
9
+\arguments{
10
+\item{tsneData}{X and y coordinates of points in the plot.}
11
+
12
+\item{Class}{Object returned by TME_classification function.}
13
+
14
+\item{cexPoint}{Set the point size.}
15
+}
16
+\value{
17
+None
18
+}
19
+\description{
20
+Generates a plot on the t-SNE coordinates, labeling cells by TME
21
+classification.
22
+}
23
+\examples{
24
+data(scExample)
25
+result <-  scTHI_score(scExample,
26
+           cellCusterA = colnames(scExample)[1:30],
27
+           cellCusterB = colnames(scExample)[31:100],
28
+           cellCusterAName = "ClusterA",
29
+           cellCusterBName = "ClusterB", filterCutoff = 0,
30
+           pvalueCutoff = 1, nPermu = 100, ncore = 8)
31
+result <- scTHI_runTsne(result)
32
+Class <- TME_classification(scExample)
33
+TME_plot(tsneData = result$tsneData, Class)
34
+}
0 35
new file mode 100644
... ...
@@ -0,0 +1,25 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/scTHI.R
3
+\docType{package}
4
+\name{scTHI}
5
+\alias{scTHI}
6
+\alias{scTHI-package}
7
+\title{single cell Tumor Hist Interaction (scTHI)}
8
+\description{
9
+Indentification of  significantly activated ligand-receptor 
10
+interactions across clusters of cells from single-cell RNA sequencing data.
11
+Single-cell RNA sequencing is the reference technique to characterize the heterogeneity of tumor microenvironment. 
12
+The composition of the various cell types making up the microenvironment can significantly affect the way 
13
+in which the immune system activates cancer rejection mechanisms. Understanding the cross-talk signals
14
+ between immune cells and cancer cells is a fundamental for the identification immuno-oncology therapeutic targets. 
15
+scTHI is  a novel method, single cell Tumor-Host Interaction tool (scTHI), 
16
+to identify significantly activated ligand-receptor interactions across clusters of cells from 
17
+single-cell RNA sequencing data. 
18
+cTHI is based on the hypothesis that when patterns of interaction are active, 
19
+they are also simultaneously and highly expressed in
20
+homogeneous cell populations. We also model the autocrine and paracrine 
21
+signalling effects of L-R partners
22
+}
23
+\details{
24
+Please have a look at the vignette for a in-depth introduction to the package.
25
+}
0 26
new file mode 100644
... ...
@@ -0,0 +1,33 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/scTHI_plotCluster.R
3
+\name{scTHI_plotCluster}
4
+\alias{scTHI_plotCluster}
5
+\title{scTHI_plotCluster}
6
+\usage{
7
+scTHI_plotCluster(scTHIresult, cexPoint = 0.8, legendPos = c("topleft",
8
+  "topright", "bottomright", "bottomleft"))
9
+}
10
+\arguments{
11
+\item{scTHIresult}{scTHI object.}
12
+
13
+\item{cexPoint}{Set the point size.}
14
+
15
+\item{legendPos}{Character string to custom the legend position.}
16
+}
17
+\value{
18
+None
19
+}
20
+\description{
21
+Graphs the output of scTHI_runTsne, labeling cells by clusters.
22
+}
23
+\examples{
24
+data(scExample)
25
+result <-  scTHI_score(scExample,
26
+                       cellCusterA = colnames(scExample)[1:30],
27
+                       cellCusterB = colnames(scExample)[31:100],
28
+                       cellCusterAName = "ClusterA",
29
+                       cellCusterBName = "ClusterB", filterCutoff = 0,
30
+                       pvalueCutoff = 1, nPermu = 100, ncore = 8)
31
+result <- scTHI_runTsne(result)
32
+scTHI_plotCluster(result)
33
+}
0 34
new file mode 100644
... ...
@@ -0,0 +1,37 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/scTHI_plotPairs.R
3
+\name{scTHI_plotPairs}
4
+\alias{scTHI_plotPairs}
5
+\title{scTHI_plotPairs}
6
+\usage{
7
+scTHI_plotPairs(scTHIresult, cexPoint = 0.8, interactionToplot)
8
+}
9
+\arguments{
10
+\item{scTHIresult}{scTHI object.}
11
+
12
+\item{cexPoint}{Set the point size.}
13
+
14
+\item{interactionToplot}{Interaction pair to plot.}
15
+}
16
+\value{
17
+None
18
+}
19
+\description{
20
+Generates a plot on the t-SNE coordinates to
21
+show the expression levels
22
+of an interaction pair of interest. Each cell is
23
+colored according to the
24
+corresponding
25
+gene expression value.
26
+}
27
+\examples{
28
+data(scExample)
29
+result <-  scTHI_score(scExample,
30
+                 cellCusterA = colnames(scExample)[1:30],
31
+                 cellCusterB = colnames(scExample)[31:100],
32
+                 cellCusterAName = "ClusterA",
33
+                 cellCusterBName = "ClusterB", filterCutoff = 0,
34
+                 pvalueCutoff = 1, nPermu = 100, ncore = 8)
35
+result <- scTHI_runTsne(result)
36
+scTHI_plotPairs(result,interactionToplot = "CXCL12_CD4")
37
+}
0 38
new file mode 100644
... ...
@@ -0,0 +1,43 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/scTHI_plotResult.R
3
+\name{scTHI_plotResult}
4
+\alias{scTHI_plotResult}
5
+\title{scTHI_plotResult}
6
+\usage{
7
+scTHI_plotResult(scTHIresult, cexNames = 0.8, plotType = c("score",
8
+  "pair"), nRes = NULL)
9
+}
10
+\arguments{
11
+\item{scTHIresult}{scTHI object.}
12
+
13
+\item{cexNames}{Size of names in barplot.}
14
+
15
+\item{plotType}{Type of plot to be generated. Default is
16
+"score", can be  also "pair". The "score" option will generate
17
+a barplot for each resulted interaction pair, representing
18
+the calculated interaction score
19
+  and the related p-Value.The "pair" option will generate
20
+  two barplot for each resulted interaction pair, representing
21
+  the percentage of cells of each cluster expressing partnerA
22
+  and partnerB gene, respectively.}
23
+
24
+\item{nRes}{Number of pairs to plot (all if NULL).}
25
+}
26
+\value{
27
+None
28
+}
29
+\description{
30
+Creates barplots of scTHI_score results.
31
+}
32
+\examples{
33
+data(scExample)
34
+result <-  scTHI_score(scExample,
35
+                       cellCusterA = colnames(scExample)[1:30],
36
+                       cellCusterB = colnames(scExample)[31:100],
37
+                       cellCusterAName = "ClusterA",
38
+                       cellCusterBName = "ClusterB", filterCutoff = 0,
39
+                       pvalueCutoff = 1, nPermu = 100, ncore = 8)
40
+
41
+scTHI_plotResult(result, plotType = "score")
42
+scTHI_plotResult(result, plotType = "pair")
43
+}
0 44
new file mode 100644
... ...
@@ -0,0 +1,29 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/scTHI_runTsne.R
3
+\name{scTHI_runTsne}
4
+\alias{scTHI_runTsne}
5
+\title{scTHI_runTsne}
6
+\usage{
7
+scTHI_runTsne(scTHIresult)
8
+}
9
+\arguments{
10
+\item{scTHIresult}{scTHI object.}
11
+}
12
+\value{
13
+The same object as scTHI_score with a fifth item tsneData
14
+  (data.frame)
15
+}
16
+\description{
17
+Runs t-SNE dimensionality reduction on selected features based on Rtsne
18
+package.
19
+}
20
+\examples{
21
+data(scExample)
22
+result <-  scTHI_score(scExample,
23
+                       cellCusterA = colnames(scExample)[1:30],
24
+                       cellCusterB = colnames(scExample)[31:100],
25
+                       cellCusterAName = "ClusterA",
26
+                       cellCusterBName = "ClusterB", filterCutoff = 0,
27
+                       pvalueCutoff = 1, nPermu = 100, ncore = 8)
28
+result <- scTHI_runTsne(result)
29
+}
0 30
new file mode 100644
... ...
@@ -0,0 +1,73 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/scTHI_score.R
3
+\name{scTHI_score}
4
+\alias{scTHI_score}
5
+\title{scTHI_score}
6
+\usage{
7
+scTHI_score(expMat, cellCusterA, cellCusterB, cellCusterAName,
8
+  cellCusterBName, topRank = 10, autocrineEffect = TRUE,
9
+  fileNameBase = "scTHI", filterCutoff = 0.5, PValue = TRUE,
10
+  pvalueCutoff = 0.05, nPermu = 1000, ncore = 8)
11
+}
12
+\arguments{
13
+\item{expMat}{ScRNA-seq gene expression matrix where rows are genes
14
+presented with Hugo Symbols and columns are cells. Gene expression
15
+values should be counts or normalized counts.}
16
+
17
+\item{cellCusterA}{Vector of columns of expMat that belong to the
18
+first cluster.}
19
+
20
+\item{cellCusterB}{Vector of columns of expMat that belong to the
21
+second cluster.}
22
+
23
+\item{cellCusterAName}{A character string labeling the clusterA.}
24
+
25
+\item{cellCusterBName}{A character string labeling the clusterB.}
26
+
27
+\item{topRank}{Filter threshold. Set to 10 (default) means that
28
+each gene of the interaction pair will be considered as expressed
29
+in a cell if it's in the top rank 10 percent.}
30
+
31
+\item{autocrineEffect}{if TRUE remove the paracrine filter}
32
+
33
+\item{fileNameBase}{Project name.}
34
+
35
+\item{filterCutoff}{Score threshold (default is 0.50). For each
36
+interaction pair, if the score calculated (for the partnerA
37
+or partnerB)
38
+will be less than filterCutoff the interaction pair will be
39
+discarded.}
40
+
41
+\item{PValue}{Logical, set to TRUE (default) compute statistical
42
+iterations. If p.value < 0.05, the value will be returned.}
43
+
44
+\item{pvalueCutoff}{cutoff of the p-value}
45
+
46
+\item{nPermu}{Number of iterations to perform (default is 1000).}
47
+
48
+\item{ncore}{Number of processors to use.}
49
+}
50
+\value{
51
+A list of results, with four items: result (data.frame),
52
+  expMat (matrix), clusterA (character),  clusterA (character)
53
+}
54
+\description{
55
+This function allows the user to compute a score for a set of
56
+ligand-receptor pairs, from a single cell gene expression matrix,
57
+and detect specific Tumor-Host interactions. You must specify at
58
+least two clusters of cells (for example tumor cells and immune
59
+cells).
60
+}
61
+\examples{
62
+
63
+####################### example of scTHI_score
64
+data(scExample)
65
+result <-  scTHI_score(scExample,
66
+      cellCusterA = colnames(scExample)[1:30],
67
+      cellCusterB = colnames(scExample)[31:100],
68
+      cellCusterAName = "ClusterA",
69
+      cellCusterBName = "ClusterB", filterCutoff = 0,
70
+     pvalueCutoff = 1, nPermu = 100, ncore = 8)
71
+
72
+}
73
+\keyword{interaction}