TCGA CNV全攻略
明白什么是CNV
对正常人来说，基因组应该是二倍体的，所以凡是测到非2倍体的地方都是CNV。但是CNV本身就是人群遗传物质多样性的体现，所以对癌症样本来说，是需要过滤掉正常人体内的germline的CNV，得到somatic的CNV。
TCGA的CNV测量及计算
TCGA里面主要是通过Affymetrix SNP6.0 array这款芯片来测拷贝数变异！
值得注意的是，并不是只有TCGA利用了SNP6.这个芯片数据，著名的CCLE计划也对一千多细胞系处理了SNP6.0芯片，数据也是可以下载的。
对SNP6.0的拷贝数芯片来说，通常是用PICNIC等软件处理原始数据，就可以得到的segment记录文件，每个样本一个结果，下面是示例结果：
Chromosome Start End Num_Probes Segment_Mean 1 61735 1510801 226 -0.0397 1 1627918 1672603 17 -0.92 1 1687587 16153497 8176 0.0077 1 16153536 16153925 5 -2.7441 1 16154201 16155010 4 -0.8711 1 16165661 72768498 34630 0.0048 1 72768916 72811148 46 -1.7394 1 72811904 95674710 14901 0.0026 1 95676511 95676518 2 -1.6636
表明了某条染色体的某个区域内，SNP6.0芯片设计了多少个探针，芯片结果的拷贝数值是多少(这个区域的拷贝数用Segment_Mean)。
通常二倍体的Segment_Mean值为0，可以用-0.2和0.2来作为该区域是否缺失或者扩增。
具体数据处理流程见NIH的TCGA官网：https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/
参考文献：http://mcr.aacrjournals.org/content/12/4/485.long
TCGA的CNV数据下载
众所周知，TCGA的数据的开放程度分成了4个等级，一般人都是下载level 3 的数据，对CNV数据也是如此。
我比较喜欢去broad institute下载TCGA的数据，所有的文件都以目录的形式存放着：
https://gdac.broadinstitute.org/runs/stddata__latest/ https://gdac.broadinstitute.org/runs/analyses__latest/
如果要下载level3的数据，就用stddata__latest 这个url即可，打开可以看到里面列出了所有的癌症种类，假如我们感兴趣的是BRCA，就直接点击进入，用下面的url即可。
https://gdac.broadinstitute.org/runs/analyses__latest/data/BRCA/20160128/
打开url可以看到非常多的文件，这里我们感兴趣的是snp6芯片的拷贝数结果，而且一般是基于hg19版本的。
wget -c -r -np -nH -k -L --cut-dirs 6 -p -A "*snp_6*hg19*Level_3*" http://gdac.broadinstitute.org/runs/stddata__2016_01_28/data/BRCA/20160128/
如果要下载其它癌症种类，只需要改变url里面的BRCA即可。如果要下载其它类型的数据，只需要改变-A 后面的匹配规则即可，其实就是打开上面url看到的几十个文件的文件名的规律。
"*snp_6*hg19*Level_3*"
几分钟就下载完数据啦，然后你就会看到下面两个截然不同的：
Merge_snp__genome_wide_snp_6__broad_mit_edu__Level_3__segmented_scna_hg19__seg Merge_snp__genome_wide_snp_6__broad_mit_edu__Level_3__segmented_scna_minus_germline_cnv_hg19__seg
其中minus了germline的CNV的就是我们想要的癌症相关的somatic CNV咯!
拿到CNV做什么？
首先两个segment文本文件已经可以直接载入IGV查看所有BRCA样本的CNV情况啦，如下所示：
CNV深度分析
注释基因
前面我们下载的CNV都是基于基因组区域的，比如1号染色体的61735起始坐标到1510801终止坐标。在IGV里面倒是可以看出一些pattern，但是人们感兴趣的往往是这些位置上面到底有哪些基因。接下来就可以对基因进行各种下游分析。
既然是对基因组片段做基因注释，那么首先就需要拿到基因的坐标信息咯，我是在gencode数据库里面下载，然后解析成下面的bed格式的，如下：
head ~/reference/gtf/gencode/protein_coding.hg19.position chr1 69091 70008 OR4F5 chr1 367640 368634 OR4F29 chr1 621096 622034 OR4F16 chr1 859308 879961 SAMD11 chr1 879584 894689 NOC2L chr1 895967 901095 KLHL17 chr1 901877 911245 PLEKHN1 chr1 910584 917473 PERM1 chr1 934342 935552 HES4 chr1 936518 949921 ISG15
然后要把我们下载的CNV文本文件，转为bed格式的，就是把列的顺序调换一下：
head Features.bed chr1 3218610 95674710 TCGA-3C-AAAU-10A-01D-A41E-01 53225 0.0055 chr1 95676511 95676518 TCGA-3C-AAAU-10A-01D-A41E-01 2 -1.6636 chr1 95680124 167057183 TCGA-3C-AAAU-10A-01D-A41E-01 24886 0.0053 chr1 167057495 167059336 TCGA-3C-AAAU-10A-01D-A41E-01 3 -1.0999 chr1 167059760 181602002 TCGA-3C-AAAU-10A-01D-A41E-01 9213 -8e-04 chr1 181603120 181609567 TCGA-3C-AAAU-10A-01D-A41E-01 6 -1.2009 chr1 181610685 201473647 TCGA-3C-AAAU-10A-01D-A41E-01 12002 0.0055 chr1 201474400 201474544 TCGA-3C-AAAU-10A-01D-A41E-01 2 -1.4235 chr1 201475220 247813706 TCGA-3C-AAAU-10A-01D-A41E-01 29781 -4e-04
避免重复造轮子，我就用我擅长的bedtools解决这个需求吧，命令很简单，如下：
bedtools intersect -a Features.bed -b ~/reference/gtf/gencode/protein_coding.hg19.position -wa -wb \ | bedtools groupby -i - -g 1-4 -c 10 -o collapse
注释结果，我挑了几个可以看的给大家，可以看到，每个CNV片段都注释到了对应的基因，有些特别大的片段，会被注释到非常多的基因。
chr8 42584924 42783715 TCGA-5T-A9QA-01A-11D-A41E-01 CHRNB3,CHRNA6,THAP1,RNF170,HOOK3 chr8 42789728 42793594 TCGA-5T-A9QA-01A-11D-A41E-01 HOOK3 chr8 42797957 42933372 TCGA-5T-A9QA-01A-11D-A41E-01 RP11-598P20.5,FNTA,HOOK3 chr8 70952673 70964372 TCGA-5T-A9QA-01A-11D-A41E-01 PRDM14 chr10 42947970 43833200 TCGA-5T-A9QA-01A-11D-A41E-01 BMS1,RET,RASGEF1A,ZNF33B,CSGALNACT2 chr10 106384615 106473355 TCGA-5T-A9QA-01A-11D-A41E-01 SORCS3 chr10 106478366 107298256 TCGA-5T-A9QA-01A-11D-A41E-01 SORCS3 chr10 117457285 117457859 TCGA-5T-A9QA-01A-11D-A41E-01 ATRNL1 chr11 68990173 69277078 TCGA-5T-A9QA-01A-11D-A41E-01 MYEOV chr11 76378708 76926535 TCGA-5T-A9QA-01A-11D-A41E-01 LRRC32,B3GNT6,OMP,TSKU,MYO7A,ACER3,CAPN5
找somatic CNVs
仔细看上面IGV的pattern你会发现某些染色体的某些片段经常会扩增或者缺失，这个现象就是人们想研究是recurrent CNV regions，当然不会用肉眼看咯，这时候需要用GISTIC这个软件。找到了recurrent CNV regions同样是需要进行基因注释，才能进行下游分析咯。
PICNIC软件用法：http://www.bio-info-trainee.com/1299.html
GISTIC软件用法：http://www.bio-info-trainee.com/1648.html