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包的安装基本就几种方法：
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 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 0X = 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 43WM = 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.000Wm = 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.000epsilonM = 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; Cchr7_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 1Y = MDA231$Y; Y ## whether SNAs are affected by CNAsnon-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
- 关于运行问题
这个软件包要想获得最后的结果，需要分四个步骤：
Markov chain Monte Carlo (MCMC)确定克隆个数；
Bayesian information criterion (BIC)确定最优个数；
样本树的后验评估；
克隆进化树的输出和绘图。
#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)确定克隆个数，如下图：
- 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.
另外补充一下，需要学习生信，欢迎进群交流！！
关注公众号，每日有更新！
Topic 6. 克隆进化之 Canopy相关推荐
- Topic 3. 克隆进化之 fishplot
上面两篇肿瘤克隆进化之后第三篇就是大家非常感兴趣的鱼图也就是sankey图,之前有发过一遍,但是好多学者反应不是很好使用,现在重新发一次,即fishplot R包 **背景:**深度大规模并行测序技术 ...
- Topic 8. 克隆进化之 RobustClone
**上一期介绍了 Cardelino 的使用,这期介绍RobustClone,我们知道单细胞数据的SNV和CNV检出率都是非常低的,那么怎么有效的获得可使用的突变位点是重中之重,那么该软件包Robus ...
- Topic 5. 克隆进化之 CITUP
这期介绍CITUP 全称 Conality Inference in Tumors Using Phylogeny,是一种使用系统发育论进行肿瘤的克隆推断生物信息学工具,可用于从单个患者获得的多个样本 ...
- Topic 9. 克隆进化之 TimeScape
我们利用Pyclone和CITUP得到了三个文件即cellfreq.txt和tree.txt 和sample_id,下面我们就利用TimeScape搞一下可视化,在这里不会出现具体的基因或突变位点,但 ...
- IF：伴FLT3-ITD突变的急性髓系白血病在米哚妥林治疗下的克隆进化
点击关注,桓峰基因 桓峰基因的教程不但教您怎么使用,还会定期分析一些相关的文章,学会教程只是基础,但是如果把分析结果整合到文章里面才是目的,觉得我们这些教程还不错,并且您按照我们的教程分析出来不错的结 ...
- SCS【13】单细胞转录组之识别细胞对“基因集”的响应 (AUCell)
点击关注,桓峰基因 桓峰基因公众号推出单细胞系列教程,有需要生信分析的老师可以联系我们!首选看下转录分析教程整理如下: Topic 6. 克隆进化之 Canopy Topic 7. 克隆进化之 Car ...
- SCS【12】单细胞转录组之评估不同单细胞亚群的分化潜能 (Cytotrace)
点击关注,桓峰基因 Topic 6. 克隆进化之 Canopy Topic 7. 克隆进化之 Cardelino Topic 8. 克隆进化之 RobustClone SCS[1]今天开启单细胞之旅, ...
- IF: 8+ 基于单细胞 RNA-seq 构建非小细胞肺癌免疫反应的中性粒细胞预后模型
. 单细胞生信分析教程 桓峰基因公众号推出单细胞生信分析教程并配有视频在线教程,目前整理出来的相关教程目录如下: Topic 6. 克隆进化之 Canopy Topic 7. 克隆进化之 Cardel ...
- SCS【14】单细胞调节网络推理和聚类 (SCENIC)
点击关注,桓峰基因 桓峰基因公众号推出单细胞系列教程,有需要生信分析的老师可以联系我们!单细胞系列分析教程整理如下: Topic 6. 克隆进化之 Canopy Topic 7. 克隆进化之 Card ...
最新文章
- c++switch实现猜拳_策略模式+简单工厂+注解消除 if-else/switch-case
- Jmeter性能测试之Switch控制器使用
- navicat 批量修改列数据
- [转] 学习，思维三部曲：WHAT、HOW、WHY（通过现象看本质）
- python在电脑上的图标_在python scrip中嵌入图标
- CORBA GIOP消息格式学习
- PTA 栈 (20分)(全网首发)(实现一个栈Stack，要求实现Push（出栈）、Pop（入栈）、Min（返回最小值的操作）的时间复杂度为O(1))
- NLP-Beginner：自然语言处理入门练习-任务一
- MySQL month()函数
- guava 并发_使用Guava对并发应用程序进行基于对象的微锁定
- Python time mktime()方法
- Android MotionEvent中getX()、getRawX()和getTop()的区别
- PHP付费资源下载交易平台网站源码
- 送给初学.net兄弟们的一些话
- 玩抖音，你喜欢的，都是对自身没好处的
- 计算机软硬件逻辑等价性是指,南航计算机组成原理复习ppt.ppt
- 英文的pdf文件怎么翻译成中文
- DigiCert SSL证书支持中文域名申请吗？
- Python 爬取豆瓣电影Top250
- JenKins添加Git报错Error performing git command: git ls-remote -h