前言

接下來我們要介紹的是 RNA-seq 數(shù)據(jù)的處理分析流程吕粹，根據(jù) RNA-seq 測序技術(shù)的不同，可以分為三種：

Stark et al. Nat Rev Genet(2019)

1. short-read
2. long-read
3. direct RNA-seq

而我們一般的 RNA-seq 測序數(shù)據(jù)分析流程算法啰扛，基本上都是基于 short-read（短讀長）技術(shù)所產(chǎn)生的數(shù)據(jù)文件

目前幔烛，我們可以從 Short Read Archive(SRA) 數(shù)據(jù)庫獲取的 RNA-seq 數(shù)據(jù)中管行，有超過 95% 的數(shù)據(jù)是由 Illumina 公司的 short read 測序技術(shù)所產(chǎn)生的

其分析過程可以用下面的路線圖表示

Conesa et al. Genome Biology (2016)

該路線圖大致分為三個部分：

1. 數(shù)據(jù)獲取：
   • 包括實(shí)驗(yàn)設(shè)計(jì)准脂、測序設(shè)計(jì)以及數(shù)據(jù)下機(jī)后的 raw reads 數(shù)據(jù)的質(zhì)控
2. 數(shù)據(jù)分析
   • 在獲取到干凈的數(shù)據(jù)之后劫扒，可以進(jìn)行 reads 的比對，然后進(jìn)行基因表達(dá)的量化意狠、差異表達(dá)分析粟关、功能富集分析等
3. 高級分析
   • 包括數(shù)據(jù)的可視化，其他小分子 RNA 分析环戈、融合分析以及與其他類型的數(shù)據(jù)進(jìn)行整合分析等

而我們分析的起始點(diǎn)闷板，是從原始數(shù)據(jù)開始的，也就是獲取 raw reads 數(shù)據(jù)院塞。通常這種高通量測序數(shù)據(jù)會保存為 FASTQ 格式的文件遮晚。

FASTQ 格式是一種以 ASCII 碼字符的形式保存生物序列及其對應(yīng)的每個堿基的質(zhì)量的文本文件。

FASTQ 文件中每條序列（通常是一條 read）是由 4 行組成拦止，其中：

• 第一行以 @ 字符開頭县遣，之后的字符為序列的標(biāo)識符和描述信息
• 第二行為具體的序列
• 第三行以 + 符號開頭糜颠，之后可以可選地加上與第一行一樣的序列標(biāo)識或描述信息
• 第四行為堿基質(zhì)量分?jǐn)?shù)（Phred），其字符數(shù)量與第二行相等萧求，每個字符表示對應(yīng)堿基的質(zhì)量得分其兴，例如

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

其中，堿基質(zhì)量值的編碼方式為

1. 先將堿基錯誤率 P 進(jìn)行負(fù)對數(shù)轉(zhuǎn)換夸政，得到 Q 值
   

2. 然后將 Q 值加上 33 或 64 得到的值所對應(yīng)的 ASCII 碼即為堿基質(zhì)量分?jǐn)?shù)

例如元旬，錯誤率 P = 0.01，則 Q = 20守问，如果是 Phred33 則對應(yīng)的質(zhì)量為字符 5（53）匀归，如果是 Phred64 則對應(yīng)的字符為 T(84)

分析流程

1. 數(shù)據(jù)獲取

一般情況下，如果自己有送樣檢測數(shù)據(jù)的話耗帕，測序公司會提供原始的 FASTQ 格式的數(shù)據(jù)穆端。如果我們要使用別人文章中發(fā)表的公開數(shù)據(jù)，還需要從數(shù)據(jù)庫中下載對應(yīng)的數(shù)據(jù)

例如仿便，我們從 SRA 數(shù)據(jù)中下載的原始測序文件是 sra 格式体啰，我們需要先使用工具將其轉(zhuǎn)換為 FASTQ 格式

2. 質(zhì)量控制

主要在三個地方需要對數(shù)據(jù)的質(zhì)量進(jìn)行監(jiān)控

• 獲取原始數(shù)據(jù)之后
• 比對完之后
• 表達(dá)定量之后

2.1 raw read

對 raw reads 數(shù)據(jù)進(jìn)行質(zhì)量控制，需要分析序列的質(zhì)量探越、GC 含量狡赐、是否存在接頭、短重復(fù)序列的分布钦幔、測序錯誤以及 PCR 重復(fù)和污染

質(zhì)控軟件：

• FastQC：用于分析 Illumina 測序平臺的數(shù)據(jù)
• NGSQC：可應(yīng)用于所有平臺

一般來說，reads 的質(zhì)量會朝著 3' 端遞減常柄，如果堿基的質(zhì)量太低鲤氢，我們需要刪除它以提高比對率

FASTX-Toolkit 和 Trimmomatic 兩個軟件可以用于切除低質(zhì)量的堿基和接頭序列

2.2 比對后

reads 通常需要比對到一個參考基因組或轉(zhuǎn)錄組，而比對的質(zhì)量是評估測序準(zhǔn)確率和是否存在 DNA 污染的一個重要指標(biāo)

比對質(zhì)量通常為比對到的 reads 數(shù)占總 reads 數(shù)的比例西潘。

例如卷玉，比對到人類參考基因組的比對質(zhì)量通常需要在 70-90%，且有大量的 reads 映射到一個相同的區(qū)間內(nèi)喷市。如果是比對到轉(zhuǎn)錄本上相种，由于可變剪切的影響，可以適當(dāng)放寬比對質(zhì)量

在外顯子和比對方向上的 read 覆蓋率的均一性品姓，也是評估質(zhì)量的重要指標(biāo)寝并。如果 reads 主要聚集在轉(zhuǎn)錄本的 3' 端，可能表明原始樣本的 RNA 質(zhì)量較低

比對上的 reads 的 GC 含量腹备，可能揭示了 PCR 的錯誤率

主要軟件有：Picard衬潦、RSeQC 和 Qualimap

2.3 定量后

在計(jì)算完表達(dá)的量化值之后，可以計(jì)算 GC 含量和基因長度的誤差植酥，在必要時可以使用標(biāo)準(zhǔn)化方法來進(jìn)行校正

如果參考轉(zhuǎn)錄組注釋得很好镀岛，則可以分析樣本的生物構(gòu)成弦牡，來評估 RNA 純化步驟的質(zhì)量。例如漂羊，rRNA 和 small RNA 不能出現(xiàn)在 polyA longRNA 的制備中

NOISeq 和 EDASeq 等 R 包可以使用圖形來展示 count 數(shù)據(jù)的質(zhì)量控制

2.4 可重復(fù)性

上面的質(zhì)量控制都只是針對單個樣本的驾锰，此外，不同樣本之間的可重復(fù)性評估走越，對于評價整個數(shù)據(jù)集的質(zhì)量也是至關(guān)重要的

技術(shù)重復(fù)樣本的可重現(xiàn)性一般很高（spearman ）椭豫，但是生物學(xué)重復(fù)樣本之間并沒有明確的標(biāo)準(zhǔn)，取決于實(shí)驗(yàn)系統(tǒng)的異質(zhì)性买喧。如果不同實(shí)驗(yàn)系統(tǒng)之間存在差異基因捻悯，則同一條件下的生物學(xué)重復(fù)在主成分分析（PCA）中會被聚類在一起。

3. 序列比對

在對樣本的 raw reads 進(jìn)行質(zhì)控之后淤毛，就可以進(jìn)行序列比對了今缚，序列比對主要有三種策略，如下圖

Conesa et al. Genome Biology (2016)

如果有參考序列低淡，根據(jù)參考序列的不同姓言，可以分為

• 比對到基因組：使用間隔比對算法，如 TopHat蔗蹋、STAR 等何荚，然后根據(jù)是否提供了注釋文件（GFF 格式文件，包含轉(zhuǎn)錄本位置信息）猪杭，又可以分為轉(zhuǎn)錄本識別和轉(zhuǎn)錄本發(fā)現(xiàn)并進(jìn)行定量分析
• 比對到轉(zhuǎn)錄組：使用非間隔比對算法餐塘，如 Bowtie 等，然后使用 RSEM 或 Kallisto 方法識別轉(zhuǎn)錄本并計(jì)算定量信息

如果沒有參考序列皂吮，則需要先把序列組裝成轉(zhuǎn)錄本戒傻，再將 reads 比對到組裝后的參考轉(zhuǎn)錄本上，然后使用 HTseq-count 等算法對轉(zhuǎn)錄本進(jìn)行定量

3.1 轉(zhuǎn)錄本發(fā)現(xiàn)

使用 Illumina 技術(shù)檢測的 short reads 來發(fā)現(xiàn)新的轉(zhuǎn)錄本是 RNA-seq 分析中的一個挑戰(zhàn)蜂筹。通常來說需纳，短 reads 很少會跨越多個剪切位點(diǎn)，這就很難直接推斷出一個轉(zhuǎn)錄本的整體長度艺挪。

此外不翩，轉(zhuǎn)錄的起始和終止位置也比較難識別，一些像 GRIT 的工具麻裳，通過合并 5' 端的信息可以提高異構(gòu)體識別的準(zhǔn)確性口蝠。其他如 Cufflinks、iReckon掂器、SLIDE 和 StringTie 等方法亚皂，通過結(jié)合現(xiàn)有的注釋信息，作為一個可能的異構(gòu)體列表

一些尋找基因的工具国瓮，如 Augustus灭必，結(jié)合 RNA-seq 數(shù)據(jù)狞谱，可以更好的注釋蛋白編碼轉(zhuǎn)錄本，但是對非編碼轉(zhuǎn)錄本的性能更差禁漓。

3.2 De novo 轉(zhuǎn)錄本重構(gòu)

在沒有轉(zhuǎn)錄本或轉(zhuǎn)錄本不全的情況下跟衅，可以對 reads 進(jìn)行組裝來重構(gòu)一份轉(zhuǎn)錄本〔ゼ撸可選的方法很多伶跷，如 SOAPdenovoTrans、Oases秘狞、Trans-ABySS 或 Trinity

通常來說叭莫，使用雙端鏈特異性測序和 long reads 測序包含更多的信息，會有更好的效果

雖然烁试，對于低表達(dá)的轉(zhuǎn)錄本進(jìn)行組裝的可靠性較低雇初，但是 reads 太多也會導(dǎo)致潛在的組裝錯誤和較長的時間消耗等問題。因此减响，在深度測序的樣本中靖诗，可以適當(dāng)減少 reads 的數(shù)量

對于多樣本的比較，可以將所有樣本作為一個輸入來構(gòu)建參考轉(zhuǎn)錄本支示，然后分別對每個樣本的 reads 進(jìn)行比對

無論是使用參考序列還是從頭開始組裝刊橘，使用短 reads 的 Illumina 技術(shù)來完全重構(gòu)轉(zhuǎn)錄組仍然是一個具有挑戰(zhàn)性的問題

4. 轉(zhuǎn)錄組定量

RNA-seq 最廣泛的應(yīng)用就是用來評估基因和轉(zhuǎn)錄本的表達(dá)，這一應(yīng)用主要是基于比對到轉(zhuǎn)錄組區(qū)間內(nèi)的 reads 的數(shù)量

最簡單的方法是颂鸿，使用 HTSeq-count 或 featureCounts 計(jì)算區(qū)間內(nèi)的 reads 數(shù)來量化基因的表達(dá)促绵。這種基因水平的（不是轉(zhuǎn)錄本水平）的量化方法使用的是 GTF 文件，這種文件包含外顯子和基因在基因組上的坐標(biāo)嘴纺。

但一般不能直接使用 read count 來比較基因的表達(dá)水平绞愚，因?yàn)樵撝禃艿睫D(zhuǎn)錄本長度、reads 總數(shù)以及測序偏差等因素的影響颖医。所以需要先進(jìn)行標(biāo)準(zhǔn)化，標(biāo)準(zhǔn)化方法有

1. RPKM/FPKM: 每百萬 reads 每一千堿基對中包含的 reads 數(shù)

該方法先計(jì)算測序深度系數(shù)裆蒸，即總 reads 數(shù)除以 一百萬熔萧，然后計(jì)算基因或轉(zhuǎn)錄本的長度（單位為 kb），標(biāo)準(zhǔn)化順序?yàn)橄认郎y序深度的影響僚祷，再消除長度的影響：

其中

• x 表示一個基因或轉(zhuǎn)錄本佛致，或基因組上一段特定的區(qū)域
• 表示比對到 x 外顯子區(qū)域的 reads 數(shù)；
• R 表示當(dāng)前樣本中包含的全部 reads 數(shù)
• 表示 x 外顯子區(qū)域包含的堿基數(shù)（長度辙谜，bp）

FPKM 與 RPKM 的計(jì)算公式一樣俺榆，只是 RPKM 用于單端測序，FPKM 用于雙端測序

2. TPM: 其與 RPKM 最大的區(qū)別是装哆，標(biāo)準(zhǔn)化順序?yàn)橄认蜷L度的影響罐脊，再消除測序深度的影響

首先定嗓，將 reads count 除以基因或轉(zhuǎn)錄本的長度（kb）得到 RPK(reads per kilobase)，然后將樣本中所有的 RPK 加起來除以 萍桌，得到標(biāo)準(zhǔn)化系數(shù)宵溅，最后使用 RPK 除以標(biāo)注化系數(shù)

其中

• x 表示一個基因或轉(zhuǎn)錄本，或基因組上一段特定的區(qū)域
• 表示比對到 x 外顯子區(qū)域的 reads 數(shù)
• 表示 x 外顯子區(qū)域包含的堿基數(shù)（kp）
• N 表示基因或轉(zhuǎn)錄本總數(shù)

這樣上炎，每個樣本的 TPM 總和是一樣的恃逻，便于比較樣本間的差異

目前，也有許多復(fù)雜的算法通過解決相關(guān)轉(zhuǎn)錄本共享 reads 的問題來評估轉(zhuǎn)錄本水平的表達(dá)藕施，例如寇损，Cufflinks 使用 TopHat 的比對結(jié)果，應(yīng)用期望最大化算法來評估轉(zhuǎn)錄本的豐度裳食。這一方法考慮到長度不同的基因的 reads 分布并不均勻等因素的影響矛市。

還有其他算法也可以量化轉(zhuǎn)錄組的表達(dá)，例如 RSEM胞谈、eXpress尘盼、Sailfish 和 kallisto 等。這些方法允許轉(zhuǎn)錄本之間存在多比對的 reads烦绳，并輸出經(jīng)測序偏差校正的樣本內(nèi)歸一化值卿捎。

5. 差異表達(dá)分析

差異表達(dá)分析是對樣本間基因的表達(dá)值進(jìn)行比較，雖然 RPKM径密、FPKM 和 TPM 標(biāo)準(zhǔn)化方法消除了測序深度和基因或轉(zhuǎn)錄本的長度因素的影響午阵，但這些方法依賴于總的或有效的 reads 數(shù)，當(dāng)樣本的具有異質(zhì)性轉(zhuǎn)錄本分布或當(dāng)高表達(dá)或差異表達(dá)的特征扭曲了 count 分布時享扔，表現(xiàn)欠佳

而像 TMM底桂、DESeq、PoissonSeq 和 UpperQuartile 等方法會忽略高變異或高表達(dá)的特征。

干擾樣本內(nèi)比較的其他因素包括不同樣本的轉(zhuǎn)錄本長度變化、轉(zhuǎn)錄本覆蓋位置的偏差规辱、平均片段大小以及基因的 GC 含量等

NOISeq 這個 R 包提供了多種繪圖稻艰，來識別 RNA-seq 數(shù)據(jù)中的誤差來源，并應(yīng)用相應(yīng)的方法來標(biāo)準(zhǔn)化這些誤差

除了這些樣本內(nèi)特異的標(biāo)準(zhǔn)化方法，還需要解決數(shù)據(jù)集之間的批次效應(yīng)（不同實(shí)驗(yàn)條件下產(chǎn)生的數(shù)據(jù)之間存在的差異），批次矯正方法有 COMBAT 和 ARSyN 等，雖然這些方法是針對芯片數(shù)據(jù)設(shè)計(jì)的捶码，但是在 RNA-seq 數(shù)據(jù)中也有很好的效果

計(jì)算差異表達(dá)的方法有很多，有些方法或链，如 edgeR 將原始的 read counts 作為輸入惫恼，并在統(tǒng)計(jì)模型中加入了標(biāo)準(zhǔn)化，另一些方法澳盐，需要先對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化祈纯，如 DESeq2 使用的是負(fù)二項(xiàng)分布作為參考分布令宿，并提供了自己的標(biāo)準(zhǔn)化方法。

baySeq 和 EBSeq 是貝葉斯方法盆繁，還有一些基于線性模型的方法掀淘。最后，一些非參數(shù)方法油昂，如 NOISeq 和 SAMseq

對于小樣本量的研究革娄，負(fù)二項(xiàng)分布會存在噪音污染，這種情況下冕碟，一些簡單點(diǎn)方法拦惋，如基于 Poisson 分布的 DEGseq，或者基于經(jīng)驗(yàn)分布的 NOISeq 可能會更好些安寺。

但是需要強(qiáng)調(diào)的是厕妖，在沒有足夠生物學(xué)重復(fù)的情況下，無法進(jìn)行總體的推斷挑庶，因此任何 p 值計(jì)算都是無效的言秸。

許多獨(dú)立的研究都已經(jīng)證實(shí)，選擇不同的方法會對結(jié)果有一定的影響迎捺，而且沒有哪一種方法能夠適用于所有的數(shù)據(jù)举畸，所以，推薦在分析的時候使用多個軟件進(jìn)行相互驗(yàn)證凳枝。

6. 可變剪切分析

可變剪接(Alternative Splicing) 是指轉(zhuǎn)錄形成的前體 RNA 通過去除內(nèi)含子抄沮、連接外顯子而形成成熟 RNA 的過程，從而實(shí)現(xiàn)一個基因同時編碼多種蛋白質(zhì)岖瑰，實(shí)現(xiàn)生物功能多樣性

在不同組織或者發(fā)育的不同階段叛买，可變剪接不是一成不變的，在特定的組織或條件下蹋订，通過連接不同的外顯子率挣，會產(chǎn)生特定的剪接異構(gòu)體（isoform）。有大量的研究發(fā)現(xiàn)露戒，可變剪接的變化與癌癥等多種疾病相關(guān)难礼，所以研究可變剪接在不同組織中的作用是非常有意義的。

轉(zhuǎn)錄本水平的差異表達(dá)分析可以潛在地檢測同一基因的轉(zhuǎn)錄異構(gòu)體表達(dá)的變化玫锋，已經(jīng)有一些算法應(yīng)用于 RNA-seq 數(shù)據(jù)的中進(jìn)行可變剪切分析

這些方法主要分為兩大類：

1. 異構(gòu)體表達(dá)估計(jì)與差異表達(dá)相結(jié)合，來揭示總基因表達(dá)中每種異構(gòu)體的比例變化

例如讼呢，BASIS 方法使用分層貝葉斯模型來直接推斷轉(zhuǎn)錄異構(gòu)體的差異表達(dá)撩鹿；CuffDiff2 方法先評估異構(gòu)體的表達(dá)，然后比較它們之間的差異悦屏；rSeqDiff 方法使用分層似然率檢驗(yàn)同時檢測無剪接變化的差異基因表達(dá)和差異異構(gòu)體表達(dá)节沦。

所有這些方法通常都受限于短讀長測序的內(nèi)在局限性键思，無法在異構(gòu)體水平上進(jìn)行準(zhǔn)確識別

2. 一種所謂的 exon-based 的方法，它跳過了對異構(gòu)體表達(dá)的估計(jì)甫贯，通過比較樣本之間基因外顯子和連接點(diǎn)上的 reads 分布來檢測可變剪接的信號

其基本假設(shè)為：可以在外顯子及其連接點(diǎn)的信號中追蹤異構(gòu)體表達(dá)的差異吼鳞。

DEXseq 和 DSGSeq 采用類似的思路，通過檢測基因的外顯子（和連接點(diǎn)）上 read counts 的差異顯著性來識別不同的異構(gòu)體叫搁。

rMATS 是通過比較用連接點(diǎn)的 reads 定義的外顯子 inclusion levels 表達(dá)水平的差異

7. 融合分析

基因融合是指兩個基因的全部或一部分的序列相互融合為一個新的基因的過程赔桌。其有可能是染色體易位、中間缺失或染色體倒置所導(dǎo)致的渴逻，可在 DNA 或 RNA 層面上表達(dá)疾党。

融合基因通過基因失調(diào)、融合產(chǎn)生嵌合體蛋白這兩種機(jī)制引發(fā)癌癥的發(fā)生惨奕。

目前雪位，RNA-seq 融合算法 100 多種，有人對常用的 15 中融合檢測算法進(jìn)行了比較

Liu et al Nucleic Acids Research, 2016

沒有哪一個算法具有明顯的優(yōu)勢梨撞，整體來看雹洗，SOAPfuse 可能會好一些，FusionCatcher 和 JAFFA 其次卧波。

8. 功能注釋

標(biāo)準(zhǔn)的轉(zhuǎn)錄組分析的最后一步时肿，是使用差異表達(dá)基因來進(jìn)行功能或通路的注釋。最常用的兩類方法是：

• 基于超幾何分布的過表達(dá)富集分析
• GSEA 富集分析

一些工具幽勒，如 GOseq 考慮了基因長度等因素對差異表達(dá)結(jié)果的影響嗜侮，并使用超幾何分布進(jìn)行富集分析，GSVA 和 SeqGSEA 使用類似 GSEA 的方法進(jìn)行功能富集

功能富集需要預(yù)先定義的基因集合或通路啥容，包括 GO锈颗、KEGG、Reactome 等數(shù)據(jù)庫咪惠。

通過在蛋白質(zhì)數(shù)據(jù)庫（例如 SwissProt）和包含保守蛋白質(zhì)結(jié)構(gòu)域（例如 Pfam 和 InterPro）的數(shù)據(jù)庫中搜索相似序列击吱，使用直系同源分析對蛋白質(zhì)編碼的轉(zhuǎn)錄本進(jìn)行功能注釋。而 Rfam 數(shù)據(jù)庫包含許多特征明確的 RNA 家族遥昧，例如 rRNA 或 tRNA覆醇，而 mirBase 和 Miranda 是專門研究 miRNA 的

9. 整合分析

將 RNA-seq 數(shù)據(jù)與其他類型的全基因組數(shù)據(jù)進(jìn)行整合分析，使我們能夠?qū)⒒虮磉_(dá)的調(diào)控與分子生理學(xué)和功能基因組學(xué)的特定方面聯(lián)系起來炭臭。

1. DNA 測序

將 RNA 測序和 DNA 測序相結(jié)合永脓，可以進(jìn)行 SNP、RNA 編輯和表達(dá)數(shù)量性狀基因座（eQTL）比對等分析鞋仍。

2. DNA 甲基化

將 DNA 甲基化和 RNA-seq 整合常摧，可以分析差異表達(dá)基因和甲基化模式之間的相關(guān)性。使用的算法包括：廣義線性模型、logistic 回歸模型和經(jīng)驗(yàn)貝葉斯模型等

3. 染色體特征

通過整合 RNA-seq 和 Chip-seq 數(shù)據(jù)落午，可以降低 Chip-seq 分析的假陽性谎懦，并展示 TF 對其靶基因的激活或抑制作用。

4. MicroRNA

整合 RNA-seq 和 miRNA-seq 有可能揭示 miRNAs 對轉(zhuǎn)錄穩(wěn)態(tài)水平的調(diào)節(jié)作用

5. 蛋白質(zhì)組和甲基化組

RNA-seq 與蛋白質(zhì)組學(xué)的整合是有爭議的溃斋，因?yàn)檫@兩種測量結(jié)果通常顯示出很低的相關(guān)性（~0.4）界拦。盡管如此，蛋白質(zhì)組學(xué)和 RNA-seq 的配對分析可用于識別新的異構(gòu)體

轉(zhuǎn)錄組與代謝組數(shù)據(jù)的整合已被用于識別在基因表達(dá)和代謝物水平上受調(diào)控的通路梗劫，并且可以使用工具來可視化通路上下文的結(jié)果享甸，如 MassTRIX、Paintomics在跳、VANTED v2 和 SteinerNet