1. 個(gè)體相同，位點(diǎn)累加，相當(dāng)于cat文件

vcf-concat A.vcf.gz B.vcf.gz C.vcf.gz | gzip -c > out.vcf.gz

vcf-concat *.vcf.gz | gzip -c > out.vcf.gz

2. 位點(diǎn)相同，個(gè)體累加，相當(dāng)于paste文件

bcftools merge file1.vcf.gz fle2.vcf.gz file3.vcf.gz > out.vcf

bcftools merge file1.vcf.gz fle2.vcf.gz file3.vcf.gz -o out.vcf

3. 位點(diǎn)不同，個(gè)體也不同，取兩個(gè)VCF文件的交集

3.1 使用bedtools進(jìn)行操作

grep "#" A.vcf > header.txt grep -v "#" A.vcf | sed 's/Chr1/1/g' > temp.txt cat header.txt temp.txt > A_new.vcf bcftools isec -p isec_output -Oz A_new.vcf.gz B.vcf.gz

結(jié)果在isec_output這個(gè)文件夾中，這里面有4個(gè)文件

1.isec_output/0000.vcf.gz would be variants unique to 1.vcf.gz 2.isec_output/0001.vcf.gz would be variants unique to 2.vcf.gz 3.isec_output/0002.vcf.gz would be variants shared by 1.vcf.gz and 2.vcf.gz as represented in 1.vcf.gz 4.isec_output/0003.vcf.gz would be variants shared by 1.vcf.gz and 2.vcf.gz as represented in 2.vcf.gz

之后得到兩個(gè)VCF共有位點(diǎn)的SNP

bcftools merge --merge all 0002.vcf.gz 0003.vcf.gz > merged.vcf

3.2 使用bedops進(jìn)行操作

這個(gè)是得到兩個(gè)文件交集的bed文件，region文件

bedops --intersect <(vcf2bed < A.vcf) <(vcf2bed < B.vcf) > answer.bed

這個(gè)是得到兩個(gè)文件交集的bed文件，位點(diǎn)文件

bedops --intersect <(vcf2bed < A.vcf) <(vcf2bed < B.vcf) > common-regions.bed bedops --everything <(vcf2bed < A.vcf) <(vcf2bed < B.vcf) > all-elements.bed bedops --element-of 1 all-elements.bed common-regions.bed > common-elements.bed

這個(gè)是為了得到兩個(gè)文件的交集在A中的變異和在B中的變異，這是個(gè)體文件

bedops --element-of 1 <(vcf2bed < A.vcf) <(vcf2bed < B.vcf) > answer1.bed bedops --element-of 1 <(vcf2bed <A.vcf) <(vcf2bed < B.vcf) > answer2.bed

4. 位點(diǎn)不同，個(gè)體也不同，取兩個(gè)VCF文件的并集

針對(duì)這個(gè)情況，并集文件的獲取需要有兩個(gè)文件的bam文件

grep -v "#" A.vcf | cut -f 1,2 > pos1.txt

grep -v "#" B.vcf | cut -f 1,2 > pos2.txt

cat pos1.txt pos2.txt > posAll.txt

得到這兩個(gè)文件之后，之后再進(jìn)行mpileup，得到了每個(gè)文件的mpileup格式的文件之后再合并即可

samtools mpileup -A -B -q 20 -Q 20 -f ref.fa bamfile.bam -l posAll.txt -r -o

#! /bin/bash #使用mpileup命令生成vcf文件 #這個(gè)示例中只對(duì)1號(hào)染色體進(jìn)行了處理 echo "SamtoolsMpileupByChr Begin: " `date` && \ samtools mpileup \ -l chr1Region.bed \ -r 1 \ -q 1 \ -C 50 \ -t DP,DV \ -m 2 \ -F 0.002 \ -f \ human.fasta \ test_3.bam \ --output test.chr1.raw.vcf && \ echo "SamtoolsMpileupByChr End: " `date`

mpileup的命令解釋

-C --adjust-MQ INT 用于降低比對(duì)質(zhì)量的系數(shù)，如果reads中含有過多的錯(cuò)配。不能設(shè)置為零。BWA推薦值為50。

-A --count-orphans 在檢測(cè)變異中，不忽略異常的reads對(duì)。

-I –positions FILE BED文件或者包含區(qū)域位點(diǎn)的位置列表文件。位置文件包含兩列，染色體和位置，從1開始計(jì)數(shù)。BED文件至少包含3列，染色體、開始位置和結(jié)束位置，開始端從0開始計(jì)數(shù)。

-r –region STR 只在指定區(qū)域產(chǎn)生pileup，需要已建立索引的bam文件。通常和-l參數(shù)一起使用。

-q --min-MQ 要使用的對(duì)齊的最小映射質(zhì)量。

-f --fasta-ref FASTA格式的fadix索引引用文件。文件可以選擇使用bgzip壓縮。

-o –output FILE 生成pileup格式文件或者VCF、BCF文件而不是默認(rèn)的標(biāo)準(zhǔn)輸出。

-g –BCF 計(jì)算基因型的似然值和輸出文件格式為BCF。

-v –VCF 計(jì)算基因型的似然值和輸出文件格式為VCF。

-D 輸出每個(gè)樣本的reads深度。

-V 輸出每個(gè)樣本未比對(duì)到參考基因組的reads數(shù)量。

-t –output-tags LIST設(shè)置FORMAT和INFO的列表內(nèi)容，以逗號(hào)分割。

-u –uncompressed 生成未壓縮的VCF和BCF文件。

-I –skip-indel 不檢測(cè)INDEL。

-m –min-ireads INT 候選INDEL的最小間隔的reads。

-F –gap-frac FLOAT 含有間隔reads的最小片段。

總結(jié)

以上是生活随笔為你收集整理的C2-关于VCF文件合并的几种方法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。