从WGS测序得到的VCF文件里面提取位于外显子区域的【直播】我的基因组84

首先要下载并且得到人类基因组的外显子坐标记录文件

这里我用的参考基因组版本仍然是hg19，所以去CCDS数据库里面下载对应版本，并且格式化成BED文件。

1. wget  ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/archive/Hs37.3/CCDS.20110907.txt

2. cat CCDS.20110907.txt |perl -alne '{/[(.*?)]/;next unless $1;$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$" foreach split/,/,$exons;}' >hg19exon.bed

制作好的bed格式的人类全部的exon区域坐标文件如下：

1. 1    801942  802433

2. 1    861321  861392

3. 1    865534  865715

4. 1    866418  866468

5. 1    871151  871275

6. 1    874419  874508

7. 1    874654  874839

8. 1    876523  876685

9. 1    877515  877630

10. 1    877789  877867

从VCF文件里面根据BED文件进行抽提

这里就不自己造轮子了，用现成的工具，而且是我们用过很多次的SnpEff套件，代码如下

1. cat snp.vcf | java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar intervals hg19exon.bed >hg19exon.snp.vcf

2. cat indel.vcf | java -jar  ~/biosoft/SnpEff/snpEff/SnpSift.jar intervals hg19exon.bed >hg19exon.indel.vcf

可以把我经由GATK best practice流程得到的SNP/INDEL记录的VCF文件都进行提取，用代码 wc -l *vcf简单统计一下提取的效果，如下：

1.      1042 hg19_exon.indel.vcf

2.     25067 hg19_exon.snp.vcf

3.    754755 indel.vcf

4.   3784343 snp.vcf

很明显可以看到，位于外显子区域的mutation毕竟是少数，这时候还可以继续看看那些在外显子上面却没有被dbSNP数据库记录的mutation还有多少:

cat hg19_exon.snp.vcf |grep -v "^#" |cut -f 3 |grep '\.' |wc

仍然有2315个SNV在外显子区域，却没有被dbSNP数据库记录，可能是我的家族特异性的位点，属于正常的基因型多样性，也有极小的可能性这些位点是后发突变，也就是通常癌症研究领域的somatic mutation 。

用下面的代码可以简单浏览一下这些在外显子上面的未知突变的情况。

1. cat hg19_exon.snp.vcf |perl -alne '{print if $F[2] eq "."}' |less -S

2. cat hg19_exon.indel.vcf |perl -alne '{print if $F[2] eq "."}' |less -S

下一讲我们再进一步注释。

猜你喜欢

菜鸟入门

数据分析

编程实践

直播基因组分析

生信技能树