怎么获得共有序列（consensus sequence）

在网上大家会反复提到的一个问题是“我怎么从我的bam文件中获得共有序列（consensus)”，这个问题的翻版是“我怎么从我的bam文件中提取fasta序列”。当然，这个问题背后实际想表达的通常是获得样品中感兴趣区间的序列，然后进行多样本或多物种的序列比对。那些做过系统发育分析的人应该对这个问题很感兴趣？

上面这些问题可以通过简单的两步来解决：

进行变异检测，获得VCF文件
VCF文件与参考基因组进行比较，替换变异的位点

变异检测的软件和方法多如牛毛，后面的公众号也会分享一些目前流行的工具和流程。这里假设你用bcftools 的mpileup进行变异的检测。

$ bcftools mpileup -Ou -f ref.fa input.bam | bcftools call -Oz -mv -o output.vcf.gz
获得的VCF文件需要排序过的，并且是经过bgzip压缩和tabix 构建索引的步骤。排序可以下载vcftools，使用vcf-sort your.vcf.gz|bgzip >sort.vcf.gz命令。

压缩参考命令如下：

$ bgzip -c output.vcf > output.vcf.gz
tabix 命令如下：

$ tabix output.vcf.gz
现在可以开始进行consensus序列的获取了，这里使用bcftools工具中的consensus命令。如果你想创建一个新的基于你的变异文件的参考基因组，参考下面这条命令：

$ bcftools consensus -f ref.fa output.vcf.gz > out.fa
假如你只是对某个区间感兴趣的话，参考下面这条命令:

$ samtools faidx ref.fa 8:11870-11890 | bcftools consensus output.vcf.gz -o out.fa
当然多个区间也是可以的:

$ samtools faidx -r regions.bed ref.fa | bcftools consensus output.vcf.gz -o out.fa
还有不懂的可以后台回复。