RNA-seq 检测变异之 GATK 最佳实践流程

RNA-seq 序列比对

对 RNA-seq 产出的数据进行变异检测分析，与常规重测序的主要区别就在序列比对这一步，因为 RNA-seq 的数据是来自转录本的，比对到参考基因组需要跨越转录剪切位点，所以 RNA-seq 进行变异检测的重点就在于跨剪切位点的精确序列比对。

文献  中对 RNA-seq 数据常用的 11 款比对软件进行了详细测试，包括 STAR 2-pass，而 GATK 对 RNA-seq 数据变异检测的最佳实践流程中选用了 STAR 2-pass 这一方法进行比对，STAR 发表的文章至今已被引用 1900 余次，这款软件的比对速度很快，也是  项目的御用比对软件。

STAR 2-pass 模式需要进行两次序列比对，建立两次参考基因组索引。它的思路是第一次建参考基因组索引之后进行初步的序列比对，根据初步比对结果得到该样本所有的剪切位点信息，包括参考基因组注释 GTF 中已知的剪切位点和比对时新发现的剪切位点，然后利用第一次比对得到的剪切位点信息重新对参考基因组建立索引，然后进行第二次的序列比对，这样可以得到更精确的比对结果。

---

这里使用了一个测试数据演示流程

• 第一次对参考基因组建索引：

1. # star 1-pass index

2. STAR --runThreadN 8 --runMode genomeGenerate \

3.        --genomeDir ./star_index/ \

4.        --genomeFastaFiles ./genome/chrX.fa \

5.        --sjdbGTFfile ./genes/chrX.gtf

• 然后进行第一次序列比对：

1. #star 1-pass align

2. STAR --runThreadN 8 --genomeDir ./star_index/ \

3.        --readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \

4.        --readFilesCommand zcat \

5.        --outFileNamePrefix ./star_1pass/ERR188044

• 之后根据第一次比对得到的所有剪切位点，重新对参考基因组建立索引：

1. # star 2-pass index

2. STAR --runThreadN 8 --runMode genomeGenerate \

3.        --genomeDir ./star_index_2pass/ \

4.        --genomeFastaFiles ./genome/chrX.fa \

5.        --sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab

• 再进行 STAR 二次序列比对：

1. # star 2-pass align

2. STAR --runThreadN 8 --genomeDir ./star_index_2pass/ \

3.        --readFilesIn ./samples/ERR188044_chrX_1.fastq.gz ./samples/ERR188044_chrX_2.fastq.gz \

4.        --readFilesCommand zcat \

5.        --outFileNamePrefix ./star_2pass/ERR188044

由于后面要用 GATK 进行 call 变异，还需要对比对结果 SAM 文件进行一些处理，这些都可以用 picard 来做，包括 SAM 头文件添加 @RG 标签，SAM 文件排序并转 BAM 格式，然后标记 duplicate：

1. # picard Add read groups, sort, mark duplicates, and create index

2. java -jar picard.jar AddOrReplaceReadGroups \

3.        I=./star_2pass/ERR188044Aligned.out.sam \

4.        O=./star_2pass/ERR188044_rg_added_sorted.bam \

5.        SO=coordinate \

6.        RGID=ERR188044 \

7.        RGLB=rna \

8.        RGPL=illumina \

9.        RGPU=hiseq \

10.        RGSM=ERR188044

11. java -jar picard.jar MarkDuplicates \

12.        I=./star_2pass/ERR188044_rg_added_sorted.bam \

13.        O=./star_2pass/ERR188044_dedup.bam  \

14.        CREATE_INDEX=true \

15.        VALIDATION_STRINGENCY=SILENT \

16.        M=./star_2pass/ERR188044_dedup.metrics

到此序列比对就完成了。

---

使用 GATK 进行变异检测

感觉 GATK 里面的工具都很慢（相对于其他的软件特别慢！），都是单线程在跑，有的虽然可以设置为多线程但是感觉没啥速度上的提升，所以要有点耐心……

由于 STAR 软件使用的 MAPQ 标准与 GATK 不同，而且比对时会有 reads 的片段落到内含子区间，需要进行一步 MAPQ 同步和 reads 剪切，使用 GATK 专为 RNA-seq 应用开发的工具 SplitNCigarReads 进行操作，它会将落在内含子区间的 reads 片段直接切除，并对 MAPQ 进行调整。

DNA 测序的重测序应用中也有序列比对软件的 MAPQ 与 GATK 无法直接对接的情况，需要进行调整。

1. # samtools faidx chrX.fa

2. # samtools dict chrX.fa

3. java -jar GenomeAnalysisTK.jar -T SplitNCigarReads \

4.        -R ./genome/chrX.fa \

5.        -I ./star_2pass/ERR188044_dedup.bam \

6.        -o ./star_2pass/ERR188044_dedup_split.bam \

7.        -rf ReassignOneMappingQuality \

8.        -RMQF 255 \

9.        -RMQT 60 \

10.        -U ALLOW_N_CIGAR_READS

之后就是可选的 Indel Realignment，对已知的 indel 区域附近的 reads 重新比对，可以稍微提高 indel 检测的真阳性率，如果时间紧张也可以不做，这一步影响很轻微 :)

1. # 可选步骤 IndelRealign

2. java -jar GenomeAnalysisTK.jar -T RealignerTargetCreator \

3.        -R ./genome/chrX.fa \

4.        -I ./star_2pass/ERR188044_dedup_split.bam \

5.        -o ./star_2pass/ERR188044_realign_interval.list \

6.        -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

7. java -jar GenomeAnalysisTK.jar -T IndelRealigner \

8.        -R ./genome/chrX.fa \

9.        -I ./star_2pass/ERR188044_dedup_split.bam \

10.        -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \

11.        -o ./star_2pass/ERR188044_realign.bam \

12.        -targetIntervals ./star_2pass/ERR188044_realign_interval.list

然后还是可选的 BQSR，这一步操作主要是针对测序质量不太好的数据，进行碱基质量再校准，如果对自己的测序数据质量足够自信可以省略，2500 之后 Hiseq 仪器的数据质量都挺不错的，可以根据 FastQC 结果来决定。这一步省了又能节省时间 :)

1. # 可选步骤 BQSR

2. java -jar GenomeAnalysisTK.jar \

3.        -T BaseRecalibrator \

4.        -R ./genome/chrX.fa \

5.        -I ./star_2pass/ERR188044_realign.bam \

6.        -knownSites 1000G_phase1.snps.high_confidence.hg19.sites.vcf \

7.        -knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \

8.        -o ./star_2pass/ERR188044_recal_data.table

9. java -jar GenomeAnalysisTK.jar  \

10.        -T PrintReads \

11.        -R ./genome/chrX.fa \

12.        -I ./star_2pass/ERR188044_realign.bam \

13.        -BQSR ./star_2pass/ERR188044_recal_data.table \

14.        -o ./star_2pass/ERR188044_BQSR.bam

比如下面的数据就可以放心的省略这两步了：

现在终于可以进行变异检测了，GATK 官网说 HC 表现比 UC 好，所以这里用 HC 进行变异检测：

1. java -jar GenomeAnalysisTK.jar -T HaplotypeCaller \

2.        -R ./genome/chrX.fa \

3.        -I ./star_2pass/ERR188044_dedup_split.bam \

4.        -dontUseSoftClippedBases \

5.        -stand_call_conf 20.0 \

6.        -o ./star_2pass/ERR188044.vcf

call 完变异之后再进行过滤：

1. java -jar GenomeAnalysisTK.jar \

2.        -T VariantFiltration \

3.        -R ./genome/chrX.fa \

4.        -V ./star_2pass/ERR188044.vcf \

5.        -window 35 \

6.        -cluster 3 \

7.        -filterName FS -filter "FS > 30.0" \

8.        -filterName QD -filter "QD < 2.0" \

9.        -o ./star_2pass/ERR188044_filtered.vcf

然后就拿到变异检测结果了，可以用 ANNOVAR 或 SnpEff 或 VEP 进行注释，根据自己的需要进行筛选了。

---

欢迎点击文末阅读原文，阅读作者相关其它文章。