CAMTHC用户手册
2019-10-29
内容
1简介
混合凸分析(CAM)是一种完全无监督的计算方法，用于分析由未知数量和不同比例的不同亚群组成的组织样本(Wang et al. 2016)．CAM假设测量的表达水平是每个亚群体表达的加权和，其中单个亚群体的贡献与该亚群体的丰度和特定表达成正比。这个线性混合模型可以表述为\ (\ mathbf {X ' = '} \)．CAM可以直接从原始的混合表达矩阵中识别分子标记，\ (\ mathbf {X} \)，进一步估计构成比矩阵，\ (\ mathbf {} \)，以及亚群体特异性表达谱矩阵，\ (\ mathbf{年代}\)．
CAMTHC是利用CAM算法对组织异质性进行表征的R包。它为CAM对混合表达谱进行无监督反褶积提供了基本功能，为理解亚群体特异性结果提供了辅助功能。CAMTHC还实现了基于分子标记S矩阵或A矩阵的先验知识进行监督反褶积的函数。半监督反褶积也可以通过结合来自CAM和先验知识的分子标记来实现，以分析混合物的表达。
2快速启动
这个函数凸轮()包括分解混合表达式配置文件的矩阵的所有必要步骤。上面有一些可选的步骤凸轮()减少矩阵的采样以减少运行时间。每一步凸轮()如果您喜欢更灵活的工作流，也可以单独执行。更多细节将在下面的章节中介绍。
通过以下方式开始分析凸轮()，您需要指定可能的亚种群数量的范围和要去除的低/高表达分子的百分比。一般情况下，可以从基因表达数据中去除30% ~ 50%的低表达基因。由于蛋白质组学数据中蛋白质数量有限，低表达蛋白被去除的较少，如0% ~ 10%。去除高表达分子对结果的影响要小得多，通常设置为0% ~ 10%。
rCAM <- CAM(数据，K = 2:5, thres。低= 0.30，三。高= 0.95)3.无监督反褶积
3．1数据类型和输入格式
从理论上讲,CAMTHC接受任何分子表达数据类型，只要表达式遵循线性混合模型。我们在基因表达数据(microarray, RNAseq)，蛋白质组学数据和DNA甲基化数据中验证了CAM的可行性。对输入表达式数据的要求:
- 在非负数值的非对数线性空间中(即>=0)。
- 已经过归一化和批效应去除处理。
- 没有缺失值;含有缺失值的分子应在CAM之前去除。
- 没有全零表达的分子;否则，全零表达的分子将在内部被移除。
输入表达式数据应该存储在一个矩阵中。数据帧，summarizeexperiment或ExpressionSet对象也被接受，并将在分析之前被内部强制转换为矩阵格式。矩阵中的每一列都应该是一个组织样本。每一行应该是一个探针/基因/蛋白质等。应该提供行名，以便CAM可以返回检测到的标记的名称。否则，行将自动命名为1,2,3，…
3.2凸轮工作流程
我们使用的数据集是从GSE19830以CAM为例，展示CAM工作流程。该数据集提供了纯组织(脑、肝、肺)及其不同比例的生物混合物的基因表达谱。
警告:“CAMTHC”包已弃用，将从Bioconductor #>版本3.12中删除。包重命名为debCAM数据(ratMix3) #ratMix3$X: X矩阵，包含待分析的混合表达谱#ratMix3$A: ground truth A矩阵，包含比例#ratMix3$S: ground truth S矩阵，包含亚种群特异性表达谱数据<- ratMix3$X #10000个基因* 21个组织#满足输入数据要求3.2.1之上CAM分析
该函数可以实现无监督反褶积凸轮()简单的设置为Section2介绍了。其他关键参数包括dim.rdc(降低数据维数)和cluster.num(集群数)。增加它们将带来更多的时间复杂性。我们也可以指定核用于配置的并行计算BiocParallel．核数= 0将禁用并行计算。没有核中的每个元素将调用一个核K．在CAM之前设置随机数生成器的种子，可以生成可重复的结果。
set.seed(111) #设置种子用于内部聚类以生成可重复的结果rCAM <- CAM(data, K = 2:5, thres. set.seed(111) #低= 0.30，三。高= 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.2MDL曲线的模型选择
我们使用广泛采用且一致的信息理论准则MDL来指导模型选择。底层子种群数量可以通过最小化总描述码长度来决定:
情节(MDL (rCAM)数据。term = TRUE)
3.2.3估计A和S矩阵
对于固定的子总体数，如K = 3, CAM估计的A和S矩阵可由
Aest <- Amat(rCAM, 3)本<- Smat(rCAM, 3)3.2.4检测到的标记基因
CAM检测到的用于A矩阵估计的标记基因可以通过
MGlist <- MGsforA(rCAM, K = 3) #为三个亚种群数据预处理过滤了很多基因，其中也有一些生物学上有意义的标记基因。所以我们需要再次检查每个基因以找到所有可能的标记。基于亚群体特异性表达的两个统计数据被利用来识别具有一定阈值的标记基因。第一个是OVE-FC(一个人对所有人-折叠变化)(Yu et al. 2010)。二是自举OVE-FC的下置信界\α(\ \)的水平。
MGstat <- MGstatistic(data, Aest, boot. exe)alpha = 0.05, nboot = 1000) mlist。FC <- lapply(seq_len(3)， function(x) rownames(MGstat)[MGstat$idx == x & MGstat$OVE.]FC > 10]) mlist。FCboot <- lapply(seq_len(3)， function(x) rownames(MGstat)[MGstat$idx == x & MGstat$ overe . fc . FCboot <- lapply(seq_len(3)， function(x) rownames(MGstat)]Alpha > 10])上述阈值是任意设置的，将显著影响结果标记的数量。每个亚种群也可以有不同的阈值。为了使阈值设置更容易，我们还可以将阈值设置为一个亚群体输入标记的折叠变化和误差边际的分位数。误差边际是由A矩阵的列向量定义的一个基因与单纯形之间的距离(Wang et al. 2016)．可以将误差阈值放宽到大于最大值的值。
MGlist。re <- reselectMG(data, MGlist, fc.thres='q0.2'， err.thres='q1.2') #q0.2: 0.2-quantile #q1.2: 1-quantile (maximum) * 1.23.2.5(可选)Re-estimation
可选择基于新的标记列表重新估计A和S矩阵和/或应用交替最小二乘(ALS)方法来进一步减小均方误差。注意，允许过多的ALS迭代可能会带来与初始值显著偏离的风险。甲基化数据的约束，\ (\ mathbf{在[0,1]}\ \)，将在重新估算时施加。
rre <- redoASest(data, MGlist。re, maxIter = 2, methy = FALSE) #rre$Aest:重新估计A矩阵#rre$最强:重新估计S矩阵3.2.6单形图- 2D
CAM成功的基础是混合表达式的散射单形是纯表达式的散射单形的旋转和压缩版本，其中标记基因位于每个顶点。simplexplot ()函数可以在二维图中显示散点单纯形和检测到的标记基因。高维单纯形中的顶点仍然位于低维单纯形的极值点。
layout(matrix(c(1,2,3,4)， 2,2, byrow = TRUE)) simplexplot(data, Aest, MGlist, main=c('初始检测到的标记'))simplexplot(data, Aest, MGlist。FC, main=c(' FC > 10')) simplexplot(data, Aest, MGlist。FCboot, main=c(expression(bold(paste('fc(bootstrap，'，alpha，'=0.05) > 10'))))) simplexplot(data, Aest, MGlist。Re, main=c("fc >= 'q0.2'， error <= 'q1.2'"))
在2D图中显示的颜色和顶点顺序可以更改为
simplexplot(数据，Aest, MGlist。FCboot, data.extra=rbind(t(ratMix3$A)，t(Aest))，角。Order = c(2,1,3)， col = "blue"， mg。坳= c(“红”、“橙色”,“绿色”),ex.col =“黑人”,ex.pch = c(17日,19日,19日,19日,17日,17))传说(“bottomright cex = 1.2,插图=。01, c("Ground Truth"，"CAM Estimation")， pch=c(19,17)， col="black")
3.2.7单形图- 3D
通过主成分分析还可以在三维空间中观察到混合表达式的凸锥和单纯形。注意，PCA不能保证顶点仍然保留为降维单纯形的极值点。
显示凸锥的代码:
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)}显示simplex的代码:
library(rgl) clear3d() Xproj <- XWProj(data, PCAmat(rCAM)) Xp <- Xproj[，-1] plot3d(Xp[， 1:3]， col='gray')， 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.8基于事实的验证
如果地面真理A和S矩阵可用，则可以对CAM的估计进行评估:
cor(Aest, ratMix3$A) #>肝脑肺#> [1，]-0.3963187 0.9839210 -0.4898373 #> [2，]0.9869012 -0.5991347 -0.4865204 #> [3，]-0.4769864 -0.5245692 0.9856413 cor(本，ratMix3$S) #>肝脑肺#> [1，]0.5047310 0.9752287 0.5648873 #> [2，]0.9732946 0.2578350 0.3488423 #> [3，]0.4680238 0.5741921 0.9935207考虑到许多共表达基因(管家基因)的存在可能主导了基本事实与估计表达谱之间的相关系数，因此最好只评估相关系数而不是标记基因。
取消list(lapply(seq_len(3)， function(k) {k.match <-马克斯(软木(东部[k], ratMix3一美元));mgk <- MGlist.FCboot[[k]];corr <- cor(本[mgk, k]， ratMix3$S[mgk, k.match]);names(corr) <- colnames(ratMix3$A)[k.match];corr})) #>脑肝肺#> 0.9982171 0.9950575 0.99845333．3备选CAM工作流
CAM中的主要步骤(数据预处理、标记基因聚类检测和矩阵分解)也可以单独执行，这样可以更加灵活。
set.seed(111) #数据预处理rPrep <- CAMPrep(Data, thres. seed)低= 0.30，三。高= 0.95)#> outlier cluster number: 0 #> convex hull cluster number: 48 #Marker gene cluster detection with a fixed K rMGC <- CAMMGCluster(3, rPrep) #A and S matrix estimation rASest <- CAMASest(rMGC, rPrep, data) #Obtain A and S matrix Aest <- Amat(rASest) Sest <- Smat(rASest) #obtain marker gene list detected by CAM and used for A estimation MGlist <- MGsforA(PrepResult = rPrep, MGResult = rMGC) #obtain a full list of marker genes 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]) MGlist.re <- reselectMG(data, MGlist, fc.thres='q0.2', err.thres='q1.2')3.4可选样本聚类
我们实现了PCA凸轮()/CAMPrep ()减少数据维数。作为另一种数据降维方法，样本聚类是可选的凸轮()/CAMPrep ()．
#集群库(apcluster) apres <- apclusterK(negDistMat(r=2)， t(data)， K = 10) #> Trying p = -2427050 #>集群数量:10 #> Trying p = -15246196 #>集群数量:7 #> Trying p = -8124448 (bisection step no. 2)1) #>集群数量:7个#>尝试p = -4563575(等分步骤no.;2) #>集群数量:8 #>尝试p = -2783138(等分步骤no.;3) #>集群数量:9 #> #>集群数量:9 for p = -2783138 #你也可以使用apcluster()，但需要确保集群数量大于潜在的子种群数量。数据。clusterMean <- lapply(slot(apres，"clusters")， function(x) rowMeans(data[， x, drop = FALSE]))数据。clusterMean <- do。调用(cbind, data. clustermean) set.seed(111) rCAM <-clusterMean, K = 2:5, thres。低= 0.30，三。高= 0.95)# or rPrep <- CAMPrep(data.clusterMean, thres.low = 0.30, thres.high = 0.95)我们仍然可以遵循工作流程3.2或3．3得到标记基因列表和估计的S矩阵。但是，估计的A矩阵是由聚类中心组成的新数据。原始数据的A矩阵可由
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *3．5案例研究
3.5.1GSE11058
GSE11058在微阵列上运行四种免疫细胞系以及它们不同比例的混合物。下面的代码块显示了如何使用CAM将四种混合物盲分离成四种纯细胞系。注意，这个数据集包含54613个探针/探针集。我们可以通过降低探针/探针集采样或去除更多低表达基因(例如70%)来减少运行时间。
#下载数据和表型库(GEOquery) gsm<- getGEO('GSE11058') #>警告:62个解析失败。> 54617 SPOT_ID 1/0/T/F/TRUE/FALSE——控制文字数据#> ............ .................. ......... ............。> 54618 SPOT_ID 1/0/T/F/TRUE/FALSE——控制文字数据#> ............ .................. ......... ............参见problems(…)了解更多细节。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. aggre)。低= 0.70,thres。高= 0.95)Aest <- Amat(rCAM, 4) Aest #> [,1] [,2] [,3] [,4] #> [1,] 0.3883538 0.2389400 0.1039654 0.26874071 #> [2,] 0.2171408 0.4881505 0.2075359 0.08717281 #> [3,] 0.3480010 0.1810097 0.4342633 0.03672609 #> [4,] 0.3795312 0.3204192 0.2752069 0.02484273 #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.2958875 -0.2006075 -0.7497038 0.98630126 #> [2,] -0.1976556 -0.1508022 0.9935143 -0.88675357 #> [3,] -0.7918492 0.9725408 -0.4427615 0.03629895 #> [4,] 0.9981488 -0.8746620 -0.1050113 0.311454183.5.2GSE41826(甲基化)
GSE41826通过将单个个体的神经元和神经胶质提取的DNA从10% ~ 90%混合，进行了重组实验。下面的代码块展示了如何使用CAM将9个混合物盲分离为神经元和胶质特异性CpG甲基化量化。运行CAM对甲基化数据的附加要求:
- 甲基化定量使用β值，而不是m值。
- 如果混合组织来自男性和女性，则必须去除X或Y染色体上的CpG位点。
- CpG站点需要降低采样，因为大量的CpG站点会减慢聚类过程。
# getGEO('GSE41826') mixtureId <- unlist(lapply(paste0('Mix'，seq_len(9)))，函数(x) gsm[[1]]美元geo_accession (gsm [[1]] $ title = = x]))数据< - gsm [[1]] [, mixtureId] #删除CpG网站性染色体从雄性和雌性#如果组织没有必要在这个例子中是混合物从同一话题# gpl < - getGEO (GPL13534) # annot < datatable (gpl) @ table (c(“名字”,“科”)]# rownames (annot) < - annot名字#美元annot < - annot [rownames(数据)]# < -数据(annot $空空! = ' x ' & annot $空空! = Y) # downsample CpG网站featureId < -样本(seq_len (nrow(数据),50000) #运行CAM #当输入数据探针太多时，由于空间限制，lof必须关闭。#当集群。Num很大，很快。rCAM <- CAM(data[featureId，]， K = 2, thres. select可以使速度增加rCAM <- CAM(data[featureId，]， K = 2, thres.)低= 0.10,thres。高= 0.60，集群。num = 100, MG.num.thres = 10, lof。Thres = 0，快。Select = 20)MGlist <- MGsforA(rCAM, K = 2) #Identify markers from all CpG sites MGlist.re <- reselectMG(data, MGlist, fc.thres='q0.2', err.thres='q1.2') #re-esitmation with methylation constraint rre <- redoASest(data, MGlist.re, maxIter = 20, methy = TRUE) #Validation using ground truth A matrix Atrue <- cbind(seq(0.1, 0.9, 0.1), seq(0.9, 0.1, -0.1)) cor(rre$Aest, Atrue)我们也可以在未甲基化的定量上运行CAM, 1 - beta，以获得非常相似的结果，因为潜在的线性混合模型也适用于未甲基化的探针强度。
rCAM <- CAM(1 - exprs(data[featureId，])， K = 2, thres。低= 0.10,thres。高= 0.60，集群。num = 100, MG.num.thres = 10, lof。Thres = 0，快。Select = 20)4监督/ Semi-supervised反褶积
4．1分子标记已知
CAM算法可以基于盲检标记估计A和S矩阵。因此，我们也可以使用部分CAM算法，基于已知标记估计A和S矩阵。
该包提供AfromMarkers ()从标记表估计矩阵。
Aest <- AfromMarkers(data, MGlist) #MGlist是一个向量列表，每个向量包含一个子种群的已知标记或者我们可以用redoASest ()从标记列表中估计A和S矩阵，或者用ALS重新估计两个矩阵以进一步减小均方误差。注意，允许过多的ALS迭代可能会带来与初始值显著偏离的风险。甲基化数据的约束，\ (\ mathbf{在[0,1]}\ \)，可在重新估算时施加。
rre <- redoASest(data, MGlist, maxIter = 10) #MGlist是一个向量列表，每个向量包含一个子种群的已知标记#maxIter = 0: ALS没有重新估计#rre$Aest:估计a矩阵#rre$ gest:估计S矩阵4.2S矩阵已知
许多数据集提供了纯化细胞系甚至每一个细胞的表达谱，可以作为S矩阵的参考。一些方法使用最小二乘技术或支持向量回归(Newman et al. 2015)根据已知的S矩阵估计A矩阵。CAMTHC将从已知的S矩阵中通过识别的标记估计出A矩阵，在与ground truth A矩阵的相关系数方面具有更好的性能。
data <- ratMix3$X S <- ratMix3$ s#已知S矩阵#识别标记pMGstat <- MGstatistic(S, c("肝"，"脑"，"肺"))pMGlist。FC <- lapply(c("Liver"，"Brain"，"Lung")， function(x) rownames(pMGstat)[pMGstat$idx == x & pMGstat$OVE.]FC > 10]) #估计矩阵Aest <- AfromMarkers(data, pMGlist.FC) #(可选)替代重估计rre <- redoASest(data, pMGlist.FC) #(可选)FC, maxIter = 10)CAMTHC还支持用最小二乘法直接从已知的S矩阵估计A矩阵。仍然需要从S矩阵中识别标记，因为只计算标记的平方误差，这有利于(1)更快的计算运行时间;(2)由于标记的鉴别能力，更大的信噪比(Newman et al. 2015)．
Aest <- redoASest(数据，pMGlist。FC, S=S, maxIter = 0)$Aest4.３矩阵是已知的
对于已知的A矩阵，CAMTHC用非负最小二乘(NNLS)估计S矩阵NMF，并进一步识别标记。
data <- ratMix3$X A <- ratMix3$A #已知A矩阵#估计S矩阵本<- t(NMF::。fcnnls(A, t(data))$coef) #标识标记pMGstat <- MGstatistic(data, A) pMGlist. fcnnls(A, t(data))FC <- lapply(unique(pMGstat$idx)， function(x) rownames(pMGstat)[pMGstat$idx == x & pMGstat$OVE.]FC > 10]) #(可选)可选重估计rre <- redoASest(数据，pMGlist。FC, A=A, maxIter = 10)4.4Semi-supervised反褶积
当已知标记、S矩阵和/或A矩阵的先验信息时，还可以将先验信息中的标记与CAM识别的标记相结合，进行半监督反褶积。监督反褶积不能在没有先验信息的情况下处理底层子种群，而无监督反褶积可能会在没有足够辨别能力的情况下错过子种群。因此，半监督反褶积可以利用这两种方法。
Aest <- AfromMarkers(data, MGlist) #MGlist是一个向量列表，每个向量包含一个子种群的已知标记和/或cam检测到的标记或
rre <- redoASest(data, MGlist, maxIter = 10) #MGlist是一个向量列表，每个向量包含一个子种群的已知标记和/或cam检测到的标记#maxIter = 0: ALS没有重新估计#rre$Aest:估计的a矩阵#rre$ gest:估计的S矩阵参考文献
纽曼，Aaron M.，刘志龙，Michael R. Green, Andrew J. Gentles，冯卫国，徐岳，Chuong D. huang, Maximilian Diehn, Ash A. Alizadeh. 2015。组织表达谱中细胞亚群的稳健枚举Nat冰毒12(5): 453 - 57。http://dx.doi.org/10.1038/nmeth.3337．
王妮娅，埃里克·p·霍夫曼，陈露露，陈莉，张震，刘春雨，于国强，大卫·m·赫林顿，罗伯特·克拉克，王岳。2016。转录异质性的数学建模识别复杂组织中的新标记和亚群科学报告6:18909。http://dx.doi.org/10.1038/srep18909．