A survey of best practices for RNA-seq data analysis
RNA-seq數據分析指南
內容
前言
各位同學/老師，大家好，現在由我給大家講講我的文獻閱讀報告！
A survey of best practices for RNA-seq data analysis ,我把它叫做RNA-seq數據分析指南。這篇文章是由佛羅里達大學等單位的研究人員在1月26日發表在Genome Biology上的，該期刊的影響因子有10.8分。這是這篇文章的通訊作者，應該挺靠譜的。
新一代測序技術在爆炸式發展的同時，也衍生出許多其他技術創新。RNA-Seq就是其中之一，這項技術使我們對細胞發育及其調控機制的理解，達到了前所未有的深度和廣度。RNA-seq可以獲得相當驚人的數據量，而這恰恰是一柄雙刃劍。豐富的數據量蘊含着大量的寶貴信息，但這樣的數據需要復雜的生物信息學分析，才能從中提取到有意義的結果。
正因如此，數據分析可以說是RNA-seq的重中之重。RNA-seq有非常廣泛的應用，但沒有哪個分析軟件是萬能的。科學家們一般會根據自己的研究對象和研究目標，采用不同的數據分析策略。現在人們已經發表了大量的RNA-seq和數據分析方案，對於剛入門的新手來說難免有些無所適從。這篇文章概述了RNA-seq生物信息學分析的現行標准和現有資源，為人們提供了一份RNA-seq數據分析指南，可以作為開展RNA-seq研究的寶貴參考資料。
這份指南覆蓋了RNA-seq數據分析的所有主要步驟，比如質量控制、讀段比對、基因和轉錄本定量、差異性基因表達、功能分析、基因融合檢測、eQTL圖譜分析等等。研究人員繪制的RNA-seq分析通用路線圖(標准Illumina測序)，將主要分析步驟分為前期分析、核心分析和高級分析三類。前期預處理包括實驗設計、測序設計和質量控制。核心分析包括轉錄組圖譜分析、差異基因表達和功能分析。高級分析包括可視化、其他RNA-seq技術和數據整合。研究人員在文章中探討了每個步驟所面臨的挑戰，也評估了一些數據處理方法的潛力和局限性。此外，他們還介紹了RNA-seq數據與其他數據類型的整合，將基因表達調控與分子生理學和功能基因組學關聯起來，這種研究方式如今越來越受到研究者的歡迎。這篇文章在結尾處介紹了一些為轉錄組領域帶來改變的新技術，特別是單細胞RNA-seq和長讀段測序技術帶來的機遇和挑戰。
背景
高通量測序平台
從這張圖，橫軸是年份，縱軸是測序儀的通量，圈里面的數字代表測序讀長。我們可以看到測序儀的通量和讀長都在增加了，其中Pacfic Biosciences 的三代測序儀讀長最長，可以達到14K，illumina的測序儀通量最大，HiSeq X Ten的通量可以達到1.8T。現在Illumina生產的測序儀占主要的市場份額，我們來單獨看看他的測序儀參數。
從左到右，測序儀的通量逐漸增大，它們適合不同的樣品和測序目的。Miseq通量比較低，適合宏基因組等微生物測序；Hiseq通量太高了，如果你送去公司測序，他們一般要20~30天才能返回數據，是因為公司要累積很多樣品，加標簽后一起測序。
高通量測序技術的應用
這幅圖的橫軸是年份，縱軸是高通量技術應用的代表性文章的引用量。不同的應用技術用顏色進行分類，數據點的大小跟發表率（引用率/月）成正比。可以看出RNA-Seq測序技術的應用最為廣泛。
實驗設計
RNA-seq到底測的是什么？**
mRNA在生物個體內RNA的組分中只占很小的一部分，rRNA占絕大多數。**一般我們說
RNA-seq指的都是mRNA-seq，后面的流程也都是主要針對mRNA-seq數據分析的。在科學家們的努力下，可以把那些非編碼RNA提取出來建庫，進行測序。
一個成功的RNA-seq研究，起決定性因素的是一個好的實驗設計。還依賴於建庫的類型、測序深度和設置適於的生物重復。並且盡量減少測序本身以外帶來的數據誤差。
文庫構建
1.一般生物體中的的RNA中，rRNA占絕大多數，含量超過90%，而mRNA的含量在1-2%左右。對於真核生物，一般使用加poly(A)選擇性富集mRNA或者而原核生物則是通過去除rRNA；
2.是否建stand-preserving庫；
3.對於Illumina，測序插入片段一般小於500bp。確定合適長度的插入片段是后續測序和分析的關鍵；
4.單端還是雙端測序毫無疑問的是，單端測序更便宜一些，如果你研究的某個物種的基因表達水平，並且它的轉錄組已經被注釋很好了，單端測序產生的數據量一般是足夠的了。
雙端測序呢，它的讀長更長，更適合於那些沒有被注釋的轉錄組物種的研究，便於其轉錄本的從頭拼接。
測序深度和重復數
這里的測序深度和重復數的設置呢。
- 測序深度和數據量不能一味加大。測序深度雖然可以提高基因定量和檢測的敏感性（低豐度的基因），但是同時也會增加一些噪音和一些無用的轉錄本。
2.重復數，很重要。（舉例子：比較兩個班的語文成績，抽取一個人代表一個班來比較，顯然是不科學的。）增加重復數可以減少實驗誤差，對提高結果的可靠性，是非常有意義的。
誤差分為技術誤差和生物學差異。
技術誤差-可以通過選擇最優化的實驗測序程序；生物誤差-三個生物學重復是最基本的啦。
然后呢，設定生物學重復對差異基因的檢出率（真陽性率 TPR）的提高具有明顯效果。上面說增加測序深度可以檢測到低豐度基因，但是對任何樣品來說的當測序深度增加再增加，它就會到達平台期。
由於科研經費有限，無法無限制地增加樣本數或數據量。
所以在生物學重復數和單個樣本測序量上必須找到平衡點。在總數據量不變的情況下，將總數據量分配到更多的生物學重復樣本中，差異分析結果的可靠性在不斷提升。
對於RNA-seq，生物學重復數的價值要大於單個樣本測序量。但增加生物學重復的樣本數，意味着要增加建庫費用。因此，即使總數據不變，設置過多的生物學重復也是不合理的。
我們最終確定設置多少生物學重復還是需要看樣本個體之間的差異大不大，這點我們一般都很清楚，在測序之前，如果你所研究的現象在兩個實驗樣本之間差異很穩定的話，就可以少設置一些重復，差異不穩定的話有時候設置10個/20個都不夠。具體問題具體分析！！
3. 測序深度（Sequencing depth），也叫乘數，指每個鹼基被測序的平均次數，是用來衡量測序量的首要參數。研究表明，增加測序深度，測序量從1.6M條reads增加到20M條reads,（75bp）但到10M條reads時就已經達到平衡了，80%的雞轉錄本被檢測到。在此基礎上增加測序量，它們會比對到已經存在的轉錄本上。
因此即使提高測序深度，低表達水平的基因的檢測是比較困難的。並且提高測序深度確實能夠增加基因差異表達的敏感度，但是並不能保證檢測到的差異具有生物學意義。
這篇2012年BMC Genomics的方法學文章，影響因子3.98至今被引用次數 87次。
直觀一些說，如果某個基因在RNA-seq結果顯示差異表達，但QPCR結果表明這個基因表達差異不顯著，
可以認為這個基因RNA-seq結果為假陽性；反之，這個結果就是真陽性。
生物學重復對差異表達分析的影響
如圖所示，在單樣本測序量保持不變的情況下，隨着生物學重復(n)的提高，差異分析的假陽性率（FPR）基本穩定，但真陽性率（TPR）在不斷提高。也就是說提高生物學重復數，實驗對差異基因的檢測更加敏感，那些差異倍數較小或差異量較低的差異表達基因（此類基因的差異檢測難度較大）能夠更加容易被檢測到。
如表2、3所示，在一定的生物學重復數（ n）的情況下，隨着單樣本測序量（Depth）的提高（25% → 100%），真陽性率（TPR）都只有有限的提高。例如在n=3的情況下，單個樣本的測序量從25%提高到100%，TPR僅僅從6.24%提高到8.95%。在表3中，如果Depth等於25%不變，當n從2提高到12，TPR的提高則是非常明顯的。因此測序深度對結果改善效果不如增加生物學重復。
總數據量不變，生物學重復數與單樣品測序量最佳組合
如果保持總測序量不變（即如果生物量重復數為n,則單個樣品的測序量降低為1/n,總數據量為n*1/n=1,保持不變）。如圖A,灰色實線代表不同的生物學重復數（n）和單樣本數據量（1/n）組合的情況下，真陽性率（TPR）的變化。結果表明，隨着n的提高，TPR率不斷提高。例如n=2,TPR約為3%，如果n=6,TPR則提高到22%。
同時我們也可以對“單樣本測序量對差異表達分析的影響”再進行深入觀察。
如果n保持不變，但單個樣本的數據量不斷降低，TPR的降低十分緩慢。例如，n=3,單個樣本的數據量從100%降低到15%，TPR的值一直處於平台期，僅僅從9%降低到5%。 但是不同的生物學重復數和單樣本測序量的組合，對假陽性率（ FPR）的影響卻較小。如圖 B，灰色實線代表不同生物學重復數（n）和單樣本數據量（1/n）組合的情況下，真陽性率（FPR）的變化。雖然 n 從2 變化到 96，FPR 基本沒有太大變化。
從圖中我們很容易發現，基於負二項分布的差異分析檢驗（P value），FPR 對生物學重復數和單個樣本數據量均不敏感，始終保持低於 0.1%水平。或者說，這個算法對 FPR 的控制還是非常理想的。
討論
隨着測序單價的下降，目前市場上 RNA-seq 類項目的單樣本測序量正在不斷提高。以 2G，PE100 測序的表達譜項目為例，其對應的測序量為 20M 條 reads。如果一條長度為 1kbp 的低表達基因的表達量為 RPKM=0.5，其理論上可以檢測到的 reads 數為 20×0.5=10。所以低豐度基因的檢測，對 RNA-seq 這個技術來說並非最大問題。
第二個問題“轉錄本表達量的高低變化”比“轉錄本的有無”更具有普遍的生物學意義。雖然個別基因的表達量變化程度，可以使用 Qpcr 來驗證。但我們往往也使用所有差異基因來統計某些規律。例如使用差異基因的 pathway 富集分析來尋找與性狀相關的 pathway。如果在全局水平的差異基因集並不可靠，那么 pathway富集分析得出的結論的可靠性自然也受到影響。而全局水平的差異基因數量巨大，是難以使用 Qpcr 驗證的。因此，定量以及差異分析的准確性是在 RNA-seq 中更值得關心的問題。
測序設計
RNA-seq文庫的制備和測序過程：RNA碎裂，cDNA合成，接頭連接，PCR擴增，加標簽（多樣品混合測序），上泳池測序;
如何減少誤差：
1. 使用末端帶隨機核酸的接頭或者使用化學碎裂法代替RNAse III碎裂法；
2. 不同批次實驗或者不同runs。
a. 如果樣品太多在一個批次或者一個run跑不完，為了避免技術誤差造成太大的實驗誤差，要把樣品隨機分配到每個批次或runs中;(到底怎么設計，我們要討論一下！！)
b. 如果你的樣品是多樣品混合測序，每個樣品要單獨加上標簽，每個lanes要保證足夠的測序深度，為了保證所有的樣品在每個lane中都有。如果送給公司去做的話，我們要選擇建庫水平好些的，並且要求他們這么去做，應該會更好。
質量控制
重復數：
技術重復（spearman秩相關系數R²>0.9）;生物重復（主成分分析PCA）
原始數據的質量控制：
原始數據回來后，你做完備份以后，做的第一件事情就是看看數據質量如何，一般來自llumina測序平台用軟件FastQC看；其他平台的數據用軟件NGSQC。一般會有原始數據的序列質量，GC含量，存在的接頭以及K-mers子串圖並且重復序列太多的reads。
並且reads 3‘末端的質量低於前段，原因是隨着測序讀長的增加，酶活性下降，熒光強度也在下降，因此測序數據質量逐漸降低乃是自然趨勢。常用的數據過濾的軟件有FASTX-Toolkit and Trimmomatic，其他還有許多，你也可以自己寫代碼處理數據。
Reads比對后的質量控制(評估比對質量的指標)：**比對上的reads占總reads的百分比； Reads比對到外顯子和參考鏈上的覆蓋度是否一致；比對到基因組序列：多重比對reads？比對到轉錄組序列：來自未被注釋的轉錄本的reads會丟失； 產生更多的多重比對reads; 轉錄本被定量以后，應該看一下GC含量和基因長度偏差，確定定量的方法是否適用。
轉錄本分析
把所有樣本的reads混合用於轉錄本的拼接。二代測序的轉錄組reads用於拼接還是存在一些問題的，最終拼接結果不太理想。一個轉錄本的拼接結果會是10~100contigs。三代測序的讀長直接可以把一個轉錄本讀完了，完全不需要拼接。
RPKM/FPKM/TPM用來表示RNA-seq基因表達水平的值；對於單端測序RPKM和FPKM值是一樣的，FPKM可以轉換成TPM。Cufflinks(支持雙端測序數據，並且需要GTF格式的注釋文件)定量算法有……
提一個問題，有那么多軟件到底怎樣才是好的，選哪個軟件好呢？
功能分析
功能分析是標准轉錄組分析流程的最后一步，分析差異表達基因的分子功能和代謝通路。
其他RNA-seq應用
小RNA：
1.小RNA的長度通常在18~34個鹼基，包含了miRNAs, short-interfering RNAs (siRNAs)，PIWI-interactingRNAs (piRNAs)以及其他種類的**。
- sRNA-seq libraries are rarely sequenced as deeply as regular RNA-seq libraries because of a lack of complexity, with a typical range of 2–10 million reads.
3. 小RNA的數據分析流程跟常規RNA的分析流程不同。
- miRTools 2.0 ,
a tool for prediction and profiling of sRNA species, uses by default reads that are 18–30 bases long
5. 比對到參考基因組上，比對軟件有：
Bowtie2 ,STAR , or Burrows-Wheeler Aligner (BWA) PatMaN and MicroRazerS map short sequences
多種數據整合分析
1. Moreover, the combination of RNA-seq and re-sequencing can be used both to remove false positives when inferring fusion genes and to analyze copy number alterations.
2. The statistically significant correlations that were observed, however, accounted for relatively small effects. (DNA methylation)
3. ….
4. 一些分析軟件：CORNA, MMIA,, MAGIA, and SePIA;
5.代謝組和轉錄組數據結合進行通路分析，有一些軟件：MassTRIX, Paintomics, VANTED v2, and SteinerNet
整合多種組學數據分析還不是很成熟，但是仍有一些軟件可以用。
展望
RNA-seq技術已經成為轉錄組分析的標准方法。其相對應的技術和數據分析工具還在不斷地發展。
對低表達的基因的定量仍是一個等待解決的問題; 三代測序技術，Smart-seq和Smart-seq2應用於轉錄組測序，所需要的樣品量少，並且可以測定單細胞內的RNA表達水平; Pacbio 技術可以直接測得接近全長的轉錄本，可以有效解決二代測序技術拼接較為零碎以及潛在嵌合拼接的問題；
目前的瓶頸：價格高（建庫價格和測序價格）;
（1）需要多種長度的文庫；
（2）測序通量有限；Pacbio新推出的sequel測序儀，比舊版本測序儀，通量提高了7倍（測序芯片的波導空數量從15萬，提升到100萬。所以，有望進一步提高Pacbio在轉錄組De Novo中的應用面。
PS: 綜述類的文章主要是講解原理和優缺點，以及當前趨勢和建議，不可能說看個綜述就學會了RNA-seq的全套分析，但是不看綜述，不了解原理，所有的分析卻都是無源之水無本之木。
還有，綜述給的建議，包括樣本數，數據量這些東西都是實驗設計之初就得考慮的，但是大多數生信工程師拿到數據很多東西就沒得改變了，所以總是有人問如果沒有重復的轉錄組數據用什么R包來找差異這樣的問題。請大家分清楚理想和現實的差別。
摘要：
佛羅里達大學、加州大學Irvine分校等單位的研究人員在一月二十六日的Genome Biology雜志上發表文章，概述了RNA-seq生物信息學分析的現行標准和現有資源，為人們提供了一份帶有注釋的RNA-seq數據分析指南。這將成為開展RNA-seq研究的寶貴參考資料。
生物通報道：新一代測序技術在爆炸式發展的同時，也衍生出許多其他技術創新。RNA深度測序（RNA-Seq）就是其中之一，這項技術使我們對細胞發育及其調控機制的理解，達到了前所未有的深度和廣度。盡管研究細胞RNA並不是什么新鮮事，但RNA-Seq的出現大大拓展了轉錄組研究的規模，取得了累累碩果，這些是傳統技術難以企及的。
RNA-seq可以獲得相當驚人的數據量，而這恰恰是一柄雙刃劍。豐富的數據量蘊含着大量的寶貴信息，但這樣的數據需要復雜的生物信息學分析，才能從中提取到有意義的結果。正因如此，數據分析可以說是RNA-seq的重中之重。
RNA-seq有非常廣泛的應用，但沒有哪個分析軟件是萬能的。科學家們一般會根據自己的研究對象和研究目標，采用不同的數據分析策略。現在人們已經發表了大量的RNA-seq和數據分析方案，對於剛入門的新手來說難免有些無所適從。
佛羅里達大學、加州大學Irvine(****爾灣)分校等單位的研究人員在一月二十六日的Genome Biology雜志上發表文章，概述了RNA-seq生物信息學分析的現行標准和現有資源，為人們提供了一份帶有注釋的RNA-seq數據分析指南。這將成為開展RNA-seq研究的寶貴參考資料。
這份指南覆蓋了RNA-seq數據分析的所有主要步驟，比如質量控制、讀段比對、基因和轉錄本定量、差異性基因表達、功能分析、基因融合檢測、eQTL圖譜分析等等。研究人員繪制的RNA-seq分析通用路線圖（標准Illumina測序），將主要分析步驟分為前期分析、核心分析和高級分析三類。前期預處理包括實驗設計、測序設計和質量控制。核心分析包括轉錄組圖譜分析、差異基因表達和功能分析。高級分析包括可視化、其他RNA-seq技術和數據整合。
研究人員在文章中探討了每個步驟所面臨的挑戰，也評估了一些數據處理方法的潛力和局限。此外，他們還介紹了RNA-seq數據與其他數據類型的整合。這種數據整合可以將基因表達調控與分子生理學和功能基因組學關聯起來，如今越來越受到研究者的歡迎。
這篇文章在結尾處介紹了一些為轉錄組領域帶來改變的新技術，特別是單細胞RNA-seq和長讀取測序技術帶來的機遇和挑戰。
2015年年初，RNA-Seq的數據分析方法如雨后春筍般涌現。三月份，Nature集團旗下刊物發表了三篇介紹RNA-Seq數據分析新方法的文章，一篇發表在《Nature Methods》上，另外兩篇發表在《Nature Biotechnology》上。這三篇文章有一位共同的作者，那就是約翰霍普金斯大學計算生物學中心的Steven Salzberg，生物信息學和計算生物學領域的傑出科學家。Salzberg通過這些文章中分別介紹了三種新工具：HISAT、StringTie和Ballgown。這些工具可以取代之前開發的早期工具，為RNA-Seq提供了全新的數據分析方法，從原始數據讀取到差異表達分析。（更多詳細信息參見：三篇文章介紹RNA-Seq數據分析的新工具）
RNA測序究竟有多可靠呢？由美國FDA牽頭的測序質量控制（SEQC）項目對RNA測序的准確性、可重現性和信息含量進行了綜合性評估。其初步調查結果發表在2014年09月的Nature Biotechnology雜志上，石樂明教授是這篇文章的通訊作者之一。研究人員用RNA參照樣本在全球多個實驗室的Illumina HiSeq、Life Technologies SOLiD、Roche 454平台上進行檢測，主要評估RNA測序在接頭區域和差異性表達譜中的表現，並將其與芯片和定量PCR（qPCR）進行比較。研究表明，數據分析的算法會對RNA測序產生很大影響，不同算法生成的轉錄本數據存在很大差異。（更多詳細信息參見：石樂明教授Nature子刊：RNA測序到底可不可靠）
前幾天，浙江大學和哈佛大學的研究人員在Cell Reports雜志上發表了一項單細胞mRNA-seq研究。基因表達變異是小鼠胚胎干細胞（ESC）的一個重要特征，但人們一直不清楚這背后的具體原因。研究人員通過分析小鼠胚胎干細胞發現，這些細胞表現出的異質性是血清培養造成的。他們在其中鑒定了高度變異的基因簇，以及獨特的染色質狀態。研究顯示，雙價基因（bivalent gene）更容易出現表達變異。進一步研究表明，無血清培養可以減少小鼠ESC的異質性和轉錄組變異。這意味着，細胞內的網絡變異大多是細胞外的培養環境造成的。（更多詳細信息參見：浙大80后教授發表單細胞RNA測序研究）