变异检测：vcf文件合并
vcf文件儲存的是樣本的變異信息文件，在同一批次分析中，如果不是采用joint calling的方式進行分析，最終會獲得單個樣本的變異數(shù)據(jù)。這種文件很難對同組不同樣本進行差異SNP分析，此處就需要對文件進行合并。vcf文件的合并有很多的軟件可以做，主要的就是GATK、vcftools和bcftools三種，但是具體的合并方法需要根據(jù)不同vcf文件中的信息來判斷。
1. 樣本相同、位點獨立的vcf文件合并
樣本相同、位點獨立指的是在不同文件中包含的樣本一致，但是位點批次質檢沒有交集，分染色體call變異的結果就是這一類的典型。這類最簡單的方法可以直接使用cat命令進行合并，但是未免不太專業(yè)，以下是三款軟件的合并方法。
1.1 gatk：GatherVcfs|MergeVcfs
gatk4提供了兩種合并vcf文件的方法，分別是GatherVcfs和MergeVcfs，兩個方法都是對相同樣本數(shù)據(jù)集的變異結果進行合并，命令示例如下。
# GatherVcfs gatk GatherVcfs -I concat-a.vcf -I concat-b.vcf -O combine_a_b_samesample_diffsites.vcf # MergeVcfs gatk MergeVcfs -I concat-a.vcf -I concat-b.vcf -O combine_a_b_diffsample_allsites_gatk.vcf左右滑動查看全部內(nèi)容。
兩個命令執(zhí)行結果完全一致，結果如下。
##fileformat=VCFv4.2 ##FILTER=<ID=q10,Description="Quality below 10"> ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth"> ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth"> ##contig=<ID=1,length=17540695> ##contig=<ID=2,length=14896646> ##contig=<ID=3,length=12399606> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT A # concat-a.vcf文件位點 1 100 . GTTT G 1806 q10 DP=35 GT:GQ:DP 0/1:409:35 1 110 . C T,G 1792 PASS DP=32 GT:GQ:DP 0/1:245:32 1 120 . GA G 628 q10 DP=21 GT:GQ:DP 1/1:21:21 1 130 . GAA GG 1016 PASS DP=22 GT:GQ:DP 0/1:212:22 1 140 . GT G 727 PASS DP=30 GT:GQ:DP 0/1:150:30 1 150 . TAAAA TA,T 246 PASS DP=10 GT:GQ:DP 1/2:12:10 2 100 . GTTT G 1806 q10 DP=35 GT:GQ:DP 0/1:409:35 2 110 . CAAA C 1792 PASS DP=32 GT:GQ:DP 0/1:245:32 2 120 . GA G 628 q10 DP=21 GT:GQ:DP 1/1:21:21 2 130 . GAA G 1016 PASS DP=22 GT:GQ:DP 0/1:212:22 2 140 . GT G 727 PASS DP=30 GT:GQ:DP 0/1:150:30 2 150 . TAAAA TA,T 246 PASS DP=10 GT:GQ:DP 1/2:12:10 2 160 . TAAAA TA,TC,T 246 PASS DP=10 GT:GQ:DP 0/2:12:10 # concat-b.vcf文件位點 3 241 . GTTT G 1806 q10 DP=35 GT:GQ:DP 0/1:409:35 3 251 . CAAA C 1792 PASS DP=32 GT:GQ:DP 0/1:245:32 3 261 . GA G 628 q10 DP=21 GT:GQ:DP 1/1:21:21 3 271 . GAA G 1016 PASS DP=22 GT:GQ:DP 0/1:212:22左右滑動查看全部內(nèi)容。
1.2 vcftools：vcf-concat
vcf-concat并不是vcftools的一個子命令，而是軟件安裝目錄下附帶的功能模塊，vcftools安裝完成后可以直接使用，命令如下。
/vcftools/bin/vcf-concat concat-a.vcf concat-b.vcf > combine_a_b_samesample_diffsites_vcftools.vc左右滑動查看全部內(nèi)容。
在進行較大文件合并時，最好將vcf文件進行壓縮并創(chuàng)建索引，效率會快。如果待合并的vcf文件很多，可以將文件名寫入到一個文件，由參數(shù)-f進行操作。該命令合并完的vcf位點變異信息與gatk結果一致，只不過表頭信息順序會發(fā)生改變，不影響數(shù)據(jù)使用。
1.3 bcftools：concat
bcftools軟件在處理vcf文件時，某些情況下會優(yōu)于vcftools，合并vcf文件的命令如下。
bcftools concat concat-a.vcf concat-b.vcf -o combine_a_b_samesample_diffsites_bcftools.vcf左右滑動查看全部內(nèi)容。
合并之后的vcf文件位點信息與其余兩款軟件處理結果一致，只不過在輸出結果的header中會出現(xiàn)運行的命令，示例如下。
##bcftools_concatVersion=1.3.1+htslib-1.3.1 ##bcftools_concatCommand=concat -o combine_a_b_samesample_diffsites_bcftools.vcf concat-a.vcf concat-b.vcf左右滑動查看全部內(nèi)容。
2. 樣本不同，位點相同或不同
這種合并方式主要是對不同樣本的變異文件進行合并，合并時共有位點會進行合并統(tǒng)計，非共有位點若在某一個樣本中沒有變異，則會自動記為缺失
2.1 vcftools:vcf-merge
模塊名稱為vcf-merge，在進行merge操作時，會對文件中的位點進行重排，耗時較長。輸入文件需要壓縮后創(chuàng)建索引，示例命令如下。
bgzip merge-test-a.vcf.gz && tabix merge-test-a.vcf.gz bgzip merge-test-b.vcf.gz && tabix merge-test-b.vcf.gz /vcftools/bin/vcf-merge merge-test-a.vcf.gz merge-test-b.vcf.gz > combine_a_b_diffsamples_allsites_vcftools.vcf左右滑動查看全部內(nèi)容。
合并之后會對不同文件中的數(shù)據(jù)集進行整合，沒有變異的位點會自動標記為缺失。
2.2 bcftools:merge
使用的方法為merge，示例如下。
bcftools merge merge-test-a.vcf.gz merge-test-b.vcf.gz -o combine_a_b_diffsamples_allsites_bcftools.vcf左右滑動查看全部內(nèi)容。
該方法也需要實現(xiàn)對所有vcf文件進行壓縮并創(chuàng)建索引，否則程序無法運行。
需要注意的就是，merge合并之后，不同軟件生成的結果會存在很大的差異，主要是統(tǒng)計結果的重新計算上，示例如下。
# merge-test-a.vcf 1 3184885 . TAAAA TA,T 246 PASS DP=10 GT:GQ:DP 1/2:12:10 # merge-test-b.vcf 1 3184885 . TAAA T 598 PASS DP=16 GT:GQ:DP 0/1:435:16 # combine_a_b_diffsamples_allsites_vcftools.vcf 1 3184885 . TAAAA TA,T 422.00 PASS AC=2,1;AN=4;DP=26;SF=0,1 GT:DP:GQ 1/2:10:12 0/1:16:435 # combine_a_b_diffsamples_allsites_bcftools.vcf 1 3184885 . TAAAA TA,T 598 PASS DP=26 GT:GQ:DP 1/2:12:10 0/1:435:16左右滑動查看全部內(nèi)容。
gatk:CombineVariants
gatk3提供了一個CombineVariants可以進行變異數(shù)據(jù)的合并，而gatk4中并沒有找到功能相同的模塊。CombineVariants使用如下。
java -jar /GenomeAnalysisTK-3.8/GenomeAnalysisTK.jar -T CombineVariants -V merge-test-a.vcf -V merge-test-b.vcf -o combine_a_b_diffsample_allsites_gatk.vcf -R ref.fna左右滑動查看全部內(nèi)容。
合并結果示例如下。
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT A B 1 3062915 . GTTT G,GT 1806 q10;q20 AC=1,1;AF=0.250,0.250;AN=4;DP=49;set=FilteredInAll GT:DP:GQ 0/1:35:99 0/2:14:99 1 3106154 . CAAAA CA,C 1792 PASS AC=1,1;AF=0.250,0.250;AN=4;DP=47;set=Intersection GT:DP:GQ 0/1:32:99 0/2:15:99 1 3157410 . GA G 628 PASS AC=3;AF=0.750;AN=4;DP=32;set=filterInvariant-variant2 GT:DP:GQ 1/1:21:21 0/1:11:49 1 3162006 . GAA G 1016 PASS AC=2;AF=0.500;AN=4;DP=41;set=Intersection GT:DP:GQ 0/1:22:99 0/1:19:99 1 3177144 . GT G 727 PASS AC=2;AF=0.500;AN=4;DP=54;set=Intersection GT:DP:GQ 0/1:30:99 0/1:24:99 1 3184885 . TAAAA TA,T 246 PASS AC=2,1;AF=0.500,0.250;AN=4;DP=26;set=Intersection GT:DP:GQ 1/2:10:12 0/1:16:99 2 3188209 . GA G 162 . AC=1;AF=0.500;AN=2;DP=15;set=variant2 GT:DP:GQ ./. 0/1:15:99 2 3199812 . G GTT,GT 481 PASS AC=1,1;AF=0.500,0.500;AN=2;DP=26;set=variant GT:DP:GQ 1/2:26:99 ./. 3 3199812 . G GTT,GT 353 PASS AC=1,1;AF=0.500,0.500;AN=2;DP=19;set=variant2 GT:DP:GQ ./. 1/2:19:99 3 3199815 . G A 353 PASS AC=1;AF=0.500;AN=2;DP=19;set=variant2 GT:DP:GQ ./. 0/1:19:99 3 3212016 . CTT C,CT 565 PASS AC=1,1;AF=0.500,0.500;AN=2;DP=26;set=variant GT:DP:GQ 1/2:26:91 ./. 4 3212016 . CTT C 677 q20 AC=1;AF=0.500;AN=2;DP=15;set=FilteredInAll GT:DP:GQ ./. 0/1:15:99 4 3258448 . TACACACAC T 325 PASS AC=1;AF=0.500;AN=2;DP=31;set=variant GT:DP:GQ 0/1:31:99 ./.左右滑動查看全部內(nèi)容。
從上面的實例可以看出，在合并計算分型質量時，vcftools會對統(tǒng)計結果進行平均取值，bcftools則取其中的最大值輸出，而gatk會取最小值輸出，并且gatk和vcftools在輸出原有數(shù)據(jù)基礎之上，還會重新計算AC、AN等指標。
對于vcf數(shù)據(jù)合并，由于不同樣本發(fā)生的變異位置很難保證一致，所以在單獨合并不同樣本的數(shù)據(jù)時，合并后的結果往往具有很高的缺失率，后續(xù)的差異分型實質上只是對不同樣本的共有位點進行分析，丟失了某些樣本或群體的特異位點。為了避免這種情況，多樣本差異分析的項目最好是用joint calling進行變異檢測，能獲得更多的位點。
3. 參考資料
[1] https://gatk.broadinstitute.org/hc/en-us/articles/360046221931-GatherVcfs-Picard-
[2] http://samtools.github.io/bcftools/bcftools.html#concat
[3] https://gatk.broadinstitute.org/hc/en-us/articles/360045800732-MergeVcfs-Picard-
[4] https://vcftools.github.io/man_latest.html
總結
以上是生活随笔為你收集整理的变异检测：vcf文件合并的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： hdoj4648(Magic Pen 6
- 下一篇： android开发接收验证码