生物信息|TCGA CNV全攻略

生信圈 · 公众号 · · 2017-06-22 20:41

正文

什么是CNV

对正常人来说，基因组应该是二倍体的，所以凡是测到非2倍体的地方都是CNV。但是CNV本身就是人群遗传物质多样性的体现，所以对癌症样本来说，是需要过滤掉正常人体内的germline的CNV，得到somatic的CNV。

CNV（copy-numbervariant）是指拷贝数目变异，也称拷贝数目多态性（copy-number polymorphism，CNP），是一个大小介于1kb至3MB的DNA片段的变异，在人类及动植物基因组中广泛分布，其覆盖的核苷酸总数大大超过单核苷酸多态性（SNP）的总数，极大地丰富了基因组遗传变异的多样性。按照CNV是否致病可分为致病性CNV、非致病性CNV和不明临床意义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/BioinformaticsPipelines/CNVPipeline/

参考文献：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同样是需要进行基因注释，才能进行下游分析咯。

本文转载自生信人，如有侵权请联系删除。