從WGS測序得到的VCF檔案裡面提取位於外顯子區域的【直播】我的基因組84
阿新 • • 發佈：2022-05-03
首先要下載並且得到人類基因組的外顯子座標記錄檔案
這裡我用的參考基因組版本仍然是hg19，所以去CCDS資料庫裡面下載對應版本，並且格式化成BED檔案。
wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/archive/Hs37.3/CCDS.20110907.txtcat 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 801942 8024331 861321 8613921 865534 8657151 866418 8664681 871151 8712751 874419 8745081 874654 8748391 876523 8766851 877515 8776301 877789 877867從VCF檔案裡面根據BED檔案進行抽提
這裡就不自己造輪子了，用現成的工具，而且是我們用過很多次的SnpEff套件，程式碼如下
cat snp.vcf | java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar intervals hg19exon.bed >hg19exon.snp.vcfcat indel.vcf | java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar intervals hg19exon.bed >hg19exon.indel.vcf
可以把我經由GATK best practice流程得到的SNP/INDEL記錄的VCF檔案都進行提取，用程式碼 wc -l *vcf簡單統計一下提取的效果，如下：
1042 hg19_exon.indel.vcf 25067 hg19_exon.snp.vcf 754755 indel.vcf 3784343 snp.vcf很明顯可以看到，位於外顯子區域的mutation畢竟是少數，這時候還可以繼續看看那些在外顯子上面卻沒有被dbSNP資料庫記錄的mutation還有多少:
cat hg19_exon.snp.vcf |grep -v "^#" |cut -f 3
|grep '.' |wc
仍然有2315個SNV在外顯子區域，卻沒有被dbSNP資料庫記錄，可能是我的家族特異性的位點，屬於正常的基因型多樣性，也有極小的可能性這些位點是後發突變，也就是通常癌症研究領域的somatic mutation 。
用下面的程式碼可以簡單瀏覽一下這些在外顯子上面的未知突變的情況。
cat hg19_exon.snp.vcf |perl -alne '{print if $F[2] eq "."}' |less -Scat hg19_exon.indel.vcf |perl -alne '{print if $F[2] eq "."}' |less -S