C2-关于VCF文件合并的几种方法

1. 个体相同，位点累加，相当于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. 位点相同，个体累加，相当于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. 位点不同，个体也不同，取两个VCF文件的交集

3.1 使用bedtools进行操作

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

结果在isec_output这个文件夹中，这里面有4个文件

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

之后得到两个VCF共有位点的SNP

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

3.2 使用bedops进行操作

这个是得到两个文件交集的bed文件，region文件

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

这个是得到两个文件交集的bed文件，位点文件

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

这个是为了得到两个文件的交集在A中的变异和在B中的变异，这是个体文件

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. 位点不同，个体也不同，取两个VCF文件的并集

针对这个情况，并集文件的获取需要有两个文件的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

得到这两个文件之后，之后再进行mpileup，得到了每个文件的mpileup格式的文件之后再合并即可

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

#! /bin/bash
#使用mpileup命令生成vcf文件
#这个示例中只对1号染色体进行了处理
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 用于降低比对质量的系数，如果reads中含有过多的错配。不能设置为零。BWA推荐值为50。

-A --count-orphans 在检测变异中，不忽略异常的reads对。

-I –positions FILE BED文件或者包含区域位点的位置列表文件。位置文件包含两列，染色体和位置，从1开始计数。BED文件至少包含3列，染色体、开始位置和结束位置，开始端从0开始计数。

-r –region STR 只在指定区域产生pileup，需要已建立索引的bam文件。通常和-l参数一起使用。

-q --min-MQ 要使用的对齐的最小映射质量。

-f --fasta-ref FASTA格式的fadix索引引用文件。文件可以选择使用bgzip压缩。

-o –output FILE 生成pileup格式文件或者VCF、BCF文件而不是默认的标准输出。

-g –BCF 计算基因型的似然值和输出文件格式为BCF。

-v –VCF 计算基因型的似然值和输出文件格式为VCF。

-D 输出每个样本的reads深度。

-V 输出每个样本未比对到参考基因组的reads数量。

-t –output-tags LIST设置FORMAT和INFO的列表内容，以逗号分割。

-u –uncompressed 生成未压缩的VCF和BCF文件。

-I –skip-indel 不检测INDEL。

-m –min-ireads INT 候选INDEL的最小间隔的reads。

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