结合GATK和samtools以及picardtools call snp

Posted 2020-08-19

刚开始学生物信息学，老师给了个以snp为标记来画遗传图的课题，研究了一段时间，开始用bwa+samtools来call snp，师姐以前用这套做过，她建议我用另外的方法来做，于是准备学下用GATK来做snp calling。

call snp首先要有个比较准确的参考基因组，然后有样本，我的样本是杂交产生的F2，下面使自己的一些使用过程和心得体会，

这套流程相对于bwa和samtool来说有所不同，先需要对sample的fasta进行筛选，我自己找了下主要是用NGSQC toolkit下的perl脚本筛选，命令为:

$ IlluQC_PRLL.pl -pe read1.fq reads2.fq 2 5 -c 8 -o QC

这个是加入了环境变量的命令

然后要用bowtie2对参考基因组建立索引文件：

$ bowtie2-build genome.fasta genome

genome为索引文件前缀，此索引文件生成目录默认为当前目录而不是ref genome所在的目录

然后将sample mapping到ref genome上命令如下：

bowtie2 --rg-id sample --rg "PL:ILLUMINA" --rg "SM:samplename" -x genome -1 IlluQC_Filtered_files/sampleR1.fastq.gz_filtered -2 IlluQC_Filtered_files/sampleR2.fastq.gz_filtered -p 8 -S sample.sam

--rg-id 是设定read group ID到sample上 -rg是加上“@RG”的头文件，这个类似于bwa mem中的 -R参数

下一步将得到的sam文件转换成二进制的bam文件，这不用到samtools，命令如下：

$ samtools  view -bS sample.sam  > sample.bam

samtools view是专门处理bam文件和sam文件的，-bS命令是输入文件格式为sam 输出文件格式为bam。

接下来是标记出PCR扩增产生的大量重复片段：

这步用到了picardHome里面的SortSam.jar

命令如下：

$ java -Xmx10g -jar picardHome/SortSam.jar INPUT=sample.bam OUTPUT=sample.sorted.bam SORT_ORDER=coordinate

排序

$ java -Xmx10g -jar picardHome/MarDuplicates.jar INPUT=sample.sorted.bam OUTPUT=sample.dd.bam METRICS_FILE=sample.dd.metrics MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

再次对genome.fasta索引

$ samtools faidx genome.fasta

$ java -Xmx10g -jar picardHome/CreateSequenceDictionary.jar R=genome.fasta O=genome.dict

$ samtools index sample.dd.bam

此处有一个需要注意的地方就是genome.fasta的前缀genome必须和genome.dict的前缀genome一样，否者下一步会报错

接下来对INDEL周围进行重新比对（realignment）：

开始用刀GATK这个软件了

java -Xmx10g -jar GenomeAnalysisTK.jar -R genome.fasta -T RealignerTargetCreator  -I sample.dd.bam -o sample.realn.intervals

java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \ -R genome.fasta -T IndelRealigner \ -targetIntervals sample.realn.intervals \ -I sample.dd.bam -o sample.realn.bam

接下来是第一遍call snp：

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \ -R genome.fasta -T UnifiedGenotyper \ -I sample.realn.bam -o sample.gatk.raw1.vcf \ --read_filter BadCigar -glm BOTH \ -stand_call_conf 30.0 -stand_emit_conf 0 

$ samtools mpileup -q 20 -Q 25 -ug -t AD -t DP -f genome.fasta sample.realn.bam | bcftool call -c -V indels -v - > sample.samtools.raw1.vcf