最近做课题，要做pse-genome。什么意思呢，就是要把某个个体call出来的snp替换掉参考基因组对应位置的碱基。这个刚开始接触的时候，头皮发麻，不知道咋入手。后来发现用脚本应该挺容易实现，问了一些大神，用perl的substr功能就可以轻松做到。不过今天给大家安利用GATK软件来做，是在摸索snp calling的时候偶然发现的，网上也没有找到类似的教程。话不多说，开始！
GATK，官方认可的call snp的软件，比samtools方法要慢一些，但效果好，假阳性低。对于没有已知变异位点的种类，尤其是植物来说，结合GATK和samtools方法初步call出来的snp用作known site，然后做BQSR，绝对是降低假阳性的不二利器。
今天我们要用到的是GATK中的FastaAlternateReferenceMaker功能。到底有多简单？只需要准备snp的vcf文件即可，最好是GATK call出来的，毕竟官方认可。然后输入下面的一行命令(仅限于4.0以下的版本，4.0及以上不适用)：
java -Xms10g -Xmx10g -jar /software_path/GenomeAnalysisTK.jar -T FastaAlternateReferenceMaker -R /ref_genome_path/ref_genome.fa -o pse_genome.fa --variant GATK.SNP.vcf
是不是很简单？也可以对单个基因或一段区域进行替换，需要加个参数，如下：
java -Xms10g -Xmx10g -jar /software_path/GenomeAnalysisTK.jar -T FastaAlternateReferenceMaker -R /ref_genome_path/ref_genome.fa -L Ch12:1000-3000 -o pse_genome.fa --variant GATK.SNP.vcf
如此，就构建好了，构建目的其一可能就是为了排除突变带来的影响，在群体遗传以及一些甲基化分析的时候可能会用到，一般也要注意一些参数，比如需要纯和位点，也需要考虑下突变频率，GT，DP等等。
此外，indel并不适用，如果vcf文件二者都包含，软件会自动忽略indel。
网友评论