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的简单使用，了解如何使用等

比对结果文件预处理

在用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 
```

还是用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

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
```
Indel Realignment

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

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

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

命令同Indel Realignment可参考WGS的方法
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
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才行。。。。