資料來(lái)源https://www.bioconductor.org/packages/devel/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html
摘自原文：
在這篇工作流程文章中，我們通過(guò)分析來(lái)自小鼠乳腺的RNA測(cè)序數(shù)據(jù)，示范了如何使用流行的
edgeR包載入、整理、過(guò)濾和歸一化數(shù)據(jù)，然后用limma包的voom方法、線性模型和經(jīng)驗(yàn)貝葉斯調(diào)節(jié)（empirical Bayes moderation）來(lái)評(píng)估差異表達(dá)并進(jìn)行基因集檢驗(yàn)。通過(guò)使用Glimma包，此流程得到了增進(jìn)，實(shí)現(xiàn)了結(jié)果的互動(dòng)探索，使用戶(hù)得以查看單個(gè)樣本與基因。 這三個(gè)軟件包提供的完整分析突出了研究人員可以使用Bioconductor輕松地從RNA測(cè)序?qū)嶒?yàn)的原始計(jì)數(shù)揭示生物學(xué)意義。
1. 初始配置
library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)
2. 數(shù)據(jù)整合
2.1 讀入計(jì)數(shù)數(shù)據(jù)
- 為開(kāi)始此分析，從https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file在線下載文件GSE63310_RAW.tar，并從壓縮包中解壓出相關(guān)的文件。下方的代碼將完成此步驟，或者也可以手動(dòng)進(jìn)行這一步并繼續(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)
- 每一個(gè)文本文件均包含一個(gè)給定樣品的原始基因水平計(jì)數(shù)。需要注意的是，我們的分析僅包含了此實(shí)驗(yàn)中的basal、LP和ML樣品（請(qǐng)查看下方相關(guān)文件名）。
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)
- 盡管這九個(gè)文本文件可以分別讀入R然后合并為一個(gè)計(jì)數(shù)矩陣，edgeR提供了更方便的途徑，使用readDGE函數(shù)即可一步完成。得到的DGEList對(duì)象中包含一個(gè)計(jì)數(shù)矩陣，它的27179行分別對(duì)應(yīng)唯一的Entrez基因標(biāo)識(shí)（ID），九列分別對(duì)應(yīng)此實(shí)驗(yàn)中的每個(gè)樣品。
x <- readDGE(files, columns=c(1,3))
class(x)
dim(x)
- 如果來(lái)自所有樣品的計(jì)數(shù)存儲(chǔ)在同一個(gè)文件中，數(shù)據(jù)可以首先讀入R再使用DGEList函數(shù)轉(zhuǎn)換為一個(gè)DGEList對(duì)象。
2.2 組織樣品信息
為進(jìn)行下游分析，與實(shí)驗(yàn)設(shè)計(jì)有關(guān)的樣品水平信息需要與計(jì)數(shù)矩陣的列關(guān)聯(lián)。這里需要包括各種對(duì)表達(dá)水平有影響的實(shí)驗(yàn)變量，無(wú)論是生物變量還是技術(shù)變量。例如，細(xì)胞類(lèi)型（在這個(gè)實(shí)驗(yàn)中是basal、LP和ML），基因型（野生型、敲除），表型（疾病狀態(tài)、性別、年齡），樣品處理（用藥、對(duì)照）和批次信息（如果樣品是在不同時(shí)間點(diǎn)進(jìn)行收集和分析的，記錄進(jìn)行實(shí)驗(yàn)的時(shí)間）等。
我們的
DGEList對(duì)象中包含的samples數(shù)據(jù)框同時(shí)存儲(chǔ)了細(xì)胞類(lèi)型（group）和批次（測(cè)序泳道lane）信息，每種信息都包含三個(gè)不同的水平。需要注意的是，在x$samples中，程序會(huì)自動(dòng)計(jì)算每個(gè)樣品的文庫(kù)大小，歸一化系數(shù)會(huì)被設(shè)置為1。 為了簡(jiǎn)單起見(jiàn)，我們從我們的DGEList對(duì)象x的列名中刪去了GEO樣品ID（GSM*）。
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))
samplenames
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
2.3 組織基因注釋
我們的DGEList對(duì)象中的第二個(gè)數(shù)據(jù)框名為genes，用于存儲(chǔ)與計(jì)數(shù)矩陣的行相關(guān)聯(lián)的基因水平的信息。 為檢索這些信息，我們可以使用包含特定物種信息的包，比如小鼠的Mus.musculus (Bioconductor Core Team 2016b)（或人類(lèi)的Homo.sapiens (Bioconductor Core Team 2016a)）；或者也可以使用biomaRt 包 (Durinck et al. 2005, 2009)，它通過(guò)接入Ensembl genome數(shù)據(jù)庫(kù)來(lái)進(jìn)行基因注釋。
可以檢索的信息類(lèi)型包括基因符號(hào)（gene symbols）、基因名稱(chēng)（gene names）、染色體名稱(chēng)和位置（chromosome names and locations）、Entrez基因ID（Entrez gene IDs）、Refseq基因ID（Refseq gene IDs）和Ensembl基因ID（Ensembl gene IDs）等。biomaRt主要使用Ensembl基因ID進(jìn)行檢索，而由于Mus.musculus包含多種不同來(lái)源的信息，它允許用戶(hù)從多種不同基因ID中選取檢索鍵。
我們使用Mus.musculus包，利用我們數(shù)據(jù)集中的Entrez基因ID來(lái)檢索相關(guān)的基因符號(hào)和染色體信息。
geneid <- rownames(x)
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),
keytype="ENTREZID")
head(genes)
- 與任何基因ID一樣，Entrez基因ID可能不能一對(duì)一地匹配我們想獲得的基因信息。在處理之前，檢查重復(fù)的基因ID和弄清楚重復(fù)的來(lái)源非常重要。我們的基因注釋中包含28個(gè)匹配到不同染色體的基因（比如基因Gm1987關(guān)聯(lián)于染色體chr4和chr4_JH584294_random，小RNA Mir5098關(guān)聯(lián)于chr2，chr5，chr8，chr11和chr17）。
- 為了處理重復(fù)的基因ID，我們可以合并來(lái)自多重匹配基因的所有染色體信息，比如將基因Gm1987分配到chr4 and chr4_JH584294_random，或選取其中一條染色體來(lái)代表具有重復(fù)注釋的基因。為了簡(jiǎn)單起見(jiàn)，我們選擇后者，保留每個(gè)基因ID第一次出現(xiàn)的信息。
genes <- genes[!duplicated(genes$ENTREZID),]
- 在此例子中，注釋與數(shù)據(jù)對(duì)象中的基因順序是相同的。
- 如果由于缺失和／或重新排列基因ID導(dǎo)致其順序不一致，可以用match來(lái)正確排序基因。
- 然后將基因注釋的數(shù)據(jù)框加入數(shù)據(jù)對(duì)象，數(shù)據(jù)即被整潔地整理入一個(gè)DGEList對(duì)象中，它包含原始計(jì)數(shù)數(shù)據(jù)和相關(guān)的樣品信息和基因注釋。
x$genes <- genes
x
3. 數(shù)據(jù)預(yù)處理
3.1 原始數(shù)據(jù)尺度轉(zhuǎn)換
由于更深的測(cè)序總會(huì)產(chǎn)生更多的序列，在差異表達(dá)相關(guān)的分析中，我們很少使用原始的序列數(shù)。在實(shí)踐中，我們通常將原始的序列數(shù)進(jìn)行歸一化，來(lái)消除測(cè)序深度所導(dǎo)致的差異。通常被使用的方法有基于序列的
CPM（counts per million）、log-CPM、FPKM（fragments per kilobase of transcript per million），和基于轉(zhuǎn)錄本數(shù)目的RPKM（reads per kilobase of transcript per million）。盡管CPM和log-CPM轉(zhuǎn)換并不像RPKM和FPKM那樣考慮到基因長(zhǎng)度區(qū)別的影響，但在我們的分析中經(jīng)常會(huì)用到它們。雖然也可以使用RPKM和FPKM，
但CPM和log-CPM只使用計(jì)數(shù)矩陣即可計(jì)算，且已足以滿(mǎn)足我們所關(guān)注的比較的需要。假設(shè)不同條件之間剪接異構(gòu)體（isoform）的使用沒(méi)有差別，差異表達(dá)分析研究同一基因在不同條件下的表達(dá)差異，而不是比較多個(gè)基因之間的表達(dá)或測(cè)定絕對(duì)表達(dá)量。換而言之，基因長(zhǎng)度在我們關(guān)注的比較中始終不變，且任何觀測(cè)到的差異是來(lái)自于條件的變化而不是基因長(zhǎng)度的變化。在此處，使用
edgeR中的cpm函數(shù)將原始計(jì)數(shù)轉(zhuǎn)換為CPM和log-CPM值。如果可以提供基因長(zhǎng)度信息，可以使用edgeR中的rpkm函數(shù)計(jì)算RPKM值，就像計(jì)算CPM值那樣簡(jiǎn)單。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)
對(duì)于一個(gè)基因，CPM值為1相當(dāng)于在測(cè)序深度最低的樣品中（JMS9-P8c, 序列數(shù)量約2千萬(wàn)）有20個(gè)計(jì)數(shù)，或者在測(cè)序深度最高的樣品中（JMS8-3，序列數(shù)量約7.6千萬(wàn)）有76個(gè)計(jì)數(shù)。
log-CPM值將被用于探索性圖表中。當(dāng)設(shè)置
log=TRUE時(shí)，cpm函數(shù)會(huì)在進(jìn)行l(wèi)og2轉(zhuǎn)換前給CPM值加上一個(gè)彌補(bǔ)值。默認(rèn)的彌補(bǔ)值是2/L，其中2是“預(yù)先計(jì)數(shù)”，而L是樣本總序列數(shù)（以百萬(wàn)計(jì)）的平均值，所以log-CPM值是根據(jù)CPM值通過(guò)計(jì)算得到的。這樣的計(jì)算方式可以確保任意兩個(gè)具有相同CPM值的序列片段計(jì)數(shù)的
log-CPM值也相同。彌補(bǔ)值的使用可以避免對(duì)零取對(duì)數(shù)，并能使所有樣本間的log倍數(shù)變化（log-fold-change）向0推移而減小低表達(dá)基因間微小計(jì)數(shù)變化帶來(lái)的巨大的偽差異性，這對(duì)于繪制探索性圖表很有用。在這個(gè)數(shù)據(jù)集中，平均的樣本總序列數(shù)是4.55千萬(wàn)，所以L約等于45.5，且每個(gè)樣本中的最小log-CPM值為。換而言之，
在加上了預(yù)先計(jì)數(shù)彌補(bǔ)值后，此數(shù)據(jù)集中的零表達(dá)計(jì)數(shù)對(duì)應(yīng)的log-CPM值為-4.51：
L <- mean(x$samples$lib.size) * 1e-6
M <- median(x$samples$lib.size) * 1e-6
c(L, M)
summary(lcpm)
- 在下游的線性模型分析中，使用
limma的voom函數(shù)時(shí)也會(huì)用到log-CPM值，但voom會(huì)默認(rèn)使用更小的預(yù)先計(jì)數(shù)重新計(jì)算自己的log-CPM值。
3.2 刪除低表達(dá)基因
- 所有數(shù)據(jù)集中都混有表達(dá)的基因與不表達(dá)的基因。盡管我們想要檢測(cè)在一種條件中表達(dá)但再另一種條件中不表達(dá)的基因，也有一些基因在所有樣品中都不表達(dá)。實(shí)際上，這個(gè)數(shù)據(jù)集中19%的基因在所有九個(gè)樣品中的計(jì)數(shù)都是零。
table(rowSums(x$counts==0)==9)
對(duì)log-CPM值的分布繪制的圖表顯示每個(gè)樣本中很大一部分基因都是不表達(dá)或者表達(dá)程度相當(dāng)?shù)偷模鼈兊膌og-CPM值非常小甚至是負(fù)的（下圖A部分）。
在任何樣本中都沒(méi)有足夠多的序列片段的基因應(yīng)該從下游分析中過(guò)濾掉。這樣做的原因有好幾個(gè)。 從生物學(xué)的角度來(lái)看，在任何條件下的表達(dá)水平都不具有生物學(xué)意義的基因都不值得關(guān)注，因此最好忽略。 從統(tǒng)計(jì)學(xué)的角度來(lái)看，去除低表達(dá)計(jì)數(shù)基因使數(shù)據(jù)中的均值 - 方差關(guān)系可以得到更精確的估計(jì)，并且還減少了在觀察差異表達(dá)的下游分析中需要進(jìn)行的統(tǒng)計(jì)檢驗(yàn)的數(shù)量。
edgeR包中的filterByExpr函數(shù)提供了自動(dòng)過(guò)濾基因的方法，可保留盡可能多的有足夠表達(dá)計(jì)數(shù)的基因。
keep.exprs <- filterByExpr(x, group=group)
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
dim(x)
此函數(shù)默認(rèn)選取最小的組內(nèi)的樣本數(shù)量為最小樣本數(shù)，保留至少在這個(gè)數(shù)量的樣本中有10個(gè)或更多序列片段計(jì)數(shù)的基因。對(duì)基因表達(dá)量進(jìn)行過(guò)濾時(shí)使用CPM值而不是表達(dá)計(jì)數(shù)，以避免對(duì)總序列數(shù)大的樣本的偏向性。在這個(gè)數(shù)據(jù)集中，總序列數(shù)的中位數(shù)是5.1千萬(wàn)，且10/51約等于0.2，所以
filterByExpr函數(shù)保留在至少三個(gè)樣本中CPM值大于等于0.2的基因。此處，一個(gè)具有生物學(xué)意義的基因需要在至少三個(gè)樣本中表達(dá)，因?yàn)槿N細(xì)胞類(lèi)型組內(nèi)各有三個(gè)重復(fù)。所使用的閾值取決于測(cè)序深度和實(shí)驗(yàn)設(shè)計(jì)。如果樣本總表達(dá)計(jì)數(shù)量增大，那么可以選擇更低的CPM閾值，因?yàn)楦蟮目偙磉_(dá)計(jì)數(shù)量提供了更好的分辨率來(lái)探究更多表達(dá)水平更低的基因。使用這個(gè)標(biāo)準(zhǔn)，基因的數(shù)量減少到了16624個(gè)，約為開(kāi)始時(shí)數(shù)量的60%。過(guò)濾后的log-CPM值顯示出每個(gè)樣本的分布基本相同（下圖B部分）。需要注意的是，從整個(gè)DGEList對(duì)象中取子集時(shí)同時(shí)刪除了被過(guò)濾的基因的計(jì)數(shù)和其相關(guān)的基因信息。過(guò)濾后的DGEList對(duì)象為留下的基因保留了相對(duì)應(yīng)的基因信息和計(jì)數(shù)。
下方給出的是繪圖所用代碼。
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")
3.3 歸一化基因表達(dá)分布
在樣品制備或測(cè)序過(guò)程中，不具備生物學(xué)意義的外部因素會(huì)影響單個(gè)樣品的表達(dá)。比如說(shuō)，在實(shí)驗(yàn)中第一批制備的樣品會(huì)總體上表達(dá)高于第二批制備的樣品。假設(shè)所有樣品表達(dá)值的范圍和分布都應(yīng)當(dāng)相似，需要進(jìn)行歸一化來(lái)確保整個(gè)實(shí)驗(yàn)中每個(gè)樣本的表達(dá)分布都相似。
密度圖和箱線圖等展示每個(gè)樣品基因表達(dá)量分布的圖表可以用于判斷是否有樣品和其他樣品分布有差異。在此數(shù)據(jù)集中，所有樣品的log-CPM分布都很相似（上圖B部分）。
盡管如此，我們依然需要使用
edgeR中的calcNormFactors函數(shù)，用TMM(Robinson and Oshlack 2010)方法進(jìn)行歸一化。此處計(jì)算得到的歸一化系數(shù)被用作文庫(kù)大小的縮放系數(shù)。當(dāng)我們使用DGEList對(duì)象時(shí)，這些歸一化系數(shù)被自動(dòng)存儲(chǔ)在x$samples$norm.factors。對(duì)此數(shù)據(jù)集而言，TMM歸一化的作用比較溫和，這體現(xiàn)在所有的縮放因子都相對(duì)接近1。
x <- calcNormFactors(x, method = "TMM")
x$samples$norm.factors
- 為了更好地可視化表現(xiàn)出歸一化的影響，我們復(fù)制了數(shù)據(jù)并進(jìn)行了調(diào)整，使得第一個(gè)樣品的計(jì)數(shù)減少到了其原始值的5%，而第二個(gè)樣品增大到了5倍。
x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5
- 下圖顯示了沒(méi)有經(jīng)過(guò)歸一化的與經(jīng)過(guò)了歸一化的數(shù)據(jù)的樣本的表達(dá)分布，其中歸一化前的分布顯然不同，而歸一化后比較相似。此處，第一個(gè)樣品的
TMM縮放系數(shù)0.06非常小，而第二個(gè)樣品的縮放系數(shù)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
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")
3.4 對(duì)樣本的無(wú)監(jiān)督聚類(lèi)
在我們看來(lái)，用于檢查基因表達(dá)分析的最重要的探索性圖表之一便是MDS圖或其余類(lèi)似的圖。這種圖表使用無(wú)監(jiān)督聚類(lèi)方法展示出了樣品間的相似性和不相似性，能讓我們?cè)谶M(jìn)行正式的檢驗(yàn)之前對(duì)于能檢測(cè)到多少差異表達(dá)基因有個(gè)大致概念。理想情況下，樣本會(huì)在不同的實(shí)驗(yàn)組內(nèi)很好的聚類(lèi)，且可以鑒別出遠(yuǎn)離所屬組的樣本，并追蹤誤差或額外方差的來(lái)源。如果存在技術(shù)重復(fù)，它們應(yīng)當(dāng)互相非常接近。
這樣的圖可以用
limma中的plotMDS函數(shù)繪制。第一個(gè)維度表示能夠最好地分離樣品且解釋最大比例的方差的引導(dǎo)性的倍數(shù)變化（leading-fold-change），而后續(xù)的維度的影響更小，并與之前的維度正交。當(dāng)實(shí)驗(yàn)設(shè)計(jì)涉及到多個(gè)因子時(shí)，建議在多個(gè)維度上檢查每個(gè)因子。如果在其中一些維度上樣本可按照某因子聚類(lèi)，這說(shuō)明該因子對(duì)于表達(dá)差異有影響，在線性模型中應(yīng)當(dāng)將其包括進(jìn)去。反之，沒(méi)有或者僅有微小影響的因子在下游分析時(shí)應(yīng)當(dāng)被剔除。在這個(gè)數(shù)據(jù)集中，可以看出樣本在維度1和2能很好地按照實(shí)驗(yàn)分組聚類(lèi)，隨后在維度3按照測(cè)序道（樣品批次）分離（如下圖所示）。請(qǐng)記住，第一維度解釋了數(shù)據(jù)中最大比例的方差，需要注意到，當(dāng)我們向高維度移動(dòng)，維度上的取值范圍會(huì)變小。
盡管所有樣本都按組聚類(lèi)，在維度1上最大的轉(zhuǎn)錄差異出現(xiàn)在basal和LP以及basal和ML之間。因此，預(yù)期在basal樣品與其他之間的成對(duì)比較中能夠得到大量的DE基因，而在ML和LP之間的比較中得到的DE基因數(shù)量略少。在其他的數(shù)據(jù)集中，不按照實(shí)驗(yàn)組聚類(lèi)的樣本可能在下游分析中只表現(xiàn)出較小的或不表現(xiàn)出差異表達(dá)。
為繪制MDS圖，我們?yōu)椴煌囊蜃淤x予不同的色彩組合。維度1和2使用以細(xì)胞類(lèi)型定義的色彩組合進(jìn)行檢查。
維度3和4使用以測(cè)序泳道（批次）定義的色彩組合進(jìn)行檢查。
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")
作為另一種選擇，Glimma包也提供了便于探索多個(gè)維度的交互式MDS圖。其中的
glMDSPlot函數(shù)可生成一個(gè)html網(wǎng)頁(yè)（如果設(shè)置launch=TRUE，將會(huì)在瀏覽器中打開(kāi)），其左側(cè)面板含有一張MDS圖，而右側(cè)面板包含一張展示了各個(gè)維度所解釋的方差比例的柱形圖。點(diǎn)擊柱形圖中的柱可切換MDS圖繪制時(shí)所使用的維度，且將鼠標(biāo)懸浮于單個(gè)點(diǎn)上可顯示相應(yīng)的樣本標(biāo)簽。也可切換配色方案，以突顯不同細(xì)胞類(lèi)型或測(cè)序泳道（批次）。此數(shù)據(jù)集的交互式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圖鏈接(鏈接好像是掛了，沒(méi)找到)
4. 差異表達(dá)分析
4.1 創(chuàng)建設(shè)計(jì)矩陣和對(duì)比
- 在此研究中，我們想知道哪些基因在我們研究的三組細(xì)胞之間以不同水平表達(dá)。在我們的分析中，假設(shè)基礎(chǔ)數(shù)據(jù)是正態(tài)分布的，為其擬合一個(gè)線性模型。在此之前，需要?jiǎng)?chuàng)建一個(gè)包含細(xì)胞類(lèi)型以及測(cè)序泳道（批次）信息的設(shè)計(jì)矩陣。
design <- model.matrix(~0+group+lane)
colnames(design) <- gsub("group", "", colnames(design))
design
- 對(duì)于一個(gè)給定的實(shí)驗(yàn)，通常有幾種等價(jià)的方法可以創(chuàng)建一個(gè)合適的設(shè)計(jì)矩陣。 比如說(shuō)，
~0+group+lane去除了第一個(gè)因子group的截距，但第二個(gè)因子lane的截距被保留。 此外也可以使用~group+lane，來(lái)自group和lane的截距均被保留。 此處的關(guān)鍵是理解如何解釋給定模型中估計(jì)得到的系數(shù)。 我們?cè)诖朔治鲋羞x取第一種模型，因?yàn)樵跊](méi)有g(shù)roup的截距的情況下能更直截了當(dāng)?shù)卦O(shè)定模型中的對(duì)比。用于細(xì)胞群之間成對(duì)比較的對(duì)比可以在limma中用makeContrasts函數(shù)設(shè)定。
contr.matrix <- makeContrasts(
BasalvsLP = Basal-LP,
BasalvsML = Basal - ML,
LPvsML = LP - ML,
levels = colnames(design))
contr.matrix
-limma的線性模型方法的核心優(yōu)勢(shì)之一便是其適應(yīng)任意實(shí)驗(yàn)復(fù)雜程度的能力。簡(jiǎn)單的設(shè)計(jì)，比如此工作流程中關(guān)于細(xì)胞類(lèi)型和批次的實(shí)驗(yàn)設(shè)計(jì)，直到更復(fù)雜的因子設(shè)計(jì)和含有交互作用項(xiàng)的模型，都能夠被相對(duì)簡(jiǎn)單地處理。當(dāng)實(shí)驗(yàn)或技術(shù)效應(yīng)可被隨機(jī)效應(yīng)模型（random effect model）模擬時(shí)，limma中的另一種可能性是使用duplicateCorrelation函數(shù)來(lái)估計(jì)交互作用，這需要在此函數(shù)以及lmFit的線性建模步驟均指定一個(gè)block參數(shù)。
4.2 從表達(dá)計(jì)數(shù)數(shù)據(jù)中刪除異方差
據(jù)顯示對(duì)于RNA-seq計(jì)數(shù)數(shù)據(jù)而言，當(dāng)使用原始計(jì)數(shù)或當(dāng)其被轉(zhuǎn)換為log-CPM值時(shí)，方差并不獨(dú)立于均值(Law et al. 2014)。使用負(fù)二項(xiàng)分布來(lái)模擬計(jì)數(shù)的方法假設(shè)均值與方差間具有二次的關(guān)系。在
limma中，假設(shè)log-CPM值符合正態(tài)分布，并使用由voom函數(shù)計(jì)算得到的精確權(quán)重來(lái)調(diào)整均值與方差的關(guān)系，從而對(duì)log-CPM值進(jìn)行線性建模。當(dāng)操作DGEList對(duì)象時(shí)，voom從x中自動(dòng)提取文庫(kù)大小和歸一化因子，以此將原始計(jì)數(shù)轉(zhuǎn)換為log-CPM值。在voom中，對(duì)于log-CPM值額外的歸一化可以通過(guò)設(shè)定
normalize.method參數(shù)來(lái)進(jìn)行。下圖左側(cè)展示了這個(gè)數(shù)據(jù)集
log-CPM值的均值-方差關(guān)系。通常而言，方差是測(cè)序?qū)嶒?yàn)中的技術(shù)差異和不同細(xì)胞類(lèi)型的重復(fù)樣本之間的生物學(xué)差異的結(jié)合，而voom圖會(huì)顯示出一個(gè)在均值與方差之間遞減的趨勢(shì)。 生物學(xué)差異高的實(shí)驗(yàn)通常會(huì)有更平坦的趨勢(shì)，其方差值在高表達(dá)處穩(wěn)定。 生物學(xué)差異低的實(shí)驗(yàn)更傾向于急劇下降的趨勢(shì)。不僅如此，
voom圖也提供了對(duì)于上游所進(jìn)行的過(guò)濾水平的可視化檢測(cè)。如果對(duì)于低表達(dá)基因的過(guò)濾不夠充分，在圖上表達(dá)低的一端，受到非常低的表達(dá)計(jì)數(shù)的影響，可以觀察到方差水平的下降。如果觀察到了這種情況，應(yīng)當(dāng)回到最初的過(guò)濾步驟并提高用于該數(shù)據(jù)集的表達(dá)閾值。當(dāng)前面觀察的MDS圖中具有明顯的樣本水平的差異時(shí)，可以用
voomWithQualityWeights函數(shù)來(lái)同時(shí)合并樣本水平的權(quán)重和voom(Liu et al. 2015)估算得到的豐度相關(guān)的權(quán)重。關(guān)于此種方式的例子參見(jiàn)Liu等(2016) (Liu et al. 2016)。
par(mfrow=c(1,2))
v <- voom(x, design, plot=TRUE)
v
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
值得注意的是，DGEList對(duì)象中存儲(chǔ)的另一個(gè)數(shù)據(jù)框，即基因和樣本水平信息所存儲(chǔ)之處，保留在了voom創(chuàng)建的EList對(duì)象v中。
v$genes數(shù)據(jù)框等同于x$genes，v$targets等同于x$samples，而v$E中所儲(chǔ)存的表達(dá)值類(lèi)似于x$counts，盡管它進(jìn)行了尺度轉(zhuǎn)換。此外，voom的EList對(duì)象中還有一個(gè)精確權(quán)重的矩陣v$weights，而設(shè)計(jì)矩陣存儲(chǔ)于v$design。
4.3 擬合線性模型以進(jìn)行比較
-
limma的線性建模使用lmFit和contrasts.fit函數(shù)進(jìn)行，它們?cè)仁菫槲㈥嚵卸O(shè)計(jì)的。這些函數(shù)不僅可以用于微陣列數(shù)據(jù)，也可以用于RNA-seq數(shù)據(jù)。它們會(huì)單獨(dú)為每個(gè)基因的表達(dá)值擬合一個(gè)模型。 - 然后，通過(guò)利用所有基因的信息來(lái)進(jìn)行經(jīng)驗(yàn)貝葉斯調(diào)整，這樣可以獲得更精確的基因水平的變異程度估計(jì)(Smyth 2004)。下一圖為此模型的殘差關(guān)于平均表達(dá)值的圖。從圖中可以看出，方差不再與表達(dá)水平均值相關(guān)。
4.4 檢查DE基因數(shù)量
- 為快速查看差異表達(dá)水平，顯著上調(diào)或下調(diào)的基因可以匯總到一個(gè)表格中。
- 顯著性的判斷使用校正p值閾值的默認(rèn)值5%。
- 在basal與LP的表達(dá)水平之間的比較中，發(fā)現(xiàn)了4648個(gè)在basal中相較于LP下調(diào)的基因和4863個(gè)在basal中相較于LP上調(diào)的基因，即共9511個(gè)DE基因。
- 在basal和ML之間發(fā)現(xiàn)了一共9598個(gè)DE基因（4927個(gè)下調(diào)基因和4671個(gè)上調(diào)基因），而在LP和ML中發(fā)現(xiàn)了一共5652個(gè)DE基因（3135個(gè)下調(diào)基因和2517個(gè)上調(diào)基因）。
- 在包括basal細(xì)胞類(lèi)型的比較中皆找到了大量的DE基因，這與我們?cè)贛DS圖中觀察到的結(jié)果相吻合。
summary(decideTests(efit))
- 一些研究中不僅僅需要使用校正p值閾值，更為嚴(yán)格定義的顯著性可能需要差異倍數(shù)的對(duì)數(shù)（log-FCs）也高于某個(gè)最小值。
- treat方法(McCarthy and Smyth 2009)可以按照對(duì)最小log-FC值的要求，使用經(jīng)過(guò)經(jīng)驗(yàn)貝葉斯調(diào)整的t統(tǒng)計(jì)值計(jì)算p值。
- 當(dāng)我們的檢驗(yàn)要求基因的log-FC顯著大于1（等同于在原本的尺度上不同細(xì)胞類(lèi)型之間差兩倍）時(shí)，差異表達(dá)基因的數(shù)量得到了下降，basal與LP相比只有3684個(gè)DE基因，basal與ML相比只有3834個(gè)DE基因，LP與ML相比只有414個(gè)DE基因。
tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)
- 在多種比較中皆差異表達(dá)的基因可以從
decideTests的結(jié)果中提取，其中的0代表不差異表達(dá)的基因，1代表上調(diào)的基因，-1代表下調(diào)的基因。 - 共有2784個(gè)基因在basal和LP以及basal和ML的比較中都差異表達(dá)，其中的20個(gè)于下方列出。
write.fit函數(shù)可用于將三個(gè)比較的結(jié)果提取并寫(xiě)入到單個(gè)輸出文件。
de.common <- which(dt[,1]!=0 & dt[,2]!=0)
length(de.common)
head(tfit$genes$SYMBOL[de.common], n=20)
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
write.fit(tfit, dt, file="results.txt")
4.5 從上到下檢查單個(gè)DE基因
- 使用
topTreat函數(shù)可以列舉出使用treat得到的結(jié)果中靠前的DE基因（對(duì)于eBayes的結(jié)果可以使用topTable函數(shù)）。默認(rèn)情況下，topTreat將基因按照校正p值從小到大排列，并為每個(gè)基因給出相關(guān)的基因信息、log-FC、平均log-CPM、校正t值、原始及經(jīng)過(guò)多重假設(shè)檢驗(yàn)校正的p值。列出前多少個(gè)基因的數(shù)量可由用戶(hù)指定，如果設(shè)為n=Inf則會(huì)包括所有的基因?；駽ldn7和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)
head(basal.vs.ml)
4.6差異表達(dá)結(jié)果的實(shí)用圖形表示
- 為可視化地總結(jié)所有基因的結(jié)果，可使用
plotMD函數(shù)繪制均值-差異（MD）圖，其中展示了線性模型擬合所得到的log-FC與log-CPM平均值間的關(guān)系，而差異表達(dá)的基因會(huì)被重點(diǎn)標(biāo)出。
plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1],
xlim=c(-8,13))
-Glimma的glMDPlot函數(shù)提供了交互式的均值-差異圖，拓展了這種圖表的功能性。 此函數(shù)的輸出為一個(gè)html頁(yè)面，左側(cè)面板為結(jié)果的總結(jié)性圖表（與plotMD的輸出類(lèi)似），右側(cè)面板包含各個(gè)樣本的log-CPM值，下方為結(jié)果的表格。 這種交互式展示允許用戶(hù)使用提供的注釋?zhuān)ū热缁蛎麡?biāo)識(shí)）搜索特定基因，而這在R統(tǒng)計(jì)圖中是做不到的。
glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1],
side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)
交互式MD圖鏈接(此處鏈接已掛)
Glimma提供的交互性使得單個(gè)圖形窗口內(nèi)能夠呈現(xiàn)出額外的信息。Glimma是以R和Javascript實(shí)現(xiàn)的，使用R代碼生成數(shù)據(jù)，并在之后使用Javascript庫(kù)D3（https://d3js.org）轉(zhuǎn)換為圖形，使用Bootstrap庫(kù)處理界面并生成互動(dòng)性可搜索的表格的數(shù)據(jù)表。 這使得圖表可以在任何現(xiàn)代的瀏覽器中查看，對(duì)于從Rmarkdown分析報(bào)告中將其作為關(guān)聯(lián)文件而附加而言十分方便。前文所展示的圖表中，一些展示了在任意一個(gè)條件下表達(dá)的所有基因（比如共同DE基因的韋恩圖或均值-差異圖），而另一些展示單獨(dú)的基因（交互性均值-差異圖右邊面板中所展示的log-CPM值）。而熱圖使用戶(hù)得以查看一部分基因的表達(dá)。這對(duì)于查看單個(gè)組或樣本的表達(dá)很有用，而不至于在關(guān)注于單個(gè)基因時(shí)失去對(duì)于研究整體的注意，也不會(huì)造成由于上千個(gè)基因所取平均值而導(dǎo)致的失去分辨率。
使用
gplots包的heatmap.2函數(shù)，我們?yōu)閎asal與LP的對(duì)照中前100個(gè)DE基因（按調(diào)整p值排序）繪制了一幅熱圖。熱圖中正確地將樣本按照細(xì)胞類(lèi)型聚類(lèi)，并重新排序了基因，形成了表達(dá)相似的塊狀。從熱圖中，我們觀察到對(duì)于basal與LP之間的前100個(gè)DE基因，ML和LP樣本的表達(dá)非常相似。
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")
圖是出來(lái)了，但是報(bào)了個(gè)錯(cuò)
Error in plot.new() : figure margins too large比對(duì)了一下原圖發(fā)現(xiàn)少了點(diǎn)東西：Row Z-Score 和樣本間關(guān)系樹(shù)（不是太明白哪里的設(shè)置出了問(wèn)題，暫時(shí)放著）
5. 使用camera的基因集檢驗(yàn)
- 在此次分析的最后，我們要進(jìn)行一些基因集檢驗(yàn)。為此，我們將camera方法(Wu and Smyth 2012)應(yīng)用于Broad Institute的MSigDB c2中的(Subramanian et al. 2005)中適應(yīng)小鼠的c2基因表達(dá)特征，這可從http://bioinf.wehi.edu.au/software/MSigDB/以
RData對(duì)象格式獲取。 此外，對(duì)于人類(lèi)和小鼠，來(lái)自MSigDB的其他有用的基因集也可從此網(wǎng)站獲取，比如標(biāo)志（hallmark）基因集。C2基因集的內(nèi)容收集自在線數(shù)據(jù)庫(kù)、出版物以及該領(lǐng)域?qū)＜遥鴺?biāo)志基因集的內(nèi)容來(lái)自MSigDB，從而獲得具有明確定義的生物狀態(tài)或過(guò)程。
這一步需要用到一個(gè)包Bioconductor工作流程包RNAseq123
BiocManager::install("RNAseq123")
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)
cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3])
head(cam.LPvsML,5)
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP,
index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
camera函數(shù)通過(guò)比較假設(shè)檢驗(yàn)來(lái)評(píng)估一個(gè)給定基因集中的基因是否相對(duì)于不在集內(nèi)的基因而言在差異表達(dá)基因的排序中更靠前。 它使用limma的線性模型框架，并同時(shí)采用設(shè)計(jì)矩陣和對(duì)比矩陣（如果有的話），且在測(cè)試的過(guò)程中會(huì)使用來(lái)自voom的觀測(cè)水平權(quán)重。 在通過(guò)基因間相關(guān)性（默認(rèn)設(shè)定為0.01，但也可通過(guò)數(shù)據(jù)估計(jì)）和基因集的規(guī)模得到方差膨脹因子（variance inflation factor），并使用它調(diào)整基因集檢驗(yàn)統(tǒng)計(jì)值的方差后，將會(huì)返回根據(jù)多重假設(shè)檢驗(yàn)進(jìn)行了校正的p值。-
此實(shí)驗(yàn)是與Lim等人(2010)(Lim et al. 2010)的數(shù)據(jù)集等價(jià)的RNA-seq，而他們使用Illumina微陣列分析了相同的分選細(xì)胞群，因此該早期文獻(xiàn)中的基因表達(dá)特征出現(xiàn)在每種對(duì)比的列表頂部正符合我們的預(yù)期。在LP和ML的對(duì)比中，我們?yōu)長(zhǎng)im等人（2010）的成熟管腔基因集（上調(diào)及下調(diào)）繪制了條碼圖（barcodeplot）。需要注意的是，由于我們的對(duì)比是將LP與ML相比而不是相反，這些基因集的方向在我們的數(shù)據(jù)集中是反過(guò)來(lái)的（如果將對(duì)比反過(guò)來(lái)，基因集的方向?qū)?huì)與對(duì)比一致）。
Figure 8 limma也有其他的基因集檢驗(yàn)，比如mroast(Wu et al. 2010)的自包含檢驗(yàn)。雖然camera非常適合檢驗(yàn)基因集的大型數(shù)據(jù)庫(kù)并觀察其中哪些相對(duì)于其他的在排序上位次更高（如前文所示），自包含檢驗(yàn)更善于集中檢驗(yàn)一個(gè)或少個(gè)選中的集合是否本身差異表達(dá)。換句話說(shuō)，camera更適用于搜尋具有意義的基因集，而mroast測(cè)試的是已經(jīng)確定有意義的基因集的顯著性。
6 使用到的軟件和代碼
- 此RNA-seq工作流程使用了Bioconductor項(xiàng)目3.8版本中的多個(gè)軟件包，運(yùn)行于R 3.5.1或更高版本。除了本文中重點(diǎn)提到的軟件（limma、Glimma以及edgeR），亦需要一些其他軟件包，包括gplots和RColorBrewer還有基因注釋包Mus.musculus。 此文檔使用knitr編譯。所有用到的包的版本號(hào)如下所示。 Bioconductor工作流程包RNAseq123（可訪問(wèn)https://bioconductor.org/packages/RNAseq123查看）內(nèi)包含此文章的英文和簡(jiǎn)體中文版以及進(jìn)行整個(gè)分析流程所需要的代碼。安裝此包即可管理以上提到的所有需要的包。對(duì)于RNA-seq數(shù)據(jù)分析實(shí)踐培訓(xùn)而言，此包也是非常有用的資源。
sessionInfo()
> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] gplots_3.0.1.1 RColorBrewer_1.1-2
[3] Mus.musculus_1.3.1 TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.7
[5] org.Mm.eg.db_3.6.0 GO.db_3.8.2
[7] OrganismDbi_1.26.0 GenomicFeatures_1.36.4
[9] GenomicRanges_1.36.0 GenomeInfoDb_1.20.0
[11] AnnotationDbi_1.46.0 IRanges_2.18.1
[13] S4Vectors_0.22.0 Biobase_2.44.0
[15] BiocGenerics_0.30.0 edgeR_3.26.8
[17] Glimma_1.12.0 limma_3.40.6
loaded via a namespace (and not attached):
[1] httr_1.4.1 bit64_0.9-7 jsonlite_1.6
[4] R.utils_2.9.0 gtools_3.8.1 assertthat_0.2.1
[7] BiocManager_1.30.4 RBGL_1.60.0 blob_1.2.0
[10] GenomeInfoDbData_1.2.1 Rsamtools_2.0.0 progress_1.2.2
[13] pillar_1.4.2 RSQLite_2.1.2 backports_1.1.4
[16] lattice_0.20-38 digest_0.6.20 XVector_0.24.0
[19] Matrix_1.2-17 R.oo_1.22.0 XML_3.98-1.20
[22] pkgconfig_2.0.2 biomaRt_2.40.4 zlibbioc_1.30.0
[25] gdata_2.18.0 BiocParallel_1.18.0 tibble_2.1.3
[28] SummarizedExperiment_1.14.0 magrittr_1.5 crayon_1.3.4
[31] memoise_1.1.0 R.methodsS3_1.7.1 graph_1.62.0
[34] tools_3.6.1 prettyunits_1.0.2 hms_0.5.1
[37] matrixStats_0.54.0 stringr_1.4.0 locfit_1.5-9.1
[40] DelayedArray_0.10.0 packrat_0.5.0 Biostrings_2.52.0
[43] compiler_3.6.1 caTools_1.17.1.2 rlang_0.4.0
[46] grid_3.6.1 RCurl_1.95-4.12 rstudioapi_0.10
[49] bitops_1.0-6 DBI_1.0.0 R6_2.4.0
[52] GenomicAlignments_1.20.1 rtracklayer_1.44.4 bit_1.1-14
[55] zeallot_0.1.0 KernSmooth_2.23-15 stringi_1.4.3
[58] Rcpp_1.0.2 vctrs_0.2.0
對(duì)比了一下和原教程使用的信息，發(fā)現(xiàn)有一點(diǎn)點(diǎn)不同，所以才會(huì)有的地方出錯(cuò)吧，用不同的電腦還是需要調(diào)試一下
參考文獻(xiàn)
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. “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.” 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.
