高通量測序技術,就是二代測序,已經成為現代生物學研究的一個較為常規的實驗手段。這一技術的發展極大地推動了基因組學,表觀基因組學以及翻譯組學的研究。RNA-seq 通過測定穩定狀態下的RNA樣品的序列來對RNA樣品進行研究,從而避免了許多之前研究手段的不足,比如象基因芯片或者 PCR 就需要背景知識。而且 RNA-seq 還可以觸及以前無法研究的領域,比如復雜結構的轉錄體。
RNA-seq可以應用於以下幾個方面的研究,1. SNPs;2. novel transcripts;3. alternative splicing;4. RNA editing。但歸根結底,RNA-seq最主要的分析還是篩選差異基因。
常用的RNA-seq操作平台有Illumina GA/ HiSeq, SOLiD 還有Roche 454。它們都是提取RNA后,純化,打碎,逆轉錄成cDNA,然后測序。測序的結果被稱為short reads。通常一個reads的長度為25-300bp之間。如果測序只測一端可能會帶來比對時的困難,於是這些操作平台提供了兩端都測的辦法,這樣的結果成對出現,中間有一定的間隔,但是因為測序長度一下子提高了一倍,所以比對會精准很多。人們把這種測序結果稱為’paired-end’ reads。一般來講,測序結果會直接轉換成一行一行的由字母組成的短序列,可以是fasta,fastq等不同格式。
然而,這一技術產生的海量數據分析卻給生物學家帶來了難題。一個測序的結果文件少則幾Gb,多則幾十Gb,單獨對比拼接,就會用去幾個小時,而后再得出差異表達的結果,其耗時耗力,並非實驗生物學家可以應付得了的。於是生物信息學的研究人員努力做出一些軟件,以降低結果分析的難度。但是,即使這樣,還是必須對分析過程有個較為細致地了解,才能正確地使用這些軟件,從而得到比較接近事實的結果。
一般的來講,RNA-seq后DE的工作流程是這樣的(圖1),首先,將短序映射到基因組相應的位置上去,其次,對映射的結果進行基因水平,外顯子水平,以及轉錄水平的拼接,而后對結果進行數據統計,標准化之后生成表達水平報告文件,最后由生物學者依據系統生物學相關知識,來對數據結果進行分析。
RNA-seq分析工作流程
不同步驟涉汲的軟件和方法:
分析步驟 | 方法 | 軟件 |
mapping | General aligner | GMAP/GSNAP |
BFAST | ||
BOWTIE | ||
CloudBurst | ||
GNUmap | ||
MAQ/BWA | ||
PerM | ||
RzaerS | ||
Mrfast/mrsfast | ||
SOAP/SOAP2 | ||
SHRiMP | ||
De novo annotator | QPALMA/GenomeMapper/PALMapper | |
SpliceMap | ||
SOAPals | ||
G-Mo.R-Se | ||
TopHat | ||
SplitSeek | ||
De novo transcript assembler | Qases | |
MIRA | ||
Summarization | Isoform-based | Cufflinks |
ALEXA-seq | ||
Gene-based | Count exons only | |
Exon junction libraries | ||
Normalization | library size | |
RPKM: reads per kilobase of exon model per million mapped reads | ERANGE | |
TMM: trimmed mean of M-values | edgeR | |
Upper quartile | Myrna | |
Differential expression | Poisson GLM (generalized linear model) | DEGseq |
Myrna | ||
Negative binomial | edgeR | |
DESeq | ||
baySeq | ||
Systems biology | Gene Ontology analysis | GOseq |
Mapping:
第一步的工作是alignment。對於RNA-seq 的 alignment,從來都不是一件容易的事情。其難點如下:
- 沒有很好的比對模板。現在的比對模板都是基因組模板,而不是真正的轉錄組模板,也就是說,這對本來就不是很長的短序來說,它很有可能是界於兩個 exon之間。我們在比對junction的時候,一般還是假設它如果沒能在基因組模板中找到合適的位置的時候,才考慮它是否是界於junction上。這種人為的假設可能並不准確。
- SNPs,鹼基插入,刪除,錯配,或者質量不高的測序結果,從模板至比對序列本身,都存在着比基因比對更為復雜的問題。
- reads 可能會有多個100%的匹配位點。
- 有些基因組可能需要龐大的內存空間。
為了解決最后一個問題,人們使用了很多辦法,但基本上都會基於事先建立的引索庫。即所謂“啟發式”比對(heuristic match)。首先使用一定長度的(通常是11個鹼基)的序列做為索引用的關鍵字,在匹配這一索引字之后,就很大程度地縮小了其需要匹配的模板范圍。但是這一辦法的問題在於不容易解決問題2中的空格,錯配問題。所以在很多軟件使用時,會要求人工確認高保真區,以及最高允許2-3個錯配。
現在比較快的“啟發式”比對主要有兩種算法,一種是哈希表(hash table),一種是BW壓縮轉換(Burrows Wheeler transform, BWT)。前者速度快,但是對內存要求比后者要高。
對於問題3,一般而言,大部分軟件使用的辦法是只保留一個匹配位點,其中,有些是只保留第一個匹配位點,有些是按照概率分布選取保留的位點。當然,前面已經提到過,可以使用paired-end read來盡量避免問題3的出現。
對於問題1,可以使用外顯子庫來確定junction reads。有兩種辦法,一種是依靠已知的外顯子庫來構建,另一種辦法就是依據已經匹配好的短序來構建外顯子庫(de novo assembly of transcriptome)。后者的不足是運算量大,對測序覆蓋范圍要求高,最好是使用paired-end reads。
還有人發現,對於 ploy(A) 的處理會減少不能 mapping 的 reads 數。比如,Pickrell et al.就發現,對於46bp的Illumina reads,87%的短序可以映射至模板,7%可以映射至junction library。如果對那些不能 mapping 的 reads,將在頭或者尾含有的超過連續4個的A或者T去除,就可以得到約0.005%的 mapping reads。
Summarizing mapped reads:
這一步,主要是基本於不同水平(外顯子水平,轉錄水平,或者基因水平)進行統計。最簡單的辦法就是統計落在每個外顯上的 reads 數。但是有研究表明,很多(可能超過15%)的 reads 會落在外顯子兩側,這會影響統計的結果。另一種辦法就是統會落在內顯子區域的 reads 數。
無論如何,即使是基因水平的綜合評價,也還是有其它的一些問題。比如overlapping的基因的統計。比如junction的統計。
Normalization:
標准化對於樣品內及樣品間的比較而言是非常重要的。標准化被分為兩類,樣品內及樣品間(between- and within-library)。
樣品內標准化使得在同一樣品內不得基因之間的表達差異變得有意義。最常用到的一個辦法就是使用落在同一基因內的 reads 數除以單位基因長度。比較常用的單位是RPKM (reads per kilobase of exon model per million mapped reads)。但是這一方法也受到樣品制備和測序方法的干擾。
而對於樣品間標准化,最簡單而直接的辦法使用 reads 總數來平衡表達量。然而 reads 總數受測序深度的干擾,而且單個基因的短序數與實際的表達量並不一定會呈線性比較關系。一些研究者推薦使用 quantile normlization,但是有研究說這一辦法並沒有實際的價值。還有提出使用對數分布法則(power law distributions)來進行樣品間標准化。但沒有研究對這一處理方式進行驗證。
Differential expression:
差異表達分析的最終目的是篩選差異表達的基因(外顯子等等)。最終的結果顯示一般來說是表格化的,這一表格按照一定的規則排序,讓人們能夠盡可能簡單地拿到想要的結果。
由於 RNA-seq 結果的離散性,人們一般都會使用統計模型來擬合實驗得到的結果。一般而言,RNA-seq的結果是比較附合伯松分布 (poisson distribution)的。這一結果得到了單通道Illumina GA測序結果的實驗驗證。但是,伯松分布分析結果常常在多組重復的樣品間帶來較高的假陽性,因為它低估了生物取樣的樣品間誤差。所以RNA-seq如何設置重復是一個很重要的問題。為了平衡重復樣品所帶來的誤差,人們使用了serial analysis of gene expression (SAGE) data。
現有的軟件一般都是針對較為簡單的實驗設計的。而對於復雜的實驗設計,比如說成對樣品,時間依賴樣品等等,還沒有專門的,較好的解決方案。大多數都使用edgeR的線性模型來進行分析。
Subsequence analysis:
簡單地講,前景是廣闊的,但目前為止手段還是比較有限的,基本上就是GO分析
From RNA-seq reads to differential expression, Oshlack et al. Genome Biology 2010.
Cited from: http://blog.sina.com.cn/s/blog_4ac87c050100xfdn.html