|
轉(zhuǎn)自 bioconductor論壇
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 2018年12月17日Contents1摘要簡(jiǎn)單且高效地分析RNA測(cè)序數(shù)據(jù)的能力正是Bioconductor的核心優(yōu)勢(shì),。 RNA-seq分析通常從基因水平的序列計(jì)數(shù)開(kāi)始,，涉及到數(shù)據(jù)預(yù)處理，探索性數(shù)據(jù)分析,，差異表達(dá)檢驗(yàn)以及通路分析,，得到的結(jié)果可用于指導(dǎo)進(jìn)一步實(shí)驗(yàn)和驗(yàn)證研究。 在這篇工作流程文章中,，我們通過(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)探索,，使用戶得以查看單個(gè)樣本與基因,。 這三個(gè)軟件包提供的完整分析突出了研究人員可以使用Bioconductor輕松地從RNA測(cè)序?qū)嶒?yàn)的原始計(jì)數(shù)揭示生物學(xué)意義。 2背景介紹RNA測(cè)序（RNA-seq）已成為基因表達(dá)研究的主要技術(shù),。其中,，基因組規(guī)模的多條件基因差異表達(dá)檢測(cè)是研究者最常探究的問(wèn)題之一。對(duì)于RNA-seq數(shù)據(jù),，來(lái)自Bioconductor項(xiàng)目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用于處理此問(wèn)題的完善的統(tǒng)計(jì)學(xué)方法,。 在這篇文章中，我們描述了一個(gè)用于分析RNA-seq數(shù)據(jù)的edgeR - limma工作流程,，使用基因水平計(jì)數(shù)作為輸入,，經(jīng)過(guò)預(yù)處理和探索性數(shù)據(jù)處理，然后得到差異表達(dá)（DE）基因和基因表達(dá)特征（gene signatures）的列表,。Glimma包(Su et al. 2017)提供的交互式圖表可以同時(shí)呈現(xiàn)整體樣本和單個(gè)基因水平的數(shù)據(jù)，使得我們相對(duì)靜態(tài)的R圖表而言,，可以探索更多的細(xì)節(jié),。 此工作流程中所分析的實(shí)驗(yàn)來(lái)自Sheridan等（2015）(Sheridan et al. 2015)，由三個(gè)細(xì)胞群組成,，即基底（basal）,、管腔祖細(xì)胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。細(xì)胞群皆分選自雌性處女小鼠的乳腺,，每種都設(shè)三個(gè)生物學(xué)重復(fù),。RNA樣品分三個(gè)批次使用Illumina HiSeq 2000進(jìn)行測(cè)序，得到長(zhǎng)為100堿基對(duì)的單端序列片段,。 本文所描述的分析假設(shè)從RNA-seq實(shí)驗(yàn)獲得的序列片段已經(jīng)與適當(dāng)?shù)膮⒖蓟蚪M比對(duì),，并已經(jīng)在基因水平上對(duì)序列進(jìn)行了統(tǒng)計(jì)計(jì)數(shù),。在本文條件下，使用Rsubread包提供的基于R的流程將序列片段與小鼠參考基因組（mm10）比對(duì)（具體而言,，先使用align函數(shù)(Liao, Smyth, and Shi 2013),，然后使用featureCounts (Liao, Smyth, and Shi 2014)進(jìn)行基因水平的總結(jié)，利用其內(nèi)置的mm10基于RefSeq的注釋）,。 這些樣本的計(jì)數(shù)數(shù)據(jù)可以從Gene Expression Omnibus (GEO)數(shù)據(jù)庫(kù) http://www.ncbi.nlm./geo/使用GEO序列登記號(hào)GSE63310下載,。更多關(guān)于實(shí)驗(yàn)設(shè)計(jì)和樣品制備的信息也可以在GEO使用該登記號(hào)查看。 3初始配置library(limma)library(Glimma)library(edgeR)library(Mus.musculus)
4數(shù)據(jù)整合4.1讀入計(jì)數(shù)數(shù)據(jù)為開(kāi)始此分析,，從https://www.ncbi.nlm./geo/download/?acc=GSE63310&format=file在線下載文件GSE63310_RAW.tar,，并從壓縮包中解壓出相關(guān)的文件。下方的代碼將完成此步驟,，或者您也可以手動(dòng)進(jìn)行這一步并繼續(xù)后續(xù)分析,。 url <- "https://www.ncbi.nlm./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)
## EntrezID GeneLength Count
## 1 497097 3634 1
## 2 100503874 3259 0
## 3 100038431 1634 0
## 4 19888 9747 0
## 5 20671 3130 1
盡管這九個(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)
## [1] "DGEList"
## attr(,"package")
## [1] "edgeR"
dim(x)
## [1] 27179 9
如果來(lái)自所有樣品的計(jì)數(shù)存儲(chǔ)在同一個(gè)文件中,，數(shù)據(jù)可以首先讀入R再使用DGEList函數(shù)轉(zhuǎn)換為一個(gè)DGEList對(duì)象,。 4.2組織樣品信息為進(jìn)行下游分析，與實(shí)驗(yàn)設(shè)計(jì)有關(guān)的樣品水平信息需要與計(jì)數(shù)矩陣的列關(guān)聯(lián),。這里需要包括各種對(duì)表達(dá)水平有影響的實(shí)驗(yàn)變量,，無(wú)論是生物變量還是技術(shù)變量。例如,，細(xì)胞類型（在這個(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ì)胞類型（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
## [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對(duì)象中的第二個(gè)數(shù)據(jù)框名為genes,，用于存儲(chǔ)與計(jì)數(shù)矩陣的行相關(guān)聯(lián)的基因水平的信息,。 為檢索這些信息，我們可以使用包含特定物種信息的包,，比如小鼠的Mus.musculus(Bioconductor Core Team 2016b)（或人類的Homo.sapiens (Bioconductor Core Team 2016a)）,；或者也可以使用biomaRt 包 (Durinck et al. 2005, 2009)，它通過(guò)接入Ensembl genome數(shù)據(jù)庫(kù)來(lái)進(jìn)行基因注釋,。 可以檢索的信息類型包括基因符號(hào)（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進(jìn)行檢索,，而由于Mus.musculus包含多種不同來(lái)源的信息,，它允許用戶從多種不同基因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)
## 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可能不能一對(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
## 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數(shù)據(jù)預(yù)處理5.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ì)算,，且已足以滿足我們所關(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ò)log2(CPM + 2/L)計(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值為log2(2/45.5) = -4.51,。換而言之,，在加上了預(yù)先計(jì)數(shù)彌補(bǔ)值后，此數(shù)據(jù)集中的零表達(dá)計(jì)數(shù)對(duì)應(yīng)的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.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函數(shù)時(shí)也會(huì)用到log-CPM值,，但voom會(huì)默認(rèn)使用更小的預(yù)先計(jì)數(shù)重新計(jì)算自己的log-CPM值。 5.2刪除低表達(dá)基因所有數(shù)據(jù)集中都混有表達(dá)的基因與不表達(dá)的基因,。盡管我們想要檢測(cè)在一種條件中表達(dá)但再另一種條件中不表達(dá)的基因,，也有一些基因在所有樣品中都不表達(dá)。實(shí)際上,，這個(gè)數(shù)據(jù)集中19%的基因在所有九個(gè)樣品中的計(jì)數(shù)都是零,。 table(rowSums(x$counts==0)==9)
##
## FALSE TRUE
## 22026 5153
對(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)
## [1] 16624 9
此函數(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ì)胞類型組內(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")
Figure 1: 每個(gè)樣本過(guò)濾前的原始數(shù)據(jù)（A）和過(guò)濾后（B）的數(shù)據(jù)的log-CPM值密度,。豎直虛線標(biāo)出了過(guò)濾步驟中所用閾值（相當(dāng)于CPM值為約0
2）。 5.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
## [1] 0.894 1.025 1.046 1.046 1.016 0.922 0.996 1.086 0.984
為了更好地可視化表現(xiàn)出歸一化的影響，我們復(fù)制了數(shù)據(jù)并進(jìn)行了調(diào)整,，使得第一個(gè)樣品的計(jì)數(shù)減少到了其原始值的5%,，而第二個(gè)樣品增大到了5倍。 x2 <- x
x2$samples$norm.factors <- 1x2$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
## [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: 樣例數(shù)據(jù)：log-CPM值的箱線圖展示了未經(jīng)歸一化的數(shù)據(jù)（A）及歸一化后的數(shù)據(jù)（B）中每個(gè)樣本的表達(dá)分布,。數(shù)據(jù)集經(jīng)過(guò)調(diào)整，樣本1和2中的表達(dá)計(jì)數(shù)分別被縮放到其原始值的5%和500%,。 5.4對(duì)樣本的無(wú)監(jiān)督聚類在我們看來(lái),，用于檢查基因表達(dá)分析的最重要的探索性圖表之一便是MDS圖或其余類似的圖。這種圖表使用無(wú)監(jiān)督聚類方法展示出了樣品間的相似性和不相似性,，能讓我們?cè)谶M(jìn)行正式的檢驗(yàn)之前對(duì)于能檢測(cè)到多少差異表達(dá)基因有個(gè)大致概念,。理想情況下，樣本會(huì)在不同的實(shí)驗(yàn)組內(nè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è)因子,。如果在其中一些維度上樣本可按照某因子聚類，這說(shuō)明該因子對(duì)于表達(dá)差異有影響,，在線性模型中應(yīng)當(dāng)將其包括進(jìn)去,。反之，沒(méi)有或者僅有微小影響的因子在下游分析時(shí)應(yīng)當(dāng)被剔除,。 在這個(gè)數(shù)據(jù)集中,，可以看出樣本在維度1和2能很好地按照實(shí)驗(yàn)分組聚類,，隨后在維度3按照測(cè)序道（樣品批次）分離（如下圖所示）。請(qǐng)記住,，第一維度解釋了數(shù)據(jù)中最大比例的方差,，需要注意到，當(dāng)我們向高維度移動(dòng),，維度上的取值范圍會(huì)變小,。 盡管所有樣本都按組聚類，在維度1上最大的轉(zhuǎn)錄差異出現(xiàn)在basal和LP以及basal和ML之間,。因此，預(yù)期在basal樣品與其他之間的成對(duì)比較中能夠得到大量的DE基因,，而在ML和LP之間的比較中得到的DE基因數(shù)量略少,。在其他的數(shù)據(jù)集中，不按照實(shí)驗(yàn)組聚類的樣本可能在下游分析中只表現(xiàn)出較小的或不表現(xiàn)出差異表達(dá),。 為繪制MDS圖,，我們?yōu)椴煌囊蜃淤x予不同的色彩組合。維度1和2使用以細(xì)胞類型定義的色彩組合進(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")
Figure 3: log-CPM值在維度1和2的MDS圖，以樣品分組上色并標(biāo)記（A）和維度3和4的MDS圖,，以測(cè)序道上色并標(biāo)記（B）,。圖中的距離對(duì)應(yīng)于最主要的倍數(shù)變化（fold change），默認(rèn)情況下也就是前500個(gè)在每對(duì)樣品之間差異最大的基因的平均（均方根）log2倍數(shù)變化,。 作為另一種選擇,，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ì)胞類型或測(cè)序泳道（批次）,。此數(shù)據(jù)集的交互式MDS圖可以從http://bioinf./folders/limmaWorkflow/glimma-plots/MDS-Plot.html看到。 glMDSPlot(lcpm, labels=paste(group, lane, sep="_"),
groups=x$samples[,c(2,5)], launch=FALSE)
交互式MDS圖鏈接 6差異表達(dá)分析6.1創(chuàng)建設(shè)計(jì)矩陣和對(duì)比在此研究中,，我們想知道哪些基因在我們研究的三組細(xì)胞之間以不同水平表達(dá),。在我們的分析中，假設(shè)基礎(chǔ)數(shù)據(jù)是正態(tài)分布的,，為其擬合一個(gè)線性模型,。在此之前,，需要?jiǎng)?chuàng)建一個(gè)包含細(xì)胞類型以及測(cè)序泳道（批次）信息的設(shè)計(jì)矩陣。 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"
對(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)有group的截距的情況下能更直截了當(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
## Contrasts
## Levels BasalvsLP BasalvsML LPvsML
## Basal 1 1 0
## LP -1 0 1
## ML 0 -1 -1
## laneL006 0 0 0
## laneL008 0 0 0
limma的線性模型方法的核心優(yōu)勢(shì)之一便是其適應(yīng)任意實(shí)驗(yàn)復(fù)雜程度的能力。簡(jiǎn)單的設(shè)計(jì),，比如此工作流程中關(guān)于細(xì)胞類型和批次的實(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ù),。 6.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ì)胞類型的重復(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
## 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: 圖中繪制了每個(gè)基因的均值（x軸）和方差（y軸）,，顯示了在該數(shù)據(jù)上使用voom前它們之間的相關(guān)性（左），以及當(dāng)運(yùn)用voom的精確權(quán)重后這種趨勢(shì)是如何消除的（右）,。左側(cè)的圖是使用voom函數(shù)繪制的，它為進(jìn)行l(wèi)og-CPM轉(zhuǎn)換后的數(shù)據(jù)擬合線性模型從而提取殘差方差,。然后,，對(duì)方差取平方根（或?qū)?biāo)準(zhǔn)差取平方根），并相對(duì)每個(gè)基因的平均表達(dá)作圖,。均值通過(guò)平均計(jì)數(shù)加上2再進(jìn)行l(wèi)og2轉(zhuǎn)換計(jì)算得到,。右側(cè)的圖使用plotSA繪制了log2殘差標(biāo)準(zhǔn)差與log-CPM均值的關(guān)系。平均log2殘差標(biāo)準(zhǔn)差由水平藍(lán)線標(biāo)出,。在這兩幅圖中,，每個(gè)黑點(diǎn)表示一個(gè)基因，紅線為對(duì)這些點(diǎn)的擬合,。 值得注意的是,，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á)值類似于x$counts,，盡管它進(jìn)行了尺度轉(zhuǎn)換。此外,，voom的EList對(duì)象中還有一個(gè)精確權(quán)重的矩陣v$weights,，而設(shè)計(jì)矩陣存儲(chǔ)于v$design。 6.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),。 6.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ì)胞類型的比較中皆找到了大量的DE基因,，這與我們?cè)贛DS圖中觀察到的結(jié)果相吻合,。 summary(decideTests(efit))
## BasalvsLP BasalvsML LPvsML
## Down 4648 4927 3135
## NotSig 7113 7026 10972
## Up 4863 4671 2517
一些研究中不僅僅需要使用校正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ì)胞類型之間差兩倍）時(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)
## BasalvsLP BasalvsML LPvsML
## Down 1632 1777 224
## NotSig 12976 12790 16210
## Up 2016 2057 190
在多種比較中皆差異表達(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)
## [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（右）的對(duì)比的DE基因數(shù)量,，還有兩種對(duì)比中共同的DE基因數(shù)量（中）。在任何對(duì)比中均不差異表達(dá)的基因數(shù)量標(biāo)于右下,。 write.fit(tfit, dt, file="results.txt")
6.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ù)量可由用戶指定，如果設(shè)為n=Inf則會(huì)包括所有的基因,?；?em style="box-sizing: border-box;">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差異表達(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的輸出類似）,，右側(cè)面板包含各個(gè)樣本的log-CPM值，下方為結(jié)果的表格,。 這種交互式展示允許用戶使用提供的注釋（比如基因名標(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生成的均值-差異圖,。左側(cè)面板顯示了總結(jié)性數(shù)據(jù)（log-FC與log-CPM值的關(guān)系），其中選中的基因在每個(gè)樣本中的數(shù)值顯示于右側(cè)面板,。下方為結(jié)果的表格,，其搜索框使用戶得以使用可行的注釋信息查找某個(gè)特定基因，如基因符號(hào)Clu,。 上方指令生成的均值-差異圖可以在線訪問(wèn)（詳見(jiàn)http://bioinf./folders/limmaWorkflow/glimma-plots/MD-Plot.html）,。 Glimma提供的交互性使得單個(gè)圖形窗口內(nèi)能夠呈現(xiàn)出額外的信息。 Glimma是以R和Javascript實(shí)現(xiàn)的,，使用R代碼生成數(shù)據(jù),，并在之后使用Javascript庫(kù)D3（https://）轉(zhuǎn)換為圖形，使用Bootstrap庫(kù)處理界面并生成互動(dòng)性可搜索的表格的數(shù)據(jù)表,。 這使得圖表可以在任何現(xiàn)代的瀏覽器中查看,，對(duì)于從Rmarkdown分析報(bào)告中將其作為關(guān)聯(lián)文件而附加而言十分方便。 前文所展示的圖表中,，一些展示了在任意一個(gè)條件下表達(dá)的所有基因（比如共同DE基因的韋恩圖或均值-差異圖）,，而另一些展示單獨(dú)的基因（交互性均值-差異圖右邊面板中所展示的log-CPM值）,。而熱圖使用戶得以查看一部分基因的表達(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ì)胞類型聚類，并重新排序了基因,，形成了表達(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")
Figure 6: 在basal和LP的對(duì)比中前100個(gè)DE基因log-CPM值的熱圖。經(jīng)過(guò)縮放調(diào)整后,，每個(gè)基因（每行）的表達(dá)均值為0,，并且標(biāo)準(zhǔn)差為1。給定基因相對(duì)高表達(dá)的樣本被標(biāo)記為紅色,，相對(duì)低表達(dá)的樣本被標(biāo)記為藍(lán)色,。淺色和白色代表中等表達(dá)水平的基因。樣本和基因已通過(guò)分層聚類的方法重新排序,。圖中顯示有樣本聚類的樹(shù)狀圖,。 7使用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./software/MSigDB/以RData對(duì)象格式獲取,。 此外,，對(duì)于人類和小鼠，來(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ò)程。 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-05
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ì)比一致）,。 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（藍(lán)色條形，圖表下方）基因集在LP和ML的對(duì)比中的條碼圖,，每個(gè)基因集都有一條富集線展示了豎直條形在圖表每部分的相對(duì)富集程度,。Lim等人的實(shí)驗(yàn)(Lim et al
2010)非常類似于我們的，用了相同的分選方式來(lái)獲取不同的細(xì)胞群,，只是他們使用的是微陣列而不是RNA-seq來(lái)測(cè)定基因表達(dá),。需要注意的是，上調(diào)基因集發(fā)生下調(diào)而下調(diào)基因集發(fā)生上調(diào)的逆相關(guān)性來(lái)自于對(duì)比的設(shè)定方式（LP相比于ML）,，如果將其對(duì)調(diào),，方向性將會(huì)吻合。 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)確定有意義的基因集的顯著性,。 8使用到的軟件和代碼此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:///packages/RNAseq123查看）內(nèi)包含此文章的英文和簡(jiǎn)體中文版以及進(jìn)行整個(gè)分析流程所需要的代碼,。安裝此包即可管理以上提到的所有需要的包,。對(duì)于RNA-seq數(shù)據(jù)分析實(shí)踐培訓(xùn)而言，此包也是非常有用的資源,。 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.6
參考文獻(xiàn)Bioconductor Core Team. 2016a. Homo.sapiens: Annotation package for the Homo.sapiens object. https:///packages/release/data/annotation/html/Homo.sapiens.html. ———. 2016b. Mus.musculus: Annotation package for the Mus.musculus object. https:///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./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. |
