肿瘤变异数据分析和可视化工具maftools:CNV的可视化
Maftools系列文章:
上篇文章《肿瘤变异数据分析和可视化工具maftools:突变的数据分析》介绍了使用maftools分析MAF格式的突变数据,但maftools本身功能不限于此,它也可以针对拷贝数变异进行一些分析。
CNV数据下载和处理
还是继续之前TCGA-LUAD的例子。Maftools接受的CNV数据需要是GISTIC输出的结果,包括4个文件——all_lesions.conf_XX.txt
、 amp_genes.conf_XX.txt
、del_genes.conf_XX.txt
以及scores.gistic
(其中XX代表置信水平)。但是从TCGA官网能下载到的GISTIC的输出结果只有LUAD.focal_score_by_genes.txt
,因此需要下载更上游的DNAcopy的输出结果,然后自己再根据TCGA的流程和参数跑一遍GISTIC。
但这样有点麻烦,更简单的办法就是下载其他数据库已经处理好的。这里我直接从Firehose下载GISTIC2的输出结果,把其中all_lesions.conf_99.txt
、amp_genes.conf_99.txt
、del_genes.conf_99.txt
、scores.gistic
四个文件挑出来。需要注意的是,Firehose的分析流程用的是hg19,而TCGA官方目前已经用hg38了,所以参考基因组的坐标是不一致的,因为后面画oncoplot需要用到突变数据,所以索性在Firehose也下载对应的临床数据和MAF文件。
不过怕麻烦就得付出代价,Firehose的数据更新速度实在太慢,直到今天9012年了还在用2016_01_28版的TCGA数据,导致TCGA-LUAD这个项目只有230个case的MAF文件(占总585个case的39%),而官网和UCSC Xena的数据虽然比较新,但是都下不了maftools需要的GISTIC2输出的4个文件,所以最好还是自己用DNAcopy的结果跑一下GISTIC2(如果有更好的办法可以在评论区留言!)。这里只是做演示,所以我还是偷下懒算了。
通过简单的数据整理后文件如下(包括各case的MAF文件合并以及《肿瘤变异数据分析和可视化工具maftools:突变数据下载和可视化》一文中提到的两个需要处理的地方):
$ ls -lhrt
total 237M
-rw-r--r-- 1 xiaofei xiaofei 459K Jun 7 19:26 all_lesions.conf_99.txt
-rw-r--r-- 1 xiaofei xiaofei 14K Jun 7 19:26 amp_genes.conf_99.txt
-rw-r--r-- 1 xiaofei xiaofei 63K Jun 7 19:26 del_genes.conf_99.txt
-rw-r--r-- 1 xiaofei xiaofei 2.5M Jun 7 19:26 scores.gistic
-rw-r--r-- 1 xiaofei xiaofei 223M Jun 7 19:34 TCGA.LUAD.maf
-rw-r--r-- 1 xiaofei xiaofei 11M Jun 7 19:40 LUAD-TP.samplefeatures.txt
使用maftools读取GISTIC输出的CNV数据并统计
使用函数readGistic
读入GISTIC2输出的4个文件:
> library(maftools)
> luad.gistic <- readGistic(gisticAllLesionsFile="all_lesions.conf_99.txt", gisticAmpGenesFile="amp_genes.conf_99.txt", gisticDelGenesFile="del_genes.conf_99.txt", gisticScoresFile="scores.gistic", isTCGA=TRUE)
## -Processing Gistic files..
## --Processing amp_genes.conf_99.txt
## --Processing del_genes.conf_99.txt
## --Processing scores.gistic
## --Summarizing by samples
直接键入GISTIC对象luad.gistic
就能得到一些简单的统计:
> luad.gistic
## An object of class GISTIC
## ID summary
## 1: Samples 516
## 2: nGenes 6057
## 3: cytoBands 75
## 4: Amp 177966
## 5: Del 896197
## 6: total 1074163
和之前介绍的MAF文件的统计方法类似,GISTIC对象可以使用getSampleSummary
、getGeneSummary
、getCytobandSummary
进行统计,并用write.GisticSummary
保存统计结果:
getSampleSummary(luad.gistic)
getGeneSummary(luad.gistic)
getCytobandSummary(luad.gistic)
write.GisticSummary(gistic=luad.gistic, basename="luad_gistic2")
CNV数据的可视化
GISTIC的输出结果可以用染色体图、气泡图(点图)、oncoplot进行可视化。
1. 染色体图
使用gisticChromPlot
可以绘制染色体各位点G-score:
gisticChromPlot(gistic=luad.gistic, markBands="all")
2. 气泡图
gisticBubblePlot
函数可以将展示显著变化的cytobands对应的样本数(X轴)和包含的基因数(Y轴),圆点的大小为-log10(qvalue):
gisticBubblePlot(gistic=luad.gistic)
3. oncoplot
可以使用gisticOncoPlot
函数以oncoplot的形式展示CNV。是否根据临床数据进行分类是可选的,如果要标注临床信息,需要先用read.maf
建立MAF对象,这里选择性别:
luad <- read.maf(maf="TCGA.LUAD.maf", clinicalData="LUAD-TP.samplefeatures.txt")
gisticOncoPlot(gistic=luad.gistic, clinicalData=getClinicalData(x=luad), clinicalFeatures="CLI_gender", sortByAnnotation=TRUE, top=10)
也可以用oncoplot
同时展示CNV和突变:
luad.plus.gistic <- read.maf(maf="TCGA.LUAD.maf", gisticAllLesionsFile="all_lesions.conf_99.tx
t", gisticAmpGenesFile="amp_genes.conf_99.txt", gisticDelGenesFile="del_genes.conf_99.txt", gist
icScoresFile="scores.gistic")
# 因为样本数太多,使用borderCol=NULL不显示样本边界。可惜的是目前gisticOncoPlot还没有加上这个参数
oncoplot(maf=luad.plus.gistic, borderCol=NULL, top=10)
4. 可视化segment文件
Maftools除了可以可视化GISTIC的输出结果以外,还可以可视化DNAcopy输出的segment文件。这里以《肿瘤变异数据分析和可视化工具maftools:突变的数据分析》一文中已经处理过的文件TCGA-38-4625.seg
为例:
plotCBSsegments("TCGA-38-4625.seg")
要注意的是,目前这个函数感觉还不太成熟,连help都没有写,而且对于女性,segment文件中没有Y染色体,会引发报错,因此需要在文件最后补上一行Y染色体才能正常绘图:
TCGA-38-4625Y 0 0 0 0
写的真好,可以写一个gistic2的应用的博客么,谢谢
学习了!
写的很好,谢谢分享
很好