RNA-seq 分析流程 —— 概述

前言

接下來我們要介紹的是 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è)計、測序設(shè)計以及數(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)行整合分析等

而我們分析的起始點，是從原始數(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 值
    Q = -10log_{10}(P)

  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 定量后

在計算完表達(dá)的量化值之后，可以計算 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 R^2 > 0.9），但是生物學(xué)重復(fù)樣本之間并沒有明確的標(biāo)準(zhǔn)，取決于實驗系統(tǒng)的異質(zhì)性。如果不同實驗系統(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)錄本并計算定量信息

如果沒有參考序列，則需要先把序列組裝成轉(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 很少會跨越多個剪切位點，這就很難直接推斷出一個轉(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)錄本?？蛇x的方法很多，如 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 計算區(qū)間內(nèi)的 reads 數(shù)來量化基因的表達(dá)。這種基因水平的（不是轉(zhuǎn)錄本水平）的量化方法使用的是 GTF 文件，這種文件包含外顯子和基因在基因組上的坐標(biāo)。

但一般不能直接使用 read count 來比較基因的表達(dá)水平，因為該值會受到轉(zhuǎn)錄本長度、reads 總數(shù)以及測序偏差等因素的影響。所以需要先進(jìn)行標(biāo)準(zhǔn)化，標(biāo)準(zhǔn)化方法有

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

該方法先計算測序深度系數(shù)，即總 reads 數(shù)除以 一百萬，然后計算基因或轉(zhuǎn)錄本的長度（單位為 kb），標(biāo)準(zhǔn)化順序為先消除測序深度的影響，再消除長度的影響：
RPKM(x)=\frac{Reads \ per \ transcript}{\frac{total \ reads}{10^6} \times \frac{transcript \ length}{1000}} = \frac{Reads \ per \ transcript}{million \ reads \times transcript \ length (kb)} = \frac{10^9 \times C_x}{R \times L_x}

其中

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

FPKM 與 RPKM 的計算公式一樣，只是 RPKM 用于單端測序，FPKM 用于雙端測序

  1. TPM: 其與 RPKM 最大的區(qū)別是，標(biāo)準(zhǔn)化順序為先消除基因長度的影響，再消除測序深度的影響

首先，將 reads count 除以基因或轉(zhuǎn)錄本的長度（kb）得到 RPK(reads per kilobase)，然后將樣本中所有的 RPK 加起來除以 10^6，得到標(biāo)準(zhǔn)化系數(shù)，最后使用 RPK 除以標(biāo)注化系數(shù)
TPM(x) = \frac{C_x / L_x \times 10^6}{\sum_{i=1}^{N} C_i / L_i}

其中

  • 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)（不同實驗條件下產(chǎn)生的數(shù)據(jù)之間存在的差異），批次矯正方法有 COMBAT 和 ARSyN 等，雖然這些方法是針對芯片數(shù)據(jù)設(shè)計的，但是在 RNA-seq 數(shù)據(jù)中也有很好的效果

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

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

對于小樣本量的研究，負(fù)二項分布會存在噪音污染，這種情況下，一些簡單點方法，如基于 Poisson 分布的 DEGseq，或者基于經(jīng)驗分布的 NOISeq 可能會更好些。

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

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

6. 可變剪切分析

可變剪接(Alternative Splicing) 是指轉(zhuǎn)錄形成的前體 RNA 通過去除內(nèi)含子、連接外顯子而形成成熟 RNA 的過程，從而實現(xiàn)一個基因同時編碼多種蛋白質(zhì)，實現(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á)估計與差異表達(dá)相結(jié)合，來揭示總基因表達(dá)中每種異構(gòu)體的比例變化

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

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

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

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

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

rMATS 是通過比較用連接點的 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）比對等分析。

  1. DNA 甲基化

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

  1. 染色體特征

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

  1. MicroRNA

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

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

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

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

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容