完整轉錄組分析（整理）

轉載來自https://zhuanlan.zhihu.com/p/393674599

寫的非常好 怕找不到留着自己看！如果作者不同意我會刪除。

前言

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

Stark et al. Nat Rev Genet(2019)

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

而我們一般的 RNA-seq 測序數據分析流程算法，基本上都是基於 short-read（短讀長）技術所產生的數據文件

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

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

Conesa et al. Genome Biology (2016)

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

1. 數據獲取：

• 包括實驗設計、測序設計以及數據下機后的 raw reads 數據的質控

1. 數據分析

• 在獲取到干凈的數據之后，可以進行 reads 的比對，然后進行基因表達的量化、差異表達分析、功能富集分析等

1. 高級分析

• 包括數據的可視化，其他小分子 RNA 分析、融合分析以及與其他類型的數據進行整合分析等

而我們分析的起始點，是從原始數據開始的，也就是獲取 raw reads 數據。通常這種高通量測序數據會保存為 FASTQ 格式的文件。

FASTQ 格式是一種以 ASCII 碼字符的形式保存生物序列及其對應的每個鹼基的質量的文本文件。

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

• 第一行以 @ 字符開頭，之后的字符為序列的標識符和描述信息
• 第二行為具體的序列
• 第三行以 + 符號開頭，之后可以可選地加上與第一行一樣的序列標識或描述信息
• 第四行為鹼基質量分數（Phred），其字符數量與第二行相等，每個字符表示對應鹼基的質量得分，例如

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

其中，鹼基質量值的編碼方式為

1. 先將鹼基錯誤率 P 進行負對數轉換，得到 Q 值 
2. 然后將 Q 值加上 33 或 64 得到的值所對應的 ASCII 碼即為鹼基質量分數

例如，錯誤率 P = 0.01，則 Q = 20，如果是 Phred33 則對應的質量為字符 5（53），如果是 Phred64 則對應的字符為 T(84)

分析流程

1. 數據獲取

一般情況下，如果自己有送樣檢測數據的話，測序公司會提供原始的 FASTQ 格式的數據。如果我們要使用別人文章中發表的公開數據，還需要從數據庫中下載對應的數據

例如，我們從 SRA 數據中下載的原始測序文件是 sra 格式，我們需要先使用工具將其轉換為 FASTQ 格式

2. 質量控制

主要在三個地方需要對數據的質量進行監控

• 獲取原始數據之后
• 比對完之后
• 表達定量之后

2.1 raw read

對 raw reads 數據進行質量控制，需要分析序列的質量、GC 含量、是否存在接頭、短重復序列的分布、測序錯誤以及 PCR 重復和污染

質控軟件：

• FastQC：用於分析 Illumina 測序平台的數據
• NGSQC：可應用於所有平台

一般來說，reads 的質量會朝着 3' 端遞減，如果鹼基的質量太低，我們需要刪除它以提高比對率

FASTX-Toolkit 和 Trimmomatic 兩個軟件可以用於切除低質量的鹼基和接頭序列

2.2 比對后

reads 通常需要比對到一個參考基因組或轉錄組，而比對的質量是評估測序准確率和是否存在 DNA 污染的一個重要指標

比對質量通常為比對到的 reads 數占總 reads 數的比例。

例如，比對到人類參考基因組的比對質量通常需要在 70-90%，且有大量的 reads 映射到一個相同的區間內。如果是比對到轉錄本上，由於可變剪切的影響，可以適當放寬比對質量

在外顯子和比對方向上的 read 覆蓋率的均一性，也是評估質量的重要指標。如果 reads 主要聚集在轉錄本的 3' 端，可能表明原始樣本的 RNA 質量較低

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

主要軟件有：Picard、RSeQC 和 Qualimap

2.3 定量后

在計算完表達的量化值之后，可以計算 GC 含量和基因長度的誤差，在必要時可以使用標准化方法來進行校正

如果參考轉錄組注釋得很好，則可以分析樣本的生物構成，來評估 RNA 純化步驟的質量。例如，rRNA 和 small RNA 不能出現在 polyA longRNA 的制備中

NOISeq 和 EDASeq 等 R 包可以使用圖形來展示 count 數據的質量控制

2.4 可重復性

上面的質量控制都只是針對單個樣本的，此外，不同樣本之間的可重復性評估，對於評價整個數據集的質量也是至關重要的

技術重復樣本的可重現性一般很高（spearman ），但是生物學重復樣本之間並沒有明確的標准，取決於實驗系統的異質性。如果不同實驗系統之間存在差異基因，則同一條件下的生物學重復在主成分分析（PCA）中會被聚類在一起。

3. 序列比對

在對樣本的 raw reads 進行質控之后，就可以進行序列比對了，序列比對主要有三種策略，如下圖

Conesa et al. Genome Biology (2016)

如果有參考序列，根據參考序列的不同，可以分為

• 比對到基因組：使用間隔比對算法，如 TopHat、STAR 等，然后根據是否提供了注釋文件（GFF 格式文件，包含轉錄本位置信息），又可以分為轉錄本識別和轉錄本發現並進行定量分析
• 比對到轉錄組：使用非間隔比對算法，如 Bowtie 等，然后使用 RSEM 或 Kallisto 方法識別轉錄本並計算定量信息

如果沒有參考序列，則需要先把序列組裝成轉錄本，再將 reads 比對到組裝后的參考轉錄本上，然后使用 HTseq-count 等算法對轉錄本進行定量

3.1 轉錄本發現

使用 Illumina 技術檢測的 short reads 來發現新的轉錄本是 RNA-seq 分析中的一個挑戰。通常來說，短 reads 很少會跨越多個剪切位點，這就很難直接推斷出一個轉錄本的整體長度。

此外，轉錄的起始和終止位置也比較難識別，一些像 GRIT 的工具，通過合並 5' 端的信息可以提高異構體識別的准確性。其他如 Cufflinks、iReckon、SLIDE 和 StringTie 等方法，通過結合現有的注釋信息，作為一個可能的異構體列表

一些尋找基因的工具，如 Augustus，結合 RNA-seq 數據，可以更好的注釋蛋白編碼轉錄本，但是對非編碼轉錄本的性能更差。

3.2 De novo 轉錄本重構

在沒有轉錄本或轉錄本不全的情況下，可以對 reads 進行組裝來重構一份轉錄本。可選的方法很多，如 SOAPdenovoTrans、Oases、Trans-ABySS 或 Trinity

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

雖然，對於低表達的轉錄本進行組裝的可靠性較低，但是 reads 太多也會導致潛在的組裝錯誤和較長的時間消耗等問題。因此，在深度測序的樣本中，可以適當減少 reads 的數量

對於多樣本的比較，可以將所有樣本作為一個輸入來構建參考轉錄本，然后分別對每個樣本的 reads 進行比對

無論是使用參考序列還是從頭開始組裝，使用短 reads 的 Illumina 技術來完全重構轉錄組仍然是一個具有挑戰性的問題

4. 轉錄組定量

RNA-seq 最廣泛的應用就是用來評估基因和轉錄本的表達，這一應用主要是基於比對到轉錄組區間內的 reads 的數量

最簡單的方法是，使用 HTSeq-count 或 featureCounts 計算區間內的 reads 數來量化基因的表達。這種基因水平的（不是轉錄本水平）的量化方法使用的是 GTF 文件，這種文件包含外顯子和基因在基因組上的坐標。

但一般不能直接使用 read count 來比較基因的表達水平，因為該值會受到轉錄本長度、reads 總數以及測序偏差等因素的影響。所以需要先進行標准化，標准化方法有

1. RPKM/FPKM: 每百萬 reads 每一千鹼基對中包含的 reads 數

該方法先計算測序深度系數，即總 reads 數除以 一百萬，然后計算基因或轉錄本的長度（單位為 kb），標准化順序為先消除測序深度的影響，再消除長度的影響：

其中

• x 表示一個基因或轉錄本，或基因組上一段特定的區域
•  表示比對到 x 外顯子區域的 reads 數；
• R 表示當前樣本中包含的全部 reads 數
•  表示 x 外顯子區域包含的鹼基數（長度，bp）

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

1. TPM: 其與 RPKM 最大的區別是，標准化順序為先消除基因長度的影響，再消除測序深度的影響

首先，將 reads count 除以基因或轉錄本的長度（kb）得到 RPK(reads per kilobase)，然后將樣本中所有的 RPK 加起來除以 ，得到標准化系數，最后使用 RPK 除以標注化系數

其中

• x 表示一個基因或轉錄本，或基因組上一段特定的區域
•  表示比對到 x 外顯子區域的 reads 數
•  表示 x 外顯子區域包含的鹼基數（kp）
• N 表示基因或轉錄本總數

這樣，每個樣本的 TPM 總和是一樣的，便於比較樣本間的差異

目前，也有許多復雜的算法通過解決相關轉錄本共享 reads 的問題來評估轉錄本水平的表達，例如，Cufflinks 使用 TopHat 的比對結果，應用期望最大化算法來評估轉錄本的豐度。這一方法考慮到長度不同的基因的 reads 分布並不均勻等因素的影響。

還有其他算法也可以量化轉錄組的表達，例如 RSEM、eXpress、Sailfish 和 kallisto 等。這些方法允許轉錄本之間存在多比對的 reads，並輸出經測序偏差校正的樣本內歸一化值。

5. 差異表達分析

差異表達分析是對樣本間基因的表達值進行比較，雖然 RPKM、FPKM 和 TPM 標准化方法消除了測序深度和基因或轉錄本的長度因素的影響，但這些方法依賴於總的或有效的 reads 數，當樣本的具有異質性轉錄本分布或當高表達或差異表達的特征扭曲了 count 分布時，表現欠佳

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

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

NOISeq 這個 R 包提供了多種繪圖，來識別 RNA-seq 數據中的誤差來源，並應用相應的方法來標准化這些誤差

除了這些樣本內特異的標准化方法，還需要解決數據集之間的批次效應（不同實驗條件下產生的數據之間存在的差異），批次矯正方法有 COMBAT 和 ARSyN 等，雖然這些方法是針對芯片數據設計的，但是在 RNA-seq 數據中也有很好的效果

計算差異表達的方法有很多，有些方法，如 edgeR 將原始的 read counts 作為輸入，並在統計模型中加入了標准化，另一些方法，需要先對數據進行標准化，如 DESeq2 使用的是負二項分布作為參考分布，並提供了自己的標准化方法。

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

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

但是需要強調的是，在沒有足夠生物學重復的情況下，無法進行總體的推斷，因此任何 p 值計算都是無效的。

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

6. 可變剪切分析

可變剪接(Alternative Splicing) 是指轉錄形成的前體 RNA 通過去除內含子、連接外顯子而形成成熟 RNA 的過程，從而實現一個基因同時編碼多種蛋白質，實現生物功能多樣性

在不同組織或者發育的不同階段，可變剪接不是一成不變的，在特定的組織或條件下，通過連接不同的外顯子，會產生特定的剪接異構體（isoform）。有大量的研究發現，可變剪接的變化與癌症等多種疾病相關，所以研究可變剪接在不同組織中的作用是非常有意義的。

轉錄本水平的差異表達分析可以潛在地檢測同一基因的轉錄異構體表達的變化，已經有一些算法應用於 RNA-seq 數據的中進行可變剪切分析

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

1. 異構體表達估計與差異表達相結合，來揭示總基因表達中每種異構體的比例變化

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

所有這些方法通常都受限於短讀長測序的內在局限性，無法在異構體水平上進行准確識別

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

其基本假設為：可以在外顯子及其連接點的信號中追蹤異構體表達的差異。

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

rMATS 是通過比較用連接點的 reads 定義的外顯子 inclusion levels 表達水平的差異

可變剪接作為轉錄后修飾的重要環節，對細胞的活性和疾病的發生發展都有廣泛的作用。盡管目前軟件很多，但是軟件間的比較卻比較少。作者系統地比較了10款軟件，利用了一致性（consistency）、可重復性（reproducibility）、精確度（precision）、召回率（recall）和錯誤發現率（false discovery rate）、報道的可變剪接基因的一致性（agreement upon reported differentially spliced genes）與功能富集分析一致性等方面來評價不同的可變剪接軟件。

軟件主要是三大類：

• exon-base ：DEXseq、edgeR、JunctionSeq、limma
• isoform-base：cuffdiff2、diffSplice2
• event-base：dSpliceType、MAJIQ、rMATS、SUPPA

可變剪接的檢測策略
isoform-base : 先構建全長轉錄本和定量，再來做差異分析
count-base：將基因用單個統計單元來表示，並且進行差異分析

• exon-base ：利用exon來作為統計單元
• event-base ：利用junction region作為統計單元

各個軟件說明

isoform base

• cufflinks：首先構建overlap graphs，之后估計轉錄本豐度，最后檢查差異基因和isoform；
• DiffSplice：基於圖論的方法來進行分析，首先基於比對構建轉錄本，然后對不同的path來定量豐度，最終鑒定出可變剪接體。

count base

• exon-base：將序列分配到不同的特征上，例如exon或者junction；這個只能分析已知的特征，而不能推斷出新的可變剪接事件；
• isoform-base：通過計算每個事件(PSI)值中的拼接百分比來量化剪接事件，PSI表示該值測量從包含該事件特定形式的基因中表達的mRNA的分數。

7. 融合分析

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

融合基因通過基因失調、融合產生嵌合體蛋白這兩種機制引發癌症的發生。

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

Liu et al Nucleic Acids Research, 2016

沒有哪一個算法具有明顯的優勢，整體來看，SOAPfuse 可能會好一些，FusionCatcher 和 JAFFA 其次。

8. 功能注釋

標准的轉錄組分析的最后一步，是使用差異表達基因來進行功能或通路的注釋。最常用的兩類方法是：

• 基於超幾何分布的過表達富集分析
• GSEA 富集分析

一些工具，如 GOseq 考慮了基因長度等因素對差異表達結果的影響，並使用超幾何分布進行富集分析，GSVA 和 SeqGSEA 使用類似 GSEA 的方法進行功能富集

功能富集需要預先定義的基因集合或通路，包括 GO、KEGG、Reactome 等數據庫。

通過在蛋白質數據庫（例如 SwissProt）和包含保守蛋白質結構域（例如 Pfam 和 InterPro）的數據庫中搜索相似序列，使用直系同源分析對蛋白質編碼的轉錄本進行功能注釋。而 Rfam 數據庫包含許多特征明確的 RNA 家族，例如 rRNA 或 tRNA，而 mirBase 和 Miranda 是專門研究 miRNA 的

9. 整合分析

將 RNA-seq 數據與其他類型的全基因組數據進行整合分析，使我們能夠將基因表達的調控與分子生理學和功能基因組學的特定方面聯系起來。

1. DNA 測序

將 RNA 測序和 DNA 測序相結合，可以進行 SNP、RNA 編輯和表達數量性狀基因座（eQTL）比對等分析。

1. DNA 甲基化

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

1. 染色體特征

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

1. MicroRNA

整合 RNA-seq 和 miRNA-seq 有可能揭示 miRNAs 對轉錄穩態水平的調節作用

1. 蛋白質組和甲基化組

RNA-seq 與蛋白質組學的整合是有爭議的，因為這兩種測量結果通常顯示出很低的相關性（~0.4）。盡管如此，蛋白質組學和 RNA-seq 的配對分析可用於識別新的異構體

轉錄組與代謝組數據的整合已被用於識別在基因表達和代謝物水平上受調控的通路，並且可以使用工具來可視化通路上下文的結果，如 MassTRIX、Paintomics、VANTED v2 和 SteinerNet

- END -

 

 

更新：更詳細了一點：https://www.jianshu.com/p/2055db183907