GATK calling variants in RNA-seq

如果在WGS上call variants的话,有不少软件以及相关流程,比如有名的GATK best practices。其实GATK也有一套对于RNA-Seq相关的call variants的流程方法,粗略一看其实跟WGS的差不多,但是有一些地方还是有差别的,我以一个小鼠的公共数据为例尝试一下,参考文章Calling variants in RNAseq

数据来源

其实是之前转录组小实战的原始数据,可在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916可在下方ftp下载测序数据,小鼠的4个样本的编号分别为SRR3589959到SRR3589962,然后用sratoolkit工具将SRR文件转换为fastq格式的测序数据做接下来的分析

比对至参考基因组

对于WGS数据,GATK建议使用BWA做比对,但是RNA-Seq数据,则建议使用STAR以便对call SNP和INDEL有最佳的灵敏度。因此使用STAR的2-pass mode作为比对的首选方法,所以我在此之前对STAR做了个笔记比对软件STAR的简单使用,了解如何使用等

比对结果文件预处理

  1. 在用STAR的2-pass mode比对时,由于考虑到后续还要给bam文件添加RG标签(GATK要求的),所以就没有在比对输出时就对其先排序,反正用picard在添加RG标签时也能顺便排序(后来发现picard运行真慢),以STAR输出的一个sam文件(SRR3589959_2-passAligned.out.sam)为例:

    java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/picard/picard.jar AddOrReplaceReadGroups \
    I=SRR3589959_2-passAligned.out.sam \
    O=SRR3589959_sorted.bam \
    SO=coordinate RGID=SRR3589959 RGLB=mRNA RGPL=illumina RGPU=HiSeq2500 RGSM=SRR3589959 
    
  2. 还是用picard套件对duplicate reads进行标记,并建index;第一次才知道picard的CREATE_INDEX参数可以这样使用。。。以前都是用samtools index来建bam文件的索引

    java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/picard/picard.jar MarkDuplicates \
    I=SRR3589959_sorted.bam \
    O=SRR3589959_sorted_maked.bam \
    CREATE_INDEX=true \
    VALIDATION_STRINGENCY=SILENT \
    METRICS_FILE=SRR3589959.metrics
    

GATK call variants

  1. Split’N’Trim and reassign mapping qualities

    这一步是跟WGS有点不一样,GATK使用专门给RNA-Seq开发的SplitNCigarReads工具,将落在外显子上的reads分离出来同时去除N错误碱基(getting rid of Ns but maintaining grouping information),并去除掉落在内含子区域的reads,以减少假阳性变异;这工具还使用ReassignOneMappingQuality将STAR软件产生的比对质量标准MAPQ转化为GATK设定的标准(由于这个标准GATK是不识别的),比如将MAPQ 255转化为GATK的默认值60;最后还指定了允许reads含有N,至少GATK现在是这么设定的,这个我不太理解为何,原文(Finally, be sure to specify that reads with N cigars should be allowed. This is currently still classified as an “unsafe” option, but this classification will change to reflect the fact that this is now a supported option for RNAseq processing)

    java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/gatk/GATK/GenomeAnalysisTK.jar \
    -T SplitNCigarReads \
    -R ~/reference/genome/mm10/GRCm38.p5.genome.fa \
    -I SRR3589959_sorted_maked.bam \
    -o SRR3589959_sorted_maked_split.bam \
    -rf ReassignOneMappingQuality \
    -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS
    

    注:参考基因组序列需要.dict文件和.fai文件,可参考https://software.broadinstitute.org/gatk/documentation/article?id=1601

    java -jar ~/biosoft/picard/picard.jar CreateSequenceDictionary R=GRCm38.p5.genome.fa O=GRCm38.p5.genome.dict
    samtools faidx GRCm38.p5.genome.fa
    
  2. Indel Realignment

    这步Indel Realignment(indel局部区域重比对),就是根据你提供的indel信息,在indel区域进行重新校正,官网解释是为了防止错失一些indel变异;可能是由于比对过程的中对gap与错配偏好性造成的(简单的说明明是indel,却被误认为是snp,这样的错误),当然也是由于indel周围的比对结果本来就不太准确。

    但是GATK也说了,这步对最后的结果影响比较少(WGS的话确实如此),但是GATK还是建议做这步的,特别是如果有对应物种的可信度较高snp和Indel的参考变异文件,具体命令可参考之前写的WGS的GATK best practices流程

  3. Base Recalibration

    这步是为了重新校正碱基质量值(BQSR),其通过机器学习方法构建了测序碱基错误率模型,根据模型对测序的碱基进行相应的调整。GATK也是建议做的,但是如果你的测序数据质量较好,其实做这步的话效果并不大,而且这步可选提供对应物种已知的snp和indel变异文件,如果没有的话,我觉得是不是更没必要做这步了?

    命令同Indel Realignment可参考WGS的方法

  4. Variant calling

    这步就是用来call variants的,用的是GATK的HaplotypeCaller模块,由于这里是单个样本,所以直接进行HaplotypeCaller即可

    java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/gatk/GATK/GenomeAnalysisTK.jar \
    -T HaplotypeCaller \
    -R ~/reference/genome/mm10/GRCm38.p5.genome.fa \
    -I SRR3589959_sorted_maked_split.bam -dontUseSoftClippedBases -stand_call_conf 20 \
    -o SRR3589959.vcf
    

    -dontUseSoftClippedBases表示GATK不考虑比对时产生的soft clipped的碱基,以减少假阳性和假阴性的变异
    -stand_call_conf相当于一个可信度打分,转录的默认是20,全基因组会考虑用30

  5. Variant filtering

    这步就是用来对上面call出来的variants进行过滤的,GATK建议过滤掉在35bp范围内出现3个以上的SNP的情况(-window 35 -cluster 3),这个是针对RNA-Seq数据的;还有其他的过滤则跟WGS相同了,过滤掉Fisher Strand values(使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值,值越小越好)大于30的,过滤掉Qual By Depth values(经过序列深度标准化的SNP质量值)小于2的;其实这个过滤标准还是根据自己的情况需求来定,GATK只是给了个建议的标准。而且我比较好奇的是,这里对snp和indel都是同一个标准吗?暂时是这样了。。

    java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/gatk/GATK/GenomeAnalysisTK.jar \
    -T VariantFiltration \
    -R ~/reference/genome/mm10/GRCm38.p5.genome.fa \
    -V SRR3589959.vcf -window 35 -cluster 3 \
    -filterName FS -filter "FS > 30.0" -filterName QD -filter "QD < 2.0" \
    -o SRR3589959_filtered.vcf
    

    注:遇到一个之前没有的报错。。结果发现是在设定阈值的时候,比如FS>30就会报错。。必须是FS>30.0才行。。。。

变异注释

这个有不少软件可以做,比如ANNOVAR,snpEff,VEP等,前2个做了个小结,可翻看前面的博文;总之就是对变异在基于基因组的做了注释,然后可能还会看看在编码区是否造成了一些影响等

如果有异的地方,可以查看最开头的官网说明文档,那里是最原始最纯正的哈~

推荐2个公众号对GATK的对WGS以及RNA-Seq做call variants的流程以及讲解:

第一个是http://www.huangshujia.me/2017/09/19/2017-09-19-Begining-WGS-Data-Analysis-The-pipeline.html

还有一个是生信技能树的一篇软件(RNA-seq检测变异之GATK最佳实践流程),由于微信链接会失效,所以就不放链接了

本文出自于http://www.bioinfo-scrounger.com转载请注明出处