翻譯自：From RNA-seq reads to differential expression, Oshlack et al. Genome Biology 2010, 11:220
高通量測(cè)序技術(shù)，也就是下一代測(cè)序技術(shù)已經(jīng)成為現(xiàn)代生物學(xué)研究的一個(gè)較為常規(guī)的實(shí)驗(yàn)手段了。這一技術(shù)的發(fā)展極大地推動(dòng)了基因組學(xué)，表觀基因組學(xué)以及翻譯組學(xué)的研究。RNA-seq通過(guò)測(cè)定穩(wěn)定狀態(tài)下的RNA樣品的序列來(lái)對(duì)RNA樣品進(jìn)行研究，從而避免了許多之前研究手段的不足，比如象基因芯片或者PCR就需要背景知識(shí)。而且RNA-seq還可以觸及以前無(wú)法研究的領(lǐng)域，比如復(fù)雜結(jié)構(gòu)的轉(zhuǎn)錄體。RNA-seq可以應(yīng)用于以下幾個(gè)方面的研究，1. SNPs；2. novel transcripts；3. alternative splicing；4. RNA editing。無(wú)論如何，使用RNA-seq最多的還是比較兩組樣品基因水平表達(dá)差異，比如野生型與突變型，用藥組與對(duì)照組，不同組織之間，癌細(xì)胞與正常細(xì)胞，等等。我們把這種基因水平差異表達(dá)，簡(jiǎn)稱為DE (differential expression，注，不是ED啊???)。
常用的RNA-seq操作平臺(tái)有Illumina GA/ HiSeq, SOLiD 還有Roche 454。它們都是提取RNA后，純化，打碎，逆轉(zhuǎn)錄成cDNA，然后測(cè)序。測(cè)序的結(jié)果被稱為short reads，短序。通常一個(gè)短序的長(zhǎng)度為25-300bp之間。如果測(cè)序只測(cè)一端可能會(huì)帶來(lái)比對(duì)時(shí)的困難，于是這些操作平臺(tái)提供了兩端都測(cè)的辦法，這樣的結(jié)果成對(duì)出現(xiàn)，中間有一定的間隔，但是因?yàn)闇y(cè)序長(zhǎng)度一下子提高了一倍，所以比對(duì)會(huì)精準(zhǔn)很多。人們把這種測(cè)序結(jié)果稱為’paired-end’ reads，成對(duì)短序。一般來(lái)講，測(cè)序結(jié)果會(huì)直接轉(zhuǎn)換成一行一行的由字母組成的短序列，可能是fasta,fastq等等不同格式。
然而，這一技術(shù)產(chǎn)生的海量數(shù)據(jù)分析卻給生物學(xué)家?guī)?lái)了難題。一個(gè)測(cè)序的結(jié)果文件少則幾Gb，多則幾十Gb，單獨(dú)對(duì)比拼接，就會(huì)用去幾個(gè)小時(shí)，而后再得出差異表達(dá)的結(jié)果，其耗時(shí)耗力，并非實(shí)驗(yàn)生物學(xué)家可以應(yīng)付得了的。于是生物信息學(xué)的研究人員努力做出一些軟件，以降低結(jié)果分析的難度。但是，即使這樣，還是必須對(duì)分析過(guò)程有個(gè)較為細(xì)致地了解，才能正確地使用這些軟件，從而得到比較接近事實(shí)的結(jié)果。
一般的來(lái)講，RNA-seq后DE的工作流程是這樣的（圖1），首先，將短序映射到基因組相應(yīng)的位置上去，其次，對(duì)映射的結(jié)果進(jìn)行基因水平，外顯子水平，以及轉(zhuǎn)錄水平的拼接，而后對(duì)結(jié)果進(jìn)行數(shù)據(jù)統(tǒng)計(jì)，標(biāo)準(zhǔn)化之后生成表達(dá)水平報(bào)告文件，最后由生物學(xué)者依據(jù)系統(tǒng)生物學(xué)相關(guān)知識(shí)，來(lái)對(duì)數(shù)據(jù)結(jié)果進(jìn)行分析。
RNA-seq分析工作流程
不同步驟涉汲的軟件和方法：
分析步驟方法
軟件mappingGeneral alignerGMAP/GSNAP
BFAST
BOWTIE
CloudBurst
GNUmap
MAQ/BWA
PerM
RzaerS
Mrfast/mrsfast
SOAP/SOAP2
SHRiMP
De novo annotatorQPALMA/GenomeMapper/PALMapper
SpliceMap
SOAPals
G-Mo.R-Se
TopHat
SplitSeek
De novo transcript assemblerQases
MIRA
SummarizationIsoform-basedCufflinks
ALEXA-seq
Gene-basedCount exons only
Exon junction libraries
Normalizationlibrary size
RPKM: reads per kilobase of exon model per million mapped readsERANGE
TMM: trimmed mean of M-valuesedgeR
Upper quartileMyrna
Differential expressionPoisson GLM (generalized linear model)DEGseq
Myrna
Negative binomialedgeR
DESeq
baySeq
Systems biologyGene Ontology analysisGOseq
映射至基因組（Mapping）
第一步的工作是比對(duì)(alignment)。對(duì)于RNA-seq的比對(duì)，從來(lái)都不是一件容易的事情。其難點(diǎn)如下：
沒(méi)有很好的比對(duì)模板?，F(xiàn)在的比對(duì)模板都是基因組模板，而不是真正的轉(zhuǎn)錄組模板，也就是說(shuō)，這對(duì)本來(lái)就不是很長(zhǎng)的短序來(lái)說(shuō)，它很有可能是界于兩個(gè)exon之間。我們?cè)诒葘?duì)junction的時(shí)候，一般還是假設(shè)它如果沒(méi)能在基因組模板中找到合適的位置的時(shí)候，才考慮它是否是界于junction上。這種人為的假設(shè)可能并不準(zhǔn)確。
SNPs，堿基插入，刪除，錯(cuò)配，或者質(zhì)量不高的測(cè)序結(jié)果，從模板至比對(duì)序列本身，都存在著比基因比對(duì)更為復(fù)雜的問(wèn)題。
短序可能會(huì)有多個(gè)100％的匹配位點(diǎn)。
有些基因組可能需要龐大的內(nèi)存空間。
為了解決最后一個(gè)問(wèn)題，人們使用了很多辦法，但基本上都會(huì)基于事先建立的引索庫(kù)。即所謂“啟發(fā)式”比對(duì)(heuristic match)。首先使用一定長(zhǎng)度的（通常是11個(gè)堿基）的序列做為索引用的關(guān)鍵字，在匹配這一索引字之后，就很大程度地縮小了其需要匹配的模板范圍。但是這一辦法的問(wèn)題在于不容易解決問(wèn)題2中的空格，錯(cuò)配問(wèn)題。所以在很多軟件使用時(shí)，會(huì)要求人工確認(rèn)高保真區(qū)，以及最高允許2?3個(gè)錯(cuò)配。
現(xiàn)在比較快的“啟發(fā)式”比對(duì)主要有兩種算法，一種是哈希表(hash table)，一種是BW壓縮轉(zhuǎn)換(Burrows Wheeler transform, BWT)。前者速度快，但是對(duì)內(nèi)存要求比后者要高。
對(duì)于問(wèn)題3，一般而言，大部分軟件使用的辦法是只保留一個(gè)匹配位點(diǎn)，其中，有些是只保留第一個(gè)匹配位點(diǎn)，有些是按照概率分布選取保留的位點(diǎn)。當(dāng)然，前面已經(jīng)提到過(guò)，可以使用paired-end read來(lái)盡量避免問(wèn)題3的出現(xiàn)。
對(duì)于問(wèn)題1，可以使用外顯子庫(kù)來(lái)確定junction reads。有兩種辦法，一種是依靠已知的外顯子庫(kù)來(lái)構(gòu)建，另一種辦法就是依據(jù)已經(jīng)匹配好的短序來(lái)構(gòu)建外顯子庫(kù)(de novo assembly of transcriptome)。后者的不足是運(yùn)算量大，對(duì)測(cè)序覆蓋范圍要求高，最好是使用paired-end reads。
還有人發(fā)現(xiàn)，對(duì)于ploy(A)的處理會(huì)減少不能映身的短序數(shù)。比如，Pickrell et al.就發(fā)現(xiàn)，對(duì)于46bp的Illumina reads，87％的短序可以映射至模板，7％可以映射至junction library。如果對(duì)那些不能映射的短序，將在頭或者尾含有的超過(guò)連續(xù)4個(gè)的A或者T去除，就可以得到約0.005%的映射。
綜合評(píng)價(jià)（Summarizing mapped reads）
這一步，主要是基本于不同水平（外顯子水平，轉(zhuǎn)錄水平，或者基因水平）進(jìn)行統(tǒng)計(jì)。最簡(jiǎn)單的辦法就是統(tǒng)計(jì)落在每個(gè)外顯上的短序數(shù)。但是有研究表明，很多（可能超過(guò)15％）的短序會(huì)落在外顯子兩側(cè)，這會(huì)影響統(tǒng)計(jì)的結(jié)果。另一種辦法就是統(tǒng)會(huì)落在內(nèi)顯子區(qū)域的短序數(shù)。
無(wú)論如何，即使是基因水平的綜合評(píng)價(jià)，也還是有其它的一些問(wèn)題。比如overlapping的基因的統(tǒng)計(jì)。比如junction的統(tǒng)計(jì)。
標(biāo)準(zhǔn)化（Normalization）
標(biāo)準(zhǔn)化對(duì)于樣品內(nèi)及樣品間的比較而言是非常重要的。標(biāo)準(zhǔn)化被分為兩類，樣品內(nèi)及樣品間（between- and within-library）。
樣品內(nèi)標(biāo)準(zhǔn)化使得在同一樣品內(nèi)不得基因之間的表達(dá)差異變得有意義。最常用到的一個(gè)辦法就是使用落在同一基因內(nèi)的短序數(shù)除以單位基因長(zhǎng)度。比較常用的單位是RPKM (reads per kilobase of exon model per million mapped reads)。但是這一方法也受到樣品制備和測(cè)序方法的干擾。
而對(duì)于樣品間標(biāo)準(zhǔn)化，最簡(jiǎn)單而直接的辦法使用短序總數(shù)來(lái)平衡表達(dá)量。然而短序總數(shù)受測(cè)序深度的干擾，而且單個(gè)基因的短序數(shù)與實(shí)際的表達(dá)量并不一定會(huì)呈線性比較關(guān)系。人們又使用四分位(quantile normlization)標(biāo)準(zhǔn)化的辦法。但是有研究說(shuō)這一辦法并沒(méi)有實(shí)際的價(jià)值。還有提出使用對(duì)數(shù)分布法則(power law distributions)來(lái)進(jìn)行樣品間標(biāo)準(zhǔn)化。但沒(méi)有研究對(duì)這一處理方式進(jìn)行驗(yàn)證。
差異表達(dá)（Differential expression）
差異表達(dá)分析的最終目的是將那些差異表達(dá)的基因（外顯子等等）從海量數(shù)據(jù)中提取出來(lái)。最終的結(jié)果顯示一般來(lái)說(shuō)是表格化的，這一表格按照一定的規(guī)則排序，讓人們能夠盡可能簡(jiǎn)單地拿到想要的結(jié)果。
由于RNA-seq結(jié)果的離散性，人們一般都會(huì)使用統(tǒng)計(jì)模型來(lái)擬合實(shí)驗(yàn)得到的結(jié)果。一般而言，RNA-seq的結(jié)果是比較附合伯松分布(poisson distribution)的。這一結(jié)果得到了單通道Illumina GA測(cè)序結(jié)果的實(shí)驗(yàn)驗(yàn)證。但是，伯松分布分析結(jié)果常常在多組重復(fù)的樣品間帶來(lái)較高的假陽(yáng)性，因?yàn)樗凸懒松锶拥臉悠烽g誤差。所以RNA-seq如何設(shè)置重復(fù)是一個(gè)很重要的問(wèn)題。為了平衡重復(fù)樣品所帶來(lái)的誤差，人們使用了serial analysis of gene expression (SAGE) data。
現(xiàn)有的軟件一般都是針對(duì)較為簡(jiǎn)單的實(shí)驗(yàn)設(shè)計(jì)的。而對(duì)于復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)，比如說(shuō)成對(duì)樣品，時(shí)間依賴樣品等等，還沒(méi)有專門的，較好的解決方案。大多數(shù)都使用edgeR的線性模型來(lái)進(jìn)行分析。
后期系統(tǒng)生物學(xué)分析
簡(jiǎn)單地講，前景是廣闊的，但目前為止手段還是比較有限的，基本上就是GO分析。
