Charity Law1, Monther Alhamdoosh2, Shian Su3, Xueyi Dong3, Luyi Tian1, Gordon K. Smyth4 and Matthew E. Ritchie5
1The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia
2CSL Limited, Bio21 Institute, 30 Flemington Road, Parkville, Victoria 3010, Australia
3The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia
4The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia
5The Walter and Eliza Hall Institute of Medical Research, 1G Royal Parade, Parkville, VIC 3052, Melbourne, Australia; Department of Medical Biology, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia; School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Melbourne, Australia
簡單且高效地分析RNA測序數據的能力正是Bioconductor的核心優(yōu)勢。 RNA-seq分析通常從基因水平的序列計數開始，涉及到數據預處理，探索性數據分析，差異表達檢驗以及通路分析，得到的結果可用于指導進一步實驗和驗證研究。 在這篇工作流程文章中，我們通過分析來自小鼠乳腺的RNA測序數據，示范了如何使用流行的edgeR包載入、整理、過濾和歸一化數據，然后用limma包的voom方法、線性模型和經驗貝葉斯調節(jié)（empirical Bayes moderation）來評估差異表達并進行基因集檢驗。通過使用Glimma包，此流程得到了增進，實現了結果的互動探索，使用戶得以查看單個樣本與基因。 這三個軟件包提供的完整分析突出了研究人員可以使用Bioconductor輕松地從RNA測序實驗的原始計數揭示生物學意義。
RNA測序（RNA-seq）已成為基因表達研究的主要技術。其中，基因組規(guī)模的多條件基因差異表達檢測是研究者最常探究的問題之一。對于RNA-seq數據，來自Bioconductor項目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用于處理此問題的完善的統(tǒng)計學方法。
在這篇文章中，我們描述了一個用于分析RNA-seq數據的edgeR - limma工作流程，使用基因水平計數作為輸入，經過預處理和探索性數據處理，然后得到差異表達（DE）基因和基因表達特征（gene signatures）的列表。Glimma包(Su et al. 2017)提供的交互式圖表可以同時呈現整體樣本和單個基因水平的數據，使得我們相對靜態(tài)的R圖表而言，可以探索更多的細節(jié)。
此工作流程中所分析的實驗來自Sheridan等（2015）(Sheridan et al. 2015)，由三個細胞群組成，即基底（basal）、管腔祖細胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。細胞群皆分選自雌性處女小鼠的乳腺，每種都設三個生物學重復。RNA樣品分三個批次使用Illumina HiSeq 2000進行測序，得到長為100堿基對的單端序列片段。
本文所描述的分析假設從RNA-seq實驗獲得的序列片段已經與適當的參考基因組比對，并已經在基因水平上對序列進行了統(tǒng)計計數。在本文條件下，使用Rsubread包提供的基于R的流程將序列片段與小鼠參考基因組（mm10）比對（具體而言，先使用align函數(Liao, Smyth, and Shi 2013)，然后使用featureCounts (Liao, Smyth, and Shi 2014)進行基因水平的總結，利用其內置的mm10基于RefSeq的注釋）。
這些樣本的計數數據可以從Gene Expression Omnibus (GEO)數據庫 http://www.ncbi.nlm.nih.gov/geo/使用GEO序列登記號GSE63310下載。更多關于實驗設計和樣品制備的信息也可以在GEO使用該登記號查看。
library(limma)library(Glimma)library(edgeR)library(Mus.musculus)為開始此分析，從https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file在線下載文件GSE63310_RAW.tar，并從壓縮包中解壓出相關的文件。下方的代碼將完成此步驟，或者您也可以手動進行這一步并繼續(xù)后續(xù)分析。
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb") utils::untar("GSE63310_RAW.tar", exdir = ".")files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")for(i in paste(files, ".gz", sep="")) R.utils::gunzip(i, overwrite=TRUE)每一個文本文件均包含一個給定樣品的原始基因水平計數。需要注意的是，我們的分析僅包含了此實驗中的basal、LP和ML樣品（請查看下方相關文件名）。
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")read.delim(files[1], nrow=5)## EntrezID GeneLength Count## 1 497097 3634 1## 2 100503874 3259 0## 3 100038431 1634 0## 4 19888 9747 0## 5 20671 3130 1盡管這九個文本文件可以分別讀入R然后合并為一個計數矩陣，edgeR提供了更方便的途徑，使用readDGE函數即可一步完成。得到的DGEList對象中包含一個計數矩陣，它的27179行分別對應唯一的Entrez基因標識（ID），九列分別對應此實驗中的每個樣品。
x <- readDGE(files, columns=c(1,3))class(x)## [1] "DGEList"## attr(,"package")## [1] "edgeR"dim(x)## [1] 27179 9如果來自所有樣品的計數存儲在同一個文件中，數據可以首先讀入R再使用DGEList函數轉換為一個DGEList對象。
為進行下游分析，與實驗設計有關的樣品水平信息需要與計數矩陣的列關聯。這里需要包括各種對表達水平有影響的實驗變量，無論是生物變量還是技術變量。例如，細胞類型（在這個實驗中是basal、LP和ML），基因型（野生型、敲除），表型（疾病狀態(tài)、性別、年齡），樣品處理（用藥、對照）和批次信息（如果樣品是在不同時間點進行收集和分析的，記錄進行實驗的時間）等。
我們的DGEList對象中包含的samples數據框同時存儲了細胞類型（group）和批次（測序泳道lane）信息，每種信息都包含三個不同的水平。需要注意的是，在x$samples中，程序會自動計算每個樣品的文庫大小，歸一化系數會被設置為1。 為了簡單起見，我們從我們的DGEList對象x的列名中刪去了GEO樣品ID（GSM*）。
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))samplenames## [1] "10_6_5_11" "9_6_5_11" "purep53" "JMS8-2" "JMS8-3" "JMS8-4" "JMS8-5" ## [8] "JMS9-P7c" "JMS9-P8c"colnames(x) <- samplenamesgroup <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", "Basal", "ML", "LP"))x$samples$group <- grouplane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))x$samples$lane <- lanex$samples## files group lib.size norm.factors lane## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004## purep53 GSM1545538_purep53.txt Basal 57160817 1 L004## JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006## JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006## JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006## JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008我們的DGEList對象中的第二個數據框名為genes，用于存儲與計數矩陣的行相關聯的基因水平的信息。 為檢索這些信息，我們可以使用包含特定物種信息的包，比如小鼠的Mus.musculus(Bioconductor Core Team 2016b)（或人類的Homo.sapiens (Bioconductor Core Team 2016a)）；或者也可以使用biomaRt 包 (Durinck et al. 2005, 2009)，它通過接入Ensembl genome數據庫來進行基因注釋。
可以檢索的信息類型包括基因符號（gene symbols）、基因名稱（gene names）、染色體名稱和位置（chromosome names and locations）、Entrez基因ID（Entrez gene IDs）、Refseq基因ID（Refseq gene IDs）和Ensembl基因ID（Ensembl gene IDs）等。biomaRt主要使用Ensembl基因ID進行檢索，而由于Mus.musculus包含多種不同來源的信息，它允許用戶從多種不同基因ID中選取檢索鍵。
我們使用Mus.musculus包，利用我們數據集中的Entrez基因ID來檢索相關的基因符號和染色體信息。
geneid <- rownames(x)genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"), keytype="ENTREZID")head(genes)## ENTREZID SYMBOL TXCHROM## 1 497097 Xkr4 chr1## 2 100503874 Gm19938 <NA>## 3 100038431 Gm10568 <NA>## 4 19888 Rp1 chr1## 5 20671 Sox17 chr1## 6 27395 Mrpl15 chr1與任何基因ID一樣，Entrez基因ID可能不能一對一地匹配我們想獲得的基因信息。在處理之前，檢查重復的基因ID和弄清楚重復的來源非常重要。我們的基因注釋中包含28個匹配到不同染色體的基因（比如基因Gm1987關聯于染色體chr4和chr4_JH584294_random，小RNA Mir5098關聯于chr2，chr5，chr8，chr11和chr17）。 為了處理重復的基因ID，我們可以合并來自多重匹配基因的所有染色體信息，比如將基因Gm1987分配到chr4 and chr4_JH584294_random，或選取其中一條染色體來代表具有重復注釋的基因。為了簡單起見，我們選擇后者，保留每個基因ID第一次出現的信息。
genes <- genes[!duplicated(genes$ENTREZID),]在此例子中，注釋與數據對象中的基因順序是相同的。如果由于缺失和／或重新排列基因ID導致其順序不一致，可以用match來正確排序基因。然后將基因注釋的數據框加入數據對象，數據即被整潔地整理入一個DGEList對象中，它包含原始計數數據和相關的樣品信息和基因注釋。
x$genes <- genesx## An object of class "DGEList"## $samples## files group lib.size norm.factors lane## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004## purep53 GSM1545538_purep53.txt Basal 57160817 1 L004## JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006## JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006## JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006## JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008## ## $counts## Samples## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c## 497097 1 2 342 526 3 3 535 2 0## 100503874 0 0 5 6 0 0 5 0 0## 100038431 0 0 0 0 0 0 1 0 0## 19888 0 1 0 0 17 2 0 1 0## 20671 1 1 76 40 33 14 98 18 8## 27174 more rows ...## ## $genes## ENTREZID SYMBOL TXCHROM## 1 497097 Xkr4 chr1## 2 100503874 Gm19938 <NA>## 3 100038431 Gm10568 <NA>## 4 19888 Rp1 chr1## 5 20671 Sox17 chr1## 27174 more rows ...由于更深的測序總會產生更多的序列，在差異表達相關的分析中，我們很少使用原始的序列數。在實踐中，我們通常將原始的序列數進行歸一化，來消除測序深度所導致的差異。通常被使用的方法有基于序列的CPM（counts per million）、log-CPM、FPKM（fragments per kilobase of transcript per million），和基于轉錄本數目的RPKM（reads per kilobase of transcript per million）。
盡管CPM和log-CPM轉換并不像RPKM和FPKM那樣考慮到基因長度區(qū)別的影響，但在我們的分析中經常會用到它們。雖然也可以使用RPKM和FPKM，但CPM和log-CPM只使用計數矩陣即可計算，且已足以滿足我們所關注的比較的需要。假設不同條件之間剪接異構體（isoform）的使用沒有差別，差異表達分析研究同一基因在不同條件下的表達差異，而不是比較多個基因之間的表達或測定絕對表達量。換而言之，基因長度在我們關注的比較中始終不變，且任何觀測到的差異是來自于條件的變化而不是基因長度的變化。
在此處，使用edgeR中的cpm函數將原始計數轉換為CPM和log-CPM值。如果可以提供基因長度信息，可以使用edgeR中的rpkm函數計算RPKM值，就像計算CPM值那樣簡單。
cpm <- cpm(x)lcpm <- cpm(x, log=TRUE, prior.count=2)對于一個基因，CPM值為1相當于在測序深度最低的樣品中（JMS9-P8c, 序列數量約2千萬）有20個計數，或者在測序深度最高的樣品中（JMS8-3，序列數量約7.6千萬）有76個計數。
log-CPM值將被用于探索性圖表中。當設置log=TRUE時，cpm函數會在進行l(wèi)og2轉換前給CPM值加上一個彌補值。默認的彌補值是2/L，其中2是“預先計數”，而L是樣本總序列數（以百萬計）的平均值，所以log-CPM值是根據CPM值通過log2(CPM + 2/L)計算得到的。這樣的計算方式可以確保任意兩個具有相同CPM值的序列片段計數的log-CPM值也相同。彌補值的使用可以避免對零取對數，并能使所有樣本間的log倍數變化（log-fold-change）向0推移而減小低表達基因間微小計數變化帶來的巨大的偽差異性，這對于繪制探索性圖表很有用。在這個數據集中，平均的樣本總序列數是4.55千萬，所以L約等于45.5，且每個樣本中的最小log-CPM值為log2(2/45.5) = -4.51。換而言之，在加上了預先計數彌補值后，此數據集中的零表達計數對應的log-CPM值為-4.51：
L <- mean(x$samples$lib.size) * 1e-6M <- median(x$samples$lib.size) * 1e-6c(L, M)## [1] 45.5 51.4summary(lcpm)## 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 ## Min. :-4.51 Min. :-4.51 Min. :-4.51 Min. :-4.51 Min. :-4.51 ## 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 ## Median :-0.68 Median :-0.36 Median :-0.10 Median :-0.09 Median :-0.43 ## Mean : 0.17 Mean : 0.33 Mean : 0.44 Mean : 0.41 Mean : 0.32 ## 3rd Qu.: 4.29 3rd Qu.: 4.56 3rd Qu.: 4.60 3rd Qu.: 4.55 3rd Qu.: 4.58 ## Max. :14.76 Max. :13.50 Max. :12.96 Max. :12.85 Max. :12.96 ## JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c ## Min. :-4.51 Min. :-4.51 Min. :-4.51 Min. :-4.51 ## 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 1st Qu.:-4.51 ## Median :-0.41 Median :-0.07 Median :-0.17 Median :-0.33 ## Mean : 0.25 Mean : 0.40 Mean : 0.37 Mean : 0.27 ## 3rd Qu.: 4.32 3rd Qu.: 4.43 3rd Qu.: 4.60 3rd Qu.: 4.44 ## Max. :14.85 Max. :13.19 Max. :12.94 Max. :14.01在下游的線性模型分析中，使用limma的voom函數時也會用到log-CPM值，但voom會默認使用更小的預先計數重新計算自己的log-CPM值。
所有數據集中都混有表達的基因與不表達的基因。盡管我們想要檢測在一種條件中表達但再另一種條件中不表達的基因，也有一些基因在所有樣品中都不表達。實際上，這個數據集中19%的基因在所有九個樣品中的計數都是零。
table(rowSums(x$counts==0)==9)## ## FALSE TRUE ## 22026 5153對log-CPM值的分布繪制的圖表顯示每個樣本中很大一部分基因都是不表達或者表達程度相當低的，它們的log-CPM值非常小甚至是負的（下圖A部分）。
在任何樣本中都沒有足夠多的序列片段的基因應該從下游分析中過濾掉。這樣做的原因有好幾個。 從生物學的角度來看，在任何條件下的表達水平都不具有生物學意義的基因都不值得關注，因此最好忽略。 從統(tǒng)計學的角度來看，去除低表達計數基因使數據中的均值 - 方差關系可以得到更精確的估計，并且還減少了在觀察差異表達的下游分析中需要進行的統(tǒng)計檢驗的數量。
edgeR包中的filterByExpr函數提供了自動過濾基因的方法，可保留盡可能多的有足夠表達計數的基因。
keep.exprs <- filterByExpr(x, group=group)x <- x[keep.exprs,, keep.lib.sizes=FALSE]dim(x)## [1] 16624 9此函數默認選取最小的組內的樣本數量為最小樣本數，保留至少在這個數量的樣本中有10個或更多序列片段計數的基因。對基因表達量進行過濾時使用CPM值而不是表達計數，以避免對總序列數大的樣本的偏向性。在這個數據集中，總序列數的中位數是5.1千萬，且10/51約等于0.2，所以filterByExpr函數保留在至少三個樣本中CPM值大于等于0.2的基因。此處，一個具有生物學意義的基因需要在至少三個樣本中表達，因為三種細胞類型組內各有三個重復。所使用的閾值取決于測序深度和實驗設計。如果樣本總表達計數量增大，那么可以選擇更低的CPM閾值，因為更大的總表達計數量提供了更好的分辨率來探究更多表達水平更低的基因。
使用這個標準，基因的數量減少到了16624個，約為開始時數量的60%。過濾后的log-CPM值顯示出每個樣本的分布基本相同（下圖B部分）。需要注意的是，從整個DGEList對象中取子集時同時刪除了被過濾的基因的計數和其相關的基因信息。過濾后的DGEList對象為留下的基因保留了相對應的基因信息和計數。
下方給出的是繪圖所用代碼。
lcpm.cutoff <- log2(10/M + 2/L)library(RColorBrewer)nsamples <- ncol(x)col <- brewer.pal(nsamples, "Paired")par(mfrow=c(1,2))plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")title(main="A. Raw data", xlab="Log-cpm")abline(v=lcpm.cutoff, lty=3)for (i in 2:nsamples){den <- density(lcpm[,i])lines(den$x, den$y, col=col[i], lwd=2)}legend("topright", samplenames, text.col=col, bty="n")lcpm <- cpm(x, log=TRUE)plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.26), las=2, main="", xlab="")title(main="B. Filtered data", xlab="Log-cpm")abline(v=lcpm.cutoff, lty=3)for (i in 2:nsamples){den <- density(lcpm[,i])lines(den$x, den$y, col=col[i], lwd=2)}legend("topright", samplenames, text.col=col, bty="n")Figure 1: 每個樣本過濾前的原始數據（A）和過濾后（B）的數據的log-CPM值密度。豎直虛線標出了過濾步驟中所用閾值（相當于CPM值為約0
2）。
在樣品制備或測序過程中，不具備生物學意義的外部因素會影響單個樣品的表達。比如說，在實驗中第一批制備的樣品會總體上表達高于第二批制備的樣品。假設所有樣品表達值的范圍和分布都應當相似，需要進行歸一化來確保整個實驗中每個樣本的表達分布都相似。
密度圖和箱線圖等展示每個樣品基因表達量分布的圖表可以用于判斷是否有樣品和其他樣品分布有差異。在此數據集中，所有樣品的log-CPM分布都很相似（上圖B部分）。
盡管如此，我們依然需要使用edgeR中的calcNormFactors函數，用TMM(Robinson and Oshlack 2010)方法進行歸一化。此處計算得到的歸一化系數被用作文庫大小的縮放系數。當我們使用DGEList對象時，這些歸一化系數被自動存儲在x$samples$norm.factors。對此數據集而言，TMM歸一化的作用比較溫和，這體現在所有的縮放因子都相對接近1。
x <- calcNormFactors(x, method = "TMM")x$samples$norm.factors## [1] 0.894 1.025 1.046 1.046 1.016 0.922 0.996 1.086 0.984為了更好地可視化表現出歸一化的影響，我們復制了數據并進行了調整，使得第一個樣品的計數減少到了其原始值的5%，而第二個樣品增大到了5倍。
x2 <- xx2$samples$norm.factors <- 1x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)x2$counts[,2] <- x2$counts[,2]*5下圖顯示了沒有經過歸一化的與經過了歸一化的數據的樣本的表達分布，其中歸一化前的分布顯然不同，而歸一化后比較相似。此處，第一個樣品的TMM縮放系數0.06非常小，而第二個樣品的縮放系數6.08很大，它們都并不接近1。
par(mfrow=c(1,2))lcpm <- cpm(x2, log=TRUE)boxplot(lcpm, las=2, col=col, main="")title(main="A. Example: Unnormalised data",ylab="Log-cpm")x2 <- calcNormFactors(x2) x2$samples$norm.factors## [1] 0.0577 6.0829 1.2202 1.1648 1.1966 1.0466 1.1505 1.2543 1.1090lcpm <- cpm(x2, log=TRUE)boxplot(lcpm, las=2, col=col, main="")title(main="B. Example: Normalised data",ylab="Log-cpm")Figure 2: 樣例數據：log-CPM值的箱線圖展示了未經歸一化的數據（A）及歸一化后的數據（B）中每個樣本的表達分布。數據集經過調整，樣本1和2中的表達計數分別被縮放到其原始值的5%和500%。
在我們看來，用于檢查基因表達分析的最重要的探索性圖表之一便是MDS圖或其余類似的圖。這種圖表使用無監(jiān)督聚類方法展示出了樣品間的相似性和不相似性，能讓我們在進行正式的檢驗之前對于能檢測到多少差異表達基因有個大致概念。理想情況下，樣本會在不同的實驗組內很好的聚類，且可以鑒別出遠離所屬組的樣本，并追蹤誤差或額外方差的來源。如果存在技術重復，它們應當互相非常接近。
這樣的圖可以用limma中的plotMDS函數繪制。第一個維度表示能夠最好地分離樣品且解釋最大比例的方差的引導性的倍數變化（leading-fold-change），而后續(xù)的維度的影響更小，并與之前的維度正交。當實驗設計涉及到多個因子時，建議在多個維度上檢查每個因子。如果在其中一些維度上樣本可按照某因子聚類，這說明該因子對于表達差異有影響，在線性模型中應當將其包括進去。反之，沒有或者僅有微小影響的因子在下游分析時應當被剔除。
在這個數據集中，可以看出樣本在維度1和2能很好地按照實驗分組聚類，隨后在維度3按照測序道（樣品批次）分離（如下圖所示）。請記住，第一維度解釋了數據中最大比例的方差，需要注意到，當我們向高維度移動，維度上的取值范圍會變小。
盡管所有樣本都按組聚類，在維度1上最大的轉錄差異出現在basal和LP以及basal和ML之間。因此，預期在basal樣品與其他之間的成對比較中能夠得到大量的DE基因，而在ML和LP之間的比較中得到的DE基因數量略少。在其他的數據集中，不按照實驗組聚類的樣本可能在下游分析中只表現出較小的或不表現出差異表達。
為繪制MDS圖，我們?yōu)椴煌囊蜃淤x予不同的色彩組合。維度1和2使用以細胞類型定義的色彩組合進行檢查。
維度3和4使用以測序泳道（批次）定義的色彩組合進行檢查。
lcpm <- cpm(x, log=TRUE)par(mfrow=c(1,2))col.group <- grouplevels(col.group) <- brewer.pal(nlevels(col.group), "Set1")col.group <- as.character(col.group)col.lane <- lanelevels(col.lane) <- brewer.pal(nlevels(col.lane), "Set2")col.lane <- as.character(col.lane)plotMDS(lcpm, labels=group, col=col.group)title(main="A. Sample groups")plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))title(main="B. Sequencing lanes")Figure 3: log-CPM值在維度1和2的MDS圖，以樣品分組上色并標記（A）和維度3和4的MDS圖，以測序道上色并標記（B）。圖中的距離對應于最主要的倍數變化（fold change），默認情況下也就是前500個在每對樣品之間差異最大的基因的平均（均方根）log2倍數變化。
作為另一種選擇，Glimma包也提供了便于探索多個維度的交互式MDS圖。其中的glMDSPlot函數可生成一個html網頁（如果設置launch=TRUE，將會在瀏覽器中打開），其左側面板含有一張MDS圖，而右側面板包含一張展示了各個維度所解釋的方差比例的柱形圖。點擊柱形圖中的柱可切換MDS圖繪制時所使用的維度，且將鼠標懸浮于單個點上可顯示相應的樣本標簽。也可切換配色方案，以突顯不同細胞類型或測序泳道（批次）。此數據集的交互式MDS圖可以從http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MDS-Plot.html看到。
glMDSPlot(lcpm, labels=paste(group, lane, sep="_"), groups=x$samples[,c(2,5)], launch=FALSE)在此研究中，我們想知道哪些基因在我們研究的三組細胞之間以不同水平表達。在我們的分析中，假設基礎數據是正態(tài)分布的，為其擬合一個線性模型。在此之前，需要創(chuàng)建一個包含細胞類型以及測序泳道（批次）信息的設計矩陣。
design <- model.matrix(~0+group+lane)colnames(design) <- gsub("group", "", colnames(design))design## Basal LP ML laneL006 laneL008## 1 0 1 0 0 0## 2 0 0 1 0 0## 3 1 0 0 0 0## 4 1 0 0 1 0## 5 0 0 1 1 0## 6 0 1 0 1 0## 7 1 0 0 1 0## 8 0 0 1 0 1## 9 0 1 0 0 1## attr(,"assign")## [1] 1 1 1 2 2## attr(,"contrasts")## attr(,"contrasts")$group## [1] "contr.treatment"## ## attr(,"contrasts")$lane## [1] "contr.treatment"對于一個給定的實驗，通常有幾種等價的方法可以創(chuàng)建一個合適的設計矩陣。 比如說，~0+group+lane去除了第一個因子group的截距，但第二個因子lane的截距被保留。 此外也可以使用~group+lane，來自group和lane的截距均被保留。 此處的關鍵是理解如何解釋給定模型中估計得到的系數。 我們在此分析中選取第一種模型，因為在沒有group的截距的情況下能更直截了當地設定模型中的對比。用于細胞群之間成對比較的對比可以在limma中用makeContrasts函數設定。
contr.matrix <- makeContrasts( BasalvsLP = Basal-LP, BasalvsML = Basal - ML, LPvsML = LP - ML, levels = colnames(design))contr.matrix## Contrasts## Levels BasalvsLP BasalvsML LPvsML## Basal 1 1 0## LP -1 0 1## ML 0 -1 -1## laneL006 0 0 0## laneL008 0 0 0limma的線性模型方法的核心優(yōu)勢之一便是其適應任意實驗復雜程度的能力。簡單的設計，比如此工作流程中關于細胞類型和批次的實驗設計，直到更復雜的因子設計和含有交互作用項的模型，都能夠被相對簡單地處理。當實驗或技術效應可被隨機效應模型（random effect model）模擬時，limma中的另一種可能性是使用duplicateCorrelation函數來估計交互作用，這需要在此函數以及lmFit的線性建模步驟均指定一個block參數。
據顯示對于RNA-seq計數數據而言，當使用原始計數或當其被轉換為log-CPM值時，方差并不獨立于均值(Law et al. 2014)。使用負二項分布來模擬計數的方法假設均值與方差間具有二次的關系。在limma中，假設log-CPM值符合正態(tài)分布，并使用由voom函數計算得到的精確權重來調整均值與方差的關系，從而對log-CPM值進行線性建模。
當操作DGEList對象時，voom從x中自動提取文庫大小和歸一化因子，以此將原始計數轉換為log-CPM值。在voom中，對于log-CPM值額外的歸一化可以通過設定normalize.method參數來進行。
下圖左側展示了這個數據集log-CPM值的均值-方差關系。通常而言，方差是測序實驗中的技術差異和不同細胞類型的重復樣本之間的生物學差異的結合，而voom圖會顯示出一個在均值與方差之間遞減的趨勢。 生物學差異高的實驗通常會有更平坦的趨勢，其方差值在高表達處穩(wěn)定。 生物學差異低的實驗更傾向于急劇下降的趨勢。
不僅如此，voom圖也提供了對于上游所進行的過濾水平的可視化檢測。如果對于低表達基因的過濾不夠充分，在圖上表達低的一端，受到非常低的表達計數的影響，可以觀察到方差水平的下降。如果觀察到了這種情況，應當回到最初的過濾步驟并提高用于該數據集的表達閾值。
當前面觀察的MDS圖中具有明顯的樣本水平的差異時，可以用voomWithQualityWeights函數來同時合并樣本水平的權重和voom(Liu et al. 2015)估算得到的豐度相關的權重。關于此種方式的例子參見Liu等(2016) (Liu et al. 2016)。
par(mfrow=c(1,2))v <- voom(x, design, plot=TRUE)v## An object of class "EList"## $genes## ENTREZID SYMBOL TXCHROM## 1 497097 Xkr4 chr1## 5 20671 Sox17 chr1## 6 27395 Mrpl15 chr1## 7 18777 Lypla1 chr1## 9 21399 Tcea1 chr1## 16619 more rows ...## ## $targets## files group lib.size norm.factors lane## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 29387429 0.894 L004## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 36212498 1.025 L004## purep53 GSM1545538_purep53.txt Basal 59771061 1.046 L004## JMS8-2 GSM1545539_JMS8-2.txt Basal 53711278 1.046 L006## JMS8-3 GSM1545540_JMS8-3.txt ML 77015912 1.016 L006## JMS8-4 GSM1545541_JMS8-4.txt LP 55769890 0.922 L006## JMS8-5 GSM1545542_JMS8-5.txt Basal 54863512 0.996 L006## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 23139691 1.086 L008## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19634459 0.984 L008## ## $E## Samples## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c## 497097 -4.29 -3.86 2.519 3.293 -4.46 -3.99 3.287 -3.210 -5.30## 20671 -4.29 -4.59 0.356 -0.407 -1.20 -1.94 0.844 -0.323 -1.21## 27395 3.88 4.41 4.517 4.562 4.34 3.79 3.899 4.340 4.12## 18777 4.71 5.57 5.396 5.162 5.65 5.08 5.060 5.751 5.14## 21399 4.79 4.75 5.370 5.122 4.87 4.94 5.138 5.031 4.98## 16619 more rows ...## ## $weights## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]## [1,] 1.08 1.33 19.8 20.27 1.99 1.40 20.49 1.11 1.08## [2,] 1.17 1.46 4.8 8.66 3.61 2.63 8.76 3.21 2.54## [3,] 20.22 25.57 30.4 28.53 31.35 25.74 28.72 21.20 16.66## [4,] 26.95 32.51 33.6 33.23 34.23 32.35 33.33 30.35 24.26## [5,] 26.61 28.50 33.6 33.21 33.57 32.00 33.31 25.17 23.57## 16619 more rows ...## ## $design## Basal LP ML laneL006 laneL008## 1 0 1 0 0 0## 2 0 0 1 0 0## 3 1 0 0 0 0## 4 1 0 0 1 0## 5 0 0 1 1 0## 6 0 1 0 1 0## 7 1 0 0 1 0## 8 0 0 1 0 1## 9 0 1 0 0 1## attr(,"assign")## [1] 1 1 1 2 2## attr(,"contrasts")## attr(,"contrasts")$group## [1] "contr.treatment"## ## attr(,"contrasts")$lane## [1] "contr.treatment"vfit <- lmFit(v, design)vfit <- contrasts.fit(vfit, contrasts=contr.matrix)efit <- eBayes(vfit)plotSA(efit, main="Final model: Mean-variance trend")Figure 4: 圖中繪制了每個基因的均值（x軸）和方差（y軸），顯示了在該數據上使用voom前它們之間的相關性（左），以及當運用voom的精確權重后這種趨勢是如何消除的（右）。左側的圖是使用voom函數繪制的，它為進行l(wèi)og-CPM轉換后的數據擬合線性模型從而提取殘差方差。然后，對方差取平方根（或對標準差取平方根），并相對每個基因的平均表達作圖。均值通過平均計數加上2再進行l(wèi)og2轉換計算得到。右側的圖使用plotSA繪制了log2殘差標準差與log-CPM均值的關系。平均log2殘差標準差由水平藍線標出。在這兩幅圖中，每個黑點表示一個基因，紅線為對這些點的擬合。
值得注意的是，DGEList對象中存儲的另一個數據框，即基因和樣本水平信息所存儲之處，保留在了voom創(chuàng)建的EList對象v中。v$genes數據框等同于x$genes，v$targets等同于x$samples，而v$E中所儲存的表達值類似于x$counts，盡管它進行了尺度轉換。此外，voom的EList對象中還有一個精確權重的矩陣v$weights，而設計矩陣存儲于v$design。
limma的線性建模使用lmFit和contrasts.fit函數進行，它們原先是為微陣列而設計的。這些函數不僅可以用于微陣列數據，也可以用于RNA-seq數據。它們會單獨為每個基因的表達值擬合一個模型。然后，通過利用所有基因的信息來進行經驗貝葉斯調整，這樣可以獲得更精確的基因水平的變異程度估計(Smyth 2004)。下一圖為此模型的殘差關于平均表達值的圖。從圖中可以看出，方差不再與表達水平均值相關。
為快速查看差異表達水平，顯著上調或下調的基因可以匯總到一個表格中。顯著性的判斷使用校正p值閾值的默認值5%。在basal與LP的表達水平之間的比較中，發(fā)現了4648個在basal中相較于LP下調的基因和4863個在basal中相較于LP上調的基因，即共9511個DE基因。在basal和ML之間發(fā)現了一共9598個DE基因（4927個下調基因和4671個上調基因），而在LP和ML中發(fā)現了一共5652個DE基因（3135個下調基因和2517個上調基因）。在包括basal細胞類型的比較中皆找到了大量的DE基因，這與我們在MDS圖中觀察到的結果相吻合。
summary(decideTests(efit))## BasalvsLP BasalvsML LPvsML## Down 4648 4927 3135## NotSig 7113 7026 10972## Up 4863 4671 2517一些研究中不僅僅需要使用校正p值閾值，更為嚴格定義的顯著性可能需要差異倍數的對數（log-FCs）也高于某個最小值。treat方法(McCarthy and Smyth 2009)可以按照對最小log-FC值的要求，使用經過經驗貝葉斯調整的t統(tǒng)計值計算p值。當我們的檢驗要求基因的log-FC顯著大于1（等同于在原本的尺度上不同細胞類型之間差兩倍）時，差異表達基因的數量得到了下降，basal與LP相比只有3684個DE基因，basal與ML相比只有3834個DE基因，LP與ML相比只有414個DE基因。
tfit <- treat(vfit, lfc=1)dt <- decideTests(tfit)summary(dt)## BasalvsLP BasalvsML LPvsML## Down 1632 1777 224## NotSig 12976 12790 16210## Up 2016 2057 190在多種比較中皆差異表達的基因可以從decideTests的結果中提取，其中的0代表不差異表達的基因，1代表上調的基因，-1代表下調的基因。共有2784個基因在basal和LP以及basal和ML的比較中都差異表達，其中的20個于下方列出。write.fit函數可用于將三個比較的結果提取并寫入到單個輸出文件。
de.common <- which(dt[,1]!=0 & dt[,2]!=0)length(de.common)## [1] 2784head(tfit$genes$SYMBOL[de.common], n=20)## [1] "Xkr4" "Rgs20" "Cpa6" "A830018L16Rik" "Sulf1" ## [6] "Eya1" "Msc" "Sbspon" "Pi15" "Crispld1" ## [11] "Kcnq5" "Rims1" "Khdrbs2" "Ptpn18" "Prss39" ## [16] "Arhgef4" "Cnga3" "2010300C02Rik" "Aff3" "Npas2"vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))Figure 5: 韋恩圖展示了僅basal和LP（左）、僅basal和ML（右）的對比的DE基因數量，還有兩種對比中共同的DE基因數量（中）。在任何對比中均不差異表達的基因數量標于右下。
write.fit(tfit, dt, file="results.txt")使用topTreat函數可以列舉出使用treat得到的結果中靠前的DE基因（對于eBayes的結果可以使用topTable函數）。默認情況下，topTreat將基因按照校正p值從小到大排列，并為每個基因給出相關的基因信息、log-FC、平均log-CPM、校正t值、原始及經過多重假設檢驗校正的p值。列出前多少個基因的數量可由用戶指定，如果設為n=Inf則會包括所有的基因。基因Cldn7和Rasef在basal與LP和basal于ML的比較中都位于DE基因的前幾名。
basal.vs.lp <- topTreat(tfit, coef=1, n=Inf)basal.vs.ml <- topTreat(tfit, coef=2, n=Inf)head(basal.vs.lp)## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val## 12759 12759 Clu chr14 -5.46 8.86 -33.6 1.72e-10 1.71e-06## 53624 53624 Cldn7 chr11 -5.53 6.30 -32.0 2.58e-10 1.71e-06## 242505 242505 Rasef chr4 -5.94 5.12 -31.3 3.08e-10 1.71e-06## 67451 67451 Pkp2 chr16 -5.74 4.42 -29.9 4.58e-10 1.74e-06## 228543 228543 Rhov chr2 -6.26 5.49 -29.1 5.78e-10 1.74e-06## 70350 70350 Basp1 chr15 -6.08 5.25 -28.3 7.27e-10 1.74e-06head(basal.vs.ml)## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val## 242505 242505 Rasef chr4 -6.53 5.12 -35.1 1.23e-10 1.24e-06## 53624 53624 Cldn7 chr11 -5.50 6.30 -31.7 2.77e-10 1.24e-06## 12521 12521 Cd82 chr2 -4.69 7.07 -31.4 2.91e-10 1.24e-06## 20661 20661 Sort1 chr3 -4.93 6.70 -30.7 3.56e-10 1.24e-06## 71740 71740 Nectin4 chr1 -5.58 5.16 -30.6 3.72e-10 1.24e-06## 12759 12759 Clu chr14 -4.69 8.86 -28.0 7.69e-10 1.48e-06為可視化地總結所有基因的結果，可使用plotMD函數繪制均值-差異（MD）圖，其中展示了線性模型擬合所得到的log-FC與log-CPM平均值間的關系，而差異表達的基因會被重點標出。
plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1], xlim=c(-8,13))Glimma的glMDPlot函數提供了交互式的均值-差異圖，拓展了這種圖表的功能性。 此函數的輸出為一個html頁面，左側面板為結果的總結性圖表（與plotMD的輸出類似），右側面板包含各個樣本的log-CPM值，下方為結果的表格。 這種交互式展示允許用戶使用提供的注釋（比如基因名標識）搜索特定基因，而這在R統(tǒng)計圖中是做不到的。
glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1], side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)使用Glimma生成的均值-差異圖。左側面板顯示了總結性數據（log-FC與log-CPM值的關系），其中選中的基因在每個樣本中的數值顯示于右側面板。下方為結果的表格，其搜索框使用戶得以使用可行的注釋信息查找某個特定基因，如基因符號Clu。
上方指令生成的均值-差異圖可以在線訪問（詳見http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MD-Plot.html）。 Glimma提供的交互性使得單個圖形窗口內能夠呈現出額外的信息。 Glimma是以R和Javascript實現的，使用R代碼生成數據，并在之后使用Javascript庫D3（https://d3js.org）轉換為圖形，使用Bootstrap庫處理界面并生成互動性可搜索的表格的數據表。 這使得圖表可以在任何現代的瀏覽器中查看，對于從Rmarkdown分析報告中將其作為關聯文件而附加而言十分方便。
前文所展示的圖表中，一些展示了在任意一個條件下表達的所有基因（比如共同DE基因的韋恩圖或均值-差異圖），而另一些展示單獨的基因（交互性均值-差異圖右邊面板中所展示的log-CPM值）。而熱圖使用戶得以查看一部分基因的表達。這對于查看單個組或樣本的表達很有用，而不至于在關注于單個基因時失去對于研究整體的注意，也不會造成由于上千個基因所取平均值而導致的失去分辨率。
使用gplots包的heatmap.2函數，我們?yōu)閎asal與LP的對照中前100個DE基因（按調整p值排序）繪制了一幅熱圖。熱圖中正確地將樣本按照細胞類型聚類，并重新排序了基因，形成了表達相似的塊狀。從熱圖中，我們觀察到對于basal與LP之間的前100個DE基因，ML和LP樣本的表達非常相似。
library(gplots)basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100]i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes)mycol <- colorpanel(1000,"blue","white","red")heatmap.2(lcpm[i,], scale="row", labRow=v$genes$SYMBOL[i], labCol=group, col=mycol, trace="none", density.info="none", margin=c(8,6), lhei=c(2,10), dendrogram="column")Figure 6: 在basal和LP的對比中前100個DE基因log-CPM值的熱圖。經過縮放調整后，每個基因（每行）的表達均值為0，并且標準差為1。給定基因相對高表達的樣本被標記為紅色，相對低表達的樣本被標記為藍色。淺色和白色代表中等表達水平的基因。樣本和基因已通過分層聚類的方法重新排序。圖中顯示有樣本聚類的樹狀圖。
在此次分析的最后，我們要進行一些基因集檢驗。為此，我們將camera方法(Wu and Smyth 2012)應用于Broad Institute的MSigDB c2中的(Subramanian et al. 2005)中適應小鼠的c2基因表達特征，這可從http://bioinf.wehi.edu.au/software/MSigDB/以RData對象格式獲取。 此外，對于人類和小鼠，來自MSigDB的其他有用的基因集也可從此網站獲取，比如標志（hallmark）基因集。C2基因集的內容收集自在線數據庫、出版物以及該領域專家，而標志基因集的內容來自MSigDB，從而獲得具有明確定義的生物狀態(tài)或過程。
load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123"))idx <- ids2indices(Mm.c2,id=rownames(v))cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1])head(cam.BasalvsLP,5)## NGenes Direction PValue FDR## LIM_MAMMARY_STEM_CELL_UP 791 Up 1.77e-18 8.36e-15## LIM_MAMMARY_STEM_CELL_DN 683 Down 4.03e-14 8.69e-11## ROSTY_CERVICAL_CANCER_PROLIFERATION_CLUSTER 170 Up 5.52e-14 8.69e-11## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 94 Down 2.74e-13 3.23e-10## SOTIRIOU_BREAST_CANCER_GRADE_1_VS_3_UP 190 Up 5.16e-13 4.87e-10cam.BasalvsML <- camera(v,idx,design,contrast=contr.matrix[,2])head(cam.BasalvsML,5)## NGenes Direction PValue FDR## LIM_MAMMARY_STEM_CELL_UP 791 Up 1.68e-22 7.92e-19## LIM_MAMMARY_STEM_CELL_DN 683 Down 7.79e-18 1.84e-14## LIM_MAMMARY_LUMINAL_MATURE_DN 172 Up 9.74e-16 1.53e-12## LIM_MAMMARY_LUMINAL_MATURE_UP 204 Down 1.15e-12 1.36e-09## NAKAYAMA_SOFT_TISSUE_TUMORS_PCA2_UP 137 Up 2.24e-12 1.88e-09cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3])head(cam.LPvsML,5)## NGenes Direction PValue FDR## LIM_MAMMARY_LUMINAL_MATURE_DN 172 Up 6.73e-14 2.35e-10## LIM_MAMMARY_LUMINAL_MATURE_UP 204 Down 9.97e-14 2.35e-10## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 94 Up 1.32e-11 2.08e-08## REACTOME_RESPIRATORY_ELECTRON_TRANSPORT 94 Down 7.01e-09 8.28e-06## REACTOME_RNA_POL_I_PROMOTER_OPENING 46 Down 2.04e-08 1.93e-05camera函數通過比較假設檢驗來評估一個給定基因集中的基因是否相對于不在集內的基因而言在差異表達基因的排序中更靠前。 它使用limma的線性模型框架，并同時采用設計矩陣和對比矩陣（如果有的話），且在測試的過程中會使用來自voom的觀測水平權重。 在通過基因間相關性（默認設定為0.01，但也可通過數據估計）和基因集的規(guī)模得到方差膨脹因子（variance inflation factor），并使用它調整基因集檢驗統(tǒng)計值的方差后，將會返回根據多重假設檢驗進行了校正的p值。
此實驗是與Lim等人(2010)(Lim et al. 2010)的數據集等價的RNA-seq，而他們使用Illumina微陣列分析了相同的分選細胞群，因此該早期文獻中的基因表達特征出現在每種對比的列表頂部正符合我們的預期。在LP和ML的對比中，我們?yōu)長im等人（2010）的成熟管腔基因集（上調及下調）繪制了條碼圖（barcodeplot）。需要注意的是，由于我們的對比是將LP與ML相比而不是相反，這些基因集的方向在我們的數據集中是反過來的（如果將對比反過來，基因集的方向將會與對比一致）。
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP, index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")Figure 7: LIM_MAMMARY_LUMINAL_MATURE_UP （紅色條形，圖表上方）和LIM_MAMMARY_LUMINAL_MATURE_DN（藍色條形，圖表下方）基因集在LP和ML的對比中的條碼圖，每個基因集都有一條富集線展示了豎直條形在圖表每部分的相對富集程度。Lim等人的實驗(Lim et al
2010)非常類似于我們的，用了相同的分選方式來獲取不同的細胞群，只是他們使用的是微陣列而不是RNA-seq來測定基因表達。需要注意的是，上調基因集發(fā)生下調而下調基因集發(fā)生上調的逆相關性來自于對比的設定方式（LP相比于ML），如果將其對調，方向性將會吻合。
limma也有其他的基因集檢驗，比如mroast(Wu et al. 2010)的自包含檢驗。雖然camera非常適合檢驗基因集的大型數據庫并觀察其中哪些相對于其他的在排序上位次更高（如前文所示），自包含檢驗更善于集中檢驗一個或少個選中的集合是否本身差異表達。換句話說，camera更適用于搜尋具有意義的基因集，而mroast測試的是已經確定有意義的基因集的顯著性。
此RNA-seq工作流程使用了Bioconductor項目3.8版本中的多個軟件包，運行于R 3.5.1或更高版本。除了本文中重點提到的軟件（limma、Glimma以及edgeR），亦需要一些其他軟件包，包括gplots和RColorBrewer還有基因注釋包Mus.musculus。 此文檔使用knitr編譯。所有用到的包的版本號如下所示。 Bioconductor工作流程包RNAseq123（可訪問https://bioconductor.org/packages/RNAseq123查看）內包含此文章的英文和簡體中文版以及進行整個分析流程所需要的代碼。安裝此包即可管理以上提到的所有需要的包。對于RNA-seq數據分析實踐培訓而言，此包也是非常有用的資源。
sessionInfo()## R version 3.6.0 (2019-04-26)## Platform: x86_64-pc-linux-gnu (64-bit)## Running under: Ubuntu 18.04.2 LTS## ## Matrix products: default## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so## ## locale:## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 ## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C ## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages:## [1] parallel stats4 stats graphics grDevices utils datasets methods ## [9] base ## ## other attached packages:## [1] knitr_1.22 gplots_3.0.1.1 ## [3] RColorBrewer_1.1-2 Mus.musculus_1.3.1 ## [5] TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.7 org.Mm.eg.db_3.8.2 ## [7] GO.db_3.8.2 OrganismDbi_1.26.0 ## [9] GenomicFeatures_1.36.0 GenomicRanges_1.36.0 ## [11] GenomeInfoDb_1.20.0 AnnotationDbi_1.46.0 ## [13] IRanges_2.18.0 S4Vectors_0.22.0 ## [15] Biobase_2.44.0 BiocGenerics_0.30.0 ## [17] edgeR_3.26.0 Glimma_1.12.0 ## [19] limma_3.40.0 BiocStyle_2.12.0 ## ## loaded via a namespace (and not attached):## [1] httr_1.4.0 bit64_0.9-7 jsonlite_1.6 ## [4] R.utils_2.8.0 gtools_3.8.1 assertthat_0.2.1 ## [7] BiocManager_1.30.4 highr_0.8 RBGL_1.60.0 ## [10] blob_1.1.1 GenomeInfoDbData_1.2.1 Rsamtools_2.0.0 ## [13] yaml_2.2.0 progress_1.2.0 RSQLite_2.1.1 ## [16] lattice_0.20-38 digest_0.6.18 XVector_0.24.0 ## [19] htmltools_0.3.6 Matrix_1.2-17 R.oo_1.22.0 ## [22] XML_3.98-1.19 pkgconfig_2.0.2 biomaRt_2.40.0 ## [25] bookdown_0.9 zlibbioc_1.30.0 gdata_2.18.0 ## [28] BiocParallel_1.18.0 SummarizedExperiment_1.14.0 magrittr_1.5 ## [31] crayon_1.3.4 memoise_1.1.0 evaluate_0.13 ## [34] R.methodsS3_1.7.1 graph_1.62.0 tools_3.6.0 ## [37] prettyunits_1.0.2 hms_0.4.2 matrixStats_0.54.0 ## [40] stringr_1.4.0 locfit_1.5-9.1 DelayedArray_0.10.0 ## [43] Biostrings_2.52.0 compiler_3.6.0 caTools_1.17.1.2 ## [46] rlang_0.3.4 grid_3.6.0 RCurl_1.95-4.12 ## [49] bitops_1.0-6 rmarkdown_1.12 DBI_1.0.0 ## [52] R6_2.4.0 GenomicAlignments_1.20.0 rtracklayer_1.44.0 ## [55] bit_1.1-14 KernSmooth_2.23-15 stringi_1.4.3 ## [58] Rcpp_1.0.1 xfun_0.6Bioconductor Core Team. 2016a. Homo.sapiens: Annotation package for the Homo.sapiens object. https://bioconductor.org/packages/release/data/annotation/html/Homo.sapiens.html.
———. 2016b. Mus.musculus: Annotation package for the Mus.musculus object. https://bioconductor.org/packages/release/data/annotation/html/Mus.musculus.html.
Durinck, S., Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber. 2005. “BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis.” Bioinformatics 21:3439–40.
Durinck, S., P. Spellman, E. Birney, and W. Huber. 2009. “Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.” Nature Protocols 4:1184–91.
Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. “Orchestrating High-Throughput Genomic Analysis with Bioconductor.” Nature Methods 12 (2):115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-seq Read Counts.” Genome Biology 15:R29.
Liao, Y., G. K. Smyth, and W. Shi. 2013. “The Subread Aligner: Fast, Accurate and Scalable Read Mapping by Seed-and-Vote.” Nucleic Acids Res 41 (10):e108.
———. 2014. “featureCounts: an Efficient General-Purpose Program for Assigning Sequence Reads to Genomic Features.” Bioinformatics 30 (7):923–30.
Lim, E., D. Wu, B. Pal, T. Bouras, M. L. Asselin-Labat, F. Vaillant, H. Yagita, G. J. Lindeman, G. K. Smyth, and J. E. Visvader. 2010. “Transcriptome analyses of mouse and human mammary cell subpopulations reveal multiple conserved genes and pathways.” Breast Cancer Research 12 (2):R21.
Liu, R., K. Chen, N. Jansz, M. E. Blewitt, and M. E. Ritchie. 2016. “Transcriptional Profiling of the Epigenetic Regulator Smchd1.” Genomics Data 7:144–7.
Liu, R., A. Z. Holik, S. Su, N. Jansz, K. Chen, H. S. Leong, M. E. Blewitt, M. L. Asselin-Labat, G. K. Smyth, and M. E. Ritchie. 2015. “Why weight? Combining voom with estimates of sample quality improves power in RNA-seq analyses.” Nucleic Acids Res 43:e97.
McCarthy, D. J., and G. K. Smyth. 2009. “Testing significance relative to a fold-change threshold is a TREAT.” Bioinformatics 25:765–71.
Ritchie, M. E., B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. 2015. “l(fā)imma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.” Nucleic Acids Res 43 (7):e47.
Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics26:139–40.
Robinson, M. D., and A. Oshlack. 2010. “A Scaling Normalization Method for Differential Expression Analysis of RNA-seq data.” Genome Biology 11:R25.
Sheridan, J. M., M. E. Ritchie, S. A. Best, K. Jiang, T. J. Beck, F. Vaillant, K. Liu, et al. 2015. “A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1.” BMC Cancer 15 (1). BioMed Central:221.
Smyth, G. K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Stat Appl Genet Mol Biol 3 (1):Article 3.
Su, S., C. W. Law, C. Ah-Cann, M. L. Asselin-Labat, M. E. Blewitt, and M. E. Ritchie. 2017. “Glimma: Interactive Graphics for Gene Expression Analysis.” Bioinformatics 33:2050–52.
Subramanian, A., P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc Natl Acad Sci U S A 102 (43):15545–50.
Wu, D., E. Lim, F. Vaillant, M. L. Asselin-Labat, J. E. Visvader, and G. K. Smyth. 2010. “ROAST: rotation gene set tests for complex microarray experiments.” Bioinformatics 26 (17):2176–82.
Wu, D., and G. K. Smyth. 2012. “Camera: a competitive gene set test accounting for inter-gene correlation.” Nucleic Acids Res 40 (17):e133.
聯系客服
