聊聊轉錄組測序——2.數據分析與解讀（上）

01-28

前言

「你把表達量和差異表達的表格給我就行了」

花那麼多錢就換兩張表格，你的良心不會痛么！（當年表達晶元都這麼做的也沒見誰良心痛啊）

當然，畢竟別人不一定靠這個吃飯，但（找不到工作的）學生物的你，稍微了解一下測序的分析流程還是值得的，畢竟技多不壓身嘛。

所以這一部分主要介紹轉錄組測序的分析流程和原理，從拿到原始數據開始，講到KEGG/Gene Ontology等功能注釋，順便推薦一下常用軟體。字數所限，這一篇先講（不生成文章能用的圖表的） Data Cleaning和比對，如果只想知道怎麼看懂文獻裡面的結果，可以直接等下一篇了。

流程概覽

轉錄組測序的分析流程大致可以分成三類，包括基因組比對（Genome mapping）、轉錄組比對（Transcriptome mapping）、轉錄組組裝（Reference-free assembly），見下圖。其中第三種主要是用於分析沒有參考基因組和基因注釋的物種，應用場合較少且不適合新手入門。對於人、小鼠、大鼠等模式物種，通常用前兩種方法進行分析。雖然轉錄組比對相關軟體和流程同樣層出不窮，但對於基因組信息較為完善的模式物種，推薦使用基因組比對的方式進行分析，具體原因下文的「比對」部分會有說明。我們下面也主要對基因組比對的方法進行介紹。

圖片來源：https://www.ncbi.nlm.nih.gov/pubmed/26813401

1. Data Cleaning

從原始數據（Raw Data）到乾淨數據（CleannData）的過程，有人翻譯成「數據清洗」，實在叫不習慣，那我就不翻譯了。

Illumina測序儀下機的數據通常為Bcl格式，是將同一個測序通道（Lane）所有樣品的數據混雜在一起的，所以公司一般不會提供Bcl文件。測序公司使用Illumina官方出品的Bcl2FastQ軟體，根據Index序列分割轉換成每個樣品的FastQ文件，打開長這樣：

每一條序列（read）包含四行，第一行是read的ID，第二行是序列，第四行是序列中每個鹼基的測序質量（更具體的細節可參考FASTQ format - Wikipedia）。

原始數據沒法直接分析，是因為部分reads測序質量較低，可能會誤導後續結果，因此需要對低質量鹼基太多或N（未能識別的鹼基）太多的reads進行去除；此外，部分測序文庫的插入片段太短，導致測到兩側的接頭序列（請參考上一篇的測序文庫結構圖理解），這些序列接頭也需要從reads中去除。最後，我們也會對清洗前後的Raw Data和Clean Data進行評估，評估內容包括鹼基質量、序列長度、鹼基比例、GC含量、重複序列、Kmers等（詳情請參考FastQC說明文檔FastQC A Quality Control tool for High Throughput Sequence Data）。

最後說一句，其實大多數測序公司都會提供CleannData的。

#常用軟體#

我以前都是用cutadapt + FASTX-Toolkit的組合，直到同事們給我推薦了Trim Galore，質量評估使用FastQC。

2. 比對

由於二代測序的reads長度通常介於50~300個鹼基，因此即便使用雙端測序，也基本不可能覆蓋完整的mRNA轉錄本，因此想直接用FastQ文件從頭分析測到了哪些轉錄本需要非常複雜的分析和計算。好在通常情況下，公共資料庫已經提供了測序樣品的基因組和轉錄本的序列。因此我們只需要知道，每一條reads來自哪一條轉錄本就可以了，這個將reads與參考（Reference）基因組/轉錄組的序列進行比較和匹配的過程，我們通常稱之為「比對」（文獻中提到的read alignment和mapping通常說的都是這個）。

正如前文所述，轉錄組測序的比對通常分為基因組比對和轉錄組比對兩種，顧名思義，基因組比對就是把reads比對到完整的基因組序列上，而轉錄組比對則是把reads比對到所有已知的轉錄本序列上。如果不是很急或者只想知道已知轉錄本表達量，個人建議使用基因組比對的方法進行分析，理由如下：

① 轉錄組比對需要準確的已知轉錄本的序列，對於來自未知轉錄本（比如一些未被資料庫收錄的lncRNA）或序列不準確的reads無法正確比對；

② 與上一條類似，轉錄組比對不能對轉錄本的可變剪接進行分析，資料庫中未收錄的剪接位點會被直接丟棄；

③ 由於同一個基因存在不同的轉錄本，因此很多reads可以同時完美比對到多個轉錄本，reads的比對評分會偏低，可能被後續計算表達量的軟體捨棄，影響後續分析（有部分軟體解決了這個問題）；

④ 由於與DNA測序使用的參考序列不同，因此不利於RNA和DNA數據的整合分析。

而上面的問題使用基因組比對都可以解決。

此外，值得注意的是，RNA測序並不能直接使用DNA測序常用的BWA、Bowtie等比對軟體，這是由於真核生物內含子的存在，導致測到的reads並不與基因組序列完全一致（如下圖所示），因此需要使用Tophat/HISAT/STAR等專門為RNA測序設計的軟體進行比對。

圖片來源：GATK | Best Practices

比對結果會展示為BAM/SAM文件，其中BAM格式是SAM格式的二進位版本（請理解為壓縮後的版本，用Samtools可以打開），打開之後長這樣：

BAM文件中每行代表一條reads的比對信息，其中第一列是read的ID，第二列為FLAG（包括是否雙端比對，比對位點是否唯一等信息），第三列為比對的染色體，第四列為比對的起始位置，第六列為CIGAR值，代表比對的具體方式（例60M2D80M代表60個鹼基完美匹配+2個鹼基缺失+80個鹼基完美匹配）等等，BAM文件的具體內容可參考SAM - Genome Analysis Wiki和http://samtools.github.io/hts-specs/SAMv1.pdf（後面這個要翻牆）。

#常用軟體#

基因組比對：

Tophat2：可以說是最被公認的RNA測序比對軟體（實際上是在DNA比對軟體Bowtie的基礎上做了一個殼），相信很多做RNA測序的同學都是看著Tophat發表在Nature Protocol上的步驟一步步入門RNA測序的；

HISAT2：Tophat2的非正式升級版本（因為據說還會有Tophat3），在Tophat的演算法基礎了上做了大量的改進，而且克服了Tophat最大的缺點——速度慢，Nature Protocol上同樣發表了操作流程；

STAR：ENCODE計劃御用比對軟體，權威程度可以與Tophat平起平坐，並且比對速度極快；

MapSplice：TCGA使用的比對軟體，我自己沒用過；

RSEM：RSEM更像一個軟體包而不是一個比對軟體，能夠提供從比對到計算差異表達的所有步驟，由於不需要自己寫代碼串聯不同軟體生成的數據格式，因此用起來比較省時省力，值得注意的是，TCGA使用MapSplice比對後再用RSEM計算表達量，並沒有直接只用RSEM原裝的Bowtie的比對結果。

轉錄組比對：這類型的軟體我用的不多，最近嘗試過Nature Methods上面發表的Salmon，能從Clean reads直接算到表達量，優點是，快，非常快。然而這個軟體連BAM文件都沒生成，雖然只是定量的話BAM文件的確沒什麼用就是了…

1.5 #可選步驟# 核糖體RNA（rRNA）去除

嗯，寫完2再寫1.5是我不對。

如果對上一篇還有印象的話，我們曾提到，轉錄組測序有一種偏貴的使用核糖體RNA去除技術構建文庫的測序。但是經常做實驗的你一定知道，這種去除是沒法做到100%去除rRNA的，更糟糕的是，同一批測序的每個樣品，rRNA的去除效率也會有一定差別的！

由於rRNA都是非編碼RNA序列，因此如果我們後續分析需要使用轉錄本組裝的方法鑒別新的lncRNA（long non-coding RNA，長非編碼RNA），這些rRNA的序列特徵很容易對lncRNA的鑒定造成干擾，因此我們必須對這些rRNA序列進行去除。

當然，如果不涉及組裝新lncRNA的話，rRNA的存在對分析結果的影響並不大。但如果樣品間rRNA殘留率差別較大，對定量的準確性會有較大影響，因此有能力的話還是建議去除rRNA序列。就算不去除，用比對軟體算算rRNA序列佔總數據量的比例也是好的，一旦不小心發現12G的數據裡面6G都來自於rRNA…（嗯我不是教你們跟公司撕X…）

#常用軟體#

核糖體去除實際上也是通過比對來進行，我在Rfam上下載rRNA的序列後，直接使用Bowtie2進行比對。

至於比對核糖體之後怎麼拿到沒有rRNA的FastQ文件，我不太清楚別人是怎麼做的，我是用Python把沒比對上的Reads的ID提取出來存成一個表格，再用Seqtk提取FastQ文件。

#########################未完待續，我們下回分解#########################

最後關於常用軟體：

首先推薦OMIC TOOLS網站（https://omictools.com/rna-seq-category），上面收集了大量高通量測序相關的軟體以及軟體之間的對比評測文獻。

其次，我寫的常用軟體基本基於我自己和身邊同事的使用情況，以及在文獻中看到的情況，並不涉及對軟體自身性能的評判。比如華大基因開發的比對用軟體SOAPsplice，沒有提到不是說它不好，而是我的確沒用過，而且用的人也比較少。

至於如何判斷哪個軟體更好，可以參考軟體評測的文獻，但沒必要以此為絕對標準，例如我見到有評測說STAR完爆Tophat的，也有說Tophat完爆STAR的，而且這兩篇評測都發在Nature Method上。我個人認為能一直存活下來也被廣泛使用的軟體大多各有千秋，對於一個沒精力看源代碼（也看不懂java和C++）的使用者來說，想要評估一個軟體的好壞，直接參考高分雜誌的使用情況來推算，可能是更靠譜的選擇。

另外，長期歡迎推薦各類軟體。