TCGA数据库copy number variation数据分析

在上一篇笔记里，学习了利用maftools分析TCGA数据库里的simple nucleotide variation数据，现在来学习一下copy number variation数据的分析。不同于snv数据的是，cnv的数据从TCGA下载后需要一些数据的预处理、以及下载其他的一些必需文件。这篇笔记就是记录从文件下载、数据预处理、以及使用GISTIC在线软件生成maftools的输入文件。

（一）下载并整理TCGA的CNV数据

这里使用TCGAbiolinks包下载cnv数据：

> library(TCGAbiolinks)
#下载临床信息
> cancer_type <- "TCGA-HNSC"
> clinical <- GDCquery_clinic(project= cancer_type,type = "clinical")
> library(SummarizedExperiment)
#下载cnv数据的seg文件
> data_type <- "Masked Copy Number Segment"#选择数据类型
> data_category <- "Copy Number Variation" #选择数据类别
> workflow_type <- "DNAcopy"
> query_hnsc_cnv <- GDCquery(project = cancer_type, 
                                        data.category = data_category, 
                                        data.type =  data_type, 
                                        workflow.type = workflow_type)
> GDCdownload(query_hnsc_cnv, method = "api")
#保存
> CNV_files <- GDCprepare(query = query_hnsc_cnv,save = TRUE, save.filename = "CNV_TCGA_HNSC.rda") 

上面data_type选择“Masked”版本，是否是masked的区别请见文章CNV拷贝数变异分析（GISTIC、maftools），masked版本是去掉了gemline里已知的mutations。因为我们想研究的是肿瘤样品里特异的copy number，所以体细胞里已经包含的就不需要了。

现在这个seg文件是一个7列173575行的表，需要把这个表的列重新排列一下，把第一列删掉，然后把sample放到第一列

#读取数据并重新排列
> cnv <- load("D:/yanfang/TCGA_maf/CNV_TCGA_HNSC.rda")
> hnsc_seg <- eval(parse(text = cnv))
> hnsc_seg <- hnsc_seg[,-1]
> hnsc_seg <- hnsc_seg[,c('Sample','Chromosome','Start','End','Num_Probes','Segment_Mean')]

整理完是这样的

提取肿瘤样品（sample列里第14,15位是01的样品）：

> tumor_seg <- hnsc_seg[substr(hnsc_seg$Sample,14,15) == "01",]

（二）Marker文件的下载以及格式格式调整

运行GISTIC除了上面我们已经整理好的seg文件，还需要marker file。

Markers File数据下载：
在TCGA官网下载，地址：(https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files

选择：

下载后，在R里读取：

> hg_marker_file <- read.delim("snp6.na35.remap.hg38.subset.txt.gz")
> View(hg_marker_file)

NOTE：如果上面你下载的seg文件是masked版本的，要提取marker文件里freqcnv=FALSE的行。

> hg_marker_file <- hg_marker_file[hg_marker_file$freqcnv =="FALSE",]

另外根据GenePattern官网的说明，marker file只需要3列（官网）：

所以把上面的marker文件取前三列：

> hg_marker_file <- hg_marker_file[,c(1,2,3)]
> write.table(hg_marker_file,"hg_marker_file.txt",sep = "\t",col.names = TRUE,row.names = F)

（三）运行GISTIC（online版）

拿到了seg文件和marker file文件，就可以用GISTIC软件分析了。GISTIC可以下载到你的电脑使用，也可以在线分析。这里我只使用在线软件进行分析，安装的方法可以参考：用GISTIC多个segment文件来找SCNA变异，安装过程还需要配置matlab环境，比较麻烦，我这么懒就直接使用在线的分析了。

使用的在线平台是GenePattern，这是一个非常强大的应用平台，需要用邮箱注册账号，注册的过程也非常快，2分钟内搞定。GenePattern平台里包含了150多个工具，可以进行基因表达分析、蛋白质组分析、SNP分析、流式细胞分析、RNA-seq分析等等。

GenePatter网站：https://cloud.genepattern.org/gp/pages/index.jsf

打开以后是这样的：

我这是注册后的界面了，如果你第一次使用，页面会提示里注册一个新账户

首先，点击界面左边的Browse Modules查看所有工具：

找到GISTIC2.0：

点击GISTIC2.0后，界面会变成这样：

然后选择基因组版本，我使用的数据是TCGA数据，要用human的hg38：

其他参数没有改动，直接点击Run，如果你想修改job memory，cpucount也可以，我用的都是默认的:

点击Run后，界面会变成这样（运行时间长短取决于你之前选择的内存、cpu数，我这个运行了不到一个小时吧）：

运行完成后，界面会显示：

一共是19个输出文件，单击文件，再点击save，就可以保存在电脑里了：

（四）理解主要输出文件

（1）All Lesions File (all_lesions.conf_XX.txt, XX 指的是置信水平)

这个文件是GISTIC运行结果的总结文件。里面包含所有扩增和缺失的显著区域，以及每一个样品在这些区域是扩增还是缺失的信息。从第10列开始，就是每一个样品的信息，用命令行查看一下这个文件：

前9列的信息的含义分别是：

1. Unique Name: 区域的名称（如：Amplification Peak 1）
2. Descriptor: 在基因组上对该区域的描述
3. Wide Peak Limits: “宽峰”边界，最可能包含某些靶向基因。这一列内容报刊染色体坐标以及probes。
4. Peak Limits: 最大扩增以及最大缺失区域的边界。
5. Region Limits: 完整显著的扩增区域或者缺失区域的边界。
6. q-values: 该peak区域的q值。
7. Residual q-values: 除去扩增/缺失后的peak（在同一条染色体上与其他更为显著peak区域重叠）的q值。
8. Broad or Focal: 鉴定该区域显著性是由于“broad”事件，还是“focal”事件。还是独立于这些事件（“both”）
9. Amplitude Threshold: 以键值对显示，后面的样品的值都代表什么意思：

0代表拷贝数在超出阈值的部分没有扩正或者缺失。1代表样品有低水平的拷贝数变异。2代表样品有高水平的拷贝数变异。

（2）Amplification Genes File (Amp_genes.conf_XX.txt, XX 是置信水平)

这个表是扩增的peaks，里面包含了基因名称，每一列就是一个扩增区域。前4行是：

1. cytoband
2. q-value
3. residual q-value
4. wide peak boundaries

上面4行与all lesions file里的内容是一样的，但是剩下的行里列出了每一个wide peak里包含的基因。如果该区域没有包含基因，则列出离这个区域最近的基因，并用[]括起来。

该文件长这样：

（3）Deletion Genes File (Del_genes.conf_XX.txt, XX是置信水平)

这个表是缺失peaks，格式与上面介绍的文件一样：

（4）GISTIC Scores File (scores.gistic)

这个表包含segment的q值，分数，以及样品的扩增/缺失频率。这个分数文件列出了q值（实际上以-log10(q)来表示的），G分数，在变异的样品里平均振幅，变异的频率。这个文件可以在IGV里进行可视化。

放在IGV里看看：

deletion scores as a blue line and amplification scores as a red line

（5）Segmented Copy Number (raw_copy_number.pdf and raw_copy_number.png)

这个文件生成两种格式，一种PDF，另一种PNG，用热图的方式展示copy number。这里我生成的文件名字和这个不一样，我的文件名字是raw_copy_number，但是里面内容是一样的。

（6）Amplification GISTIC plot (amp_qplot.pdf and amp_qplot.png)

扩增GISTIC图，这个文件也是生成了两种格式的版本。图上方的数字是G-score（代表了该变异的振幅，以及在所有样品里该变异发生的频率），下方的数字是q值。

（7）Deletion GISTIC plot (del_qplot.pdf and del_quplot.png)

缺失图。与上面介绍的图是一样的，只不过显示的都是缺失peak。

参考文章：
1.CNV拷贝数变异分析（GISTIC、maftools）
2.The Use of Big Cancer Genomics Data- How to Get and Interpret the Data?
3.TCGA CNV全攻略