原文來自這裡[1]，略編輯。
這方面的工作做的少，所以懂得真不多。
簡單且高效地分析RNA測序數據的能力是Bioconductor的核心優勢。RNA-seq分析通常從基因水平的序列計數開始，涉及到數據預處理，探索性數據分析，差異表達檢驗以及通路分析，得到的結果可用於指導進一步實驗和驗證研究。在這篇工作流程文章中，我們通過分析來自小鼠乳腺的RNA測序數據，示範了如何使用流行的edgeR包載入、整理、過濾和歸一化數據，然後用limma包的voom方法、線性模型和經驗貝葉斯調節（empirical Bayes moderation）來評估差異表達並進行基因集檢驗。通過使用Glimma包，此流程得到了增進，實現了結果的互動探索，使用戶得以查看單個樣本與基因。這三個軟體包提供的完整分析突出了研究人員可以使用Bioconductor輕鬆地從RNA測序實驗的原始計數揭示生物學意義。
2背景介紹RNA測序（RNA-seq）已成為基因表達研究的主要技術。其中，基因組規模的多條件基因差異表達檢測是研究者最常探究的問題之一。對於RNA-seq數據，來自Bioconductor項目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用於處理此問題的完善的統計學方法。
在這篇文章中，我們描述了一個用於分析RNA-seq數據的edgeR - limma工作流程，使用基因水平計數作為輸入，經過預處理和探索性數據處理，然後得到差異表達（DE）基因和基因表達特徵（gene signatures）的列表。Glimma包(Su et al. 2017)提供的交互式圖表可以同時呈現整體樣本和單個基因水平的數據，使得我們相對靜態的R圖表而言，可以探索更多的細節。
此工作流程中所分析的實驗來自Sheridan等（2015）(Sheridan et al. 2015)，由三個細胞群組成，即基底（basal）、管腔祖細胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。細胞群皆分選自雌性處女小鼠的乳腺，每種都設三個生物學重複。RNA樣品分三個批次使用Illumina HiSeq 2000進行測序，得到長為100鹼基對的單端序列片段。
本文所描述的分析假設從RNA-seq實驗獲得的序列片段已經與適當的參考基因組比對，並已經在基因水平上對序列進行了統計計數。在本文條件下，使用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使用該登記號查看。
3初始配置library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)
4數據整合4.1讀入計數數據為開始此分析，從https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file在線下載文件GSE63310_RAW.tar，並從壓縮包中解壓出相關的文件。下方的代碼將完成此步驟，或者您也可以手動進行這一步並繼續後續分析。
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對象。
4.2組織樣品信息為進行下遊分析，與實驗設計有關的樣品水平信息需要與計數矩陣的列關聯。這裡需要包括各種對表達水平有影響的實驗變量，無論是生物變量還是技術變量。例如，細胞類型（在這個實驗中是basal、LP和ML），基因型（野生型、敲除），表型（疾病狀態、性別、年齡），樣品處理（用藥、對照）和批次信息（如果樣品是在不同時間點進行收集和分析的，記錄進行實驗的時間）等。
我們的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) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$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
4.3組織基因注釋我們的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 <- genes
x
## 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 ...
5數據預處理5.1原始數據尺度轉換由於更深的測序總會產生更多的序列，在差異表達相關的分析中，我們很少使用原始的序列數。在實踐中，我們通常將原始的序列數進行歸一化，來消除測序深度所導致的差異。通常被使用的方法有基於序列的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那樣考慮到基因長度區別的影響，但在我們的分析中經常會用到它們。雖然也可以使用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函數會在進行log2轉換前給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-6
M <- median(x$samples$lib.size) * 1e-6
c(L, M)
## [1] 45.5 51.4
summary(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值。
5.2刪除低表達基因所有數據集中都混有表達的基因與不表達的基因。儘管我們想要檢測在一種條件中表達但再另一種條件中不表達的基因，也有一些基因在所有樣品中都不表達。實際上，這個數據集中19%的基因在所有九個樣品中的計數都是零。
table(rowSums(x$counts==0)==9)
##
## FALSE TRUE
## 22026 5153對log-CPM值的分布繪製的圖表顯示每個樣本中很大一部分基因都是不表達或者表達程度相當低的，它們的log-CPM值非常小甚至是負的（下圖A部分）。
在任何樣本中都沒有足夠多的序列片段的基因應該從下遊分析中過濾掉。這樣做的原因有好幾個。從生物學的角度來看，在任何條件下的表達水平都不具有生物學意義的基因都不值得關注，因此最好忽略。從統計學的角度來看，去除低表達計數基因使數據中的均值 - 方差關係可以得到更精確的估計，並且還減少了在觀察差異表達的下遊分析中需要進行的統計檢驗的數量。
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值為約02）。5.3歸一化基因表達分布在樣品製備或測序過程中，不具備生物學意義的外部因素會影響單個樣品的表達。比如說，在實驗中第一批製備的樣品會總體上表達高於第二批製備的樣品。假設所有樣品表達值的範圍和分布都應當相似，需要進行歸一化來確保整個實驗中每個樣本的表達分布都相似。
密度圖和箱線圖等展示每個樣品基因表達量分布的圖表可以用於判斷是否有樣品和其他樣品分布有差異。在此數據集中，所有樣品的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 <- x
x2$samples$norm.factors <- 1
x2$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.1090
lcpm <- 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%。5.4對樣本的無監督聚類在我們看來，用於檢查基因表達分析的最重要的探索性圖表之一便是MDS圖或其餘類似的圖。這種圖表使用無監督聚類方法展示出了樣品間的相似性和不相似性，能讓我們在進行正式的檢驗之前對於能檢測到多少差異表達基因有個大致概念。理想情況下，樣本會在不同的實驗組內很好的聚類，且可以鑑別出遠離所屬組的樣本，並追蹤誤差或額外方差的來源。如果存在技術重複，它們應當互相非常接近。
這樣的圖可以用limma中的plotMDS函數繪製。第一個維度表示能夠最好地分離樣品且解釋最大比例的方差的引導性的倍數變化（leading-fold-change），而後續的維度的影響更小，並與之前的維度正交。當實驗設計涉及到多個因子時，建議在多個維度上檢查每個因子。如果在其中一些維度上樣本可按照某因子聚類，這說明該因子對於表達差異有影響，在線性模型中應當將其包括進去。反之，沒有或者僅有微小影響的因子在下遊分析時應當被剔除。
在這個數據集中，可以看出樣本在維度1和2能很好地按照實驗分組聚類，隨後在維度3按照測序道（樣品批次）分離（如下圖所示）。請記住，第一維度解釋了數據中最大比例的方差，需要注意到，當我們向高維度移動，維度上的取值範圍會變小。
儘管所有樣本都按組聚類，在維度1上最大的轉錄差異出現在basal和LP以及basal和ML之間。因此，預期在basal樣品與其他之間的成對比較中能夠得到大量的DE基因，而在ML和LP之間的比較中得到的DE基因數量略少。在其他的數據集中，不按照實驗組聚類的樣本可能在下遊分析中只表現出較小的或不表現出差異表達。
為繪製MDS圖，我們為不同的因子賦予不同的色彩組合。維度1和2使用以細胞類型定義的色彩組合進行檢查。
維度3和4使用以測序泳道（批次）定義的色彩組合進行檢查。
lcpm <- cpm(x, log=TRUE)
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
col.lane <- lane
levels(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)交互式MDS圖連結[2]
6差異表達分析6.1創建設計矩陣和對比在此研究中，我們想知道哪些基因在我們研究的三組細胞之間以不同水平表達。在我們的分析中，假設基礎數據是正態分布的，為其擬合一個線性模型。在此之前，需要創建一個包含細胞類型以及測序泳道（批次）信息的設計矩陣。
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"對於一個給定的實驗，通常有幾種等價的方法可以創建一個合適的設計矩陣。比如說，~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的線性模型方法的核心優勢之一便是其適應任意實驗複雜程度的能力。簡單的設計，比如此工作流程中關於細胞類型和批次的實驗設計，直到更複雜的因子設計和含有交互作用項的模型，都能夠被相對簡單地處理。當實驗或技術效應可被隨機效應模型（random effect model）模擬時，limma中的另一種可能性是使用duplicateCorrelation函數來估計交互作用，這需要在此函數以及lmFit的線性建模步驟均指定一個block參數。
6.2從表達計數數據中刪除異方差據顯示對於RNA-seq計數數據而言，當使用原始計數或當其被轉換為log-CPM值時，方差並不獨立於均值(Law et al. 2014)。使用負二項分布來模擬計數的方法假設均值與方差間具有二次的關係。在limma中，假設log-CPM值符合正態分布，並使用由voom函數計算得到的精確權重來調整均值與方差的關係，從而對log-CPM值進行線性建模。
當操作DGEList對象時，voom從x中自動提取文庫大小和歸一化因子，以此將原始計數轉換為log-CPM值。在voom中，對於log-CPM值額外的歸一化可以通過設定normalize.method參數來進行。
下圖左側展示了這個數據集log-CPM值的均值-方差關係。通常而言，方差是測序實驗中的技術差異和不同細胞類型的重複樣本之間的生物學差異的結合，而voom圖會顯示出一個在均值與方差之間遞減的趨勢。生物學差異高的實驗通常會有更平坦的趨勢，其方差值在高表達處穩定。生物學差異低的實驗更傾向於急劇下降的趨勢。
不僅如此，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函數繪製的，它為進行log-CPM轉換後的數據擬合線性模型從而提取殘差方差。然後，對方差取平方根（或對標準差取平方根），並相對每個基因的平均表達作圖。均值通過平均計數加上2再進行log2轉換計算得到。右側的圖使用plotSA繪製了log2殘差標準差與log-CPM均值的關係。平均log2殘差標準差由水平藍線標出。在這兩幅圖中，每個黑點表示一個基因，紅線為對這些點的擬合。值得注意的是，DGEList對象中存儲的另一個數據框，即基因和樣本水平信息所存儲之處，保留在了voom創建的EList對象v中。v$genes數據框等同於x$genes，v$targets等同於x$samples，而v$E中所儲存的表達值類似於x$counts，儘管它進行了尺度轉換。此外，voom的EList對象中還有一個精確權重的矩陣v$weights，而設計矩陣存儲於v$design。
6.3擬合線性模型以進行比較limma的線性建模使用lmFit和contrasts.fit函數進行，它們原先是為微陣列而設計的。這些函數不僅可以用於微陣列數據，也可以用於RNA-seq數據。它們會單獨為每個基因的表達值擬合一個模型。然後，通過利用所有基因的信息來進行經驗貝葉斯調整，這樣可以獲得更精確的基因水平的變異程度估計(Smyth 2004)。下一圖為此模型的殘差關於平均表達值的圖。從圖中可以看出，方差不再與表達水平均值相關。
6.4檢查DE基因數量為快速查看差異表達水平，顯著上調或下調的基因可以匯總到一個表格中。顯著性的判斷使用校正p值閾值的默認值5%。在basal與LP的表達水平之間的比較中，發現了4648個在basal中相較於LP下調的基因和4863個在basal中相較於LP上調的基因，即共9511個DE基因。在basal和ML之間發現了一共9598個DE基因（4927個下調基因和4671個上調基因），而在LP和ML中發現了一共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統計值計算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] 2784
head(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")
6.5從上到下檢查單個DE基因使用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-06
head(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
6.6差異表達結果的實用圖形表示為可視化地總結所有基因的結果，可使用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統計圖中是做不到的。
glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1],
side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)交互式MD圖連結[3]
使用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[4]）轉換為圖形，使用Bootstrap庫處理界面並生成互動性可搜索的表格的數據表。這使得圖表可以在任何現代的瀏覽器中查看，對於從Rmarkdown分析報告中將其作為關聯文件而附加而言十分方便。
前文所展示的圖表中，一些展示了在任意一個條件下表達的所有基因（比如共同DE基因的韋恩圖或均值-差異圖），而另一些展示單獨的基因（交互性均值-差異圖右邊面板中所展示的log-CPM值）。而熱圖使用戶得以查看一部分基因的表達。這對於查看單個組或樣本的表達很有用，而不至於在關注於單個基因時失去對於研究整體的注意，也不會造成由於上千個基因所取平均值而導致的失去解析度。
使用gplots包的heatmap.2函數，我們為basal與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。給定基因相對高表達的樣本被標記為紅色，相對低表達的樣本被標記為藍色。淺色和白色代表中等表達水平的基因。樣本和基因已通過分層聚類的方法重新排序。圖中顯示有樣本聚類的樹狀圖。7使用camera的基因集檢驗在此次分析的最後，我們要進行一些基因集檢驗。為此，我們將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，從而獲得具有明確定義的生物狀態或過程。
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-10
cam.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-09
cam.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，但也可通過數據估計）和基因集的規模得到方差膨脹因子（variance inflation factor），並使用它調整基因集檢驗統計值的方差後，將會返回根據多重假設檢驗進行了校正的p值。
此實驗是與Lim等人(2010)(Lim et al. 2010)的數據集等價的RNA-seq，而他們使用Illumina微陣列分析了相同的分選細胞群，因此該早期文獻中的基因表達特徵出現在每種對比的列表頂部正符合我們的預期。在LP和ML的對比中，我們為Lim等人（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 al2010)非常類似於我們的，用了相同的分選方式來獲取不同的細胞群，只是他們使用的是微陣列而不是RNA-seq來測定基因表達。需要注意的是，上調基因集發生下調而下調基因集發生上調的逆相關性來自於對比的設定方式（LP相比於ML），如果將其對調，方向性將會吻合。limma也有其他的基因集檢驗，比如mroast(Wu et al. 2010)的自包含檢驗。雖然camera非常適合檢驗基因集的大型資料庫並觀察其中哪些相對於其他的在排序上位次更高（如前文所示），自包含檢驗更善於集中檢驗一個或少個選中的集合是否本身差異表達。換句話說，camera更適用於搜尋具有意義的基因集，而mroast測試的是已經確定有意義的基因集的顯著性。
8使用到的軟體和代碼此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.5.2 (2018-12-20)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.8-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.8-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] gplots_3.0.1 RColorBrewer_1.1-2
## [3] Mus.musculus_1.3.1 TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.4
## [5] org.Mm.eg.db_3.7.0 GO.db_3.7.0
## [7] OrganismDbi_1.24.0 GenomicFeatures_1.34.1
## [9] GenomicRanges_1.34.0 GenomeInfoDb_1.18.1
## [11] AnnotationDbi_1.44.0 IRanges_2.16.0
## [13] S4Vectors_0.20.1 Biobase_2.42.0
## [15] BiocGenerics_0.28.0 edgeR_3.24.2
## [17] Glimma_1.10.0 limma_3.38.3
## [19] knitr_1.21 BiocStyle_2.10.0
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.0 bit64_0.9-7 jsonlite_1.6
## [4] R.utils_2.7.0 gtools_3.8.1 assertthat_0.2.0
## [7] BiocManager_1.30.4 highr_0.7 RBGL_1.58.1
## [10] blob_1.1.1 GenomeInfoDbData_1.2.0 Rsamtools_1.34.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.22.0
## [19] htmltools_0.3.6 Matrix_1.2-15 R.oo_1.22.0
## [22] XML_3.98-1.16 pkgconfig_2.0.2 biomaRt_2.38.0
## [25] bookdown_0.9 zlibbioc_1.28.0 gdata_2.18.0
## [28] BiocParallel_1.16.2 SummarizedExperiment_1.12.0 magrittr_1.5
## [31] crayon_1.3.4 memoise_1.1.0 evaluate_0.12
## [34] R.methodsS3_1.7.1 graph_1.60.0 tools_3.5.2
## [37] prettyunits_1.0.2 hms_0.4.2 matrixStats_0.54.0
## [40] stringr_1.3.1 locfit_1.5-9.1 DelayedArray_0.8.0
## [43] Biostrings_2.50.1 compiler_3.5.2 caTools_1.17.1.1
## [46] rlang_0.3.0.1 grid_3.5.2 RCurl_1.95-4.11
## [49] bitops_1.0-6 rmarkdown_1.11 DBI_1.0.0
## [52] R6_2.3.0 GenomicAlignments_1.18.0 rtracklayer_1.42.1
## [55] bit_1.1-14 KernSmooth_2.23-15 stringi_1.2.4
## [58] Rcpp_1.0.0 xfun_0.4
參考文獻Bioconductor 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. 「limma 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.」 Bioinformatics 26: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.
參考資料[1]這裡: https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html
[2]交互式MDS圖連結: https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/glimma-plots/MDS-Plot.html
[3]交互式MD圖連結: https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/glimma-plots/MD-Plot.html
[4]https://d3js.org: https://d3js.org/