Recombinant dimeric PICK1 peptide inhibitors for long-term relief of chronic pain by AAV therapeutics.

1. #' @import ggrepel
2. #' @import gridExtra
3. #' @import leidenAlg
4. #' @import R6
5. #' @import Rtsne
6. NULL
8. ## exporting to inherit parameters below, leiden.community
9. #' @export
10. leidenAlg::leiden.community
13. #' @title Conos R6 class
14. #' @description The class encompasses sample collections, providing methods for calculating and visualizing joint graph and communities.
15. #' @import methods
16. #' @param x a named list of pagoda2 or Seurat objects (one per sample)
17. #' @param n.cores numeric Number of cores to use (default=parallel::detectCores(logical=FALSE))
18. #' @param verbose boolean Whether to provide verbose output (default=TRUE)
19. #' @param clustering string Name of the clustering to use
20. #' @param groups a factor on cells to use for coloring
21. #' @param colors a color factor (named with cell names) use for cell coloring
22. #' @param gene show expression of a gene
23. #' @export Conos
24. Conos <- R6::R6Class("Conos", lock_objects=FALSE,
25. public = list(
26. #' @field samples list of samples (Pagoda2 or Seurat objects)
27. samples = list(),
29. #' @field pairs pairwise alignment results
30. pairs = list(),
32. #' @field graph alignment graph
33. graph = NULL,
35. #' @field clusters list of clustering results named by clustering type
36. clusters = list(),
38. #' @field expression.adj adjusted expression values
39. expression.adj = list(),
41. #' @field embeddings list of joint embeddings
42. embeddings = list(),
44. #' @field embedding joint embedding
45. embedding = NULL,
47. #' @field n.cores number of cores
48. n.cores = 1,
50. #' @field misc list with unstructured additional info
51. misc = list(),
53. #' @field override.conos.plot.theme boolean Whether to override the conos plot theme
54. override.conos.plot.theme = FALSE,
56. #' @description initialize Conos class
57. #'
58. #' @param override.conos.plot.theme boolean Whether to reset plot settings to the ggplot2 default (default=FALSE)
59. #' @param ... additional parameters upon initializing Conos
60. #' @return a new 'Conos' object
61. #' @examples
62. #' con <- Conos$new(small_panel.preprocessed, n.cores=1)
63. #'
64. initialize=function(x, ..., n.cores=parallel::detectCores(logical=FALSE), verbose=TRUE, override.conos.plot.theme=FALSE) {
65. self$n.cores <- n.cores
66. self$override.conos.plot.theme <- override.conos.plot.theme
68. if (missing(x)){
69. return()
70. }
72. if ('Conos' %in% class(x)) { # copy constructor
73. for(n in ls(x)) {
74. if (!is.function(get(n, x))) assign(n, get(n, x), self)
75. }
76. } else {
77. if (!is.list(x)) {
78. stop("x is not a list of pagoda2 or Seurat objects")
79. }
80. if (inherits(x = x[[1]], what = c('Pagoda2', 'seurat', 'Seurat'))) {
81. self$addSamples(x)
82. } else {
83. stop("Only Pagoda2 or Seurat result lists are currently supported")
84. }
85. }
86. },
88. #' @description Initialize or add a set of samples to the conos panel. Note: this will simply add samples, but will not update graph, clustering, etc.
89. #'
90. #' @param replace boolean Whether the existing samples should be purged before adding new ones (default=FALSE)
91. #' @param verbose boolean Whether to provide verbose output (default=FALSE)
92. #' @return invisible view of the full sample list
93. addSamples=function(x, replace=FALSE, verbose=FALSE) {
94. # check names
95. if(is.null(names(x))) {
96. stop("The sample list must be named")
97. }
98. if(replace || length(self$samples)==0) {
99. if(length(x)<2) {
100. stop("The provided list contains less than 2 samples; 2 required, >3 recommended")
101. }
102. }
103. if(any(duplicated(names(self$samples)))) {
104. stop("duplicate names found in the supplied samples")
105. }
106. if(any(names(x) %in% names(self$samples))) {
107. stop("Some of the names in the provided sample list are identical to already existing samples")
108. }
110. # TODO: package-independent wrapper
111. self$samples <- c(self$samples, x)
112. },
114. #' @description Build the joint graph that encompasses all the samples, establishing weighted inter-sample cell-to-cell links
115. #'
116. #' @param k integer integer Size of the inter-sample neighborhood (default=15)
117. #' @param k.self integer Size of the with-sample neighborhoods (default=10).
118. #' @param k.self.weight numeric Weight multiplier on the intra-sample edges relative to inter-sample edges (default=0.1)
119. #' @param alignment.strength numeric Alignment strength (default=NULL will result in alignment.strength=0)
120. #' @param space character Reduced expression space used to establish putative alignments between pairs of samples (default='PCA'). Currently supported spaces are:
121. #' --- "CPCA" Common principal component analysis
122. #' --- "JNMF" Joint NMF
123. #' --- "genes" Gene expression space (log2 transformed)
124. #' --- "PCA" Principal component analysis
125. #' --- "CCA" Canonical correlation analysis
126. #' --- "PMA" (Penalized Multivariate Analysis <https://cran.r-project.org/web/packages/PMA/index.html>)
127. #' @param matching.method character Matching method (default='mNN'). Currently supported methods are "NN" (nearest neighbors) or "mNN" (mututal nearest neighbors).
128. #' @param metric character Distance metric to measure similarity (default='angular'). Currenlty supported metrics are "angular" and "L2".
129. #' @param k1 numeric Neighborhood radius for identifying mutually-matching neighbors (default=k). Note that k1 must be greater than or equal to k, i.e. k1>=k. Increasing k1 beyond k will lead to more aggressive alignment of distinct subpopulations (i.e. increased alignment strengths).
130. #' @param data.type character Type of data type in the input pagoda2 objects within r.n (default='counts').
131. #' @param l2.sigma numeric L2 distances get transformed as exp(-d/sigma) using this value (default=1e5)
132. #' @param var.scale boolean Whether to use common variance scaling (default=TRUE). If TRUE, use geometric means for variance, as we're trying to focus on the common variance components. See scaledMatricesP2() code.
133. #' @param ncomps integer Number of components (default=40)
134. #' @param n.odgenes integer Number of overdispersed genes to be used in each pairwise alignment (default=2000)
135. #' @param matching.mask an optional matrix explicitly specifying which pairs of samples should be compared (a symmetrical matrix of logical values with row and column names corresponding to sample names). (default=NULL). By default, comparisons between all paris are allowed. The argument can be used to exclude comparisons across certain pairs of samples (e.g. techincal replicates, which are expected to show very high similarity).
136. #' @param exclude.samples optional list of sample names that should be excluded from the alignment and the resulting graph (default=NULL)
137. #' @param common.centering boolean When calculating reduced expression space for a given sample pair, whether the expression of genes should be centered using the mean from both samples (TRUE) or using the mean within each sample (FALSE) (default=TRUE)
138. #' @param base.groups an optional factor on cells specifying previously-obtained cell grouping to be used for adjusting the sample alignment (default: NULL). Specifically, cell clusters specfiieid by the base.groups can be used to i) calculate global expression axes which are appended to the overall set of eigenvectors, ii) adding decoy cells.
139. #' @param append.global.axes boolean Whether to project samples on global expression axes, as defined by pre-defined (typically crude) set of cell subpopulations as specified by the base.gruops parameter (default=TRUE, but works only if base.groups is specified)
140. #' @param append.decoys boolean Whether to use pre-defined cell groups (specified by base.groups) to append decoy cells to the samples which are otherwise lacking any of the pre-specified cell groups (default=TRUE, but works only if base.groups is specified). The decoy cells can reduce the number of erroneous matches in highly heterogeneous sample collections, where some of the samples lack entire cell subpopulations which are found in other samples. The approach only works if the base.groups (typically a crude clustering of top-level cell types) can be established with a reasonable confidence.
141. #' @param decoy.threshold integer Minimal number of cells of a given cell type that should exist in a given sample (according to base.groups) to avoid addition of decoy cells to that sample for the purposes of alignment (default=1)
142. #' @param n.decoys integer Number of decoy cells that should be added to a sample that had less than decoy.threshold cells of a given cell type (default=k*2)
143. #' @param score.component.variance boolean Whether to score the amount of total variance explained by different components (default=FALSE as it takes extra time to calculate)
144. #' @param snn boolean Whether to transform the joint graph by computing a shared nearest neighborhood graph (analogous to Seurat 3), further weighting the edges between two matched cells based on the similarity (measured by Jaccard coefficient) of all of their predicted neighbors (across all of the samples) (default: FALSE)
145. #' @param snn.quantile numeric Specifies how the shared neighborhood graph transformation will determine final edge weights. If snn.quantile=NULL, the edge weight will be simply equal to the Jaccard coefficient of the neighborhoods. If snn.quantile is a vector of two numeric values (p1, p2), they will be treated as quantile probabilities, and quantile values (q1,q2) on the set of all Jaccard coefficients (for all edges) will be determiend. The edge weights will then be reset, so that edges with Jaccard coefficients below or equal to q1 will be set to 0, and those with coefficients >=q2 will be set to 1. The rest of the weights will be mapped uniformly from [q1,q2]->[0,1] range. If a single numeric value is supplied, it will be treated as a symmetric quantile probability (i.e. snn.quantile=0.8 is equivalent to specifying snn.quantile=c(1-0.8,0.8)). (default: 0.9)
146. #' @param min.snn.jaccard numeric Minimum Jaccard coefficient required for a shared neighborhood graph edge (default: 0). The edges with Jaccard coefficients below this threshold will be removed (i.e. weight set to 0)
147. #' @param min.snn.weight numeric Shared nearest neighbor procedure will adjust the weights of the edges, and even eliminate some of the edges (by setting their weight to zero). The min.snn.weight parameter allows to set a minimal adjusted edge weight, so that the edge weight is never reduced beyond this level (and hence never deleted) (default: 0 - no adjustments)
148. #' @param snn.k.self integer Size of the within-sample neighorhood to be used in shared nearest neighbor calculations (default=k.self)
149. #' @param balance.edge.weights boolean Whether to balance edge weights to control for a cell- or sample- specific factor (default=FALSE)
150. #' @param balancing.factor.per.sample A covariate factor per sample that should be controlled for by adjusting edge weights in the joint graph (default=NULL)
151. #' @param balancing.factor.per.cell A per-cell factor (discrete factor, named with cell names) specifying a design difference should be controlled for by adjusting edge weights in the joint graph (default=NULL)
152. #' @param same.factor.downweight numeric Optional weighting factor for edges connecting cells with the same cell factor level per cell balancing (default=1.0)
153. #' @param k.same.factor integer An neighborhood size that should be used when aligning samples of the same balancing.factor.per.sample level. Setting a value smaller than k will lead to reduction of alingment strenth within the sample batches (default=k)
154. #' @return joint graph to be used for downstream analysis
155. #' @examples
156. #' con <- Conos$new(small_panel.preprocessed, n.cores=1)
157. #' con$buildGraph(k=10, k.self=5, space='PCA', ncomps=10, n.odgenes=20, matching.method='mNN',
158. #' metric='angular', score.component.variance=TRUE, verbose=TRUE)
159. #'
160. #'
161. buildGraph=function(k=15, k.self=10, k.self.weight=0.1, alignment.strength=NULL, space='PCA', matching.method='mNN', metric='angular', k1=k, data.type='counts', l2.sigma=1e5, var.scale=TRUE, ncomps=40,
162. n.odgenes=2000, matching.mask=NULL, exclude.samples=NULL, common.centering=TRUE, verbose=TRUE,
163. base.groups=NULL, append.global.axes=TRUE, append.decoys=TRUE, decoy.threshold=1, n.decoys=k*2, score.component.variance=FALSE,
164. snn=FALSE, snn.quantile=0.9, min.snn.jaccard=0, min.snn.weight=0, snn.k.self=k.self,
165. balance.edge.weights=FALSE, balancing.factor.per.cell=NULL, same.factor.downweight=1.0, k.same.factor=k, balancing.factor.per.sample=NULL) {
167. supported.spaces <- c("CPCA","JNMF","genes","PCA","PMA","CCA")
168. if (!space %in% supported.spaces) {
169. stop(paste0("only the following spaces are currently supported: [",paste(supported.spaces,collapse=' '),"]"))
170. }
171. supported.matching.methods <- c("mNN", "NN")
172. if (!matching.method %in% supported.matching.methods) {
173. stop(paste0("only the following matching methods are currently supported: ['",paste(supported.matching.methods,collapse="' '"),"']"))
174. }
175. supported.metrics <- c("L2","angular")
176. if (!metric %in% supported.metrics) {
177. stop(paste0("only the following distance metrics are currently supported: ['",paste(supported.metrics,collapse="' '"),"']"))
178. }
180. if (!is.null(snn.quantile) && !is.na(snn.quantile)) {
181. if(length(snn.quantile)==1) {
182. snn.quantile <- c(1-snn.quantile,snn.quantile)
183. }
184. snn.quantile <- sort(snn.quantile,decreasing=FALSE)
185. if(snn.quantile[1]<0 | snn.quantile[2]>1) {
186. stop("snn.quantile must be one or two numbers in the [0,1] range")
187. }
188. }
190. if (!is.null(alignment.strength)) {
191. alignment.strength %<>% max(0) %>% min(1)
192. k1 <- sapply(self$samples, function(sample) ncol(getCountMatrix(sample))) %>% max() %>%
193. `*`(alignment.strength ^ 2) %>% round() %>% max(k)
194. } else {
195. alignment.strength <- 0 # otherwise, estimation of cor.base uses NULL value
196. }
198. if (k1<k) { stop("k1 must be >= k") }
199. # calculate or update pairwise alignments
200. sn.pairs <- private$updatePairs(space=space, ncomps=ncomps, n.odgenes=n.odgenes, verbose=verbose, var.scale=var.scale, matching.mask=matching.mask, exclude.samples=exclude.samples, score.component.variance=score.component.variance)
201. if (ncol(sn.pairs)<1) { stop("insufficient number of comparable pairs") }
202. if (!is.null(base.groups)) {
203. samf <- lapply(self$samples,getCellNames)
204. base.groups <- as.factor(base.groups[names(base.groups) %in% unlist(samf)]) # clean up the group factor
205. if (length(base.groups) < 2) stop("provided base.groups doesn't cover enough cells")
206. # make a sample factor
207. samf <- setNames(rep(names(samf),unlist(lapply(samf,length))),unlist(samf))
208. if (append.global.axes) {
209. cms.clust <- self$getClusterCountMatrices(groups=base.groups, common.genes=FALSE)
210. global.proj <- projectSamplesOnGlobalAxes(self$samples, cms.clust, data.type, verbose, self$n.cores)
211. }
212. }
214. if (snn){
215. local.neighbors <- getLocalNeighbors(self$samples[! names(self$samples) %in% exclude.samples], snn.k.self, k.self.weight, metric, l2.sigma=l2.sigma, verbose, self$n.cores)
216. } else {
217. local.neighbors <- NULL
218. }
220. # determine inter-sample mapping
221. if (verbose) message('inter-sample links using ',matching.method,' ')
222. cached.pairs <- self$pairs[[space]]
223. cor.base <- 1 + min(1, alignment.strength * 10) ## see convertDistanceToSimilarity()
224. mnnres <- papply(1:ncol(sn.pairs), function(j) {
225. # we'll look up the pair by name (possibly reversed), not to assume for the ordering of $pairs[[space]] to be the same
226. i <- match(paste(sn.pairs[,j],collapse='.vs.'),names(cached.pairs))
227. if(is.na(i)) { i <- match(paste(rev(sn.pairs[,j]),collapse='.vs.'),names(cached.pairs)) }
228. if(is.na(i)) { stop(paste("unable to find alignment for pair",paste(sn.pairs[,j],collapse='.vs.'))) }
230. k.cur <- k
231. if (!is.null(balancing.factor.per.sample) && (balancing.factor.per.sample[sn.pairs[1,j]] == balancing.factor.per.sample[sn.pairs[2,j]])) {
232. k.cur <- min(k.same.factor, k1) # It always should be less then k1, though never supposed to be set higher
233. }
235. if(space=='JNMF') {
236. mnn <- getNeighborMatrix(cached.pairs[[i]]$rot1, cached.pairs[[i]]$rot2,
237. k=k.cur, k1=k1, matching=matching.method, metric=metric, l2.sigma=l2.sigma, cor.base=cor.base)
238. } else if (space %in% c("CPCA","GSVD","PCA")) {
239. #common.genes <- Reduce(intersect,lapply(r.ns, getGenes))
240. if(!is.null(cached.pairs[[i]]$CPC)) {
241. # CPCA or PCA
242. #od.genes <- intersect(rownames(cached.pairs[[i]]$CPC),common.genes)
243. od.genes <- rownames(cached.pairs[[i]]$CPC)
244. rot <- cached.pairs[[i]]$CPC[od.genes,]
245. } else if(!is.null(cached.pairs[[i]]$o$Q)) {
246. # GSVD
247. rot <- cached.pairs[[i]]$o$Q
248. od.genes <- rownames(rot) <- colnames(cached.pairs[[i]]$o$A)
249. } else {
250. stop("unknown reduction provided")
251. }
253. # TODO: a more careful analysis of parameters used to calculate the cached version
254. if(ncomps > ncol(rot)) {
255. warning(paste0("specified ncomps (",ncomps,") is greater than the cached version (",ncol(rot),")"))
256. } else {
257. rot <- rot[,1:ncomps,drop=FALSE]
258. }
260. mnn <- getPcaBasedNeighborMatrix(self$samples[sn.pairs[,j]], od.genes=od.genes, rot=rot, data.type=data.type,
261. k=k.cur, k1=k1, matching.method=matching.method, metric=metric, l2.sigma=l2.sigma, cor.base=cor.base,
262. var.scale=var.scale, common.centering=common.centering,
263. base.groups=base.groups, append.decoys=append.decoys, samples=self$samples, samf=samf, decoy.threshold=decoy.threshold,
264. n.decoys=n.decoys, append.global.axes=append.global.axes, global.proj=global.proj)
265. } else if (space=='genes') { ## Overdispersed Gene space
266. mnn <- getNeighborMatrix(as.matrix(cached.pairs[[i]]$genespace1), as.matrix(cached.pairs[[i]]$genespace2),
267. k=k.cur, k1=k1, matching=matching.method, metric=metric, l2.sigma=l2.sigma, cor.base=cor.base)
268. } else if(space=='PMA' || space=='CCA') {
269. mnn <- getNeighborMatrix(cached.pairs[[i]]$u, cached.pairs[[i]]$v,
270. k=k.cur, k1=k1, matching=matching.method, metric=metric, l2.sigma=l2.sigma, cor.base=cor.base)
271. } else {
272. stop("Unknown space: ", space)
273. }
275. if(snn) { # optionally, perform shared neighbor weighting, a la SeuratV3, scran
277. m1 <- cbind(mnn,t(local.neighbors[[sn.pairs[1,j] ]]))
278. m2 <- rbind(local.neighbors[[sn.pairs[2,j] ]], mnn)
280. mnn1 <- mnn
281. mnn1@x <- rep(1,length(mnn1@x))
282. m1@x <- rep(1,length(m1@x))
283. m2@x <- rep(1,length(m2@x))
285. x <- ((m1 %*% m2) * mnn1) / pmax(outer(rowSums(m1),colSums(m2),FUN=pmin),1)
287. # scale by Jaccard coefficient
289. if(min.snn.jaccard>0) {
290. x@x[x@x<min.snn.jaccard] <- 0
291. }
293. if(!is.null(snn.quantile) && !is.na(snn.quantile)) {
294. xq <- quantile(x@x,p=c(snn.quantile[1],snn.quantile[2]))
295. x@x <- pmax(0,pmin(1,(x@x-xq[1])/pmax(1,diff(xq))))
296. }
298. x <- drop0(x)
299. if (min.snn.weight>0) {
300. if (packageVersion("Matrix") >= '1.4.2'){
301. mnn <- as(drop0(mnn*min.snn.weight + mnn*x),'TsparseMatrix')
302. } else {
303. mnn <- as(drop0(mnn*min.snn.weight + mnn*x),'dgTMatrix')
304. }
305. } else {
306. if (packageVersion("Matrix") >= '1.4.2'){
307. mnn <- as(drop0(mnn*x),'TsparseMatrix')
308. } else {
309. mnn <- as(drop0(mnn*x),'dgTMatrix')
310. }
311. }
313. }
315. if(verbose) cat(".")
316. return(data.frame('mA.lab'=rownames(mnn)[mnn@i+1],'mB.lab'=colnames(mnn)[mnn@j+1],'w'=mnn@x, stringsAsFactors=FALSE))
317. }, n.cores=self$n.cores,mc.preschedule=TRUE)
319. if (verbose) message(" done")
320. ## Merge the results into a edge table
321. el <- do.call(rbind,mnnres)
322. if (nrow(el)==0) {
323. el = data.frame('mA.lab'=0,'mB.lab'=0,'w'=0, 'type'=1, stringsAsFactors=FALSE)
324. } else {
325. el$type <- 1 # encode connection type 1- intersample, 0- intrasample
326. }
327. # append local edges
328. if(k.self>0) {
329. if(is.null(local.neighbors) || snn.k.self != k.self) { # recalculate local neighbors
330. local.neighbors <- getLocalNeighbors(self$samples[! names(self$samples) %in% exclude.samples], k.self, k.self.weight, metric, l2.sigma=l2.sigma, verbose, self$n.cores)
331. }
332. el <- rbind(el,getLocalEdges(local.neighbors))
333. }
334. if(verbose) message('building graph .')
335. el <- el[el[,3]>0,]
336. g <- graph_from_edgelist(as.matrix(el[,c(1,2)]), directed =FALSE)
337. E(g)$weight <- el[,3]
338. E(g)$type <- el[,4]
340. if(verbose) cat(".")
342. # collapse duplicate edges
343. g <- simplify(g, edge.attr.comb=list(weight="sum", type = "first"))
344. if(verbose) message('done')
346. if (!is.null(balancing.factor.per.sample)) {
347. if (is.null(balancing.factor.per.cell)) {
348. sf <- self$getDatasetPerCell()
349. balancing.factor.per.cell <- setNames(balancing.factor.per.sample[as.character(sf)], names(sf))
350. } else {
351. warning("Both balancing.factor.per.cell and balancing.factor.per.sample are provided. Used the former for balancing edge weights")
352. }
353. }
355. if (balance.edge.weights || !is.null(balancing.factor.per.cell)) {
356. if(verbose) message('balancing edge weights ')
358. if (is.null(balancing.factor.per.cell)) {
359. balancing.factor.per.cell <- self$getDatasetPerCell()
360. }
362. g <- igraph::as_adjacency_matrix(g, attr="weight") %>%
363. adjustWeightsByCellBalancing(factor.per.cell=balancing.factor.per.cell, balance.weights=balance.edge.weights,
364. same.factor.downweight=same.factor.downweight) %>%
365. igraph::graph_from_adjacency_matrix(mode="undirected", weighted=TRUE)
367. if(verbose) message('done')
368. }
369. self$graph <- g
370. return(invisible(g))
371. },
374. #' @description Calculate genes differentially expressed between cell clusters. Estimates base mean, z-score, p-values, specificity, precision, expressionFraction, AUC (if append.auc=TRUE)
375. #'
376. #' @param groups a cell factor (a factor named with cell names) specifying clusters of cells to be compared (one against all). To compare two cell clusters against each other, simply pass a factor containing only two levels (default: NULL, see clustering)
377. #' @param clustering character Name of the clustering to use (see names(con$clusters)) for the value of the groups factor (default: NULL - if groups are not specified, the first clustering will be used)
378. #' @param z.threshold numeric Minimum absolute value of a Z score for which the genes should be reported (default=3.0).
379. #' @param upregulated.only boolean If TRUE, will report only genes significantly upregulated in each cluster; otherwise both up- and down-regulated genes will be reported (default=FALSE)
380. #' @param append.specificity.metrics boolean Whether to append specificity metrics (default=TRUE)
381. #' @param append.auc boolean Whether to append AUC scores (default=TRUE)
382. #' @return list of DE results; each is a data frame with rows corresponding to the differentially expressed genes, and columns listing log2 fold change (M), signed Z scores (both raw and adjusted for mulitple hypothesis using BH correction), optional specificty/sensitivity and AUC metrics.
383. #'
384. getDifferentialGenes=function(clustering=NULL, groups=NULL, z.threshold=3.0, upregulated.only=FALSE, verbose=TRUE, append.specificity.metrics=TRUE, append.auc=TRUE) {
386. groups <- parseCellGroups(self, clustering, groups)
388. groups %<>% as.factor() %>% droplevels()
389. # TODO: add Seurat
390. '%ni%' <- Negate('%in%')
391. if ('Pagoda2' %ni% class(self$samples[[1]])) {
392. stop("Only Pagoda2 objects are supported for marker genes")
393. }
395. de.genes <- getDifferentialGenesP2(self$samples, groups=groups, z.threshold=z.threshold, upregulated.only=upregulated.only, verbose=verbose, n.cores=self$n.cores)
396. de.genes <- de.genes[levels(groups)]
397. if (append.specificity.metrics) {
398. if (verbose) message("Estimating specificity metrics")
400. cm.merged <- self$getJointCountMatrix(raw=TRUE)
401. groups.clean <- groups %>% .[!is.na(.)] %>% .[names(.) %in% rownames(cm.merged)]
402. de.genes %<>% lapply(function(x) if ((length(x) > 0) && (nrow(x) > 0)) subset(x, complete.cases(x)) else x)
403. de.genes %<>% names() %>% setNames(., .) %>%
404. sccore::plapply(function(n) appendSpecificityMetricsToDE(de.genes[[n]], groups.clean, n, p2.counts=cm.merged, append.auc=append.auc),
405. progress=verbose, n.cores=self$n.cores, fail.on.error=TRUE)
406. }
408. if (verbose) message("All done!")
410. return(de.genes)
411. },
413. #' @description Find cell clusters (as communities on the joint graph)
414. #'
415. #' @param method community detection method (igraph syntax) (default=leiden.community)
416. #' @param min.group.size numeric Minimal allowed community size (default=0)
417. #' @param name character Optional name of the clustering result (will default to the algorithm name) (default=NULL will try to obtain the name from the community detection method, or will use 'community' as a default)
418. #' @param test.stability boolean Whether to test stability of community detection (default=FALSE)
419. #' @param stability.subsampling.fraction numeric Fraction of clusters to subset (default=0.95). Must be within range [0, 1].
420. #' @param stability.subsamples integer Number of subsampling iterations (default=100)
421. #' @param cls optional pre-calculated community result (may be useful for stability testing) (default: NULL)
422. #' @param sr optional pre-calculated subsampled community results (useful for stability testing) (default: NULL)
423. #' @param ... extra parameters are passed to the specified community detection method
424. #' @return invisible list containing identified communities (groups) and the full community detection result (result); The results are stored in $clusters$name slot in the conos object. Each such slot contains an object with elements: $results which stores the raw output of the community detection method, and $groups which is a factor on cells describing the resulting clustering. The later can be used, for instance, in plotting: con$plotGraph(groups=con$clusters$leiden$groups). If test.stability==TRUE, then the result object will also contain a $stability slot.
425. #' @examples
426. #' con <- Conos$new(small_panel.preprocessed, n.cores=1)
427. #' con$buildGraph(k=10, k.self=5, space='PCA', ncomps=10, n.odgenes=20, matching.method='mNN',
428. #' metric='angular', score.component.variance=TRUE, verbose=TRUE)
429. #' con$findCommunities(method = igraph::walktrap.community, steps=5)
430. #'
431. findCommunities=function(method=leiden.community, min.group.size=0, name=NULL, test.stability=FALSE, stability.subsampling.fraction=0.95, stability.subsamples=100, verbose=TRUE, cls=NULL, sr=NULL, ...) {
433. if (is.null(cls)) {
434. cls <- method(self$graph, ...)
435. }
436. if (is.null(name)) {
437. name <- cls$algorithm
438. if(is.null(name)) {
439. name <- "community"
440. }
441. }
443. ## Extract groups from this graph
444. cls.mem <- membership(cls)
445. if(suppressWarnings(any(is.na(as.numeric(cls.mem))))) {
446. cls.groups <- as.factor(cls.mem)
447. } else {
448. cls.groups <- factor(setNames(as.character(cls.mem),names(cls.mem)),levels=sort(as.numeric(unique(as.character(cls.mem))),decreasing=FALSE))
449. }
450. cls.levs <- levels(cls.groups)
451. res <- list(groups=cls.groups,result=cls)
453. # test stability
454. if(test.stability) {
456. subset.clustering <- function(g,f=stability.subsampling.fraction,seed=NULL, ...) {
457. if(!is.null(seed)) { set.seed(seed) }
458. vi <- sample(1:length(V(g)),ceiling(length(V(g))*(f)))
459. sg <- induced_subgraph(g,vi)
460. method(sg,...)
461. }
462. if(verbose) { message("running ",stability.subsamples," subsampling iterations ... ")}
463. if (is.null(sr)) {
464. sr <- papply(1:stability.subsamples,function(i) subset.clustering(self$graph,f=stability.subsampling.fraction,seed=i),n.cores=self$n.cores)
465. }
467. if(verbose) { message("done")}
469. if(verbose) message("calculating flat stability stats ... ")
470. # Jaccard coefficient for each cluster against all, plus random expectation
471. jc.stats <- do.call(rbind, papply(sr,function(o) {
472. p1 <- membership(o)
473. p2 <- cls.groups[names(p1)]
474. p1 <- as.character(p1)
475. #x <- tapply(1:length(p2),factor(p2,levels=cls.levs),function(i1) {
476. x <- tapply(1:length(p2),p2,function(i1) {
477. i2 <- which(p1==p1[i1[[1]]])
478. length(intersect(i1,i2))/length(unique(c(i1,i2)))
479. })
481. }, n.cores=self$n.cores, mc.preschedule=TRUE))
484. # based on the following C code from clues:
485. # v0.2.4 on Feb. 3, 2009 by Weiliang Qiu, <https://github.com/cran/clues/blob/master/R/adjustedRand.R>
486. # (1) moved some code out of the for loop
487. #
488. # cl1 --- partition 1 of the data set
489. # cl2 --- partition 2 of the data set
490. #
491. # flag = 1 --- Rand index
492. # flag = 2 --- Hubert and Arabie's adjusted Rand index
493. # flag = 3 --- Morey and Agresti's adjusted Rand index
494. # flag = 4 --- Fowlkes and Mallows's index
495. # flag = 5 --- Jaccard index
497. # modified for C++, 19 January 2022
498. adjustedRand <- function(cl1, cl2, randMethod = c("Rand","HA", "MA", "FM", "Jaccard")){
499. if(!is.vector(cl1)){
500. stop("cl1 is not a vector!\n")
501. }
502. if(!is.vector(cl2)){
503. stop("cl2 is not a vector!\n")
504. }
505. if(length(cl1) != length(cl2)){
506. stop("Two vectors have different lengths!\n")
507. }
509. len <- length(randMethod)
510. if(len == 0){
511. stop("The argument 'randMethod' is empty!\n")
512. }
514. # unique values of elements in 'cl1'
515. cl1u <- unique(cl1)
516. # number of clusters in partition 1
517. m1 <- length(cl1u)
519. # unique values of elements in 'cl2'
520. cl2u <- unique(cl2)
521. # number of clusters in partition 2
522. m2 <- length(cl2u)
524. n <- length(cl1)
525. randVec <- rep(0, len)
526. names(randVec) <- randMethod
527. for(i in 1:len){
528. randMethod[i] <- match.arg(arg = randMethod[i],
529. choices = c("Rand","HA", "MA", "FM", "Jaccard"))
531. flag <- match(randMethod[i],
532. c("Rand","HA", "MA", "FM", "Jaccard"))
534. c.res <- adjustedRandcpp(as.integer(cl1),
535. as.integer(cl1u),
536. as.integer(cl2),
537. as.integer(cl2u),
538. as.integer(m1),
539. as.integer(m2),
540. as.integer(n),
541. as.integer(flag))
542. randVec[i] <- c.res
543. }
544. return(randVec)
545. }
547. # Adjusted rand index
548. if(verbose) message("adjusted Rand ... ")
549. ari <- unlist(papply(sr,function(o) { ol <- membership(o); adjustedRand(as.integer(ol),as.integer(cls.groups[names(ol)]),randMethod='HA') }, n.cores=self$n.cores))
550. if(verbose) message("done")
552. res$stability <- list(flat=list(jc=jc.stats,ari=ari))
554. # hierarchical measures
555. if(verbose) message("calculating hierarchical stability stats ... ")
556. if(is.hierarchical(cls)) {
557. # hierarchical to hierarchical stability analysis - cut reference
558. # determine hierarchy of clusters (above the cut)
559. t.get.walktrap.upper.merges <- function(res,n=length(unique(membership(res)))) {
560. clm <- complete.dend(res,FALSE)
561. x <- tail(clm,n-1)
562. x <- x - 2*nrow(res$merges) + nrow(x)-1
563. # now all >=0 ids are cut leafs and need to be reassigned ids according to their rank
564. xp <- x+nrow(x)+1
565. xp[x<=0] <- rank(-x[x<=0])
566. xp
567. }
569. clm <- t.get.walktrap.upper.merges(cls)
570. res$stability$upper.tree <- clm
572. if(verbose) message("tree Jaccard ... ")
573. jc.hstats <- do.call(rbind, papply(sr,function(z) bestClusterThresholds(z,cls.groups,clm)$threshold, n.cores=self$n.cores))
574. } else {
575. # compute cluster hierarchy based on cell mixing (and then something)
576. # assess stability for that hierarchy (to visualize internal node stability)
577. # for the original clustering and every subsample clustering,
578. if(verbose) message("upper clustering ... ")
579. cgraph <- getClusterGraph(self$graph,cls.groups,plot=FALSE,normalize=FALSE)
580. chwt <- walktrap.community(cgraph,steps=9)
581. clm <- complete.dend(chwt,FALSE)
583. if(verbose) message("clusterTree Jaccard ... ")
584. jc.hstats <- do.call(rbind, papply(sr,function(st1) {
585. mf <- membership(st1)
586. mf <- as.factor(setNames(as.character(mf),names(mf)))
588. st1g <- getClusterGraph(self$graph,mf, plot=FALSE, normalize=TRUE)
590. st1w <- walktrap.community(st1g, steps=8)
592. #merges <- st1w$merge; leaf.factor <- mf; clusters <- cls.groups
593. x <- bestClusterTreeThresholds(st1w,mf,cls.groups,clm)
594. x$threshold
595. }, n.cores=self$n.cores))
597. }
598. res$stability$upper.tree <- clm
599. res$stability$sr <- sr
600. res$stability$hierarchical <- list(jc=jc.hstats)
601. if(verbose) message("done")
603. }
605. ## Filter groups
606. if(min.group.size>0) {
607. lvls.keep <- names(which(table(cls.groups) > min.group.size))
608. cls.groups[! as.character(cls.groups) %in% as.character(lvls.keep)] <- NA
609. cls.groups <- as.factor(cls.groups)
610. }
611. res$groups <- cls.groups
612. self$clusters[[name]] <- res
613. return(invisible(res))
615. },
617. #' @description Plot panel of individual embeddings per sample with joint coloring
618. #'
619. #' @param groups a cell factor (a factor named with cell names) specifying clusters of cells to be compared (one against all). To compare two cell clusters against each other, simply pass a factor containing only two levels (default=NULL, see clustering)
620. #' @param clustering character Name of the clustering to use (see names(con$clusters)) for the value of the groups factor (default=NULL - if groups are not specified, the first clustering will be used)
621. #' @param use.local.clusters boolean Whether clusters should be taken from the individual samples; otherwise joint clusters in the conos object will be used (see clustering) (default=FALSE).
622. #' @param plot.theme string Theme for the plot, passed to plotSamples() (default=NULL)
623. #' @param use.common.embedding boolean Whether a joint embedding in the conos object should be used (or embeddings determined for the individual samples) (default=FALSE)
624. #' @param embedding (default=NULL) If a character value is passed, it is interpreted as an embedding name (a name of a joint embedding in conos when use.commmon.embedding=TRUE, or a name of an embedding within the individual objects when use.common.embedding=FALSE).
625. #' If a matrix is passed, it is interpreted as an actual embedding (then first two columns are interpreted as x/y coordinates, row names must be cell names). If NULL, the default embedding will be used.
626. #' @param adj.list an optional list of additional ggplot2 directions to apply (default=NULL)
627. #' @param ... Additional parameters passed to plotSamples(), plotEmbeddings(), sccore::embeddingPlot().
628. #' @return cowplot grid object with the panel of plots
629. #'
630. plotPanel=function(clustering=NULL, groups=NULL, colors=NULL, gene=NULL, use.local.clusters=FALSE, plot.theme=NULL, use.common.embedding=FALSE, embedding=NULL, adj.list=NULL, ...) {
632. if (use.local.clusters) {
633. if (is.null(clustering) && !(inherits(x = self$samples[[1]], what = c('seurat', 'Seurat')))) {
634. stop("You have to provide 'clustering' parameter to be able to use local clusters")
635. }
637. groups <- lapply(self$samples, getClustering, clustering) %>%
638. lapply(function(cls) setNames(as.character(cls), names(cls))) %>% Reduce(c, .)
640. if (is.null(groups)) {
641. stop(paste0("No clustering '", clustering, "' presented in the samples"))
642. }
643. } else if (is.null(groups) && is.null(colors) && is.null(gene)) {
644. groups <- getClusteringGroups(self$clusters, clustering)
645. }
647. if (use.common.embedding) { # look up the embedding within the conos object
648. ## if use.common.embedding, pass the Conos embedding to plotSamples
649. embedding <- private$getEmbedding(embedding)
650. adj.list <- c(ggplot2::lims(x=range(embedding[,1]), y=range(embedding[,2])), adj.list)
651. }
653. plotSamples(self$samples, groups=groups, colors=colors, gene=gene, plot.theme=private$adjustTheme(plot.theme), embedding.type=embedding, adj.list=adj.list, ...)
654. },
656. #' @description Generate an embedding of a joint graph
657. #'
658. #' @param method Embedding method (default='largeVis'). Currently 'largeVis' and 'UMAP' are supported.
659. #' @param embedding.name character Optional name of the name of the embedding set by user to store multiple embeddings (default: method name)
660. #' @param M numeric (largeVis) The number of negative edges to sample for each positive edge to be used (default=1)
661. #' @param gamma numeric (largeVis) The strength of the force pushing non-neighbor nodes apart (default=1)
662. #' @param alpha numeric (largeVis) Hyperparameter used in the default distance function, \eqn{1 / (1 + \alpha \dot ||y_i - y_j||^2)} (default=0.1). The function relates the distance
663. #' between points in the low-dimensional projection to the likelihood that the two points are nearest neighbors. Increasing \eqn{\alpha} tends
664. #' to push nodes and their neighbors closer together; decreasing \eqn{\alpha} produces a broader distribution. Setting \eqn{\alpha} to zero
665. #' enables the alternative distance function. \eqn{\alpha} below zero is meaningless.
666. #' @param perplexity (largeVis) The perplexity passed to largeVis (default=NA)
667. #' @param sgd_batches (largeVis) The number of edges to process during SGD (default=1e8). Defaults to a value set based on the size of the dataset. If the parameter given is
668. #' between \code{0} and \code{1}, the default value will be multiplied by the parameter.
669. #' @param seed numeric Random seed for the largeVis algorithm (default=1)
670. #' @param target.dims numeric Number of dimensions for the reduction (default=2). Higher dimensions can be used to generate embeddings for subsequent reductions by other methods, such as tSNE
671. #' @param ... additional arguments, passed to UMAP embedding (run ?conos:::embedGraphUmap for more info)
672. #'
673. embedGraph=function(method='largeVis', embedding.name=method, M=1, gamma=1, alpha=0.1, perplexity=NA, sgd_batches=1e8, seed=1, verbose=TRUE, target.dims=2, ...) {
674. supported.methods <- c('largeVis', 'UMAP')
675. if(!method %in% supported.methods) {
676. stop(paste0("Currently, only the following embeddings are supported: ",paste(supported.methods,collapse=' ')))
677. }
679. ## check if embedding.name already in list
680. ## if so, throw warning
681. if (length(self$embeddings)>0){
682. ## check if embedding.name already created
683. if (embedding.name %in% names(self$embeddings)){
684. warning(paste0("Already created an embedding: ", embedding.name, ". Overwriting."))
685. }
686. }
688. if (method == 'largeVis') {
689. wij <- as_adj(self$graph,attr='weight')
690. if(!is.na(perplexity)) {
691. wij <- buildWijMatrix(wij,perplexity=perplexity, threads=self$n.cores)
692. }
693. coords <- projectKNNs(wij = wij, dim=target.dims, verbose = verbose,sgd_batches = sgd_batches,gamma=gamma, M=M, seed=seed, alpha=alpha, rho=1, threads=self$n.cores)
694. colnames(coords) <- V(self$graph)$name
695. self$embedding <- t(coords)
696. embedding.result <- self$embedding
697. } else {
698. ## method == 'UMAP'
699. if (!requireNamespace("uwot", quietly=TRUE)){
700. stop("You need to install package 'uwot' to be able to use UMAP embedding. Please install it.")
701. }
703. self$embedding <- embedGraphUmap(self$graph, verbose=verbose, return.all=FALSE, n.cores=self$n.cores, target.dims=target.dims, ...)
704. embedding.result <- self$embedding
705. }
707. self$embeddings[[embedding.name]] <- embedding.result
708. self$embedding <- embedding.result # hang on to the latest embedding for backwards compatibility
709. return(invisible(embedding.result))
710. },
712. #' @description Plot cluster stability statistics.
713. #'
714. #' @param clustering string Name of the clustering result to show (default=NULL)
715. #' @param what string Show a specific plot (ari - adjusted rand index, fjc - flat Jaccard, hjc - hierarchical Jaccard, dend - cluster dendrogram, all - everything except 'dend') (default='all')
716. #' @return cluster stability statistics
717. plotClusterStability=function(clustering=NULL, what='all') {
718. if(is.null(clustering)){
719. clustering <- names(self$clusters)[[1]]
720. }
722. if(is.null(self$clusters[[clustering]])){
723. stop(paste("clustering",clustering,"doesn't exist, run findCommunity() first"))
724. }
726. if(is.null(self$clusters[[clustering]]$stability)){
727. stop(paste("clustering",clustering,"doesn't have stability info. Run findCommunity( ... , test.stability=TRUE) first"))
728. }
730. st <- self$clusters[[clustering]]$stability
731. nclusters <- ncol(st$flat$jc)
732. jitter.alpha <- 0.1
734. if (what=='all' || what=='ari') {
735. p.fai <- ggplot2::ggplot(data.frame(aRI=st$flat$ari), ggplot2::aes(x=1,y=aRI)) +
736. ggplot2::geom_boxplot(notch=TRUE, outlier.shape=NA) +
737. ggplot2::geom_point(shape=16, position = ggplot2::position_jitter(), alpha=jitter.alpha) +
738. ggplot2::guides(color=FALSE) +
739. ggplot2::geom_hline(yintercept=1, linetype="dashed", alpha=0.2) +
740. ggplot2::ylim(c(0,1)) + ggplot2::labs(x=" ", y="adjusted Rand Index") +
741. ggplot2::theme(legend.position="none", axis.ticks.x=ggplot2::element_blank(), axis.text.x=ggplot2::element_blank())
743. if(what=='ari')
744. return(p.fai)
745. }
747. if (what=='all' || what=='fjc') {
748. df <- reshape2::melt(st$flat$jc)
749. colnames(df) <- c('rep','cluster','jc')
750. df$cluster <- factor(colnames(st$flat$jc)[df$cluster],levels=levels(self$clusters[[clustering]]$groups))
752. p.fjc <- ggplot2::ggplot(df,aes(x=cluster,y=jc,color=cluster)) +
753. ggplot2::geom_boxplot(aes(color=cluster), notch=TRUE, outlier.shape=NA) +
754. ggplot2::geom_jitter(shape=16, position=position_jitter(0.2), alpha=jitter.alpha) +
755. ggplot2::guides(color=FALSE) +
756. ggplot2::geom_hline(yintercept=1, linetype="dashed", alpha=0.2) +
757. ggplot2::ylab("Jaccard coefficient (flat)") + ggplot2::ylim(c(0,1))
759. if(what=='fjc') return(p.fjc)
760. }
762. if (what=='all' || what=='hjc') {
763. # hierarchical
764. df <- reshape2::melt(st$hierarchical$jc[,1:nclusters])
765. colnames(df) <- c('rep','cluster','jc')
766. df$cluster <- factor(colnames(st$flat$jc)[df$cluster],levels=levels(self$clusters[[clustering]]$groups))
767. p.hjc <- ggplot2::ggplot(df,aes(x=cluster,y=jc,color=cluster)) +
768. ggplot2::geom_boxplot(aes(color=cluster), notch=TRUE, outlier.shape=NA) +
769. ggplot2::geom_jitter(shape=16, position=ggplot2::position_jitter(0.2), alpha=jitter.alpha) +
770. ggplot2::guides(color=FALSE) +
771. ggplot2::geom_hline(yintercept=1, linetype="dashed", alpha=0.2) +
772. ggplot2::ylab("Jaccard coefficient (hierarchical)") + ggplot2::ylim(c(0,1))
774. if(what=='hjc') return(p.hjc)
775. }
777. if (what=='dend') {
778. m <- st$upper.tree
779. nleafs <- nrow(m)+1
780. m[m<=nleafs] <- -1*m[m<=nleafs]
781. m[m>0] <- m[m>0]-nleafs
782. hc <- list(merge=m, height=1:nrow(m), labels=levels(self$clusters[[clustering]]$groups), order=c(1:nleafs))
783. class(hc) <- 'hclust'
784. # fix the ordering so that edges don't intersects
785. hc$order <- order.dendrogram(as.dendrogram(hc))
787. d <- as.dendrogram(hc) %>% dendextend::hang.dendrogram()
789. # depth-first traversal of a merge matrix
790. t.dfirst <- function(m,i=nrow(m)) {
791. rl <- m[i,1]; if(rl<0) { rl <- abs(rl) } else { rl <- t.dfirst(m,rl) }
792. rr <- m[i,2]; if(rr<0) { rr <- abs(rr) } else { rr <- t.dfirst(m,rr) }
793. c(i+nrow(m)+1,rl,rr)
794. }
795. xy <- dendextend::get_nodes_xy(d)
796. to <- t.dfirst(hc$merge)
797. plot(d,las=2,axes=FALSE)
798. # flat on the left
799. #x <- apply(st$flat$jc,2,median)
800. #text(xy,labels=round(x[to],2),col='blue',adj=c(-0.1,-1.24),cex=0.8)
801. x <- apply(st$hierarchical$jc,2,median)
802. text(xy,labels=round(x[to],2),col='red',adj=c(-0.1,-0.12),cex=0.8)
803. return(NULL)
804. }
806. cowplot::plot_grid(plotlist=list(p.fai,p.fjc,p.hjc),nrow=1,rel_widths=c(4,nclusters,nclusters))
807. },
809. #' @description Plot joint graph
810. #'
811. #' @param groups a cell factor (a factor named with cell names) specifying clusters of cells to be compared (one against all). To compare two cell clusters against each other, simply pass a factor containing only two levels (default: NULL, see clustering)
812. #' @param clustering a character name of the clustering to use (see names(con$clusters)) for the value of the groups factor (default: NULL - if groups are not specified, the first clustering will be used)
813. #' @param color.by character A shortcut to color the plot by 'cluster' or by 'sample' (default: 'cluster'). If any other string is input, an error is thrown.
814. #' @param embedding A character name of an embedding, or a matrix of the actual embedding (rownames should correspond to cells, first to columns to x/y coordinates). If NULL (default: NULL), the latest generated embedding will be used
815. #' @param colors a color factor (named with cell names) use for cell coloring (default=NULL)
816. #' @param gene Show expression of a gene (default=NULL)
817. #' @param plot.theme Theme for the plot, passed to sccore::embeddingPlot() (default=NULL)
818. #' @param subset A subset of cells to show (default: NULL - shows all the cells)
819. #' @param ... Additional parameters passed to sccore::embeddingPlot()
820. #' @return ggplot2 plot of joint graph
821. #'
822. plotGraph=function(color.by='cluster', clustering=NULL, embedding=NULL, groups=NULL, colors=NULL, gene=NULL, plot.theme=NULL, subset=NULL, ...) {
823. embedding <- private$getEmbedding(embedding)
825. if (!is.null(subset)) {
826. embedding <- embedding[rownames(embedding) %in% subset,,drop=FALSE]
827. }
829. if (!is.null(gene)) {
830. colors <- lapply(self$samples, getGeneExpression, gene) %>% Reduce(c, .)
831. if(all(is.na(colors))) stop(paste("Gene", gene,"is not found in any of the samples"))
832. }
834. if(is.null(groups) && is.null(colors)) {
835. if(color.by == 'cluster') {
836. groups <- getClusteringGroups(self$clusters, clustering)
837. } else if(color.by == 'sample') {
838. groups <- self$getDatasetPerCell()
839. } else {
840. stop('Supported values of color.by are ("cluster" and "sample"); Use groups/colors parameters to explicitly pass factor/numeric data for coloring')
841. }
842. }
844. return(embeddingPlot(embedding, groups=groups, colors=colors, plot.theme=private$adjustTheme(plot.theme), ...))
845. },
847. #' @description Smooth expression of genes to minimize the batch effect between samples
848. #' Use diffusion of expression on graph with the equation dv = exp(-a * (v + b))
849. #'
850. #' @param genes List of genes to be smooothed smoothing (default=NULL will smooth top n.od.genes overdispersed genes)
851. #' @param n.od.genes numeric If 'genes' is NULL, top n.od.genes of overdispersed genes are taken across all samples (default=500)
852. #' @param fading numeric Level of fading of expression change from distance on the graph (parameter 'a' of the equation) (default=10)
853. #' @param fading.const numeric Minimal penalty for each new edge during diffusion (parameter 'b' of the equation) (default=0.5)
854. #' @param max.iters numeric Maximal number of diffusion iterations (default=15)
855. #' @param tol numeric Tolerance after which the diffusion stops (default=5e-3)
856. #' @param name string Name to save the correction (default='diffusion')
857. #' @param verbose boolean Verbose mode (default=TRUE)
858. #' @param count.matrix Alternative gene count matrix to correct (rows: genes, columns: cells; has to be dense matrix). Default: joint count matrix for all datasets.
859. #' @param normalize boolean Whether to normalize values (default=TRUE)
860. #' @return smoothed expression of the input genes
861. #'
862. correctGenes=function(genes=NULL, n.od.genes=500, fading=10.0, fading.const=0.5, max.iters=15, tol=5e-3, name='diffusion', verbose=TRUE, count.matrix=NULL, normalize=TRUE) {
863. edges <- igraph::as_edgelist(self$graph)
864. edge.weights <- igraph::edge.attributes(self$graph)$weight
866. if (is.null(count.matrix)) {
867. if (is.null(genes)) {
868. genes <- getOdGenesUniformly(self$samples, n.genes=n.od.genes)
869. }
871. cms <- lapply(self$samples, getCountMatrix, transposed=TRUE)
872. genes <- Reduce(intersect, lapply(cms, colnames)) %>% intersect(genes)
873. count.matrix <- Reduce(rbind, lapply(cms, function(x) x[, genes])) %>% as.matrix()
874. } else {
875. count.matrix <- t(count.matrix)
876. }
877. vn <- V(self$graph)$name
878. if(!all(rownames(count.matrix)==vn)) { # subset to a common set of genes
879. if(!all(vn %in% rownames(count.matrix))) {
880. stop("count.matrix does not provide values for all the vertices in the alignment graph!")
881. }
882. count.matrix <- count.matrix[vn,]
883. }
885. ## Wrapper to make is.label.fixed optional
886. smoothMatrixOnGraph <- function(edges, edge.weights, matrix, is.label.fixed=logical(), ...) {
887. smooth_count_matrix(edges, edge.weights, matrix, is_label_fixed=is.label.fixed, ...)
888. }
890. cm <- smoothMatrixOnGraph(edges, edge.weights, count.matrix, max_n_iters=max.iters, diffusion_fading=fading,
891. diffusion_fading_const=fading.const, verbose=verbose, normalize=normalize)
892. return(invisible(self$expression.adj[[name]] <<- cm))
893. },
895. #' @description Estimate labeling distribution for each vertex, based on a partial labeling of the cells.
896. #' There are two methods used for the propagation to calculate the distribution of labels: "solver" and "diffusion".
897. #' * "diffusion" (default) will estimate the labeling distribution for each vertex, based on provided labels using a random walk.
898. #' * "solver" will propagate labels using the algorithm described by Zhu, Ghahramani, Lafferty (2003) <http://mlg.eng.cam.ac.uk/zoubin/papers/zgl.pdf>
899. #' Confidence values are then calculated by taking the maximum value from this distribution of labels, for each cell.
900. #'
901. #' @param labels Input labels
902. #' @param method type of propagation. Either 'diffusion' or 'solver'. 'solver' gives better result
903. #' but has bad asymptotics, so is inappropriate for datasets > 20k cells. (default='diffusion')
904. #' @param ... additional arguments for conos:::propagateLabels* functions
905. #' @return list with three fields:
906. #' * labels = matrix with distribution of label probabilities for each vertex by rows.
907. #' * uncertainty = 1 - confidence values
908. #' * label.distribution = the distribution of labels calculated using either the methods "diffusion" or "solver"
909. #'
910. propagateLabels=function(labels, method="diffusion", ...) {
911. if (method == "solver") {
912. label.dist <- propagateLabelsSolver(self$graph, labels, ...)
913. } else if (method == "diffusion") {
914. label.dist <- propagateLabelsDiffusion(self$graph, labels, ...)
915. } else {
916. stop("Unknown method: ", method, ". Only 'solver' and 'diffusion' are supported.")
917. }
919. labels <- colnames(label.dist)[apply(label.dist, 1, which.max)] %>%
920. setNames(rownames(label.dist))
922. confidence <- apply(label.dist, 1, max) %>% setNames(rownames(label.dist))
924. return(list(labels=labels, uncertainty=(1 - confidence), label.distribution=label.dist))
925. },
927. #' @description Calculate pseudo-bulk expression matrices for clusters (by adding up, for each gene, all of the molecules detected for all cells in a given cluster in a given sample)
928. #'
929. #' @param common.genes boolean Whether to bring individual sample matrices to a common gene list (default=TRUE)
930. #' @param omit.na.cells boolean If set to FALSE, the resulting matrices will include a first column named 'NA' that will report total molecule counts for all of the cells that were not covered by the provided factor. (default=TRUE)
931. #' @return a list of per-sample uniform dense matrices with rows being genes, and columns being clusters
932. #'
933. getClusterCountMatrices=function(clustering=NULL, groups=NULL, common.genes=TRUE, omit.na.cells=TRUE) {
934. if(is.null(groups)) {
935. groups <- getClusteringGroups(self$clusters, clustering)
936. }
938. groups <- as.factor(groups)
940. matl <- lapply(self$samples,function(s) {
941. m <- getRawCountMatrix(s,trans=TRUE) # rows are cells
942. cl <- factor(groups[match(rownames(m),names(groups))],levels=levels(groups));
943. tc <- colSumByFactor(m,cl)
944. if(omit.na.cells) { tc <- tc[-1,,drop=FALSE] }
945. t(tc)
946. })
947. # bring to a common gene space
948. if(common.genes) {
949. gs <- unique(unlist(lapply(matl,rownames)))
950. matl <- lapply(matl,function(m) {
951. nm <- matrix(0,nrow=length(gs),ncol=ncol(m))
952. colnames(nm) <- colnames(m)
953. rownames(nm) <- gs
954. mi <- match(rownames(m),gs)
955. nm[mi,] <- m
956. nm
957. })
958. }
959. matl
960. },
962. #' @description applies 'getCellNames()' on all samples
963. #' @return list of cellnames for all samples
964. #' @examples
965. #' con <- Conos$new(small_panel.preprocessed, n.cores=1)
966. #' con$getDatasetPerCell()
967. #'
968. getDatasetPerCell=function() {
969. getSampleNamePerCell(self$samples)
970. },
972. #' @description Retrieve joint count matrices
973. #'
974. #' @param raw boolean If TRUE, return merged "raw" count matrices, using function getRawCountMatrix(). Otherwise, return the merged count matrices, using getCountMatrix(). (default=FALSE)
975. #' @return list of merged count matrices
976. #' @examples
977. #' con <- Conos$new(small_panel.preprocessed, n.cores=1)
978. #' con$getJointCountMatrix()
979. #'
980. getJointCountMatrix=function(raw=FALSE) {
981. lapply(self$samples, (if (raw) getRawCountMatrix else getCountMatrix), transposed=TRUE) %>%
982. mergeCountMatrices(transposed=TRUE)
983. }
984. ),
988. ## Private functions
989. private = list(
990. adjustTheme=function(theme) {
991. if (is.null(theme)) {
992. theme <- ggplot2::theme()
993. }
994. main.theme <- ggplot2::theme_bw() + ggplot2::theme(
995. legend.background=ggplot2::element_rect(fill=ggplot2::alpha("white", 0.6)),
996. plot.margin=ggplot2::margin()
997. )
999. if (self$override.conos.plot.theme) {
1000. return(main.theme + ggplot2::theme_get() + theme)
1001. }
1003. return(main.theme + theme)
1004. },
1006. # a utility function to look up an embedding by name or accept an actual embedding data
1007. getEmbedding=function(embedding) {
1008. if (!is.null(embedding)) {
1009. if (class(embedding) %in% c('matrix')) { # actuall embedding was passed
1010. # check validity?
1011. } else if (inherits(embedding, 'character')) { # look up embedding by name
1012. ## check if embedding.name exists in list
1013. if (embedding %in% names(self$embeddings)) {
1014. ## embedding to plot
1015. embedding <- self$embeddings[[embedding]]
1016. } else {
1017. ## embedding.name not in list of self$embeddings, so the user is confused
1018. ## throw error
1019. stop(paste0("No embedding named '", embedding, "' found. Please generate this with embedGraph()."))
1020. }
1021. } else {
1022. stop('embedding must be either a character name of the embedding, or an actual matrix of embedding coordinates')
1023. }
1024. } else {
1025. if(!is.null(self$embedding)) { # use the latest
1026. embedding <- self$embedding
1027. } else { # pick one from the list
1028. if(is.null(self$embeddings) || length(self$embeddings)<1) stop("no joint embeddings have been generated; use embedGraph() first")
1029. embedding <- self$embedding[length(self$embedding)] # by default, pick last-named embedding
1030. }
1031. }
1032. return(embedding)
1033. },
1035. updatePairs=function(space='PCA', data.type='counts', ncomps=50, n.odgenes=1e3, var.scale=TRUE, matching.mask=NULL, exclude.samples=NULL, score.component.variance=FALSE, verbose=FALSE) {
1037. # make a list of all pairs
1038. sample.names <- names(self$samples)
1039. if(!is.null(exclude.samples)) {
1040. mi <- sample.names %in% exclude.samples
1041. if(verbose) { message("excluded ", sum(mi), " out of ", length(sample.names), " samples, based on supplied exclude.samples") }
1042. sample.names <- sample.names[!mi]
1043. }
1045. # TODO: add random subsampling for very large panels
1046. if(!is.null(matching.mask)) { # remove pairs that shouldn't be compared directly
1047. tryCatch(matching.mask <- matching.mask[sample.names, sample.names],
1048. error=function(e) stop("matching.mask should have the same row- and colnames as provided samples. Error:", e))
1050. matching.mask <- matching.mask | t(matching.mask)
1051. selected.ids <- which(lower.tri(matching.mask) & matching.mask) - 1
1052. sn.pairs <- sample.names[selected.ids %/% length(sample.names) + 1] %>%
1053. cbind(sample.names[selected.ids %% length(sample.names) + 1]) %>%
1054. t()
1056. if(verbose) message("Use ", ncol(sn.pairs), " pairs, based on the passed exclude.pairs")
1057. } else {
1058. sn.pairs <- combn(sample.names, 2)
1059. }
1061. # determine the pairs that need to be calculated
1062. if (is.null(self$pairs[[space]])) {
1063. self$pairs[[space]] <- list()
1064. }
1065. mi <- rep(NA,ncol(sn.pairs))
1066. nm <- match(apply(sn.pairs,2,paste,collapse='.vs.'),names(self$pairs[[space]]))
1067. mi[which(!is.na(nm))] <- na.omit(nm)
1068. # try reverse match as well
1069. nm <- match(apply(sn.pairs[c(2,1),,drop=FALSE],2,paste,collapse='.vs.'),names(self$pairs[[space]]));
1070. mi[which(!is.na(nm))] <- na.omit(nm)
1071. if(verbose) message('found ',sum(!is.na(mi)),' out of ',length(mi),' cached ',space,' space pairs ... ')
1072. if(any(is.na(mi))) { # some pairs are missing
1073. if(verbose) message('running ',sum(is.na(mi)),' additional ',space,' space pairs ')
1074. xl2 <- sccore::plapply(which(is.na(mi)), function(i) {
1075. if (space=='CPCA') {
1076. xcp <- quickCPCA(self$samples[sn.pairs[,i]],data.type=data.type,ncomps=ncomps,n.odgenes=n.odgenes,verbose=FALSE,var.scale=var.scale, score.component.variance=score.component.variance)
1077. } else if(space=='JNMF') {
1078. xcp <- quickJNMF(self$samples[sn.pairs[,i]],data.type=data.type,n.comps=ncomps,n.odgenes=n.odgenes,var.scale=var.scale,verbose=FALSE,max.iter=3e3)
1079. } else if (space == 'genes') {
1080. xcp <- quickNULL(p2.objs = self$samples[sn.pairs[,i]], data.type=data.type, n.odgenes=n.odgenes, var.scale = var.scale, verbose = FALSE)
1081. } else if (space == 'PCA') {
1082. xcp <- quickPlainPCA(self$samples[sn.pairs[,i]], data.type=data.type,ncomps=ncomps,n.odgenes=n.odgenes,verbose=FALSE,var.scale=var.scale, score.component.variance=score.component.variance)
1083. } else if (space == 'CCA' || space=='PMA') {
1084. xcp <- quickCCA(self$samples[sn.pairs[,i]],data.type=data.type,ncomps=ncomps,n.odgenes=n.odgenes,verbose=FALSE,var.scale=var.scale, score.component.variance=score.component.variance,PMA=(space=='PMA'))
1085. }
1086. if(verbose) cat(".")
1087. xcp
1088. }, n.cores=self$n.cores, mc.preschedule=(space=='PCA'), progress=FALSE, fail.on.error=TRUE)
1090. names(xl2) <- apply(sn.pairs[,which(is.na(mi)),drop=FALSE],2,paste,collapse='.vs.')
1091. xl2 <- xl2[!unlist(lapply(xl2,is.null))]
1092. self$pairs[[space]] <- c(self$pairs[[space]],xl2)
1093. }
1095. # re-do the match and order
1096. mi <- rep(NA,ncol(sn.pairs))
1097. nm <- match(apply(sn.pairs,2,paste,collapse='.vs.'),names(self$pairs[[space]]))
1098. mi[which(!is.na(nm))] <- na.omit(nm)
1099. nm <- match(apply(sn.pairs[c(2,1),,drop=FALSE],2,paste,collapse='.vs.'),names(self$pairs[[space]]))
1100. mi[which(!is.na(nm))] <- na.omit(nm)
1101. if(any(is.na(mi))) {
1102. warning("unable to get complete set of pair comparison results")
1103. sn.pairs <- sn.pairs[,!is.na(mi),drop=FALSE]
1104. }
1105. if(verbose) message(" done")
1106. return(invisible(sn.pairs))
1107. }
1108. )
1109. )

conclass.R at commit e656464, under GPL-3.0 · at the source