RNA-seq 檢測變異之 GATK 最佳實踐流程

RNA-seq 序列比對

對 RNA-seq 產出的數據進行變異檢測分析，與常規重測序的主要區別就在序列比對這一步，因為 RNA-seq 的數據是來自轉錄本的，比對到參考基因組需要跨越轉錄剪切位點，所以 RNA-seq 進行變異檢測的重點就在於跨剪切位點的精確序列比對。

文獻 systematic evaluation of spliced alignment programs for RNA-seq data 中對 RNA-seq 數據常用的 11 款比對軟體進行了詳細測試，包括 STAR 2-pass，而 GATK 對 RNA-seq 數據變異檢測的最佳實踐流程中選用了 STAR 2-pass 這一方法進行比對，STAR 發表的文章至今已被引用 1900 餘次，這款軟體的比對速度很快，也是 ENCODE 項目的御用比對軟體。

STAR 2-pass 模式需要進行兩次序列比對，建立兩次參考基因組索引。它的思路是第一次建參考基因組索引之後進行初步的序列比對，根據初步比對結果得到該樣本所有的剪切位點信息，包括參考基因組注釋 GTF 中已知的剪切位點和比對時新發現的剪切位點，然後利用第一次比對得到的剪切位點信息重新對參考基因組建立索引，然後進行第二次的序列比對，這樣可以得到更精確的比對結果。

這裡使用了一個測試數據演示流程

# star 1-pass index

STAR --runThreadN 8 --runMode genomeGenerate \

       --genomeDir ./star_index/ \

       --genomeFastaFiles ./genome/chrX.fa \

       --sjdbGTFfile ./genes/chrX.gtf

#star 1-pass align

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

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

       --readFilesCommand zcat \

       --outFileNamePrefix ./star_1pass/ERR188044

# star 2-pass index

STAR --runThreadN 8 --runMode genomeGenerate \

       --genomeDir ./star_index_2pass/ \

       --genomeFastaFiles ./genome/chrX.fa \

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

# star 2-pass align

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

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

       --readFilesCommand zcat \

       --outFileNamePrefix ./star_2pass/ERR188044

由於後面要用 GATK 進行 call 變異，還需要對比對結果 SAM 文件進行一些處理，這些都可以用 picard 來做，包括 SAM 頭文件添加 @RG 標籤，SAM 文件排序並轉 BAM 格式，然後標記 duplicate：

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

java -jar picard.jar AddOrReplaceReadGroups \

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

       O=./star_2pass/ERR188044_rg_added_sorted.bam \

       SO=coordinate \

       RGID=ERR188044 \

       RGLB=rna \

       RGPL=illumina \

       RGPU=hiseq \

       RGSM=ERR188044

java -jar picard.jar MarkDuplicates \

       I=./star_2pass/ERR188044_rg_added_sorted.bam \

       O=./star_2pass/ERR188044_dedup.bam  \

       CREATE_INDEX=true \

       VALIDATION_STRINGENCY=SILENT \

       M=./star_2pass/ERR188044_dedup.metrics

到此序列比對就完成了。

使用 GATK 進行變異檢測

感覺 GATK 裡面的工具都很慢（相對於其他的軟體特別慢！），都是單線程在跑，有的雖然可以設置為多線程但是感覺沒啥速度上的提升，所以要有點耐心……

由於 STAR 軟體使用的 MAPQ 標準與 GATK 不同，而且比對時會有 reads 的片段落到內含子區間，需要進行一步 MAPQ 同步和 reads 剪切，使用 GATK 專為 RNA-seq 應用開發的工具 SplitNCigarReads 進行操作，它會將落在內含子區間的 reads 片段直接切除，並對 MAPQ 進行調整。

DNA 測序的重測序應用中也有序列比對軟體的 MAPQ 與 GATK 無法直接對接的情況，需要進行調整。

# samtools faidx chrX.fa

# samtools dict chrX.fa

java -jar GenomeAnalysisTK.jar -T SplitNCigarReads \

       -R ./genome/chrX.fa \

       -I ./star_2pass/ERR188044_dedup.bam \

       -o ./star_2pass/ERR188044_dedup_split.bam \

       -rf ReassignOneMappingQuality \

       -RMQF 255 \

       -RMQT 60 \

       -U ALLOW_N_CIGAR_READS

之後就是可選的 Indel Realignment，對已知的 indel 區域附近的 reads 重新比對，可以稍微提高 indel 檢測的真陽性率，如果時間緊張也可以不做，這一步影響很輕微 :)

# 可選步驟 IndelRealign

java -jar GenomeAnalysisTK.jar -T RealignerTargetCreator \

       -R ./genome/chrX.fa \

       -I ./star_2pass/ERR188044_dedup_split.bam \

       -o ./star_2pass/ERR188044_realign_interval.list \

       -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

java -jar GenomeAnalysisTK.jar -T IndelRealigner \

       -R ./genome/chrX.fa \

       -I ./star_2pass/ERR188044_dedup_split.bam \

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

       -o ./star_2pass/ERR188044_realign.bam \

       -targetIntervals ./star_2pass/ERR188044_realign_interval.list

然後還是可選的 BQSR，這一步操作主要是針對測序質量不太好的數據，進行鹼基質量再校準，如果對自己的測序數據質量足夠自信可以省略，2500 之後 Hiseq 儀器的數據質量都挺不錯的，可以根據 FastQC 結果來決定。這一步省了又能節省時間 :)

# 可選步驟 BQSR

java -jar GenomeAnalysisTK.jar \

       -T BaseRecalibrator \

       -R ./genome/chrX.fa \

       -I ./star_2pass/ERR188044_realign.bam \

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

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

       -o ./star_2pass/ERR188044_recal_data.table

java -jar GenomeAnalysisTK.jar  \

       -T PrintReads \

       -R ./genome/chrX.fa \

       -I ./star_2pass/ERR188044_realign.bam \

       -BQSR ./star_2pass/ERR188044_recal_data.table \

       -o ./star_2pass/ERR188044_BQSR.bam

比如下面的數據就可以放心的省略這兩步了：

現在終於可以進行變異檢測了，GATK 官網說 HC 表現比 UC 好，所以這裡用 HC 進行變異檢測：

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller \

       -R ./genome/chrX.fa \

       -I ./star_2pass/ERR188044_dedup_split.bam \

       -dontUseSoftClippedBases \

       -stand_call_conf 20.0 \

       -o ./star_2pass/ERR188044.vcf

call 完變異之後再進行過濾：

java -jar GenomeAnalysisTK.jar \

       -T VariantFiltration \

       -R ./genome/chrX.fa \

       -V ./star_2pass/ERR188044.vcf \

       -window 35 \

       -cluster 3 \

       -filterName FS -filter "FS > 30.0" \

       -filterName QD -filter "QD < 2.0" \

       -o ./star_2pass/ERR188044_filtered.vcf

然後就拿到變異檢測結果了，可以用 ANNOVAR 或 SnpEff 或 VEP 進行注釋，根據自己的需要進行篩選了。

歡迎點擊文末閱讀原文，閱讀作者相關其它文章。