前五期的肿瘤可进化分析都是基于组织或血液检出的SNV和CNV，而这几年兴起的单细胞，越来越展露头角，自然，单细胞克隆进化也是目前的热点系列，接下来我们将介绍几款分析单细胞克隆进化的软件，也可以说是单细胞克隆进化的pepiline之一，即：
Canopy + Cardelino +RobustClone
Canopy：通过NGS测序评估肿瘤内异质性和追踪纵向和空间克隆进化历史。
癌症是一种由体细胞遗传和表观遗传改变的进化选择驱动的疾病。在这里，我们提出Canopy，这是一种通过一个或多个来自单个患者的样本的somatic variants 和表观变异来推断肿瘤进化系统发育的方法。该软件应用于纵向和空间实验设计的批量测序数据集，以及来自人类癌细胞MDA-MB-231的可移植转移模型。Canopy成功地识别了细胞种群，并推断出与现有知识和地面真相一致的系统发育。通过数值模拟研究，探讨了关键参数对反卷积精度的影响，并与现有方法进行了比较。
Canopy是一个开源的R包，https://cran.r-project.org/web/packages/Canopy 。
关于安装问题
R包的安装基本就几种方法：
CRAN直接安装install.packages()，
BiocManager::install()
devtools::install_github()，
其他直接安装的模式。
该软件的安装如下：
install.packages('Canopy')#or:install.packages(c("ape", "fields", "pheatmap", "scatterplot3d", "devtools"))devtools::install_github("yuchaojiang/Canopy/package")
关于输入文件
Canopy的输入是肿瘤和正常样本配对检出的SNA and CNA，该软件包给出来 例子是一个来自异质人乳腺癌细胞系MDA-MB-231的可移植转移模型系统重建肿瘤系统发育的例子。来自亲本系MDA-MB-231的癌细胞被移植到小鼠宿主体内，导致器官特异性转移。混合细胞群(MCPs)在体内选择从骨或肺转移，并生长成表型稳定和转移能力的癌细胞系。亲本系和MCP子系全外显子组测序，并对体细胞SNAs和CNAs进行分析。Canopy用于推断转移系统发育。
SNV是使用常规的软件GATK，注释使用ANNOVAR，当配对正常样本可用时，也可以使用MuTect和VarScan2；CNA输入是通过提炼并结合TITAN的等位基因特异性分割结果得到的。而该软件包输入在该软件包并未解决这个棘手的问题，需要我们自行挑选SNAs或CNAs在不同样本之间表现出不同的模式(来自同一患者，因为我们在观察肿瘤内的异质性)。对于SNV，这意味着观测到的VAFs是不同的,热图是一种很好的可视化方法。对于CNA，这意味着WM和Wm是不同的，IGV是一个很好的可视化工具，并建议关注大型CNA区域，这有助于消除错误调用和加快计算。该软件包自带6个测试集，但是文章补充文件中给出5个表格作为输入，我们也只能利用例子中的数据MDA231完成整个分析流程的代码，并复现一下文章中Figure 4 的内容，数据读取以及数据格式，如下：
library(Canopy)#data('MDA231_tree')#data(AML43)#data(toy)#data(toy2)#data(toy3)data("MDA231")projectname = MDA231$projectname ## name of projectR = MDA231$R; R ## mutant allele read depth (for SNAs)MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lungBRAF 155 59 136 77 49KRAS 44 21 54 19 17ALPK2 37 17 28 10 7RYR1 44 0 26 0 0X = MDA231$X; X ## total depth (for SNAs)MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lungBRAF 157 111 177 146 71KRAS 44 30 64 42 27ALPK2 63 17 65 24 7RYR1 107 56 165 55 43WM = MDA231$WM; WM ## observed major copy number (for CNA regions)MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lungchr7 2.998 2.002 2.603 2.000 2.001chr12 1.998 1.998 1.603 1.001 1.999chr18 1.000 2.992 1.000 1.002 2.996chr19 2.000 2.000 2.000 2.000 2.000Wm = MDA231$Wm; Wm ## observed minor copy number (for CNA regions)MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lungchr7 0.002 0.998 0.397 1.000 0.999chr12 0.002 0.998 0.397 1.000 0.999chr18 1.000 0.004 1.000 0.999 0.002chr19 1.000 1.000 1.000 1.000 1.000epsilonM = MDA231$epsilonM ## standard deviation of WM, pre-fixed hereepsilonm = MDA231$epsilonm ## standard deviation of Wm, pre-fixed here## Matrix C specifices whether CNA regions harbor specific CNAs## only needed if overlapping CNAs are observed, specifying which CNAs overlapC = MDA231$C; Cchr7_1 chr7_2 chr12_1 chr12_2 chr18 chr19chr7 1 1 0 0 0 0chr12 0 0 1 1 0 0chr18 0 0 0 0 1 0chr19 0 0 0 0 0 1Y = MDA231$Y; Y ## whether SNAs are affected by CNAsnon-cna_region chr7 chr12 chr18 chr19BRAF 0 1 0 0 0KRAS 0 0 1 0 0ALPK2 0 0 0 1 0RYR1 0 0 0 0 1
关于运行问题
这个软件包要想获得最后的结果，需要分四个步骤：
Markov chain Monte Carlo (MCMC)确定克隆个数；
Bayesian information criterion (BIC)确定最优个数；
样本树的后验评估；
克隆进化树的输出和绘图。
#2.4 MCMC samplingK = 3:6 # number of subclonesnumchain = 20 # number of chains with random initiationssampchain = canopy.sample(R = R, X = X, WM = WM, Wm = Wm, epsilonM = epsilonM,epsilonm = epsilonm, C = C, Y = Y, K = K, numchain = numchain,max.simrun = 50000, min.simrun = 10000,writeskip = 200, projectname = projectname, cell.line = TRUE,plot.likelihood = TRUE)save.image(file = paste(projectname, '_postmcmc_image.rda',sep=''),compress = 'xz')length(sampchain) ## number of subtree spaces (K=3:6)length(sampchain[[which(K==4)]]) ## number of chains for subtree space with 4 subcloneslength(sampchain[[which(K==4)]][[1]]) ## number of posterior trees in each chain#2.5 BIC for model selectionburnin = 100thin = 5 # If there is error in the bic and canopy.post step below, make sure# burnin and thinning parameters are wisely selected so that there are# posterior trees left.bic = canopy.BIC(sampchain = sampchain, projectname = projectname, K = K,numchain = numchain, burnin = burnin, thin = thin, pdf = FALSE)optK = K[which.max(bic)]#2.6 Posterior evaluation of sampled treespost = canopy.post(sampchain = sampchain, projectname = projectname, K = K,+ numchain = numchain, burnin = burnin, thin = thin, optK = optK,+ C = C, post.config.cutoff = 0.05)samptreethin = post[[1]] # list of all post-burnin and thinning treessamptreethin.lik = post[[2]] # likelihoods of trees in samptreeconfig = post[[3]] # configuration for each posterior treeconfig.summary = post[[4]] # configuration summaryprint(config.summary)# first column: tree configuration# second column: posterior configuration probability in the entire tree space# third column: posterior configuration likelihood in the subtree space#2.7 Tree output and plottingconfig.i = config.summary[which.max(config.summary[,3]),1]cat('Configuration', config.i, 'has the highest posterior likelihood!\n')# plot the most likely tree in the posterior tree spaceoutput.tree = canopy.output(post, config.i, C)canopy.plottree(output.tree)# plot the tree with configuration 1 in the posterior tree spaceoutput.tree = canopy.output(post, 1, C)canopy.plottree(output.tree, pdf=TRUE, pdf.name =paste(projectname,'_first_config.pdf',sep=''))
Markov chain Monte Carlo (MCMC)确定克隆个数，如下图：
Bayesian information criterion (BIC)确定最优个数，如下图所示：
克隆进化树的输出和绘图，如下图所示：
关于应用上存在的问题
该软件包使用起来还是有几个缺点需要改进的：
输入文件的获得：需要自行提前筛选出来在不同时期的突变频率的变化，但其实这个筛选过程在大量测序WES or WGS中还是很难拿捏的准的；
运行时间问题：我只用该软件包给的例子MDA231，等待的时间大概是1个多小时；
代码中参杂了许多复杂的参数，需要深入的解读，其实可以更简单一些。
不过估计该软件包主要也是在说算法多一些，也不一定非得做成那种傻瓜式的参数，只是我这人比较2，喜欢看简单参数，不看算法的过程，下期将介绍利用Canopy生产的结果做后续分析的软件包Cardelino。
Reference:
Jiang Y, Qiu Y, Minn AJ, Zhang NR. Assessing intratumor heterogeneity and tracking longitudinal and spatial clonal evolutionary history by next-generation sequencing. Proc Natl Acad Sci U S A. 2016 Sep 13;113(37):E5528-37.
网友评论