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 indexSTAR --runThreadN 8 --runMode genomeGenerate         --genomeDir ./star_index/         --genomeFastaFiles ./genome/chrX.fa         --sjdbGTFfile ./genes/chrX.gtf
 

• 然後進行第一次序列比對：

#star 1-pass alignSTAR --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 indexSTAR --runThreadN 8 --runMode genomeGenerate         --genomeDir ./star_index_2pass/         --genomeFastaFiles ./genome/chrX.fa         --sjdbFileChrStartEnd ./star_1pass/ERR188044SJ.out.tab
 

• 再進行 STAR 二次序列比對：

# star 2-pass alignSTAR --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 indexjava -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.fajava -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 檢測的真陽性率，如果時間緊張也可以不做，這一步影響很輕微 :)

# 可選步驟 IndelRealignjava -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 結果來決定。這一步省了又能節省時間 :)

# 可選步驟 BQSRjava -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.tablejava -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 進行註釋，根據自己的需要進行篩選了。