bcftools 处理vcf文件,寻找多个vcf文件中突变的交集
目标
通过varscan等找到somatic mutations之后,想看多个样本的突变之间的交集。bcftools是用来处理vcf文件,可以进行各种相关的处理和分析。
软件的安装可以使用conda,可以参考我往期的教程:https://www.jianshu.com/p/e82a8d799b13
1、bcftools reheader 重新命名vcf的头文件
input="input vcf"
out="output name"
bcftools reheader -h vcf.header.txt ${input}.vcf -o ${out}.vcf
-h vcf.header.txt:新的vcf头文件,用于替代旧的vcf头文件
2、bcftools sort 对vcf进行排序
vcf里面的突变如果没有按照染色体的位置排序,在后续的生成index的时候会报错,因此需要先将生成的突变按照染色体的位置进行排序
bcftools sort ${out}.vcf -o ${out}.vcf.gz -O z
-O, --output-type <b|u|z|v> b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]
3、bcftools index 对vcf建立索引
需要建立vcf的文件的索引,用于后续分析
bcftools index ${out}.vcf.gz
4、bcftools isec 计算不同vcf之间的突变的交集
这里假设有3个文件,分别是 Breast_1.vcf.gz, Breast_2.vcf.gz和 Breast_3.vcf.gz,我们想知道每一个突变在这3个样本中是否出现
#-n+1: 在3个样本中至少出现一次突变的样本
bcftools isec -n+1 -c all -p merge Breast_1.vcf.gz Breast_2.vcf.gz Breast_3.vcf.gz
#提取在Breast_1和Breast_2中存在,而在Breast_3中不存在的突变
bcftools isec -n~110 -c all -p merge Breast_1.vcf.gz Breast_2.vcf.gz Breast_3.vcf.gz
-n, --nfiles [+-=~]<int> output positions present in this many (=), this many or more (+), this many or fewer (-), the exact (~) files
-p, --prefix <dir> if given, subset each of the input files accordingly, see also -w
在merge这个文件夹下面,将会生成site.txt这个文件。
head site.txt
chr1 14717 G A 101
chr1 14907 A G 111
chr1 14930 A G 101
这里面的第一行中 101 表示突变在第一个和第三个文件中出现,但是不在第二个文件中出现。1和0分别表示突变有没有在对应的文件中出现。
5、bcftools merge 将所有的突变合并成一个文件
#合并所有的突变到一个文件中
bcftools merge --force-samples -o concate.vcf.gz -O z Breast_1.vcf.gz Breast_2.vcf.gz Breast_3.vcf.gz
6、bcftools view 根据步骤4中计算的交集情况,提取相应的突变
bcftools index concate.vcf.gz
#根据sites.txt,提取对应的突变
bcftools view -T merge/sites.txt concate.vcf.gz -O z -o concate.select.vcf.gz
#检查突变的数目和sites.txt中的是否一致
less -S concate.select.vcf.gz|wc -l
#731
wc -l merge/sites.txt
#731