tags： RNA-seq GATK STAR SNP INDEL

RNA-seq 序列比對(duì)

對(duì) RNA-seq 產(chǎn)出的數(shù)據(jù)進(jìn)行變異檢測(cè)分析，與常規(guī)重測(cè)序的主要區(qū)別就在序列比對(duì)這一步，因?yàn)?RNA-seq 的數(shù)據(jù)是來自轉(zhuǎn)錄本的，比對(duì)到參考基因組需要跨越轉(zhuǎn)錄剪切位點(diǎn)，所以 RNA-seq 進(jìn)行變異檢測(cè)的重點(diǎn)就在于跨剪切位點(diǎn)的精確序列比對(duì)。

有一篇文獻(xiàn) systematic evaluation of spliced alignment programs
for RNA-seq data 對(duì) RNA-seq 數(shù)據(jù)常用的 11 款比對(duì)軟件進(jìn)行了詳細(xì)的比較測(cè)試，例如 tophat2, STAR 等。 GATK 發(fā)布的 RNA-seq 數(shù)據(jù)變異檢測(cè)最佳實(shí)踐流程用了 STAR 2-pass 這一方法進(jìn)行序列比對(duì)，STAR 發(fā)表的文章至今已被引用 1900 余次，這款軟件的比對(duì)速度很快，也是 ENCODE 項(xiàng)目的御用比對(duì)軟件。

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

這里使用了一個(gè)測(cè)試數(shù)據(jù)演示流程，第一次對(duì)參考基因組建索引：

# star 1-pass index
STAR --runThreadN 8 --runMode genomeGenerate \
        --genomeDir ./star_index/ \
        --genomeFastaFiles ./genome/chrX.fa \
        --sjdbGTFfile ./genes/chrX.gtf

然后進(jìn)行第一次序列比對(duì)：

#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

之后根據(jù)第一次比對(duì)得到的所有剪切位點(diǎn)，重新對(duì)參考基因組建立索引：

# star 2-pass index
STAR --runThreadN 8 --runMode genomeGenerate \
        --genomeDir ./star_index_2pass/ \
        --genomeFastaFiles ./genome/chrX.fa \
        --sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab

再進(jìn)行 STAR 二次序列比對(duì)：

# 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 進(jìn)行 call 變異，還需要對(duì)比對(duì)結(jié)果 SAM 文件進(jìn)行一些處理，這些都可以用 picard 來做，包括 SAM 頭文件添加 @RG 標(biāo)簽，SAM 文件排序并轉(zhuǎn) BAM 格式，然后標(biāo)記 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

到此序列比對(duì)就完成了。

使用 GATK 進(jìn)行變異檢測(cè)

感覺 GATK 里面的工具都很慢（相對(duì)于其他的軟件特別慢?。际菃尉€程在跑，有的雖然可以設(shè)置為多線程但是感覺沒啥速度上的提升，所以要有點(diǎn)耐心……

由于 STAR 軟件使用的 MAPQ 標(biāo)準(zhǔn)與 GATK 不同，而且比對(duì)時(shí)會(huì)有 reads 的片段落到內(nèi)含子區(qū)間，需要進(jìn)行一步 MAPQ 同步和 reads 剪切，使用 GATK 專為 RNA-seq 應(yīng)用開發(fā)的工具 SplitNCigarReads 進(jìn)行操縱，它會(huì)將落在內(nèi)含子區(qū)間的 reads 片段直接切除，并對(duì) MAPQ 進(jìn)行調(diào)整。DNA 測(cè)序的重測(cè)序應(yīng)用中也有序列比對(duì)軟件的 MAPQ 與 GATK 無法直接對(duì)接的情況，需要進(jìn)行調(diào)整。

# 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，對(duì)已知的 indel 區(qū)域附近的 reads 重新比對(duì)，可以稍微提高 indel 檢測(cè)的真陽(yáng)性率，如果時(shí)間緊張也可以不做，這一步影響很輕微 :)

# 可選步驟 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，這一步操作主要是針對(duì)測(cè)序質(zhì)量不太好的數(shù)據(jù)，進(jìn)行堿基質(zhì)量再校準(zhǔn)，如果對(duì)自己的測(cè)序數(shù)據(jù)質(zhì)量足夠自信可以省略，2500 之后 Hiseq 儀器的數(shù)據(jù)質(zhì)量都挺不錯(cuò)的，可以根據(jù) FastQC 結(jié)果來決定。這一步省了又能節(jié)省時(shí)間 :)

# 可選步驟 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

比如下面的數(shù)據(jù)就可以放心的省略這兩步了：

R1 數(shù)據(jù)質(zhì)量夠好

R2 數(shù)據(jù)質(zhì)量也夠好

現(xiàn)在終于可以進(jìn)行變異檢測(cè)了，GATK 官網(wǎng)說 HC 表現(xiàn)比 UC 好，所以這里用 HC 進(jìn)行變異檢測(cè)：

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 完變異之后再進(jìn)行過濾：

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

然后就拿到變異檢測(cè)結(jié)果了，可以用 ANNOVAR 或 SnpEff 或 VEP 進(jìn)行注釋，根據(jù)自己的需要進(jìn)行篩選了。