直接来看例子吧。
- 左对齐标准化Indel，对于多等位基因位点进行拆分（annovar注释必须的）。
#首先需要压缩VCF并建立索引。
bgzip -f "$outDIR"_tmp/"$out_basename".vcf -@ 10
tabix -p vcf "$outDIR"_tmp/"$out_basename".vcf.gz
bcftools norm --fasta-ref "$fasta_file" --multiallelics -both --output "$outDIR"_norm/"$out_basename" --output-type z "$inputVCF"
- 根据条件进行筛选，比如：筛选出FILTER列是PASS的，DP >20 , GQ >100, QUAL >100的位点：
bcftools filter -i ' FILTER=="PASS" && DP>20 && GQ>100 && QUAL>100' "$outDIR"_norm/"$out_basename" --output "$outDIR"/"$out_basename" --output-type z
- 提取 突变位点的AD和DP，顺便可以计算 VAF（习惯使用awk，所以下面的都是awk了..）：
bcftools query -f '[%AD]\n' "$outDIR"/"$out_basename" |awk 'BEGIN{FS=","}{ print $2 }' > > "$outDIR"_tmp/"$out_basename"_AD_tmp.txt
bcftools query -f '[%DP]\n' "$outDIR"/"$out_basename" > "$outDIR"_tmp/"$out_basename"_DP_tmp.txt
# 计算VAF:
awk 'BEGIN{OFS="\t"}{ if(NR==FNR)AD[NR]=$1; if(NR>FNR)DP[FNR]=$1; }END{ for(i=1;i<=length(AD);i++){ print AD[i],DP[i],AD[i]/DP[i] } }' "$outDIR"_tmp/"$out_basename"_AD_tmp.txt "$outDIR"_tmp/"$out_basename"_DP_tmp.txt > "$outDIR"_tmp/"$out_basename"_AD_DP_VAF_tmp.txt
# 过滤VAF(获得VAF > 0.3的位点)
num_anno=`zcat "$outDIR"/"$out_basename" |awk '{if($0~"^#")print 1 }' |wc -l`
zcat "$outDIR"/"$out_basename" |awk -v num_anno="$num_anno" 'BEGIN{FS="\t"}{ if(NR==FNR && $0~"^#"){print $0}; if(NR==FNR && $0!~"^#")vcf[NR-num_anno]=$0; if(NR>FNR && $3>0.3)print vcf[FNR] }' - "$outDIR"_tmp/"$out_basename"_AD_DP_VAF_tmp.txt > "$outDIR"_tmp/"$out_basename"_VAF_filtered.vcf
- 合并多个样本的VCF文件, 注意需要每个文件的样本名唯一，如果不唯一使用--force-samples 将自动重命名。
# 定义用于存储待合并的vcf文件路径：
combine_vars=""
for((i=0;i<${#filtered_vcf[*]};i++)){
#为每个文件建立索引，如果没压缩要先压缩
tabix -p vcf "${filtered_vcf[$i]}"
echo "完成了第"$i"个"
combine_vars="${combine_vars}"" ""${filtered_vcf[$i]}"
}
# --missing-to-ref表示缺失的GT表示为0/0，-merge操作多等位基因位点，这里表示不产生多等位基因位点。
# 当然对于header使用--use-header指定，info列的合并使用--info-rules指定规则。
bcftools merge --missing-to-ref --merge none --output "$workdir"/test/combined.vcf.gz --output-type z --threads 30
- 上面是合并多个样本，如果是相同样本的位点合并呢，比如一个样本的SNP和INDEL进行合并，首先必须对待合并文件进行排序：
# 排序位点：
bcftools sort SNP_filtered.vcf -O z -o SNP_filtered_sorted.vcf.gz
bcftools sort INDEL_filtered.vcf -O z -o INDEL_filtered_sorted.vcf.gz
# 合并：
bcftools concat SNP_filtered_sorted.vcf.gzINDEL_filtered_sorted.vcf.gz -a -O z -o ALL_filtered_sorted.vcf.gz
- 提取指定染色体上的位点：
bcftools filter -t chr1,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr2,chr20,chr21,chr22,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chrM,chrX,chrY "$workdir"/test/combined_split.vcf.gz --output "$workdir"/test/combined_split_chr.vcf.gz --output-type z
- 移除INFO和FORMAT中除了GT和PL的列：
bcftools annotate -x INFO,^FORMAT/GT,FORMAT/PL file.vcf
- 使用 src.bcf来注释 dst.bcf，只注释ID，QUAL和TAG，如果TAG存在则不覆盖。
bcftools annotate -a src.bcf -c ID,QUAL,+TAG dst.bcf
- 除了FORMAT的GT列外，注释所有的INFO和FORMAT.
bcftools annotate -a src.bcf -c INFO,^FORMAT/GT dst.bcf
- 使用TAB分割的文件进行注释VCF(1-bae)：
# 需要 1-base的坐标系并且建立索引：
tabix -s1 -b2 -e2 annots.tab.gz
bcftools annotate -a annots.tab.gz -h annots.hdr -c CHROM,POS,REF,ALT,-,TAG file.vcf
bcftools annotate -a annots.tab.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf
- 使用bed文件进行注释(0-base)：
bcftools annotate -a annots.bed.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf
- 常用查询命令：
# 输出染色体、位置、REF、ALT:
bcftools query -f '%CHROM %POS %REF %ALT{0}\n' file.vcf.gz
# 还是输出上面的结果，但用\t代替空格并输出样本名和基因型
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' file.vcf.gz
# 输出GQ和GT:
bcftools query -f 'GQ:[ %GQ] \t GT:[ %GT]\n' file.vcf
# 创建bed文件: chr, pos (0-based), end pos (1-based), id
bcftools query -f'%CHROM\t%POS0\t%END\t%ID\n' file.bcf
# 输出样本的突变位点信息和GT：
bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]' -i'GT="alt"' file.bcf
更多原创精彩视频敬请关注生信杂谈:
网友评论