RNA-seq原始數據質控後，是否要合併PE和SE的比對結果|《解螺旋技術交流圈》精華第1期

05-13

思考這樣一個問題：「RNAseq原始的Pair-end測序數據質控之後，部分Pair-end的read變成了Single-end的read，分開比對後得到了PE的BAM和SE的BAM，這個時候要不要合併這兩個BAM文件？」

關於這個問題，我們在知識星球上對此進行過討論。現在我總結一下我們的觀點。

思考問題的熊：

這個問題可能需要多幾個角度考慮。

1. 為什麼質控之後雙端的reads變成了單端？一種可能是因為一端的read質控確實不合格，第二種情況（通常也是更常見的）是因為因為建庫的原因，去接頭之後兩條reads overlap，從信息量來說就變成了單端，這個時候類似trimmomatic的軟體就會默認把其中一條read當作無用read而扔掉。
2. 通常來說，因為第一種情況而過濾掉的reads是相對少的，因此扔掉並無大礙，但是第二種就出問題了。目前大多數軟體都不支持直接輸入PE和SE的fastq文件進行mapping，你就必須把它分開做，這個時候你要是把單端的敢扔掉其實就是扔掉了大量的信息。
3. 怎麼解決這個問題呢？目前看來比較簡單的方法是直接將trimmomatic的一個參數keepBothReads 改為true ，讓因為第二種原因而扔掉的reads得以保留，然後直接用pair的fastq做mapping就可以了。
4. 這個參數的影響可以看圖

解螺旋的礦工（星球中的星主：YellowTree+）：

RNAseq的數據在質控之後，有時候甚至會有多達10%的read會從Pair-end變為Single-End！！！丟掉的話，損失很大，所以要考慮留下來。我的看法和 @思考問題的熊有相似之處，不過在處理Single-End的read上有些不同。我覺得可以這樣做：
1. 過濾後，Pair-end read和Single-End Read分開各自去完成比對，得到各自的比對結果（BAM格式）；
2. 分別計算PE的比對結果和SE比對結果在各個轉錄本上的覆蓋數，然後把它們相加起來；
3. 再用想在常用的基因表達分析工具（如EdgeR、DESeq等）進行下游分析；
這裡在（2）中唯一要擔心的是，有些SE Read覆蓋的轉錄本也許不是最準確的那個（因為缺了另一半，你無法有效去判斷），但我覺得這部分是少數，對結果不會有影響，這就是我認為可以分開比對，分別計算覆蓋數再合併的原因。這樣做的另一個好處是，不用太擔心另一半低質量的Read誤導了比對結果（不過可能這種情況也不會太多）——這也就是我和 @思考問題的熊存在差異的地方，不過我覺得熊的看法的一個好處是，處理起來會更簡單一些。

CJ

@YellowTree+，支持哈。另，對於Trimmomatic處理之後產生大量單端數據，那麼必然是size-selection這些步驟出問題咯。如果硬是要保留，我同樣支持星主的做法。

@思考問題的熊，在此情況下，我覺得也同樣有必要考慮Trimmomatic的其他參數，Trimmomatic默認的剪切模式會產生正反完全 overlap 或者更過的情況。一些比對軟體，如bowtie（或者hisat）似乎是不支持的。

礦工註解：

CJ在這裡指的是ILLUMINACLIP參數的「keepBothReads」參數。這個參數很重要，它的作用是R1和R2在去除了接頭序列之後，如果剩餘的部分是完全反向互補的，其默認參數（false），會把整條與R1完全反向互補的 R2去除，當做重複去除掉，但在有些情況下，例如這裡需要用到Paired reads的時候，就要將這個參數改為 true，否則會損失一部分Paired reads。

我舉個例子，看一個 PE150 數據的測試，就知道 keepBothReads 參數的重要性了：

$ java -jar trimmomatic-0.36.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:51ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequencesInput Read Pairs: 2500 Both Surviving: 1633 (65.32%) Forward Only Surviving: 828 (33.12%) Reverse Only Surviving: 12 (0.48%) Dropped: 27 (1.08%)TrimmomaticPE: Completed successfully# 使用 ILLUMINACLIP 默認的第六個參數 false，只有 65.32% paired reads 保留下來$ java -jar trimmomatic-0.36.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:51ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequencesInput Read Pairs: 2500 Both Surviving: 2439 (97.56%) Forward Only Surviving: 22 (0.88%) Reverse Only Surviving: 16 (0.64%) Dropped: 23 (0.92%)TrimmomaticPE: Completed successfully# 將 ILLUMINACLIP 第六個參數改為 true，其餘所有參數均相同，結果有 97.56% paired reads 保留下來

推薦閱讀

【乾貨】這麼說，FPKM和RPKM真的是錯的咯？！——關於FPKM/RPKM的深度反思
GATK4.0和全基因組數據分析實踐（上）
GATK4.0和全基因組數據分析實踐（下）
該如何自學入門生物信息學

這是我的知識星球：解螺旋技術交流圈，是一個與讀者朋友們的私人朋友圈，歡迎你的加入。我有9年前沿而完整的生物信息學、NGS領域的工作經歷，在該領域發有多篇Nature級別的科學文章。

這是知識星球上第一個真正與基因組學和生物信息學強相關的圈子。我旨在營造一個高質量的組學知識圈和人脈圈，通過提問、彼此分享、交流經驗、心得等，更好地學習生物信息學知識提升基因組數據分析和解讀的能力。

在這裡你可以結識到全國優秀的基因組學和生物信息學專家，同時可以分享你的經驗、見解和思考，有問題也可以向我提問或者向圈裡的星友們提問。

知識星球邀請鏈接：「解螺旋技術交流圈」