首页 > 代码库 > 结合GATK和samtools以及picardtools call snp
结合GATK和samtools以及picardtools call snp
刚开始学生物信息学,老师给了个以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
结合GATK和samtools以及picardtools call snp