Topic 6. 克隆进化之 Canopy

前五期的肿瘤可进化分析都是基于组织或血液检出的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包的安装基本就几种方法：

1. CRAN直接安装install.packages()，

2. BiocManager::install()

3. devtools::install_github()，

4. 其他直接安装的模式。

该软件的安装如下：

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 project
R = MDA231$R; R ## mutant allele read depth (for SNAs)
   MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lung
BRAF           155           59          136                  77           49
KRAS            44           21           54                  19           17
ALPK2           37           17           28                  10            7
RYR1            44            0           26                   0            0
 
 X = MDA231$X; X ## total depth (for SNAs)
  MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lung
BRAF           157          111          177                 146           71
KRAS            44           30           64                  42           27
ALPK2           63           17           65                  24            7
RYR1           107           56          165                  55           43

WM = MDA231$WM; WM ## observed major copy number (for CNA regions)
  MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lung
chr7         2.998        2.002        2.603               2.000        2.001
chr12        1.998        1.998        1.603               1.001        1.999
chr18        1.000        2.992        1.000               1.002        2.996
chr19        2.000        2.000        2.000               2.000        2.000

Wm = MDA231$Wm; Wm ## observed minor copy number (for CNA regions)
 MCP1833_bone MCP1834_lung MCP2287_bone MDA-MB-231_parental MCP3481_lung
chr7         0.002        0.998        0.397               1.000        0.999
chr12        0.002        0.998        0.397               1.000        0.999
chr18        1.000        0.004        1.000               0.999        0.002
chr19        1.000        1.000        1.000               1.000        1.000

epsilonM = MDA231$epsilonM ## standard deviation of WM, pre-fixed here
epsilonm = 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 overlap
C = MDA231$C; C
 chr7_1 chr7_2 chr12_1 chr12_2 chr18 chr19
chr7       1      1       0       0     0     0
chr12      0      0       1       1     0     0
chr18      0      0       0       0     1     0
chr19      0      0       0       0     0     1

Y = MDA231$Y; Y ## whether SNAs are affected by CNAs
    non-cna_region chr7 chr12 chr18 chr19
BRAF               0    1     0     0     0
KRAS               0    0     1     0     0
ALPK2              0    0     0     1     0
RYR1               0    0     0     0     1

• 关于运行问题

这个软件包要想获得最后的结果，需要分四个步骤：

1. Markov chain Monte Carlo (MCMC)确定克隆个数；

2. Bayesian information criterion (BIC)确定最优个数；

3. 样本树的后验评估；

4. 克隆进化树的输出和绘图。

#2.4 MCMC sampling
K = 3:6 # number of subclones
numchain = 20 # number of chains with random initiations
sampchain = 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 subclones
length(sampchain[[which(K==4)]][[1]]) ## number of posterior trees in each chain

#2.5 BIC for model selection
burnin = 100
thin = 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 trees
post = 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 trees
samptreethin.lik = post[[2]] # likelihoods of trees in samptree
config = post[[3]] # configuration for each posterior tree
config.summary = post[[4]] # configuration summary
print(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 plotting
config.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 space
output.tree = canopy.output(post, config.i, C)
canopy.plottree(output.tree)
# plot the tree with configuration 1 in the posterior tree space
output.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)确定克隆个数，如下图：

1. Bayesian information criterion (BIC)确定最优个数，如下图所示：

• 克隆进化树的输出和绘图，如下图所示：

• 关于应用上存在的问题

该软件包使用起来还是有几个缺点需要改进的：

1. 输入文件的获得：需要自行提前筛选出来在不同时期的突变频率的变化，但其实这个筛选过程在大量测序WES or WGS中还是很难拿捏的准的；

2. 运行时间问题：我只用该软件包给的例子MDA231，等待的时间大概是1个多小时；

3. 代码中参杂了许多复杂的参数，需要深入的解读，其实可以更简单一些。

不过估计该软件包主要也是在说算法多一些，也不一定非得做成那种傻瓜式的参数，只是我这人比较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.

另外补充一下，需要学习生信，欢迎进群交流！！

关注公众号，每日有更新！