3.2.1 Analysis by CAM

Unsupervised deconvolution can be achieved by the function CAM() with a simple setting as Section 2 introduced. Other critical parameters are dim.rdc (reduced data dimension) and cluster.num (number of clusters). Increasing them will bring more time complexity. We can also specify cores for parallel computing configured by BiocParallel. cores = 0 will disable parallel computing. No cores argument will invoke one core for each element in K. Setting the seed of the random number generator prior to CAM can generate reproducible results.

set.seed(111) # set seed for internal clustering to generate reproducible results
rCAM <- CAM(data, K = 2:5, thres.low = 0.30, thres.high = 0.95)
#CAM return three sub results: 
#rCAM@PrepResult contains details corresponding to data preprocessing.
#rCAM@MGResult contains details corresponding to marker gene clusters detection.
#rCAM@ASestResult contains details corresponding to A and S matrix estimation.

3.2.2 Estimated A and S matrix

The A and S matrix estimated by CAM with a fixed subpopulation number, e.g. K = 3, can be obtained by

Aest <- Amat(rCAM, 3)
Sest <- Smat(rCAM, 3)

3.2.3 Detected Marker Genes

The marker genes detected by CAM and used for A matrix estimation can be obtained by

MGlist <- MGsforA(rCAM, K = 3) #for three subpopulations

Data preprocessing has filtered many genes, among which there are also some biologically-meaningful marker genes. So we need to check each gene again to find all the possible markers. Two statistics based on the subpopulation-specific expressions are exploited to identify marker genes with certain thresholds. The first is OVE-FC (one versus everyone - fold change) (Yu et al. 2010). The second is the lower confidence bound of bootstrapped OVE-FC at \(\alpha\) level.

MGstat <- MGstatistic(data, Aest, boot.alpha = 0.05, nboot = 1000)
MGlist.FC <- lapply(seq_len(3), function(x) 
    rownames(MGstat)[MGstat$idx == x & MGstat$OVE.FC > 10])
MGlist.FCboot <- lapply(seq_len(3), function(x) 
    rownames(MGstat)[MGstat$idx == x & MGstat$OVE.FC.alpha > 10])

3.2.4 Simplex Plot - 2D

Fundamental to the success of CAM is that the scatter simplex of mixed expressions is a rotated and compressed version of the scatter simplex of pure expressions, where the marker genes are located at each vertex. simplexplot() function can show the scatter simplex and detected marker genes in a 2D plot. The vertices in the high-dimensional simplex will still locate at extreme points of the low-dimensional simplex.

simplexplot(data, Aest, MGlist)

simplexplot(data, Aest, MGlist.FC)

simplexplot(data, Aest, MGlist.FCboot)

The colors and the vertex order displayed in 2D plot can be changed as

simplexplot(data, Aest, MGlist.FCboot, corner.order = c(2,1,3), 
            data.col = "blue", corner.col = c("red","orange","green"))

3.2.5 Simplex Plot - 3D

We can also observe the convex cone and simplex of mixture expressions in 3D space by PCA. Note that PCA cannot guarantee that vertices are still preserved as extreme points of the dimension-reduced simplex.

Code to show convex cone:

library(rgl)
Xp <- data %*% t(PCAmat(rCAM))
plot3d(Xp[, 1:3], col='gray', size=3, 
       xlim=range(Xp[,1]), ylim=range(Xp[,2:3]), zlim=range(Xp[,2:3]))
abclines3d(0,0,0, a=diag(3), col="black")
for(i in seq_along(MGlist)){
    points3d(Xp[MGlist[[i]], 1:3], col= rainbow(3)[i], size = 8)
}

Code to show simplex:

library(rgl)
clear3d()
Xproj <- XWProj(data, PCAmat(rCAM))
Xp <- Xproj[,-1]
plot3d(Xp[, 1:3], col='gray', size=3, 
       xlim=range(Xp[,1:3]), ylim=range(Xp[,1:3]), zlim=range(Xp[,1:3]))
abclines3d(0,0,0, a=diag(3), col="black")
for(i in seq_along(MGlist)){
    points3d(Xp[MGlist[[i]], 1:3], col= rainbow(3)[i], size = 8)
}

3.2.6 Model Selection by MDL Curves

We used MDL, a widely-adopted and consistent information theoretic criterion, to guide model selection. The underlying subpopulation number can be decided by minimizing the total description code length:

plot(MDL(rCAM), data.term = TRUE)

3.2.7 Validation with Ground Truth

If the ground truth A and S matrix are available, the estimation from CAM can be evaluated:

cor(Aest, ratMix3$A)
#>           Liver      Brain       Lung
#> [1,] -0.4727518 -0.5284658  0.9847760
#> [2,]  0.9871181 -0.5982854 -0.4875316
#> [3,] -0.3970947  0.9840614 -0.4891499
cor(Sest, ratMix3$S)
#>          Liver     Brain      Lung
#> [1,] 0.4663993 0.5756277 0.9933711
#> [2,] 0.9720974 0.2562309 0.3423406
#> [3,] 0.5077558 0.9748317 0.5664304

Considering the presence of many co-expressed genes (housekeeping genes) may dominate the correlation coefficients between ground truth and estimated expression profiles, it is better to assess correlation coefficients over marker genes only.

unlist(lapply(seq_len(3), function(k) {
    k.match <- which.max(cor(Aest[,k], ratMix3$A));
    mgk <- MGlist.FCboot[[k]];
    corr <- cor(Sest[mgk, k], ratMix3$S[mgk, k.match]);
    names(corr) <- colnames(ratMix3$A)[k.match];
    corr}))
#>      Lung     Liver     Brain 
#> 0.9985593 0.9951157 0.9982504

3.5.1 GSE11058

GSE11058 run four immune cell lines on microarrays as well as their mixtures of various relative proportions. The code chunk below shows how to use CAM to blindly separate the four mixtures into four pure cell lines. Note that this data set contains 54613 probe/probesets. We can reduce running time by downsampling probe/probesets or remove more low-expressed genes (e.g. 70%).

#download data and phenotypes
library(GEOquery)
gsm<- getGEO('GSE11058')
pheno <- pData(phenoData(gsm[[1]]))$characteristics_ch1
mat <- exprs(gsm[[1]])
mat <- mat[-grep("^AFFX", rownames(mat)),]
mat.aggre <- sapply(unique(pheno), function(x) rowMeans(mat[,pheno == x]))
data <- mat.aggre[,5:8]

#running CAM
set.seed(111)
rCAM <- CAM(data, K = 4, thres.low = 0.70, thres.high = 0.95)
Aest <- Amat(rCAM, 4)
Aest
#>           [,1]      [,2]      [,3]       [,4]
#> [1,] 0.4028161 0.2418041 0.1020896 0.25329022
#> [2,] 0.2213951 0.4948602 0.2021001 0.08164461
#> [3,] 0.3505424 0.1831931 0.4321002 0.03416430
#> [4,] 0.3780136 0.3249908 0.2740140 0.02298166

#Use ground truth A to validate CAM-estimated A matrix
Atrue <- matrix(c(2.50, 0.50, 0.10, 0.02,
              1.25, 3.17, 4.95, 3.33,
              2.50, 4.75, 1.65, 3.33,
              3.75, 1.58, 3.30, 3.33), 4,4,
              dimnames = list(c("MixA", "MixB", "MixC","MixD"),
                               c("Jurkat", "IM-9", "Raji", "THP-1")))
Atrue <- Atrue / rowSums(Atrue)
Atrue
#>           Jurkat      IM-9      Raji     THP-1
#> MixA 0.250000000 0.1250000 0.2500000 0.3750000
#> MixB 0.050000000 0.3170000 0.4750000 0.1580000
#> MixC 0.010000000 0.4950000 0.1650000 0.3300000
#> MixD 0.001998002 0.3326673 0.3326673 0.3326673
cor(Aest, Atrue)
#>          Jurkat       IM-9       Raji       THP-1
#> [1,]  0.3725502 -0.2613393 -0.7516873  0.99172438
#> [2,] -0.1987591 -0.1501336  0.9935634 -0.88653613
#> [3,] -0.7885817  0.9696095 -0.4515971  0.04908184
#> [4,]  0.9982001 -0.8746573 -0.1056444  0.31225365

3.5.2 GSE41826 (methylation)

GSE41826 conducted a reconstitution experiment by mixing neuron and glial derived DNA from a single individual from 10% ~ 90%. The code chunk below shows how to use CAM to blindly separate 9 mixtures into neuron and glial specific CpG methylation quantification. Additional requirements of running CAM on methylation data:

• Use beta-value, not M-value, as methylation quantification.
• CpG sites in X or Y chromosomes must be removed if mixture tissues are from both males and females.
• CpG sites needs to be downsampled since the huge number of CpG sites will slow the clustering process.

#download data
library(GEOquery)
gsm <- getGEO('GSE41826')
mixtureId <- unlist(lapply(paste0('Mix',seq_len(9)), 
    function(x) gsm[[1]]$geo_accession[gsm[[1]]$title==x]))
data <- gsm[[1]][,mixtureId]

#Remove CpG sites in sex chromosomes if tissues are from both males and females
#Not necessary in this example as mixtures are from the same subject
#gpl<- getGEO('GPL13534')
#annot<-dataTable(gpl)@table[,c('Name','CHR')]
#rownames(annot) <- annot$Name
#annot <- annot[rownames(data),]
#data <- data[annot$CHR != 'X' & annot$CHR != 'Y',]

#downsample CpG sites
featureId <- sample(seq_len(nrow(gsm[[1]])), 20000)

#running CAM
rCAM <- CAM(data[featureId,], K = 2, thres.low = 0.10, thres.high = 0.60)
Aest <- Amat(rCAM, 2)
Atrue <- cbind(seq(0.1, 0.9, 0.1), seq(0.9, 0.1, -0.1))
cor(Aest, Atrue)

#obtain a full list of marker CpG sites
MGstat <- MGstatistic(data, Aest)
MGlist.FC <- lapply(seq_len(2), function(x) 
    rownames(MGstat)[MGstat$idx == x & MGstat$OVE.FC > 10])

We can also run CAM on unmethylated quantification, 1 - beta, to obtain quite similar results since the underlying linear mixing model is also applicable for unmethylated probe intensity.

rCAM <- CAM(1 - exprs(data[featureId,]),
            K = 2, thres.low = 0.10, thres.high = 0.60)