目標(biāo)
通過varscan等找到somatic mutations之后，想看多個(gè)樣本的突變之間的交集。bcftools是用來處理vcf文件，可以進(jìn)行各種相關(guān)的處理和分析。
軟件的安裝可以使用conda，可以參考我往期的教程：http://www.itdecent.cn/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 對(duì)vcf進(jìn)行排序
vcf里面的突變?nèi)绻麤]有按照染色體的位置排序，在后續(xù)的生成index的時(shí)候會(huì)報(bào)錯(cuò)，因此需要先將生成的突變按照染色體的位置進(jìn)行排序
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 對(duì)vcf建立索引
需要建立vcf的文件的索引，用于后續(xù)分析
bcftools index ${out}.vcf.gz
4、bcftools isec 計(jì)算不同vcf之間的突變的交集
這里假設(shè)有3個(gè)文件，分別是 Breast_1.vcf.gz， Breast_2.vcf.gz和 Breast_3.vcf.gz，我們想知道每一個(gè)突變?cè)谶@3個(gè)樣本中是否出現(xiàn)
#-n+1: 在3個(gè)樣本中至少出現(xiàn)一次突變的樣本
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這個(gè)文件夾下面，將會(huì)生成site.txt這個(gè)文件。
head site.txt
chr1 14717 G A 101
chr1 14907 A G 111
chr1 14930 A G 101
這里面的第一行中 101 表示突變?cè)诘谝粋€(gè)和第三個(gè)文件中出現(xiàn)，但是不在第二個(gè)文件中出現(xiàn)。1和0分別表示突變有沒有在對(duì)應(yīng)的文件中出現(xiàn)。
5、bcftools merge 將所有的突變合并成一個(gè)文件
#合并所有的突變到一個(gè)文件中
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 根據(jù)步驟4中計(jì)算的交集情況，提取相應(yīng)的突變
bcftools index concate.vcf.gz
#根據(jù)sites.txt，提取對(duì)應(yīng)的突變
bcftools view -T merge/sites.txt concate.vcf.gz -O z -o concate.select.vcf.gz
#檢查突變的數(shù)目和sites.txt中的是否一致
less -S concate.select.vcf.gz|wc -l
#731
wc -l merge/sites.txt
#731
