Spatially resolved transcriptomics identifies intercellular signaling post-ischemic stroke that controls neural stem cell proliferation.

2. #' Compute the communication probability/strength between any interacting cell groups
3. #'
4. #' To further speed up on large-scale datasets, USER can downsample the data using the function 'subset' from Seurat package (e.g., pbmc.small <- subset(pbmc, downsample = 500)), or using the function `sketchData` from CellChat, in particular for the large cell clusters;
5. #'
6. #'
7. #' @param object CellChat object
8. #' @param type Methods for computing the average gene expression per cell group. By default = "triMean", producing fewer but stronger interactions;
9. #' When setting `type = "truncatedMean"`, a value should be assigned to 'trim', producing more interactions.
10. #' @param trim the fraction (0 to 0.25) of observations to be trimmed from each end of x before the mean is computed
11. #' @param LR.use A subset of ligand-receptor interactions used in inferring communication network
12. #' @param raw.use Whether use the raw data (i.e., `[email hidden]`) or the smoothed data (i.e., `[email hidden]`).
13. #' Set raw.use = FALSE to use the projected data when analyzing single-cell data with shallow sequencing depth because the projected data could help to reduce the dropout effects of signaling genes, in particular for possible zero expression of subunits of ligands/receptors.
14. #' @param population.size Whether consider the proportion of cells in each group across all sequenced cells.
15. #' Set population.size = FALSE if analyzing sorting-enriched single cells, to remove the potential artifact of population size.
16. #' Set population.size = TRUE if analyzing unsorted single-cell transcriptomes, with the reason that abundant cell populations tend to send collectively stronger signals than the rare cell populations.
17. #'
18. #' Parameters for spatial data analysis:
19. #' @param distance.use Whether to use distance constraints to compute communication probability. Setting `distance.use = TRUE` indicates that the cell-cell communication probability is inversely proportional to the computed distance.
20. #' Setting `distance.use = FALSE` will only filter out interactions between spatially distant regions, but not add distance constraints.
21. #' @param interaction.range The maximum interaction/diffusion length of ligands (Unit: microns). This hard threshold is used to filter out the connections between spatially distant regions
22. #' @param scale.distance A scale or normalization factor for the spatial distances when setting `distance.use = TRUE`. For example, scale.distance equals 1, 0.1, 0.01, 0.001, 0.11, or 0.011. We choose this values such that the minimum value of the scaled distances is in [1,2]. This value is not necessary when setting `distance.use = FALSE`.
23. #'
24. #' When comparing communication across different CellChat objects, the same scale factor should be used. For a single CellChat analysis, different scale factors will not affect the ranking of the signaling based on their interaction strength.
25. #'
26. #' @param k.min The minimum number of interacting cell pairs required for defining spatially proximal cell groups.
27. #' @param contact.dependent Whether using the `contact-dependent` manner for inference signaling, that is determining interacting cell pairs by requiring cells to be in direct membrane-membrane contact. By default `contact.dependent = TRUE` when inferring contact-dependent and juxtacrine signaling (that is "Cell-Cell Contact" signaling classified in CellChatDB$interaction$annotation).
28. #' If only focusing on `Secreted Signaling`, the `contact-dependent` manner will be not used except for setting `contact.dependent.forced = TRUE`.
29. #' @param contact.range The interaction range (Unit: microns) to restrict the contact-dependent signaling when `contact.dependent = TRUE`.
30. #' For spatial transcriptomics in a single-cell resolution, `contact.range` is approximately equal to the estimated cell diameter (i.e., the cell center-to-center distance), which means that contact-dependent and juxtacrine signaling can only happens when the two cells are contact to each other.
31. #'
32. #' Typically, `contact.range = 10`, which is a typical human cell size. However, for low-resolution spatial data such as 10X visium, it should be the cell center-to-center distance (i.e., `contact.range = 100` for visium data). The function `computeCellDistance` can compute the center-to-center distance.
33. #'
34. #' @param contact.knn.k Number of neighbors to restrict the contact-dependent signaling within the neatest neighbors when `contact.dependent = TRUE`. By default, CellChat uses `contact.range` to restrict the contact-dependent signaling; however, users can also provide a value of `contact.knn.k`, in order to determine interacting cell pairs based on the k-nearest neighbors (knn).
35. #' For 10X visium, contact.knn.k = 6. For other spatial technologies, this value may be hard to determine because the sequenced cells/spots are usually not regularly arranged.
36. #' @param do.symmetric Whether converting the adjacent matrix into symmetric one when determining spatially proximal cell groups. Default is TRUE, indicating that if adj(i,j) or adj(j,i) is zero, then both are zeros.
37. #'
38. #' @param contact.dependent.forced Whether forcing to use the `contact-dependent` manner for inference signaling for all L-R pairs including secreted signaling. Users can set `contact.dependent.forced = TRUE` if also preferring interactions within a contact manner for `Secreted Signaling`.
39. #'
40. #' @param nboot Threshold of p-values
41. #' @param seed.use Set a random seed. By default, set the seed to 1.
42. #' @param Kh Parameter in Hill function
43. #' @param n Parameter in Hill function
44. #'
45. #'
46. #' @importFrom future nbrOfWorkers
47. #' @importFrom future.apply future_sapply
48. #' @importFrom pbapply pbsapply
49. #' @importFrom stats aggregate
50. #' @importFrom Matrix crossprod
51. #' @importFrom utils txtProgressBar setTxtProgressBar
52. #'
53. #' @return A CellChat object with updated slot 'net':
54. #'
55. #' object@net$prob is the inferred communication probability (strength) array, where the first, second and third dimensions represent a source, target and ligand-receptor pair, respectively.
56. #'
57. #' USER can access all the inferred cell-cell communications using the function 'subsetCommunication(object)', which returns a data frame.
58. #'
59. #' object@net$pval is the corresponding p-values of each interaction
60. #'
61. #' @export
62. #'
63. computeCommunProb <- function(object, type = c("triMean", "truncatedMean","thresholdedMean", "median"), trim = 0.1, LR.use = NULL, raw.use = TRUE, population.size = FALSE,
64. distance.use = TRUE, interaction.range = 250, scale.distance = 0.01, k.min = 10, contact.dependent = TRUE, contact.range = NULL, contact.knn.k = NULL, contact.dependent.forced = FALSE, do.symmetric = TRUE,
65. nboot = 100, seed.use = 1L, Kh = 0.5, n = 1) {
66. type <- match.arg(type)
67. cat(type, "is used for calculating the average gene expression per cell group.", "\n")
68. FunMean <- switch(type,
69. triMean = triMean,
70. truncatedMean = function(x) mean(x, trim = trim, na.rm = TRUE),
71. thresholdedMean = function(x) thresholdedMean(x, trim = trim, na.rm = TRUE),
72. median = function(x) median(x, na.rm = TRUE))
74. if (raw.use) {
75. data <- as.matrix([email hidden])
76. } else {
77. data <- as.matrix([email hidden])
78. }
79. if (is.null(LR.use)) {
80. pairLR.use <- object@LR$LRsig
81. } else {
82. if (length(unique(LR.use$annotation)) > 1) {
83. LR.use$annotation <- factor(LR.use$annotation, levels = c("Secreted Signaling","ECM-Receptor", "Non-protein Signaling", "Cell-Cell Contact"))
84. LR.use <- LR.use[order(LR.use$annotation), , drop = FALSE]
85. LR.use$annotation <- as.character(LR.use$annotation)
86. }
87. pairLR.use <- LR.use
88. }
89. complex_input <- object@DB$complex
90. cofactor_input <- object@DB$cofactor
91. my.sapply <- ifelse(
92. test = future::nbrOfWorkers() == 1,
93. yes = sapply,
94. no = future.apply::future_sapply
95. )
97. ptm = Sys.time()
99. pairLRsig <- pairLR.use
100. group <- object@idents
101. geneL <- as.character(pairLRsig$ligand)
102. geneR <- as.character(pairLRsig$receptor)
103. nLR <- nrow(pairLRsig)
104. numCluster <- nlevels(group)
105. if (numCluster != length(unique(group))) {
106. stop("Please check `unique(object@idents)` and ensure that the factor levels are correct!
107. You may need to drop unused levels using 'droplevels' function. e.g.,
108. `meta$labels = droplevels(meta$labels, exclude = setdiff(levels(meta$labels),unique(meta$labels)))`")
109. }
111. data.use <- data/max(data)
112. nC <- ncol(data.use)
114. # compute the average expression per group
115. data.use.avg <- aggregate(t(data.use), list(group), FUN = FunMean)
116. data.use.avg <- t(data.use.avg[,-1])
117. colnames(data.use.avg) <- levels(group)
118. # compute the expression of ligand or receptor
119. dataLavg <- computeExpr_LR(geneL, data.use.avg, complex_input)
120. dataRavg <- computeExpr_LR(geneR, data.use.avg, complex_input)
121. # take account into the effect of co-activation and co-inhibition receptors
122. dataRavg.co.A.receptor <- computeExpr_coreceptor(cofactor_input, data.use.avg, pairLRsig, type = "A")
123. dataRavg.co.I.receptor <- computeExpr_coreceptor(cofactor_input, data.use.avg, pairLRsig, type = "I")
124. dataRavg <- dataRavg * dataRavg.co.A.receptor/dataRavg.co.I.receptor
126. dataLavg2 <- t(replicate(nrow(dataLavg), as.numeric(table(group))/nC))
127. dataRavg2 <- dataLavg2
129. # compute the expression of agonist and antagonist
130. index.agonist <- which(!is.na(pairLRsig$agonist) & pairLRsig$agonist != "")
131. index.antagonist <- which(!is.na(pairLRsig$antagonist) & pairLRsig$antagonist != "")
132. # quantify the communication probability
134. # compute the spatial constraint
135. if (object@options$datatype != "RNA") {
136. data.spatial <- object@images$coordinates
137. if ("spatial.factors" %in% names(object@images)) {
138. ratio <- object@images$spatial.factors$ratio
139. tol <- object@images$spatial.factors$tol
140. } else {
141. stop("`object@images$spatial.factors` is missing. Please update the object via `updateCellChat`! \n")
142. }
144. meta.t = data.frame(group = group, samples = object@meta$samples, row.names = rownames(object@meta))
145. res <- computeRegionDistance(coordinates = data.spatial, meta = meta.t, interaction.range = interaction.range, ratio = ratio, tol = tol, k.min = k.min, contact.dependent = contact.dependent, contact.range = contact.range, contact.knn.k = contact.knn.k)
146. d.spatial <- res$d.spatial # NaN if no nearby cell pairs
147. adj.contact <- res$adj.contact # zeros if no nearby cell pairs
148. if (distance.use) {
149. print(paste0('>>> Run CellChat on spatial transcriptomics data using distances as constraints of the computed communication probability <<< [', Sys.time(),']'))
150. d.spatial <- d.spatial * scale.distance
151. diag(d.spatial) <- NaN
152. d.min <- min(d.spatial, na.rm = TRUE)
153. if (d.min < 1) {
154. cat("The suggested minimum value of scaled distances is in [1,2], and the calculated value here is ", d.min,"\n")
155. stop("Please increase the value of `scale.distance` and use a value that is slighly smaller than ", format(1/d.min, digits = 2) ,"\n")
156. }
157. P.spatial <- 1/d.spatial
158. P.spatial[is.na(d.spatial)] <- 0
159. diag(P.spatial) <- max(P.spatial) # if this value is 1, the self-connections will have more larger weight.
160. d.spatial <- d.spatial/scale.distance # This is only for saving the data
161. } else {
162. print(paste0('>>> Run CellChat on spatial transcriptomics data without distance values as constraints of the computed communication probability <<< [', Sys.time(),']'))
163. P.spatial <- matrix(1, nrow = numCluster, ncol = numCluster)
164. P.spatial[is.na(d.spatial)] <- 0 # diagonal is 1
165. }
167. } else {
168. print(paste0('>>> Run CellChat on sc/snRNA-seq data <<< [', Sys.time(),']'))
169. d.spatial <- matrix(NaN, nrow = numCluster, ncol = numCluster)
170. P.spatial <- matrix(1, nrow = numCluster, ncol = numCluster)
171. adj.contact <- matrix(1, nrow = numCluster, ncol = numCluster)
172. contact.dependent = FALSE; contact.dependent.forced = FALSE; contact.range = NULL; contact.knn.k = NULL;
173. distance.use = NULL; interaction.range = NULL; ratio = NULL; tol = NULL; k.min = NULL;
174. }
176. if (object@options$datatype == "RNA") {
177. nLR1 <- nLR
178. } else {
179. if (contact.dependent.forced == TRUE) {
180. cat("Force to run CellChat in a `contact-dependent` manner for all L-R pairs including secreted signaling.\n")
181. P.spatial <- P.spatial * adj.contact
182. nLR1 <- nLR
183. } else { # contact.dependent.forced == F
184. if (contact.dependent == TRUE && length(unique(pairLRsig$annotation)) > 0) {
185. if (all(unique(pairLRsig$annotation) %in% c("Cell-Cell Contact"))) {
186. cat("All the input L-R pairs are `Cell-Cell Contact` signaling. Run CellChat in a contact-dependent manner. \n")
187. P.spatial <- P.spatial * adj.contact
188. nLR1 <- nLR
189. } else if (all(unique(pairLRsig$annotation) %in% c("Secreted Signaling", "ECM-Receptor", "Non-protein Signaling"))) {
190. cat("Molecules of the input L-R pairs are diffusible. Run CellChat in a diffusion manner based on the `interaction.range`.\n")
191. nLR1 <- nLR
192. } else {
193. cat("The input L-R pairs have both secreted signaling and contact-dependent signaling. Run CellChat in a contact-dependent manner for `Cell-Cell Contact` signaling, and in a diffusion manner based on the `interaction.range` for other L-R pairs. \n")
194. nLR1 <- max(which(pairLRsig$annotation %in% c("Secreted Signaling", "ECM-Receptor", "Non-protein Signaling")))
195. }
196. } else { # contact.dependent == F or there is no `annotation` column in the database
197. cat("Run CellChat in a diffusion manner based on the `interaction.range` for all L-R pairs. Setting `contact.dependent = TRUE` if preferring a contact-dependent manner for `Cell-Cell Contact` signaling. \n")
198. nLR1 <- nLR
199. }
200. }
201. }
203. Prob <- array(0, dim = c(numCluster,numCluster,nLR))
204. Pval <- array(0, dim = c(numCluster,numCluster,nLR))
206. set.seed(seed.use)
207. permutation <- replicate(nboot, sample.int(nC, size = nC))
208. data.use.avg.boot <- my.sapply(
209. X = 1:nboot,
210. FUN = function(nE) {
211. groupboot <- group[permutation[, nE]]
212. data.use.avgB <- aggregate(t(data.use), list(groupboot), FUN = FunMean)
213. data.use.avgB <- t(data.use.avgB[,-1])
214. return(data.use.avgB)
215. },
216. simplify = FALSE
217. )
218. pb <- txtProgressBar(min = 0, max = nLR, style = 3, file = stderr())
220. for (i in 1:nLR) {
221. # ligand/receptor
222. dataLR <- Matrix::crossprod(matrix(dataLavg[i,], nrow = 1), matrix(dataRavg[i,], nrow = 1))
223. P1 <- dataLR^n/(Kh^n + dataLR^n)
224. P1_Pspatial <- P1*P.spatial
225. if (sum(P1_Pspatial) == 0) {
226. Pnull = P1_Pspatial
227. Prob[ , , i] <- Pnull
228. p = 1
229. Pval[, , i] <- matrix(p, nrow = numCluster, ncol = numCluster, byrow = FALSE)
230. } else {
231. if (i > nLR1) {
232. P.spatial <- P.spatial * adj.contact
233. }
234. # agonist and antagonist
235. if (is.element(i, index.agonist)) {
236. data.agonist <- computeExpr_agonist(data.use = data.use.avg, pairLRsig, cofactor_input, index.agonist = i, Kh = Kh, n = n)
237. P2 <- Matrix::crossprod(matrix(data.agonist, nrow = 1))
238. } else {
239. P2 <- matrix(1, nrow = numCluster, ncol = numCluster)
240. }
241. if (is.element(i, index.antagonist)) {
242. data.antagonist <- computeExpr_antagonist(data.use = data.use.avg, pairLRsig, cofactor_input, index.antagonist = i, Kh = Kh, n = n)
243. P3 <- Matrix::crossprod(matrix(data.antagonist, nrow = 1))
244. } else {
245. P3 <- matrix(1, nrow = numCluster, ncol = numCluster)
246. }
247. # number of cells
248. if (population.size) {
249. P4 <- Matrix::crossprod(matrix(dataLavg2[i,], nrow = 1), matrix(dataRavg2[i,], nrow = 1))
250. } else {
251. P4 <- matrix(1, nrow = numCluster, ncol = numCluster)
252. }
254. # Pnull = P1*P2*P3*P4
255. Pnull = P1*P2*P3*P4*P.spatial
256. Prob[ , , i] <- Pnull
258. Pnull <- as.vector(Pnull)
260. #Pboot <- foreach(nE = 1:nboot) %dopar% {
261. Pboot <- my.sapply(
262. X = 1:nboot,
263. FUN = function(nE) {
264. data.use.avgB <- data.use.avg.boot[[nE]]
265. dataLavgB <- computeExpr_LR(geneL[i], data.use.avgB, complex_input)
266. dataRavgB <- computeExpr_LR(geneR[i], data.use.avgB, complex_input)
267. # take account into the effect of co-activation and co-inhibition receptors
268. dataRavgB.co.A.receptor <- computeExpr_coreceptor(cofactor_input, data.use.avgB, pairLRsig[i, , drop = FALSE], type = "A")
269. dataRavgB.co.I.receptor <- computeExpr_coreceptor(cofactor_input, data.use.avgB, pairLRsig[i, , drop = FALSE], type = "I")
270. dataRavgB <- dataRavgB * dataRavgB.co.A.receptor/dataRavgB.co.I.receptor
271. dataLRB = Matrix::crossprod(dataLavgB, dataRavgB)
272. P1.boot <- dataLRB^n/(Kh^n + dataLRB^n)
273. # agonist and antagonist
274. if (is.element(i, index.agonist)) {
275. data.agonist <- computeExpr_agonist(data.use = data.use.avgB, pairLRsig, cofactor_input, index.agonist = i, Kh = Kh, n = n)
276. P2.boot <- Matrix::crossprod(matrix(data.agonist, nrow = 1))
277. } else {
278. P2.boot <- matrix(1, nrow = numCluster, ncol = numCluster)
279. }
280. if (is.element(i, index.antagonist)) {
281. data.antagonist <- computeExpr_antagonist(data.use = data.use.avgB, pairLRsig, cofactor_input, index.antagonist = i, Kh = Kh, n= n)
282. P3.boot <- Matrix::crossprod(matrix(data.antagonist, nrow = 1))
283. } else {
284. P3.boot <- matrix(1, nrow = numCluster, ncol = numCluster)
285. }
287. if (population.size) {
288. groupboot <- group[permutation[, nE]]
289. dataLavg2B <- as.numeric(table(groupboot))/nC
290. dataLavg2B <- matrix(dataLavg2B, nrow = 1)
291. dataRavg2B <- dataLavg2B
292. P4.boot = Matrix::crossprod(dataLavg2B, dataRavg2B)
293. } else {
294. P4.boot = matrix(1, nrow = numCluster, ncol = numCluster)
295. }
297. # Pboot = P1.boot*P2.boot*P3.boot*P4.boot
298. Pboot = P1.boot*P2.boot*P3.boot*P4.boot*P.spatial
299. return(as.vector(Pboot))
300. }
301. )
302. Pboot <- matrix(unlist(Pboot), nrow=length(Pnull), ncol = nboot, byrow = FALSE)
303. nReject <- rowSums(Pboot - Pnull > 0)
304. p = nReject/nboot
305. Pval[, , i] <- matrix(p, nrow = numCluster, ncol = numCluster, byrow = FALSE)
306. }
307. setTxtProgressBar(pb = pb, value = i)
308. }
309. close(con = pb)
310. Pval[Prob == 0] <- 1
311. dimnames(Prob) <- list(levels(group), levels(group), rownames(pairLRsig))
312. dimnames(Pval) <- dimnames(Prob)
313. net <- list("prob" = Prob, "pval" = Pval)
314. execution.time = Sys.time() - ptm
315. object@options$run.time <- as.numeric(execution.time, units = "secs")
317. object@options$parameter <- list(type.mean = type, trim = trim, raw.use = raw.use, population.size = population.size, nboot = nboot, seed.use = seed.use, Kh = Kh, n = n,
318. distance.use = distance.use, interaction.range = interaction.range, ratio = ratio, tol = tol, k.min = k.min,
319. contact.dependent = contact.dependent, contact.range = contact.range, contact.knn.k = contact.knn.k, contact.dependent.forced = contact.dependent.forced
320. )
321. if (object@options$datatype != "RNA") {
322. object@images$distance <- d.spatial
323. }
324. object@net <- net
325. print(paste0('>>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [', Sys.time(),']'))
326. return(object)
327. }
330. #' Compute the communication probability on signaling pathway level by summarizing all related ligands/receptors
331. #'
332. #' @param object CellChat object
333. #' @param net A list from object@net; If net = NULL, net = object@net
334. #' @param pairLR.use A dataframe giving the ligand-receptor interactions; If pairLR.use = NULL, pairLR.use = object@LR$LRsig
335. #' @param thresh threshold of the p-value for determining significant interaction
336. #'
337. #' @return A CellChat object with updated slot 'netP':
338. #'
339. #' object@netP$prob is the communication probability array on signaling pathway level; USER can convert this array to a data frame using the function 'reshape2::melt()',
340. #'
341. #' e.g., `df.netP <- reshape2::melt(object@netP$prob, value.name = "prob"); colnames(df.netP)[1:3] <- c("source","target","pathway_name")` or access all significant interactions using the function \code{\link{subsetCommunication}}
342. #'
343. #' object@netP$pathways list all the signaling pathways with significant communications.
344. #'
345. #' From version >= 1.1.0, pathways are ordered based on the total communication probabilities. NB: pathways with small total communication probabilities might be also very important since they might be specifically activated between only few cell types.
346. #'
347. #' @export
348. #'
349. computeCommunProbPathway <- function(object = NULL, net = NULL, pairLR.use = NULL, thresh = 0.05) {
350. if (is.null(net)) {
351. net <- object@net
352. }
353. if (is.null(pairLR.use)) {
354. pairLR.use <- object@LR$LRsig
355. }
356. prob <- net$prob
357. prob[net$pval > thresh] <- 0
359. LR <- dimnames(prob)[[3]]
360. LR.sig <- LR[apply(prob, 3, sum) != 0]
362. pathways <- unique(pairLR.use$pathway_name)
363. group <- factor(pairLR.use$pathway_name, levels = pathways)
364. prob.pathways <- aperm(apply(prob, c(1, 2), by, group, sum), c(2, 3, 1))
365. pathways.sig <- pathways[apply(prob.pathways, 3, sum) != 0]
366. prob.pathways.sig <- prob.pathways[,,pathways.sig, drop = FALSE]
367. idx <- sort(apply(prob.pathways.sig, 3, sum), decreasing=TRUE, index.return = TRUE)$ix
368. pathways.sig <- pathways.sig[idx]
369. prob.pathways.sig <- prob.pathways.sig[, , idx]
371. if (is.null(object)) {
372. netP = list(pathways = pathways.sig, prob = prob.pathways.sig)
373. return(netP)
374. } else {
375. object@net$LRs <- LR.sig
376. object@netP$pathways <- pathways.sig
377. object@netP$prob <- prob.pathways.sig
378. return(object)
379. }
380. }
383. #' Calculate the aggregated network by counting the number of links or summarizing the communication probability
384. #'
385. #' @param object CellChat object
386. #' @param sources.use,targets.use,signaling,pairLR.use Please check the description in function \code{\link{subsetCommunication}}
387. #' @param remove.isolate whether removing the isolate cell groups without any interactions when applying \code{\link{subsetCommunication}}
388. #' @param thresh threshold of the p-value for determining significant interaction
389. #' @param return.object whether return an updated CellChat object
390. #' @importFrom dplyr group_by summarize groups
391. #' @importFrom stringr str_split
392. #'
393. #' @return Return an updated CellChat object:
394. #'
395. #' `object@net$count` is a matrix: rows and columns are sources and targets respectively, and elements are the number of interactions between any two cell groups. USER can convert a matrix to a data frame using the function `reshape2::melt()`
396. #'
397. #' `object@net$weight` is also a matrix containing the interaction weights between any two cell groups
398. #'
399. #' `object@net$sum` is deprecated. Use `object@net$weight`
400. #'
401. #' @export
402. #'
403. aggregateNet <- function(object, sources.use = NULL, targets.use = NULL, signaling = NULL, pairLR.use = NULL, remove.isolate = TRUE, thresh = 0.05, return.object = TRUE) {
404. net <- object@net
405. if (is.null(sources.use) & is.null(targets.use) & is.null(signaling) & is.null(pairLR.use)) {
406. prob <- net$prob
407. pval <- net$pval
408. pval[prob == 0] <- 1
409. prob[pval >= thresh] <- 0
410. net$count <- apply(prob > 0, c(1,2), sum)
411. net$weight <- apply(prob, c(1,2), sum)
412. net$weight[is.na(net$weight)] <- 0
413. net$count[is.na(net$count)] <- 0
414. } else {
415. df.net <- subsetCommunication(object, slot.name = "net",
416. sources.use = sources.use, targets.use = targets.use,
417. signaling = signaling,
418. pairLR.use = pairLR.use,
419. thresh = thresh)
420. df.net$source_target <- paste(df.net$source, df.net$target, sep = "|")
421. df.net2 <- df.net %>% group_by(source_target) %>% summarize(count = n(), .groups = 'drop')
422. df.net3 <- df.net %>% group_by(source_target) %>% summarize(prob = sum(prob), .groups = 'drop')
423. df.net2$prob <- df.net3$prob
424. a <- stringr::str_split(df.net2$source_target, "|", simplify = T)
425. df.net2$source <- as.character(a[, 1])
426. df.net2$target <- as.character(a[, 2])
427. cells.level <- levels(object@idents)
428. if (remove.isolate) {
429. message("Isolate cell groups without any interactions are removed. To block it, set `remove.isolate = FALSE`")
430. df.net2$source <- factor(df.net2$source, levels = cells.level[cells.level %in% unique(df.net2$source)])
431. df.net2$target <- factor(df.net2$target, levels = cells.level[cells.level %in% unique(df.net2$target)])
432. } else {
433. df.net2$source <- factor(df.net2$source, levels = cells.level)
434. df.net2$target <- factor(df.net2$target, levels = cells.level)
435. }
437. count <- tapply(df.net2[["count"]], list(df.net2[["source"]], df.net2[["target"]]), sum)
438. prob <- tapply(df.net2[["prob"]], list(df.net2[["source"]], df.net2[["target"]]), sum)
439. net$count <- count
440. net$weight <- prob
441. net$weight[is.na(net$weight)] <- 0
442. net$count[is.na(net$count)] <- 0
443. }
444. if (return.object) {
445. object@net <- net
446. return(object)
447. } else {
448. return(net)
449. }
451. }
454. #' Compute averaged expression values for each cell group
455. #'
456. #' @param object CellChat object
457. #' @param features a char vector giving the used features. default use all features
458. #' @param group.by cell group information; default is `object@idents` when input is a single object and `object@idents$joint` when input is a merged object; otherwise it should be one of the column names of the meta slot
459. #' @param type methods for computing the average gene expression per cell group.
460. #'
461. #' By default = "triMean", defined as a weighted average of the distribution's median and its two quartiles (https://en.wikipedia.org/wiki/Trimean);
462. #'
463. #' When setting `type = "truncatedMean"`, a value should be assigned to 'trim'. See the function `base::mean`.
464. #'
465. #' @param trim the fraction (0 to 0.25) of observations to be trimmed from each end of x before the mean is computed.
466. #' @param slot.name the data in the slot.name to use
467. #' @param data.use a customed data matrix. Default: data.use = NULL and the expression matrix in the 'slot.name' is used
468. #'
469. #' @return Returns a matrix with genes as rows, cell groups as columns.
471. #' @export
472. #'
473. computeAveExpr <- function(object, features = NULL, group.by = NULL, type = c("triMean", "truncatedMean", "median"), trim = NULL,
474. slot.name = c("data.signaling", "data"), data.use = NULL) {
475. type <- match.arg(type)
476. slot.name <- match.arg(slot.name)
477. FunMean <- switch(type,
478. triMean = triMean,
479. truncatedMean = function(x) mean(x, trim = trim, na.rm = TRUE),
480. median = function(x) median(x, na.rm = TRUE))
482. if (is.null(data.use)) {
483. data.use <- slot(object, slot.name)
484. }
485. if (is.null(features)) {
486. features.use <- row.names(data.use)
487. } else {
488. features.use <- intersect(features, row.names(data.use))
489. }
490. data.use <- data.use[features.use, , drop = FALSE]
491. data.use <- as.matrix(data.use)
493. if (is.null(group.by)) {
494. labels <- object@idents
495. if (!is.factor(labels)) {
496. message("Use the joint cell labels from the merged CellChat object")
497. labels <- object@idents$joint
498. }
499. } else {
500. labels <- object@meta[[group.by]]
501. }
502. if (!is.factor(labels)) {
503. labels <- factor(labels)
504. }
505. # compute the average expression per group
506. data.use.avg <- aggregate(t(data.use), list(labels), FUN = FunMean)
507. data.use.avg <- t(data.use.avg[,-1])
508. rownames(data.use.avg) <- features.use
509. colnames(data.use.avg) <- levels(labels)
510. return(data.use.avg)
511. }
515. #' Compute the expression of complex in individual cells using geometric mean
516. #' @param complex_input the complex_input from CellChatDB
517. #' @param data.use data matrix (row are genes and columns are cells or cell groups)
518. #' @param complex the names of complex
519. #' @return
520. #' @importFrom dplyr select starts_with
521. #' @importFrom future nbrOfWorkers
522. #' @importFrom future.apply future_sapply
523. #' @importFrom pbapply pbsapply
524. #' @export
525. computeExpr_complex <- function(complex_input, data.use, complex) {
526. Rsubunits <- complex_input[complex,] %>% dplyr::select(starts_with("subunit"))
527. my.sapply <- ifelse(
528. test = future::nbrOfWorkers() == 1,
529. yes = sapply,
530. no = future.apply::future_sapply
531. )
532. data.complex = my.sapply(
533. X = 1:nrow(Rsubunits),
534. FUN = function(x) {
535. RsubunitsV <- unlist(Rsubunits[x,], use.names = F)
536. RsubunitsV <- RsubunitsV[RsubunitsV != ""]
537. return(geometricMean(data.use[RsubunitsV, , drop = FALSE]))
538. }
539. )
540. data.complex <- t(data.complex)
541. return(data.complex)
542. }
544. # Compute the average expression of complex per cell group using geometric mean
545. # @param complex_input the complex_input from CellChatDB
546. # @param data.use data matrix (rows are genes and columns are cells)
547. # @param complex the names of complex
548. # @param group a factor defining the cell groups
549. # @param FunMean the function for computing mean expression per group
550. # @return
551. # @importFrom dplyr select starts_with
552. # @importFrom future nbrOfWorkers
553. # @importFrom future.apply future_sapply
554. # @importFrom pbapply pbsapply
555. # #' @export
556. .computeExprGroup_complex <- function(complex_input, data.use, complex, group, FunMean) {
557. Rsubunits <- complex_input[complex,] %>% dplyr::select(starts_with("subunit"))
558. my.sapply <- ifelse(
559. test = future::nbrOfWorkers() == 1,
560. yes = pbapply::pbsapply,
561. no = future.apply::future_sapply
562. )
563. data.complex = my.sapply(
564. X = 1:nrow(Rsubunits),
565. FUN = function(x) {
566. RsubunitsV <- unlist(Rsubunits[x,], use.names = F)
567. RsubunitsV <- RsubunitsV[RsubunitsV != ""]
568. RsubunitsV <- intersect(RsubunitsV, rownames(data.use))
569. if (length(RsubunitsV) > 1) {
570. data.avg <- aggregate(t(data.use[RsubunitsV, ,drop = FALSE]), list(group), FUN = FunMean)
571. data.avg <- t(data.avg[,-1])
572. } else if (length(RsubunitsV) == 1) {
573. data.avg <- aggregate(matrix(data.use[RsubunitsV,], ncol = 1), list(group), FUN = FunMean)
574. data.avg <- t(data.avg[,-1])
575. } else {
576. data.avg = matrix(0, nrow = 1, ncol = length(unique(group)))
577. }
578. return(geometricMean(data.avg))
579. }
580. )
581. data.complex <- t(data.complex)
582. return(data.complex)
583. }
585. #' Compute the expression of ligands or receptors using geometric mean
586. #' @param geneLR a char vector giving a set of ligands or receptors
587. #' @param data.use data matrix (row are genes and columns are cells or cell groups)
588. #' @param complex_input the complex_input from CellChatDB
589. # #' @param group a factor defining the cell groups; If NULL, compute the expression of ligands or receptors in individual cells; otherwise, compute the average expression of ligands or receptors per cell group
590. # #' @param FunMean the function for computing average expression per cell group
591. #' @return
592. #' @export
593. computeExpr_LR <- function(geneLR, data.use, complex_input){
594. nLR <- length(geneLR)
595. numCluster <- ncol(data.use)
596. index.singleL <- which(geneLR %in% rownames(data.use))
597. dataL1avg <- data.use[geneLR[index.singleL],]
598. dataLavg <- matrix(nrow = nLR, ncol = numCluster)
599. dataLavg[index.singleL,] <- dataL1avg
600. index.complexL <- setdiff(1:nLR, index.singleL)
601. if (length(index.complexL) > 0) {
602. complex <- geneLR[index.complexL]
603. data.complex <- computeExpr_complex(complex_input, data.use, complex)
604. dataLavg[index.complexL,] <- data.complex
605. }
606. return(dataLavg)
607. }
610. #' Modeling the effect of coreceptor on the ligand-receptor interaction
611. #'
612. #' @param data.use data matrix
613. #' @param cofactor_input the cofactor_input from CellChatDB
614. #' @param pairLRsig a data frame giving ligand-receptor interactions
615. #' @param type when type == "A", computing expression of co-activation receptor; when type == "I", computing expression of co-inhibition receptor.
616. #' @return
617. #' @importFrom future nbrOfWorkers
618. #' @importFrom future.apply future_sapply
619. #' @importFrom pbapply pbsapply
620. #' @export
621. computeExpr_coreceptor <- function(cofactor_input, data.use, pairLRsig, type = c("A", "I")) {
622. type <- match.arg(type)
623. if (type == "A") {
624. coreceptor.all = pairLRsig$co_A_receptor
625. } else if (type == "I"){
626. coreceptor.all = pairLRsig$co_I_receptor
627. }
628. index.coreceptor <- which(!is.na(coreceptor.all) & coreceptor.all != "")
629. if (length(index.coreceptor) > 0) {
630. my.sapply <- ifelse(
631. test = future::nbrOfWorkers() == 1,
632. yes = sapply,
633. no = future.apply::future_sapply
634. )
635. coreceptor <- coreceptor.all[index.coreceptor]
636. coreceptor.ind <- cofactor_input[coreceptor, grepl("cofactor" , colnames(cofactor_input) )]
637. data.coreceptor.ind = my.sapply(
638. X = 1:nrow(coreceptor.ind),
639. FUN = function(x) {
640. coreceptor.indV <- unlist(coreceptor.ind[x,], use.names = F)
641. coreceptor.indV <- coreceptor.indV[coreceptor.indV != ""]
642. coreceptor.indV <- intersect(coreceptor.indV, rownames(data.use))
643. if (length(coreceptor.indV) == 1) {
644. return(1 + data.use[coreceptor.indV, ])
645. } else if (length(coreceptor.indV) > 1) {
646. return(apply(1 + data.use[coreceptor.indV, ], 2, prod))
647. } else {
648. return(matrix(1, nrow = 1, ncol = ncol(data.use)))
649. }
650. }
651. )
652. data.coreceptor.ind <- t(data.coreceptor.ind)
653. data.coreceptor <- matrix(1, nrow = length(coreceptor.all), ncol = ncol(data.use))
654. data.coreceptor[index.coreceptor,] <- data.coreceptor.ind
655. } else {
656. data.coreceptor <- matrix(1, nrow = length(coreceptor.all), ncol = ncol(data.use))
657. }
658. return(data.coreceptor)
659. }
661. # Modeling the effect of coreceptor on the ligand-receptor interaction
662. #
663. # @param data.use data matrix
664. # @param cofactor_input the cofactor_input from CellChatDB
665. # @param pairLRsig a data frame giving ligand-receptor interactions
666. # @param type when type == "A", computing expression of co-activation receptor; when type == "I", computing expression of co-inhibition receptor.
667. # @param group a factor defining the cell groups
668. # @param FunMean the function for computing mean expression per group
669. # @return
670. # @importFrom future nbrOfWorkers
671. # @importFrom future.apply future_sapply
672. # @importFrom pbapply pbsapply
673. # #' @export
674. .computeExprGroup_coreceptor <- function(cofactor_input, data.use, pairLRsig, type = c("A", "I"), group, FunMean) {
675. type <- match.arg(type)
676. if (type == "A") {
677. coreceptor.all = pairLRsig$co_A_receptor
678. } else if (type == "I"){
679. coreceptor.all = pairLRsig$co_I_receptor
680. }
681. index.coreceptor <- which(!is.na(coreceptor.all) & coreceptor.all != "")
682. if (length(index.coreceptor) > 0) {
683. my.sapply <- ifelse(
684. test = future::nbrOfWorkers() == 1,
685. yes = pbapply::pbsapply,
686. no = future.apply::future_sapply
687. )
688. coreceptor <- coreceptor.all[index.coreceptor]
689. coreceptor.ind <- cofactor_input[coreceptor, grepl("cofactor" , colnames(cofactor_input) )]
690. data.coreceptor.ind = my.sapply(
691. X = 1:nrow(coreceptor.ind),
692. FUN = function(x) {
693. coreceptor.indV <- unlist(coreceptor.ind[x,], use.names = F)
694. coreceptor.indV <- coreceptor.indV[coreceptor.indV != ""]
695. coreceptor.indV <- intersect(coreceptor.indV, rownames(data.use))
696. if (length(coreceptor.indV) > 1) {
697. data.avg <- aggregate(t(data.use[coreceptor.indV,]), list(group), FUN = FunMean)
698. data.avg <- t(data.avg[,-1])
699. return(apply(1 + data.avg, 2, prod))
700. # return(1 + apply(data.avg, 2, mean))
701. } else if (length(coreceptor.indV) == 1) {
702. data.avg <- aggregate(matrix(data.use[coreceptor.indV,], ncol = 1), list(group), FUN = FunMean)
703. data.avg <- t(data.avg[,-1])
704. return(1 + data.avg)
705. } else {
706. return(matrix(1, nrow = 1, ncol = length(unique(group))))
707. }
708. }
709. )
710. data.coreceptor.ind <- t(data.coreceptor.ind)
711. data.coreceptor <- matrix(1, nrow = length(coreceptor.all), ncol = length(unique(group)))
712. data.coreceptor[index.coreceptor,] <- data.coreceptor.ind
713. } else {
714. data.coreceptor <- matrix(1, nrow = length(coreceptor.all), ncol = length(unique(group)))
715. }
717. return(data.coreceptor)
718. }
720. #' Modeling the effect of agonist on the ligand-receptor interaction
721. #' @param data.use data matrix
722. #' @param cofactor_input the cofactor_input from CellChatDB
723. #' @param pairLRsig the L-R interactions
724. #' @param group a factor defining the cell groups
725. #' @param index.agonist the index of agonist in the database
726. #' @param Kh a parameter in Hill function
727. #' @param FunMean the function for computing mean expression per group
728. #' @param n Hill coefficient
729. #' @return
730. #' @export
731. #' @importFrom stats aggregate
732. computeExprGroup_agonist <- function(data.use, pairLRsig, cofactor_input, group, index.agonist, Kh, FunMean, n) {
733. agonist <- pairLRsig$agonist[index.agonist]
734. agonist.ind <- cofactor_input[agonist, grepl("cofactor" , colnames(cofactor_input))]
735. agonist.indV <- unlist(agonist.ind, use.names = F)
736. agonist.indV <- agonist.indV[agonist.indV != ""]
737. agonist.indV <- intersect(agonist.indV, rownames(data.use))
738. if (length(agonist.indV) == 1) {
739. data.avg <- aggregate(matrix(data.use[agonist.indV,], ncol = 1), list(group), FUN = FunMean)
740. data.avg <- t(data.avg[,-1])
741. data.agonist <- 1 + data.avg^n/(Kh^n + data.avg^n)
742. } else if (length(agonist.indV) > 1) {
743. data.avg <- aggregate(t(data.use[agonist.indV,]), list(group), FUN = FunMean)
744. data.avg <- t(data.avg[,-1])
745. data.agonist <- apply(1 + data.avg^n/(Kh^n + data.avg^n), 2, prod)
746. } else {
747. data.agonist = matrix(1, nrow = 1, ncol = length(unique(group)))
748. }
749. return(data.agonist)
750. }
752. #' Modeling the effect of antagonist on the ligand-receptor interaction
753. #'
754. #' @param data.use data matrix
755. #' @param cofactor_input the cofactor_input from CellChatDB
756. #' @param pairLRsig the L-R interactions
757. #' @param group a factor defining the cell groups
758. #' @param index.antagonist the index of antagonist in the database
759. #' @param Kh a parameter in Hill function
760. #' @param n Hill coefficient
761. #' @param FunMean the function for computing mean expression per group
762. #' @return
763. #' @export
764. #' @importFrom stats aggregate
765. computeExprGroup_antagonist <- function(data.use, pairLRsig, cofactor_input, group, index.antagonist, Kh, FunMean, n) {
766. antagonist <- pairLRsig$antagonist[index.antagonist]
767. antagonist.ind <- cofactor_input[antagonist, grepl( "cofactor" , colnames(cofactor_input) )]
768. antagonist.indV <- unlist(antagonist.ind, use.names = F)
769. antagonist.indV <- antagonist.indV[antagonist.indV != ""]
770. antagonist.indV <- intersect(antagonist.indV, rownames(data.use))
771. if (length(antagonist.indV) == 1) {
772. data.avg <- aggregate(matrix(data.use[antagonist.indV,], ncol = 1), list(group), FUN = FunMean)
773. data.avg <- t(data.avg[,-1])
774. data.antagonist <- Kh^n/(Kh^n + data.avg^n)
775. } else if (length(antagonist.indV) > 1) {
776. data.avg <- aggregate(t(data.use[antagonist.indV,]), list(group), FUN = FunMean)
777. data.avg <- t(data.avg[,-1])
778. data.antagonist <- apply(Kh^n/(Kh^n + data.avg^n), 2, prod)
779. } else {
780. data.antagonist = matrix(1, nrow = 1, ncol = length(unique(group)))
781. }
782. return(data.antagonist)
783. }
786. #' Modeling the effect of agonist on the ligand-receptor interaction
787. #' @param data.use data matrix
788. #' @param cofactor_input the cofactor_input from CellChatDB
789. #' @param pairLRsig the L-R interactions
790. # #' @param group a factor defining the cell groups
791. #' @param index.agonist the index of agonist in the database
792. #' @param Kh a parameter in Hill function
793. # #' @param FunMean the function for computing mean expression per group
794. #' @param n Hill coefficient
795. #' @return
796. #' @export
797. #' @importFrom stats aggregate
798. computeExpr_agonist <- function(data.use, pairLRsig, cofactor_input, index.agonist, Kh, n) {
799. agonist <- pairLRsig$agonist[index.agonist]
800. agonist.ind <- cofactor_input[agonist, grepl("cofactor" , colnames(cofactor_input))]
801. agonist.indV <- unlist(agonist.ind, use.names = F)
802. agonist.indV <- agonist.indV[agonist.indV != ""]
803. agonist.indV <- intersect(agonist.indV, rownames(data.use))
804. if (length(agonist.indV) == 1) {
805. # data.avg <- aggregate(matrix(data.use[agonist.indV,], ncol = 1), list(group), FUN = FunMean)
806. # data.avg <- t(data.avg[,-1])
807. data.avg <- data.use[agonist.indV,, drop = FALSE]
808. data.agonist <- 1 + data.avg^n/(Kh^n + data.avg^n)
809. } else if (length(agonist.indV) > 1) {
810. # data.avg <- aggregate(t(data.use[agonist.indV,]), list(group), FUN = FunMean)
811. # data.avg <- t(data.avg[,-1])
812. data.avg <- data.use[agonist.indV,, drop = FALSE]
813. data.agonist <- apply(1 + data.avg^n/(Kh^n + data.avg^n), 2, prod)
814. } else {
815. # data.agonist = matrix(1, nrow = 1, ncol = length(unique(group)))
816. data.agonist = matrix(1, nrow = 1, ncol = ncol(data.use))
817. }
818. return(data.agonist)
819. }
821. #' Modeling the effect of antagonist on the ligand-receptor interaction
822. #'
823. #' @param data.use data matrix
824. #' @param cofactor_input the cofactor_input from CellChatDB
825. #' @param pairLRsig the L-R interactions
826. # #' @param group a factor defining the cell groups
827. #' @param index.antagonist the index of antagonist in the database
828. #' @param Kh a parameter in Hill function
829. #' @param n Hill coefficient
830. # #' @param FunMean the function for computing mean expression per group
831. #' @return
832. #' @export
833. #' @importFrom stats aggregate
834. computeExpr_antagonist <- function(data.use, pairLRsig, cofactor_input, index.antagonist, Kh, n) {
835. antagonist <- pairLRsig$antagonist[index.antagonist]
836. antagonist.ind <- cofactor_input[antagonist, grepl( "cofactor" , colnames(cofactor_input) )]
837. antagonist.indV <- unlist(antagonist.ind, use.names = F)
838. antagonist.indV <- antagonist.indV[antagonist.indV != ""]
839. antagonist.indV <- intersect(antagonist.indV, rownames(data.use))
840. if (length(antagonist.indV) == 1) {
841. # data.avg <- aggregate(matrix(data.use[antagonist.indV,], ncol = 1), list(group), FUN = FunMean)
842. # data.avg <- t(data.avg[,-1])
843. data.avg <- data.use[antagonist.indV,, drop = FALSE]
844. data.antagonist <- Kh^n/(Kh^n + data.avg^n)
845. } else if (length(antagonist.indV) > 1) {
846. # data.avg <- aggregate(t(data.use[antagonist.indV,]), list(group), FUN = FunMean)
847. # data.avg <- t(data.avg[,-1])
848. data.avg <- data.use[antagonist.indV,, drop = FALSE]
849. data.antagonist <- apply(Kh^n/(Kh^n + data.avg^n), 2, prod)
850. } else {
851. # data.antagonist = matrix(1, nrow = 1, ncol = length(unique(group)))
852. data.antagonist = matrix(1, nrow = 1, ncol = ncol(data.use))
853. }
854. return(data.antagonist)
855. }
858. #' Compute the geometric mean
859. #' @param x a numeric vector
860. #' @param na.rm whether remove na
861. #' @return
862. #' @export
863. geometricMean <- function(x,na.rm=TRUE){
864. if (is.null(nrow(x))) {
865. exp(mean(log(x),na.rm=na.rm))
866. } else {
867. exp(apply(log(x),2,mean,na.rm=na.rm))
868. }
869. }
872. #' Compute the Tukey's trimean
873. #' @param x a numeric vector
874. #' @param na.rm whether remove na
875. #' @return
876. #' @importFrom collapse fquantile
877. #' @export
878. triMean <- function(x, na.rm = TRUE) {
879. mean(collapse::fquantile(x, probs = c(0.25, 0.50, 0.50, 0.75), na.rm = na.rm))
880. }
882. #' Compute the average expression per cell group when the percent of expressing cells per cell group larger than a threshold
883. #' @param x a numeric vector
884. #' @param trim the percent of expressing cells per cell group to be considered as zero
885. #' @param na.rm whether remove na
886. #' @return
887. #' @importFrom Matrix nnzero
888. # #' @export
889. thresholdedMean <- function(x, trim = 0.1, na.rm = TRUE) {
890. percent <- Matrix::nnzero(x)/length(x)
891. if (percent < trim) {
892. return(0)
893. } else {
894. return(mean(x, na.rm = na.rm))
895. }
896. }
898. #' Filter cell-cell communication if there are only few number of cells in certain cell groups or inconsistent cell-cell communication across samples
899. #'
900. #' @param object CellChat object
901. #' @param min.cells The minmum number of cells required in each cell group for cell-cell communication
902. #' @param min.samples The minmum number of samples required for consistent cell-cell communication across samples (that is an interaction present in at least `min.samples` samples) when mutiple samples/replicates/batches are merged as an input for CellChat analysis.
903. #' @param rare.keep Whether to keep the interactions associated with the rare populations when min.samples >= 2. When a rare population is identified in the merged samples (say 15 cells in this rare population from two samples), it is likely to filter out the interactions associated with this rare population when setting min.samples >= 2. Setting `rare.keep = TRUE` to retain the identified interactions associated with this rare population.
904. #' @param nonFilter.keep Whether to keep the non-filtered cell-cell communication in the CellChat object. This is useful for avoiding re-running `computeCommunProb` if you want to adjust the parameters when running `filterCommunication`.
905. #' @return CellChat object with an updated slot net
906. #' @export
907. #'
908. filterCommunication <- function(object, min.cells = 10, min.samples = NULL, rare.keep = FALSE, nonFilter.keep = FALSE) {
909. net <- object@net
910. if (nonFilter.keep == TRUE) {
911. cat("The non-filtered cell-cell communication is stored in `object@net$prob.nonFilter` and `object@net$pval.nonFilter`. \n")
912. object@net$prob.nonFilter <- net$prob
913. object@net$pval.nonFilter <- net$pval
914. }
915. num.interaction0 <- sum(net$prob > 0)
916. cell.excludes <- which(as.numeric(table(object@idents)) <= min.cells)
917. if (length(cell.excludes) > 0) {
918. cat("The cell-cell communication related with the following cell groups are excluded due to the few number of cells: ", toString(levels(object@idents)[cell.excludes]), "!",'\t')
919. net$prob[cell.excludes,,] <- 0
920. net$prob[,cell.excludes,] <- 0
921. num.interaction1 <- sum(net$prob > 0)
922. pct.dicrease <- scales::percent((num.interaction0-num.interaction1)/num.interaction0, accuracy = .1)
923. cat(paste0(pct.dicrease, " interactions are removed!",'\n'))
924. } else {
925. num.interaction1 <- num.interaction0
926. }
928. sample.info <- object@meta$samples
929. sample.id <- levels(sample.info)
930. if (is.null(min.samples)) {
931. min.samples <- 1
932. } else if (min.samples > length(sample.id)) {
933. stop(paste0("There are only ", length(sample.id), " samples in the data. Please change the value of `min.samples`! "))
934. }
935. if (length(sample.id) >= 2 & min.samples >= 2) {
936. if (object@options$parameter$raw.use == TRUE) {
937. data <- as.matrix([email hidden])
938. } else {
939. if ("data.smooth" %in% methods::slotNames(object) == FALSE) {
940. stop("`[email hidden]` is missing. Please update the CellChat object via `updateCellChat`! \n")
941. }
942. data <- as.matrix([email hidden])
943. }
944. data.use <- data/max(data)
945. group <- object@idents
946. type <- object@options$parameter$type.mean
947. trim <- object@options$parameter$trim
948. FunMean <- switch(type,
949. triMean = triMean,
950. truncatedMean = function(x) mean(x, trim = trim, na.rm = TRUE),
951. thresholdedMean = function(x) thresholdedMean(x, trim = trim, na.rm = TRUE),
952. median = function(x) median(x, na.rm = TRUE))
954. LR <- dimnames(net$prob)[[3]]
955. idx.nonzero <- which(apply(net$prob, 3, sum) != 0)
956. LR.nonzero <- LR[idx.nonzero] # only examine the L-R pairs with nonzero communication probabilities.
958. interaction_input <- object@DB$interaction
959. complex_input <- object@DB$complex
960. geneIfo <- object@DB$geneInfo
961. idx <- match(LR.nonzero, interaction_input$interaction_name)
962. geneL <- as.character(interaction_input$ligand[idx])
963. geneR <- as.character(interaction_input$receptor[idx])
965. geneLR <- c(unique(geneL), unique(geneR))
966. geneLR <- extractGeneSubset(geneLR, complex_input, geneIfo)
967. data.use <- data.use[rownames(data.use) %in% geneLR, ]
969. score.LR <- array(0, dim = c(nlevels(group),nlevels(group),length(LR.nonzero), length(sample.id)))
970. LR.nonzero.all <- c()
971. cell.excludes.sample <- c()
972. for (i in 1:length(sample.id)) {
973. cell.use <- which(sample.info == sample.id[i])
974. group.use <- group[cell.use]
975. group.use <- droplevels(group.use)
976. # get the rare populations with few cells in each sample
977. cell.excludes.sample.i <- which(as.numeric(table(object@idents[cell.use])) <= min.cells)
978. cell.excludes.sample <- c(cell.excludes.sample, cell.excludes.sample.i)
979. # compute average expression per cell group
980. data.use.i <- data.use[, cell.use]
981. data.use.avg <- aggregate(t(data.use.i), list(group.use), FUN = FunMean)
982. data.use.avg <- t(data.use.avg[,-1])
983. group.exist <- which(levels(group) %in% unique(group.use))
984. if (length(group.exist) < nlevels(group)) {
985. data.use.avg.temp <- matrix(0, nrow = nrow(data.use), ncol = nlevels(group))
986. data.use.avg.temp[ , group.exist] <- data.use.avg
987. rownames(data.use.avg.temp) <- rownames(data.use.avg)
988. data.use.avg <- data.use.avg.temp
989. }
990. colnames(data.use.avg) <- levels(group)
991. # compute the average expression of ligand or receptor in each cell group
992. dataLavg <- computeExpr_LR(geneL, data.use.avg, complex_input)
993. dataRavg <- computeExpr_LR(geneR, data.use.avg, complex_input)
994. # compute the interaction scores for each ligand-receptor pair based on their expression
995. for (jj in 1:length(LR.nonzero)) { # It is not good to use parallel here because it will change the order of LR
996. score.LR[,,jj,i] <- Matrix::crossprod(matrix(dataLavg[jj, ], nrow = 1), matrix(dataRavg[jj, ], nrow = 1))
997. }
998. if (length(cell.excludes.sample.i) > 0) {
999. cat(paste0("The number of cells of the following cell groups in ", sample.id[i], " sample are less than ", min.cells, " cells: ",toString(levels(object@idents)[cell.excludes.sample.i]), "!",'\n'))
1000. score.LR[cell.excludes.sample.i, , , i] <- 0
1001. score.LR[ ,cell.excludes.sample.i, , i] <- 0
1002. }
1003. #LR.nonzero.all <- c(LR.nonzero.all, LR.nonzero[apply(score.LR[ , , , i], 3, sum) != 0])
1004. }
1005. #LR.nonzero.jointOnly <- setdiff(LR.nonzero, unique(LR.nonzero.all))
1007. # get the excluded cell groups that are not observed in the merged data, which is very possible for rare populations
1008. cell.excludes.sample <- unique(cell.excludes.sample)
1009. if (length(cell.excludes.sample) > 0) {
1010. cell.excludes.sample <- setdiff(cell.excludes.sample, cell.excludes)
1011. }
1013. score.LR[score.LR > 0] <- 1 # binarize the interaction score
1014. score.LR.consitent <- array(0, dim = c(nlevels(group),nlevels(group),length(LR.nonzero)))
1015. LR.inconsitent <- c()
1016. for (jj in 1:length(LR.nonzero)) {
1017. score.LR.sum <- apply(score.LR[ , , jj, ], c(1,2), sum) # elements 2 and 1 means consistent and inconsistent interactions across samples, respectively.
1018. # set communication probability to be zero for inconsistent interactions across samples
1019. if (sum((score.LR.sum > 0) * (score.LR.sum < min.samples)) > 0) {
1020. #LR.inconsitent <- c(LR.inconsitent, LR.nonzero[jj])
1021. score.LR.consitent <- (score.LR.sum >= min.samples) * 1
1022. if (rare.keep == TRUE & length(cell.excludes.sample) > 0) {
1023. score.LR.consitent[cell.excludes.sample, ] <- 1
1024. score.LR.consitent[ ,cell.excludes.sample] <- 1
1025. }
1026. net$prob[ , , LR.nonzero[jj]] <- net$prob[ , , LR.nonzero[jj]] * score.LR.consitent
1027. }
1028. }
1029. num.interaction2 <- sum(net$prob > 0)
1030. pct.dicrease <- scales::percent((num.interaction1-num.interaction2)/num.interaction1, accuracy = .1)
1031. cat(paste0(pct.dicrease, " interactions are removed due to their inconsistence across ", min.samples, " samples!",'\n'))
1032. }
1034. object@net <- net
1035. return(object)
1036. }
1039. #' Identify all the significant interactions (L-R pairs) from some cell groups to other cell groups
1040. #'
1041. #' @param object CellChat object
1042. #' @param from a vector giving the index or the name of source cell groups
1043. #' @param to a corresponding vector giving the index or the name of target cell groups. Note: The length of 'from' and 'to' must be the same, giving the corresponding pair of cell groups for communication.
1044. #' @param bidirection whether show the bidirectional communication, i.e., both 'from'->'to' and 'to'->'from'.
1045. #' @param pair.only whether only return ligand-receptor pairs without pathway names and communication strength
1046. #' @param pairLR.use0 ligand-receptor pairs to use; default is all the significant interactions
1047. #' @param thresh threshold of the p-value for determining significant interaction
1048. #'
1049. #' @return
1050. #' @export
1051. #'
1052. identifyEnrichedInteractions <- function(object, from, to, bidirection = FALSE, pair.only = TRUE, pairLR.use0 = NULL, thresh = 0.05){
1053. pairwiseLR <- object@net$pairwiseRank
1054. if (is.null(pairwiseLR)) {
1055. stop("The interactions between pairwise cell groups have not been extracted!
1056. Please first run `object <- rankNetPairwise(object)`")
1057. }
1058. group.names.all <- names(pairwiseLR)
1059. if (!is.numeric(from)) {
1060. from <- match(from, group.names.all)
1061. if (sum(is.na(from)) > 0) {
1062. message("Some input cell group names in 'from' do not exist!")
1063. from <- from[!is.na(from)]
1064. }
1065. }
1066. if (!is.numeric(to)) {
1067. to <- match(to, group.names.all)
1068. if (sum(is.na(to)) > 0) {
1069. message("Some input cell group names in 'to' do not exist!")
1070. to <- to[!is.na(to)]
1071. }
1072. }
1073. if (length(from) != length(to)) {
1074. stop("The length of 'from' and 'to' must be the same!")
1075. }
1076. if (bidirection) {
1077. from2 <- c(from, to)
1078. to <- c(to, from)
1079. from <- from2
1080. }
1081. if (is.null(pairLR.use0)) {
1082. k <- 0
1083. pairLR.use0 <- list()
1084. for (i in 1:length(from)){
1085. pairwiseLR_ij <- pairwiseLR[[from[i]]][[to[i]]]
1086. idx <- pairwiseLR_ij$pval < thresh
1087. if (length(idx) > 0) {
1088. k <- k +1
1089. pairLR.use0[[k]] <- pairwiseLR_ij[idx,]
1090. }
1091. }
1092. pairLR.use0 <- do.call(rbind, pairLR.use0)
1093. }
1095. k <- 0
1096. pval <- matrix(nrow = length(rownames(pairLR.use0)), ncol = length(from))
1097. prob <- pval
1098. group.names <- c()
1099. for (i in 1:length(from)) {
1100. k <- k+1
1101. pairwiseLR_ij <- pairwiseLR[[from[i]]][[to[i]]]
1102. pairwiseLR_ij <- pairwiseLR_ij[rownames(pairLR.use0),]
1103. pval_ij <- pairwiseLR_ij$pval
1104. prob_ij <- pairwiseLR_ij$prob
1105. pval_ij[pval_ij > 0.05] = 1
1106. pval_ij[pval_ij > 0.01 & pval_ij <= 0.05] = 2
1107. pval_ij[pval_ij <= 0.01] = 3
1108. prob_ij[pval_ij ==1] <- 0
1109. pval[,k] <- pval_ij
1110. prob[,k] <- prob_ij
1111. group.names <- c(group.names, paste(group.names.all[from[i]], group.names.all[to[i]], sep = " - "))
1112. }
1113. prob[which(prob == 0)] <- NA
1114. # remove rows that are entirely NA
1115. pval <- pval[rowSums(is.na(prob)) != ncol(prob), ,drop = FALSE]
1116. pairLR.use0 <- pairLR.use0[rowSums(is.na(prob)) != ncol(prob), ,drop = FALSE]
1117. prob <- prob[rowSums(is.na(prob)) != ncol(prob), ,drop = FALSE]
1118. if (pair.only) {
1119. pairLR.use0 <- dplyr::select(pairLR.use0, ligand, receptor)
1120. }
1121. return(pairLR.use0)
1122. }
1125. #' Compute the region distance based on the spatial locations of each splot/cell of the spatial transcriptomics
1126. #'
1127. #' @param coordinates a data matrix in which each row gives the spatial locations/coordinates of each cell/spot
1128. #' @param meta a data frame including at least two columns named `group` and `samples`. `meta$group` is a factor vector defining the regions/labels of each cell/spot. `meta$samples` is a factor vector defining the sample labels of each dataset.
1129. #' @param interaction.range The maximum interaction/diffusion range of ligands. This hard threshold is used to filter out the connections between spatially distant regions
1130. #' @param ratio a numerical vector giving the conversion factor when converting spatial coordinates from Pixels or other units to Micrometers (i.e.,Microns).
1131. #'
1132. #' For example, setting `ratio = 0.18` indicates that 1 pixel equals 0.18um in the coordinates.
1133. #' For 10X visium, it is the ratio of the theoretical spot size (i.e., 65um) over the number of pixels that span the diameter of a theoretical spot size in the full-resolution image (i.e., 'spot.size.fullres' in the 'scalefactors_json.json' file).
1134. #' @param tol a numerical vector giving the tolerance factor to increase the robustness when comparing the center-to-center distance against the `interaction.range`. This can be the half value of cell/spot size in the unit of um.
1135. #'
1136. #' For example, for 10X visium, `tol` can be set as `65/2`; for slide-seq, `tol` can be set as `10/2`.
1137. #' If the cell/spot size is not known, we provide a function `computeCellDistance` to compute the center-to-center distance. `tol` can be the the half value of the minimum center-to-center distance.
1138. #' @param k.min the minimum number of interacting cell pairs required for defining adjacent cell groups
1139. #' @param contact.dependent Whether determining spatially proximal cell groups based on either the contact.range or the k-nearest neighbors (knn). By default `contact.dependent = TRUE` when inferring contact-dependent and juxtacrine signaling (including ECM-Receptor and Cell-Cell Contact signaling classified in CellChatDB$interaction$annotation).
1140. #' If only focusing on `Secreted Signaling`, the `contact.dependent` will be automatically set as FALSE except for `contact.dependent.forced = TRUE`.
1141. #' @param contact.range The interaction range (Unit: microns) to restrict the contact-dependent signaling.
1142. #' For spatial transcriptomics in a single-cell resolution, `contact.range` is approximately equal to the estimated cell diameter (i.e., the cell center-to-center distance), which means that contact-dependent and juxtacrine signaling can only happens when the two cells are contact to each other.
1143. #'
1144. #' Typically, `contact.range = 10`, which is a typical human cell size. However, for low-resolution spatial data such as 10X visium, it should be the cell center-to-center distance (i.e., `contact.range = 100` for visium data). The function `computeCellDistance` can compute the center-to-center distance.
1145. #'
1146. #' @param contact.knn.k Number of neighbors to restrict the contact-dependent signaling within the neatest neighbors. By default, CellChat uses `contact.range` to restrict the contact-dependent signaling; however, users can also provide a value of `contact.knn.k`, in order to determine spatially proximal cell groups based on the k-nearest neighbors (knn).
1147. #' For 10X visium, contact.knn.k = 6. For other spatial technologies, this value may be hard to determine because the sequenced cells/spots are usually not regularly arranged.
1148. #' @param do.symmetric Whether converting the adjacent matrix into symmetric one when determining spatially proximal cell groups. Default is TRUE, indicating that if adj(i,j) or adj(j,i) is zero, then both are zeros.
1149. #'
1150. #' @importFrom BiocNeighbors queryKNN AnnoyParam
1151. #' @return A list including a square matrix giving the pairwise region distances and an adjacent matrix indicating physically contacting cell groups based on either the contact.range or the k-nearest neighbors
1152. #'
1153. #' @export
1154. computeRegionDistance <- function(coordinates, meta,
1155. interaction.range = NULL, ratio = NULL, tol = NULL, k.min = 10,
1156. contact.dependent = TRUE, contact.range = NULL, contact.knn.k = NULL, do.symmetric = TRUE
1157. ) {
1158. trim <- 0.1
1159. FunMean <- function(x) mean(x, trim = trim, na.rm = TRUE) # This is used for computing the average distance between two cell groups
1160. group <- meta$group
1161. numCluster <- nlevels(group)
1162. level.use <- levels(group)
1163. level.use <- level.use[level.use %in% unique(group)]
1164. samples <- meta$samples
1165. samples.use <- levels(samples)
1166. d.spatial <- array(NaN, dim = c(numCluster,numCluster,length(samples.use)))
1167. adj.spatial <- array(0, dim = c(numCluster,numCluster,length(samples.use)))
1168. adj.contact <- array(0, dim = c(numCluster,numCluster,length(samples.use)))
1169. adj.contact.knn <- array(0, dim = c(numCluster,numCluster,length(samples.use)))
1171. if (contact.dependent == TRUE & !is.null(contact.knn.k)) {
1172. ## find the k-nearest neighbors for each single cell
1173. # my.knn <- FNN::get.knn(coordinates, k = contact.knn.k)
1174. # nn.ranked <- my.knn$nn.index # this is a matrix with the size of nCell * contact.knn.k
1175. nn.ranked <- matrix(NA, nrow = nrow(coordinates), ncol = contact.knn.k)
1176. for (k in 1:length(samples.use)) {
1177. idx.k <- which(samples == samples.use[k])
1178. my.knn <- suppressWarnings(BiocNeighbors::findKNN(coordinates[idx.k, ], k = contact.knn.k, BNPARAM = BiocNeighbors::AnnoyParam(), get.index = TRUE))
1179. nn.ranked[idx.k, ] <- my.knn$index # this is a matrix with the size of nCell * contact.knn.k
1180. }
1181. k.min.contact <- k.min
1182. } else {
1183. nn.ranked <- matrix(1, nrow = nrow(coordinates), ncol = 1)
1184. k.min.contact <- -1 # this produces adj.contact.knn with all elements being 1
1185. }
1186. if (contact.dependent == TRUE) {
1187. if (is.null(contact.range) & is.null(contact.knn.k)) {
1188. stop("Please check the documentation of `computeCommunProb` and provide the value of either `contact.range` or `contact.knn.k`")
1189. }
1190. } else {
1191. contact.range <- 10000 # this produces adj.contact with all elements being 1
1192. }
1194. for (k in 1:length(samples.use)) {
1195. idx.k <- samples == samples.use[k]
1196. for (i in 1:numCluster) {
1197. for (j in 1:numCluster) {
1198. idx.i <- which((group == level.use[i]) & idx.k)
1199. idx.j <- which((group == level.use[j]) & idx.k)
1200. if (length(idx.i) == 0 | length(idx.j) == 0) {
1201. next # if one cell group is missing in one sample, just goes to next loop
1202. }
1203. data.spatial.i <- coordinates[idx.i, , drop = FALSE]
1204. data.spatial.j <- coordinates[idx.j, , drop = FALSE]
1205. # for each point in the i-th cell group, find its 1-nearest neighbor in the j-th cell group
1206. #qout <- suppressWarnings(BiocNeighbors::queryKNN(data.spatial.j, data.spatial.i, k = 1, BNPARAM = BiocNeighbors::KmknnParam(), get.index = TRUE))
1207. qout <- suppressWarnings(BiocNeighbors::queryKNN(data.spatial.j, data.spatial.i, k = 1, BNPARAM = BiocNeighbors::AnnoyParam(), get.index = TRUE))
1208. # qout$index is an one column matrix with length being `length(idx.i)`, which is the index of the 1-nearest neighbor in the j-th cell group defined by `idx.j`
1209. # qout$distance is an one column matrix with length being `length(idx.i)`, which is the distance to the 1-nearest neighbor in the j-th cell group defined by `idx.j`
1211. # conver the calculated distance into the distance in micrometers
1212. qout$distance <- qout$distance*ratio[k]
1213. # long-range distance
1214. idx <- qout$distance - interaction.range < tol[k]
1215. adj.spatial[i,j,k] <- (length(unique(qout$index[idx])) >= k.min) * 1
1216. # short-range distance based on contact.range
1217. idx2 <- qout$distance - contact.range < tol[k]
1218. adj.contact[i,j,k] <- (length(unique(qout$index[idx2])) >= k.min) * 1
1219. # short-range distance based on knn
1220. knn.i <- unique(as.vector(nn.ranked[idx.i, ]))
1221. #adj.contact.knn[i,j,k] <- (length(intersect(knn.i, idx.j)) >= k.min.contact) * 1
1222. adj.contact.knn[i,j,k] <- (length(intersect(knn.i, unique(qout$index[idx]))) >= k.min.contact) * 1 # knn within the long-range distance
1223. # computing the average distance between two cell groups
1224. d.spatial[i,j,k] <- FunMean(qout$distance) # since distances are positive values, different ways for computing the mean have little effects.
1226. }
1227. }
1228. }
1230. # merged spatial information from different samples
1231. d.spatial <- apply(d.spatial, c(1,2), function(x) mean(x, na.rm = TRUE))
1232. adj.spatial <- apply(adj.spatial, c(1,2), mean)
1233. adj.contact <- apply(adj.contact, c(1,2), mean)
1234. adj.contact.knn <- apply(adj.contact.knn, c(1,2), mean)
1235. # for multi-samples analysis, the following is needed
1236. adj.spatial[adj.spatial > 0] <- 1
1237. adj.contact[adj.contact > 0] <- 1
1238. adj.contact.knn[adj.contact.knn > 0] <- 1
1240. # make these adjacent matrix as symmetric
1241. if (do.symmetric) {
1242. adj.spatial <- adj.spatial * t(adj.spatial) # if one is zero, then both are zeros.
1243. adj.contact <- adj.contact * t(adj.contact) # if one is zero, then both are zeros.
1244. adj.contact.knn <- adj.contact.knn * t(adj.contact.knn) # if one is zero, then both are zeros.
1245. }
1246. d.spatial <- (d.spatial + t(d.spatial))/2
1248. # filter out the spatially distant cell groups
1249. adj.spatial[adj.spatial == 0] <- NaN
1250. d.spatial <- d.spatial * adj.spatial
1252. rownames(d.spatial) <- levels(group); colnames(d.spatial) <- levels(group)
1254. if (length(contact.knn.k) > 0) {
1255. adj.contact = adj.contact.knn
1256. }
1257. res <- list(d.spatial = d.spatial, adj.contact = adj.contact)
1258. return(res)
1260. }
1262. #' Compute cell-cell distance based on the spatial coordinates
1263. #'
1264. #' @param coordinates a data matrix in which each row gives the spatial locations/coordinates of each cell/spot
1265. #' @param interaction.range The maximum interaction/diffusion range of ligands. This hard threshold is used to filter out the connections between spatially distant cells
1266. #' @param ratio The conversion factor when converting spatial coordinates from Pixels or other units to Micrometers (i.e.,Microns).
1267. #'
1268. #' For example, setting `ratio = 0.18` indicates that 1 pixel equals 0.18um in the coordinates.
1269. #' For 10X visium, it is the ratio of the theoretical spot size (i.e., 65um) over the number of pixels that span the diameter of a theoretical spot size in the full-resolution image (i.e., 'spot.size.fullres' in the 'scalefactors_json.json' file).
1270. #' @param tol The tolerance factor to increase the robustness when comparing the center-to-center distance against the `interaction.range`. This can be the half value of cell/spot size in the unit of um.
1271. #'
1272. #' For example, for 10X visium, `tol` can be set as `65/2`; for slide-seq, `tol` can be set as `10/2`.
1273. #' If the cell/spot size is not known, we provide a function `computeCellDistance` to compute the center-to-center distance. `tol` can be the the half value of the minimum center-to-center distance.
1274. #'
1275. #' @return an object of class "dist" giving the pairwise cell-cell distance
1276. #' @export
1277. #'
1278. computeCellDistance <- function(coordinates, interaction.range = NULL, ratio = NULL, tol = NULL){
1279. if (ncol(coordinates) == 2) {
1280. colnames(coordinates) <- c("x_cent","y_cent")
1281. } else {
1282. stop("Please check the input 'coordinates' and make sure it is a two column matrix.")
1283. }
1285. d.spatial <- collapse::fdist(coordinates)
1286. if (!is.null(ratio)) {
1287. d.spatial <- d.spatial*ratio
1288. }
1290. if(!is.null(interaction.range) & !is.null(tol)){
1291. message("\n Apply a predefined spatial distance threshold based on the interaction length...")
1292. d.spatial[d.spatial > (interaction.range + tol)] <- NaN
1293. }
1294. return(d.spatial)
1295. }

modeling.R at commit 75253cd, under GPL-3.0 · at the source