克隆进化详解

Topic 6. 克隆进化之 Canopy

2022-03-26  本文已影响0人  桓峰基因

    前五期的肿瘤可进化分析都是基于组织或血液检出的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            0  X = 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; C chr7_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 CNAs    non-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.

    上一篇 下一篇

    猜你喜欢

    热点阅读