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
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.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
针对这个情况,并集文件的获取需要有两个文件的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的最小片段。