Overview

The mumosa package implements a variety of utilities for analyses of single-cell data with multiple modalities. This usually refers to single-cell RNA-seq experiments with proteomics, epigenomics or other data collected from the same cells. The aim is to investigate how different modalities relate to each other via analyses of correlations, and to combine data together from multiple modalities for integrated downstream analyses. The actual analyses of each individual modality are deferred to other packages; the scope of mumosa is limited to the sharing of information across modalities.

Setting up the dataset

To demonstrate the functionalities of this package, we will use a subset of a CITE-seq experiment from the scRNAseq package. The main Experiment contains the RNA-seq counts while the adt alternative Experiment contains the CITE-seq counts - see the SingleCellExperiment package documentation for more details.

library(scater)
library(scran)
library(scRNAseq)
sce <- KotliarovPBMCData()
sce <- sce[,1:1000] # subset for speed.
sce

## class: SingleCellExperiment 
## dim: 32738 1000 
## metadata(0):
## assays(1): counts
## rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
## rowData names(0):
## colnames(1000): AAACCTGAGAGCCCAA_H1B1ln1 AAACCTGAGGCGTACA_H1B1ln1 ...
##   ATAACGCGTTTGGGCC_H1B1ln1 ATAACGCTCAACGCTA_H1B1ln1
## colData names(24): nGene nUMI ... dmx_hto_match timepoint
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(1): ADT

We perform some cursory analyses on the RNA component:

stats <- perCellQCMetrics(sce, subsets=list(Mito=grep("^MT-", rownames(sce))))
filter <- quickPerCellQC(stats, sub.fields="subsets_Mito_percent")
sce <- sce[,!filter$discard]
sce <- logNormCounts(sce)
dec <- modelGeneVar(sce)
set.seed(10000)
sce <- runPCA(sce, ncomponents=25, subset_row=getTopHVGs(dec, n=5000))

And again on the protein component.

# TODO: sync this with the book.
sceA <- as(altExp(sce), "SingleCellExperiment")
statsA <- perCellQCMetrics(sceA)
keep <- statsA$detected > 50 
sceA <- sceA[,keep]

library(DropletUtils)
amb <- inferAmbience(assay(sceA))
sceA <- computeMedianFactors(sceA, reference=amb)
sceA <- logNormCounts(sceA)

set.seed(10000)
sceA <- runPCA(sceA, ncomponents=15)

Finally, we put everything back into the same object, only considering the cells passing QC in both modalities.

sce2 <- sce[,keep]
altExp(sce2) <- sceA

Rescaling by neighbors

The easiest way to combine data for the same set of cells is to simply cbind their matrices together prior to downstream analyses like clustering. However, this requires some rescaling to adjust for the differences in the number of features and variation of each modality; otherwise, the modality with more features or stronger (technical) variation would just dominate later calculations. rescaleByNeighbors() quantifies the “noise” of each modality using on the median distance to the k-th nearest neighbor, and then rescales each expression/PC matrix by that distance to equalize the noise across modalities.

library(mumosa)
output <- rescaleByNeighbors(list(reducedDim(sce2), reducedDim(altExp(sce2))))
dim(output)

## [1] 911  40

The result is a cbinded matrix that can be used directly in downstream analyses like clustering and dimensionality reduction.

set.seed(100001)
library(bluster)
sce2$combined.clustering <- clusterRows(output, NNGraphParam())

reducedDim(sce2, "combined") <- output
sce2 <- runTSNE(sce2, dimred="combined")
plotTSNE(sce2, colour_by="combined.clustering")

Advanced users can also adjust the weight given to each modality via the weights= argument; a modality assigned a weight of 2 will mean that the rescaled distance is twice that of a modality with a weight of 1.

Intersecting clusters

In other situations, we may already have a clustering for each modality. For example, we might have one clustering for each of the gene expression and ADT PCs:

g.rna <- buildSNNGraph(sce2, use.dimred="PCA")
rna.clusters <- igraph::cluster_walktrap(g.rna)$membership
table(rna.clusters)

## rna.clusters
##   1   2   3   4   5   6   7 
## 115 218  24  84  92 254 124

g.adt <- buildSNNGraph(altExp(sce2), use.dimred='PCA')
adt.clusters <- igraph::cluster_walktrap(g.adt)$membership
table(adt.clusters)

## adt.clusters
##   1   2   3   4   5   6   7   8 
## 111 191 122 181 122  82  20  82

We can naively “intersect” these clusters by concatenating their identities, as shown below. Cells are only assigned to the same intersected clusters if they belong to the same per-modality clusters. In practice, this often creates a large number of irrelevant clusters with small numbers of cells.

intersected <- paste0(adt.clusters, ".", rna.clusters)
table(intersected)

## intersected
## 1.2 1.4 1.5 2.2 2.5 2.6 3.7 4.2 4.4 4.6 5.1 5.3 5.6 6.2 6.6 7.3 7.7 8.4 8.5 
##  39   2  70 155  20  16 122   6   2 173 115   6   1  18  64  18   2  80   2

A more refined approach is implemented in the intersectClusters() function. This performs the intersection as described but then performs a series of pairwise merges to eliminate the small clusters. It will stop merging when the within-cluster sum of squares exceeds that of the original clustering in any modality; this ensures that any within-modality structure captured by the original clustering is preserved.

intersected2 <- intersectClusters(list(rna.clusters, adt.clusters),
    list(reducedDim(sce2), reducedDim(altExp(sce2))))
table(intersected2)

## intersected2
##   1   2   3   4   5   6   7   8   9 
## 124  92  39  82  84 180 115 171  24

Multi-metric UMAP

Another approach is to take advantage of UMAP’s ability to combine information on different scales. This is achieved by - roughly speaking - creating nearest-neighbor graphs within each modality, and then “intersecting” the graphs such that cells are only considered to be neighbors if they are neighbors in each modality. In this manner, we incorporate locality information from all modalities without explicitly comparing distances. We apply this strategy using the calculateMultiUMAP() function, increasing the number of components to avoid loss of information.

set.seed(100002)
umap.out <- calculateMultiUMAP(list(reducedDim(sce2), reducedDim(altExp(sce2))), 
    n_components=20)
dim(umap.out)

## [1] 911  20

Again, we can use this in downstream analyses like clustering:

library(bluster)
sce2$umap.clustering <- clusterRows(umap.out, NNGraphParam())

And also visualization, though perhaps it is more natural to just compute the UMAP on two dimensions for this purpose. runMultiMap() is just a wrapper around calculateMultiUMAP() that conveniently stores the output in the reducedDims of the input SingleCellExperiment.

set.seed(100002)
sce2 <- runMultiUMAP(sce2, dimred="PCA", extras=list(reducedDim(altExp(sce2))))
plotReducedDim(sce2, "MultiUMAP", colour_by="umap.clustering")

Intersecting graphs

Inspired by the UMAP approach, we can perform intersections on arbitrary graphs with the intersectGraphs() function. This is relevant in cases where we have nearest-neighbor graphs constructed separately for each modality, e.g., for graph-based clustering. We intersect these graphs by setting the weights of edges to the product of the edge weights from each individual graph.

g.com <- intersectGraphs(g.rna, g.adt)

Two cells will only be connected by a high-weighted edge if that edge is of high weight in each original graph, i.e., cells must be nearest neighbors in all modalities to be considered as nearest neighbors in the intersection. We can then use this graph for clustering, which hopefully gives us clusters where cells are similar across all modalities.

com.clusters <- igraph::cluster_walktrap(g.com)$membership
table(com.clusters)

## com.clusters
##   1   2   3   4   5   6   7   8   9  10 
## 123 124 118 169 178  70  21  22  84   2

Correlation analyses

Given two modalities, we may be interested in knowing which features in one modality are correlated to the features in another modality. The computeCorrelations() function does exactly this:

cor.all <- computeCorrelations(sce2, altExp(sce2)[1:5,])
cor.all[order(cor.all$p.value),]

## DataFrame with 163690 rows and 5 columns
##           feature1      feature2       rho     p.value         FDR
##        <character>   <character> <numeric>   <numeric>   <numeric>
## 1            CD79A     BTLA_PROT  0.443223 4.01052e-45 2.89580e-40
## 2         HLA-DRB1    CD123_PROT  0.430457 2.17126e-42 7.83878e-38
## 3          HLA-DRA    CD123_PROT  0.403383 5.75021e-37 1.38398e-32
## 4             FCN1     CD13_PROT  0.397864 6.38599e-36 1.15275e-31
## 5            TCL1A     BTLA_PROT  0.383289 2.97413e-33 4.29495e-29
## ...            ...           ...       ...         ...         ...
## 163686      snoU13 AnnexinV_PROT       NaN          NA          NA
## 163687      snoU13     BTLA_PROT       NaN          NA          NA
## 163688      snoU13    CD117_PROT       NaN          NA          NA
## 163689      snoU13    CD123_PROT       NaN          NA          NA
## 163690      snoU13     CD13_PROT       NaN          NA          NA

For multi-batch experiments, we can specify blocking factors to avoid being confounded by the batch effect. Each batch receives equal weight when averaging correlations and in computing the final combined p-value.

b <- rep(1:3, length.out=ncol(sce2))
cor.all.batch <- computeCorrelations(sce2, altExp(sce2)[1:5,], block=b)

However, this approach can be very slow when dealing with large numbers of features; indeed, one may note the subsetting to the first 5 features in the code above. Users can enable parallelization via BPPARAM= to speed up the process but this will only go so far.

An alternative approach is to, for each feature in one modality, perform an approximate search for the top most-correlated features in another modality. This assumes that only the strongest correlations (positive or negative) are of actual interest, while the bulk of weak correlations can be ignored. We use the findTopCorrelations() function with a specified number of top features to extract:

set.seed(100001) # for IRLBA.
top.cor <- findTopCorrelations(sce2[1:100,], y=altExp(sce2), number=10)
top.cor

## List of length 2
## names(2): positive negative

top.cor$positive

## DataFrame with 530 rows and 5 columns
##         feature1              feature2        rho   p.value       FDR
##      <character>           <character>  <numeric> <numeric> <numeric>
## 1   RP11-34P13.7 RatIgG2bkIsotype_PROT  0.0550420 0.0484271         1
## 2   RP11-34P13.7            CD194_PROT  0.0519342 0.0586247         1
## 3   RP11-34P13.7             CD19_PROT  0.0310772 0.1743962         1
## 4   RP11-34P13.7             CD34_PROT  0.0301584 0.1816169         1
## 5   RP11-34P13.7             CD1c_PROT  0.0300619 0.1823858         1
## ...          ...                   ...        ...       ...       ...
## 526     C1orf222             CD19_PROT 0.03895679  0.120066         1
## 527     C1orf222            CD141_PROT 0.03693348  0.132724         1
## 528     C1orf222             CD10_PROT 0.02722832  0.205866         1
## 529     C1orf222            CD184_PROT 0.00819402  0.402461         1
## 530     C1orf222            CD273_PROT 0.00572678  0.431477         1

This returns the top 10 features in the CITE-seq data (y) for each feature in the main RNA-seq data. As the search is approximate, some features may not be ranked as highly as they would be under an exact search - the approximation quality can be modified by increasing the number of PCs (d) used for data compression. Blocking is also supported with block=, in which case correlations are only computed within each level of the blocking factor.