ADAPTS Vignette #3

Building a Better Signature Matrix From Single Cell Data

Install ADAPTSdata4 using the code:

install.packages(‘devtools’)

library(devtools)

devtools::install_github(‘yxinyi2/ADAPTSdata4’)

library(ADAPTS)
library(pheatmap)
library(parallel)
library(ranger)
library(ADAPTSdata4)
doParallel::registerDoParallel(cores = parallel::detectCores())
set.seed(42)

Step 1: Build signature matrices from the scHenry data

This prostate data set comes from Henry, Gervaise H., et al. “A cellular anatomy of the normal adult human prostate and prostatic urethra.” Cell reports 25.12 (2018): 3530-3542.

scHenry<-log(ADAPTSdata4::scHenry+1)

Step 1a: Randomly split the dataset into training and test set

Around half of the cell types fall into the training set and the other half fall into the test set. scSample function helps collapse cell types with similar cell type IDs into n (groupsize) groups.

trainTestSet<-ADAPTS::splitSCdata(scHenry,numSets = 2,randomize = TRUE)

## Basal Epithelia : 18439; 9220, 9219

## Other : 7163; 3582, 3581

## Luminal Epithelia : 2238; 1119, 1119

## Fibroblast : 1168; 584, 584

## Hillock Epithelia : 1312; 656, 656

## Smooth Muscle : 945; 473, 472

## Endothelia : 1586; 793, 793

## Leukocyte : 459; 230, 229

## Club Epithelia : 2530; 1265, 1265

## Neuroenodcrine Epithelia : 25; 13, 12

trainSet<-as.matrix(trainTestSet[[1]])
testSet<-as.matrix(trainTestSet[[2]])
 
trainSet.30sam <- ADAPTS::scSample(RNAcounts = trainSet, groupSize = 30, randomize = TRUE)
trainSet.3sam <- ADAPTS::scSample(RNAcounts = trainSet, groupSize = 3, randomize = TRUE)

Step 1b: Build a deconvolution seed matrix using ranger forest and estimate the accuracy on pseudo bulk test set

ADAPTS provides the option of building a new seed matrix de novo based on the sample given, in addition to augmenting existing signature matrices, such as LM22. This is particularly helpful for single cell data sets, where the cell types present have come from their native tissue.

seedMat<-ADAPTS::buildSeed(trainSet=trainSet, trainSet.3sam = trainSet.3sam, trainSet.30sam = trainSet.30sam, genesInSeed = 100, groupSize = 30, randomize = TRUE, num.trees = 1000, plotIt = TRUE)

pseudobulk.test <- data.frame(test=rowSums(testSet))
pseudobulk.test.counts<-table(colnames(testSet))
actFrac.test <- 100 * pseudobulk.test.counts / sum(pseudobulk.test.counts)
 
estimates.test <- as.data.frame(ADAPTS::estCellPercent.DCQ(seedMat, pseudobulk.test))
colnames(estimates.test)<-'seed'
estimates.test$actFrac<-round(actFrac.test[rownames(estimates.test)],2)
 
seedAcc<-ADAPTS::calcAcc(estimates=estimates.test[,1], reference=estimates.test[,2])

Step 1c: Build a deconvolution matrix using all the genes and estimate the accuracy on pseudo bulk test set

This step tests if building signature matrix is really necessary by comparing the performance of signature matrices and all-gene matrix.

allGeneSig <- apply(trainSet.3sam, 1, function(x){tapply(x, colnames(trainSet.3sam), mean, na.rm=TRUE)})
 
estimates.allGene <- as.data.frame(ADAPTS::estCellPercent.DCQ(t(allGeneSig), pseudobulk.test))
colnames(estimates.allGene)<-'all'
estimates.test<-cbind(estimates.allGene,estimates.test)

allAcc<-ADAPTS::calcAcc(estimates=estimates.test[,1], reference=estimates.test[,3])

Step 1d: Augment the seed matrix and estimate the accuray on pseudo bulk test set

ADAPTS takes in the seed matrix, adds one additional gene from the full data at a time and records their condition number. The new augmented signature matrix is chosen based on the lowest condition number.

gList <- ADAPTS::gListFromRF(trainSet=trainSet.30sam)
augTrain <- ADAPTS::AugmentSigMatrix(origMatrix = seedMat, fullData = trainSet.3sam, gList = gList, nGenes = 1:100, newData = trainSet.3sam, plotToPDF = FALSE, pdfDir = '.')

estimates.augment <- as.data.frame(ADAPTS::estCellPercent.DCQ(augTrain, pseudobulk.test))
colnames(estimates.augment) <- 'aug'
estimates.test <- cbind(estimates.augment, estimates.test)

augAcc<-ADAPTS::calcAcc(estimates=estimates.test[,1], reference=estimates.test[,4])

Step 1e: Shrink the augmented matrix to minimize the condition number and estimate the accuray on pseudo bulk test set

This step iteratively removes the genes that would most improve the condition number by their absence. The condition number will usually decrease before increasing again when the number of genes becomes small, as shown in the plot below. Note that this plot should be read from right to left.

ADAPTS provides options to automatically select different points along the condition number curve corresponding to different numbers of genes. Typically, the best classifier has the largest number of genes where the condition number is at the bottom of the bowl shown in the plot. The fastStop option speeds up this process by halting execution when the algorithm has detected (perhaps incorrectly) that it has reached this point.

augTrain.shrink <- ADAPTS::shrinkSigMatrix(augTrain, numChunks=NULL, verbose=FALSE, plotIt=TRUE, sigGenesList=NULL, singleCore=FALSE, fastStop=FALSE)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

estimates.shrink <- as.data.frame(ADAPTS::estCellPercent.DCQ(augTrain.shrink, pseudobulk.test))
colnames(estimates.shrink)<-'shrink'
estimates.test <- cbind(estimates.shrink, estimates.test)

titleStr <- paste('Shrunk Signature Matrix,', '# Genes:',nrow(augTrain),'->',nrow(augTrain.shrink))
pheatmap(augTrain.shrink, main=titleStr)

augshrunkAcc<-ADAPTS::calcAcc(estimates=estimates.test[,1], reference=estimates.test[,5])

Step 2: Combine indistinguishible cell types for meta analysis

Cell-types for single cells are typically assigned by using unsupervised clustering techniques to group cells that are putatively all of the same cell type and then assigning cells in each group based on the expression of marker genes. While powerful, this approach is limited in that the delineations between clusters are often wrong. In other word, cell identified as different cell types by clustering may not really be different cell types or they may not really be detectable as different cell types based on gene expression data.

ADAPTS detects when cell clusters ought to be combined and this typically results in a much more accurate signature matrix.

This step shows that combining cell types that are highly correlated improves the ability to accurately deconvolve. The clusters are built based on n-pass spillover matrix. ADAPTS iteratively applying the spillover calculation until the results converge into clusters of cell types. See more details in Algorithm 2 of ADAPTS paper at https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0224693. All the signature matrices are built as in Step 1.

varClusts <- ADAPTS::clustWspillOver(sigMatrix = seedMat, geneExpr = trainSet.3sam)
metaCluster.id <- list()
  for(i in 1:length(varClusts$allClusters)) {
    for (x in varClusts$allClusters[[i]]) {
      metaCluster.id[[x]] <- paste('Meta',i,sep='_')
    }
  }
 
metaClust.LUT <-  unlist(metaCluster.id) 
 
metatrainSet<-trainSet
metatestSet<-testSet
 
colnames(metatrainSet) <- metaClust.LUT[colnames(metatrainSet)]
colnames(metatestSet) <- metaClust.LUT[colnames(metatestSet)]
 
metatrainSet.3sam <- ADAPTS::scSample(RNAcounts = metatrainSet, groupSize = 3, randomize = TRUE)
metatrainSet.30sam <- ADAPTS::scSample(RNAcounts = metatrainSet, groupSize = 30, randomize = TRUE)
 
metapseudobulk.test <- data.frame(test=rowSums(metatestSet))
metapseudobulk.test.counts <- table(colnames(metatestSet))
metaactFrac.test <- 100 * metapseudobulk.test.counts / sum(metapseudobulk.test.counts)

Step 2a: Build a deconvolution seed matrix for clustered cell types

metaseedMat <- ADAPTS::buildSeed(trainSet=metatrainSet, trainSet.3sam = metatrainSet.3sam, trainSet.30sam = metatrainSet.30sam, genesInSeed = 100, groupSize = 30, randomize = TRUE, num.trees = 1000, plotIt = TRUE)

estimates.Meta.onTest <- as.data.frame(ADAPTS::estCellPercent.DCQ(metaseedMat, metapseudobulk.test))
colnames(estimates.Meta.onTest)<-'meta.seed'
estimates.Meta.onTest$metaactFrac.test <- round(metaactFrac.test[rownames(estimates.Meta.onTest)],2)

metaseedAcc<-ADAPTS::calcAcc(estimates=estimates.Meta.onTest[,1], reference=estimates.Meta.onTest[,2])

Step 2b: Build a deconvolution matrix using all the genes for clustered cell types

metaallGeneSig <- apply(metatrainSet.3sam, 1, function(x){tapply(x, colnames(metatrainSet.3sam), mean, na.rm=TRUE)})

metaestimates.allGene <- as.data.frame(ADAPTS::estCellPercent.DCQ(t(metaallGeneSig), metapseudobulk.test))
colnames(metaestimates.allGene)<-'meta.all'
estimates.Meta.onTest <- cbind(metaestimates.allGene, estimates.Meta.onTest)

metaallAcc<-ADAPTS::calcAcc(estimates=estimates.Meta.onTest[,1], reference=estimates.Meta.onTest[,3])

Step 2c: Augment the seed matrix for clustered cell types

metagList <- ADAPTS::gListFromRF(trainSet=metatrainSet.30sam)
 
metaaugTrain <- ADAPTS::AugmentSigMatrix(origMatrix = metaseedMat, fullData = metatrainSet.3sam, gList = metagList, nGenes = 1:100, newData = metatrainSet.3sam, plotToPDF = FALSE, pdfDir = '.')

estimates.Meta.augment <- as.data.frame(ADAPTS::estCellPercent.DCQ(metaaugTrain, metapseudobulk.test))
colnames(estimates.Meta.augment) <- paste('meta.aug')
estimates.Meta.onTest <- cbind(estimates.Meta.augment, estimates.Meta.onTest)

metaaugAcc<-ADAPTS::calcAcc(estimates=estimates.Meta.onTest[,1], reference=estimates.Meta.onTest[,4])

Step 2d: Shrink the augmented matrix for clustered cell types

metaaugTrain.shrink <- ADAPTS::shrinkSigMatrix(sigMatrix=metaaugTrain, numChunks=50, verbose=FALSE, plotIt = TRUE, aggressiveMin=TRUE, fastStop=FALSE, singleCore=TRUE)

estimates.Meta.shrink <- as.data.frame(ADAPTS::estCellPercent.DCQ(metaaugTrain.shrink, metapseudobulk.test))
colnames(estimates.Meta.shrink) <- paste('meta.shrink')
estimates.Meta.onTest <- cbind(estimates.Meta.shrink, estimates.Meta.onTest)

metatitleStr <- paste('Shrunk Signature Matrix,', '# Genes:',nrow(metaaugTrain),'->',nrow(metaaugTrain.shrink))
pheatmap(metaaugTrain.shrink, main=metatitleStr)

metashrunkAcc <-ADAPTS::calcAcc(estimates=estimates.Meta.onTest[,1], reference=estimates.Meta.onTest[,5])

Final result

The most accurate signature matrix (as measured by Pearson and Spearman’s correlation) is usually generated by meta-clustering with ‘clustWspillOver’, building a seed matrix with ‘buildSeed’, augmenting it with ‘gListFromRF’ and ‘AugmentSigMatrix’, and then removing excess genes with ‘shrinkSigMatrix’. Compare the rho.cor and spear.cor for ‘meta.shrink’ to all other methods, but note that overall RMSE tends to increase due to combing the clusters.

The result varies in different runs due to the randomization in ranger forest. Codes for automatically running the entire simulation multiple times will be included in the next version of ADAPTS.

ADAPTS achieves impressive correlations on this more advanced single cell data set.

acc<-cbind(augshrunkAcc,augAcc,allAcc,seedAcc)
acc<-acc[c(1,3,6), ]
acc<-cbind(acc,rep(1,3))
colnames(acc)<-c('shrink','aug','all','seed','actFrac')
deconTable<-round(rbind(estimates.test,acc),2)
 
metacc<-cbind(metashrunkAcc,metaaugAcc,metaallAcc,metaseedAcc)
metacc<-metacc[c(1,3,6), ]
metacc<-cbind(metacc,rep(1,3))
colnames(metacc)<-c('meta.shrink','meta.aug','meta.all','meta.seed','metaactFrac.test')
metadeconTable<-round(rbind(estimates.Meta.onTest,metacc),2)
 
print(deconTable)

##                          shrink   aug   all  seed actFrac
## Basal Epithelia           20.56 13.90 11.71 18.38   51.42
## Club Epithelia            13.54 10.25  8.93 12.63    7.06
## Endothelia                 6.03  9.15  9.59  6.45    4.42
## Fibroblast                 2.03  6.68  9.25  0.76    3.26
## Hillock Epithelia         12.06  9.33  9.36  8.99    3.66
## Leukocyte                  1.25  7.62 10.13  0.00    1.28
## Luminal Epithelia         10.06 10.67 10.23 16.72    6.24
## Neuroenodcrine Epithelia  14.19 12.56  9.99 17.83    0.07
## Other                     17.02 12.23 10.77 18.24   19.97
## Smooth Muscle              3.26  7.61 10.04  0.00    2.63
## others                     0.00  0.00  0.00  0.00      NA
## rho.cor                    0.69  0.67  0.80  0.52    1.00
## spear.rho                  0.61  0.49  0.32  0.63    1.00
## rmse                      11.37 13.36 14.17 12.64    1.00

print(metadeconTable)

##           meta.shrink meta.aug meta.all meta.seed metaactFrac.test
## Meta_1          68.10    28.42    25.95     27.22            78.44
## Meta_2           0.00    23.40    24.51     24.16             5.70
## Meta_3           0.00    21.08    25.06     21.92             5.89
## Meta_4          31.90    27.10    24.49     26.70             9.97
## others           0.00     0.00     0.00      0.00               NA
## rho.cor          0.91     0.71     0.91      0.64             1.00
## spear.rho        0.95     0.80     0.40      0.80             1.00
## rmse            12.80    28.89    30.36     29.59             1.00