在群体水平分析DNA甲基化数据或转录组数据时，需要注意的一点是不同的材料/个体/基因型对于我们选定的参考基因组，在reads比对上有怎样的bias，以及这种bias对于后续差异分析有何影响。
不同基因型之间是存在大量变异的，可通过vcf文件获取。我们目前最经常考虑的变异是SNP，所以在群体水平比对DNA甲基化数据或转录组数据时，往往先利用每种基因型与所选参考基因组之间的SNP数据生成一个“伪参考基因组”，把数据比对到伪参考基因组上有助于提高reads的有效比对率。
Fig 1
在大豆群体DNA甲基化文章中（前一篇推送）作者在方法部分提到了伪参考基因组（Fig 1）。
如何生成伪参考基因组
生成伪参考基因组需要两步：
（1）拿到该个体/基因型的变异信息（vcf文件）
（2）将变异信息替换掉原参考基因组的相应位点
Fig 2
对于第一步，如果我们手中是重测序数据或者转录组数据，比对后利用gatk或者samtools/bcftools都可以拿到相应的vcf文件。
如果手中只有DNA甲基化数据，可以利用中国农大郭伟龙老师开发的DNA甲基化数据比对软件BS-Seeker2和配套的下游分析软件cgmaptools（Fig 2）拿到vcf文件。BS-Seeker2+cgmaptools也是我认为非常简单有效的DNA甲基化数据分析的上有流程，两者更丰富的操作可参考https://www.docin.com/p-2195164655.html和https://cgmaptools.github.io/
第一步：DNA甲基化测序fq数据转为Cgmap和ATCGmap格式
Fig 3
Fig 4
现在我们有A、B、C三个个体的DNA甲基化数据，利用BS-Seeker2比对到同一个ref上，可以得到A.bam、B.bam和C.bam，利用cgmaptools convert（Fig3-4）可以将bam转为Cgmap和ATCGmap格式。
需要注意的是，这一次比对的目的只是帮助我们拿到变异信息，并不用于下游的大量分析，所以可以设置较高的错配率以及更宽松的比对阈值。
（2）Call SNP
Fig 5
Fig 6
Cgmaptools 的snv选项通过输入ATCGmap格式的文件检测变异位点（Fig 5），可输出vcf文件，输出文件中Y代表T或C，R代表A或G（Fig 6）。
（3）利用bcftools替换参考基因组的SNP
无论是用重测序数据，还是转录组数据还是甲基化数据，都能拿到SNP位点的vcf文件，下一步需要利用vcf文件替换原参考基因组。
Fig 7
Fig 8
bcftools consensus（Fig 7-8）可以根据输入的vcf文件对参考基因组生成伪参考基因组。
-H 可以指定进行怎样的替换。
本文使用 文章同步助手 同步
网友评论