原文來自這里厢拭，略編輯
摘要
簡單且高效地分析RNA測序數(shù)據(jù)的能力是Bioconductor的核心優(yōu)勢它匕。 RNA-seq分析通常從基因水平的序列計(jì)數(shù)開始吩跋，涉及到數(shù)據(jù)預(yù)處理寞射，探索性數(shù)據(jù)分析，差異表達(dá)檢驗(yàn)以及通路分析钞澳，得到的結(jié)果可用于指導(dǎo)進(jìn)一步實(shí)驗(yàn)和驗(yàn)證研究怠惶。 在這篇工作流程文章中，我們通過分析來自小鼠乳腺的RNA測序數(shù)據(jù)轧粟，示范了如何使用流行的edgeR包載入策治、整理、過濾和歸一化數(shù)據(jù)兰吟，然后用limma包的voom方法通惫、線性模型和經(jīng)驗(yàn)貝葉斯調(diào)節(jié)（empirical Bayes moderation）來評(píng)估差異表達(dá)并進(jìn)行基因集檢驗(yàn)。通過使用Glimma包混蔼，此流程得到了增進(jìn)履腋，實(shí)現(xiàn)了結(jié)果的互動(dòng)探索，使用戶得以查看單個(gè)樣本與基因惭嚣。 這三個(gè)軟件包提供的完整分析突出了研究人員可以使用Bioconductor輕松地從RNA測序?qū)嶒?yàn)的原始計(jì)數(shù)揭示生物學(xué)意義遵湖。
2背景介紹
RNA測序（RNA-seq）已成為基因表達(dá)研究的主要技術(shù)。其中晚吞，基因組規(guī)模的多條件基因差異表達(dá)檢測是研究者最常探究的問題之一延旧。對(duì)于RNA-seq數(shù)據(jù)，來自Bioconductor項(xiàng)目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用于處理此問題的完善的統(tǒng)計(jì)學(xué)方法槽地。
在這篇文章中迁沫，我們描述了一個(gè)用于分析RNA-seq數(shù)據(jù)的edgeR - limma工作流程芦瘾，使用基因水平計(jì)數(shù)作為輸入，經(jīng)過預(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)來自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)行測序句占，得到長為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ù)庫 http://www.ncbi.nlm.nih.gov/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ù)
為開始此分析，從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)
## 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
如果來自所有樣品的計(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)變量，無論是生物變量還是技術(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）和批次（測序泳道lane）信息籽慢，每種信息都包含三個(gè)不同的水平浸遗。需要注意的是，在x$samples中箱亿，程序會(huì)自動(dòng)計(jì)算每個(gè)樣品的文庫大小跛锌，歸一化系數(shù)會(huì)被設(shè)置為1。 為了簡單起見届惋，我們從我們的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)译秦，它通過接入Ensembl genome數(shù)據(jù)庫來進(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包含多種不同來源的信息械拍，它允許用戶從多種不同基因ID中選取檢索鍵突勇。
我們使用Mus.musculus包装盯，利用我們數(shù)據(jù)集中的Entrez基因ID來檢索相關(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ù)的來源非常重要定躏。我們的基因注釋中包含28個(gè)匹配到不同染色體的基因（比如基因Gm1987關(guān)聯(lián)于染色體chr4和chr4_JH584294_random账磺，小RNA Mir5098關(guān)聯(lián)于chr2，chr5痊远，chr8垮抗，chr11和chr17）。 為了處理重復(fù)的基因ID碧聪，我們可以合并來自多重匹配基因的所有染色體信息冒版，比如將基因Gm1987分配到chr4 and chr4_JH584294_random，或選取其中一條染色體來代表具有重復(fù)注釋的基因逞姿。為了簡單起見壤玫，我們選擇后者，保留每個(gè)基因ID第一次出現(xiàn)的信息哼凯。
genes <- genes[!duplicated(genes$ENTREZID),]
在此例子中，注釋與數(shù)據(jù)對(duì)象中的基因順序是相同的楚里。如果由于缺失和／或重新排列基因ID導(dǎo)致其順序不一致断部，可以用match來正確排序基因。然后將基因注釋的數(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)換
由于更深的測序總會(huì)產(chǎn)生更多的序列达址，在差異表達(dá)相關(guān)的分析中蔑祟，我們很少使用原始的序列數(shù)。在實(shí)踐中沉唠，我們通常將原始的序列數(shù)進(jìn)行歸一化疆虚，來消除測序深度所導(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那樣考慮到基因長度區(qū)別的影響篇亭，但在我們的分析中經(jīng)常會(huì)用到它們。雖然也可以使用RPKM和FPKM锄贷，但CPM和log-CPM只使用計(jì)數(shù)矩陣即可計(jì)算译蒂，且已足以滿足我們所關(guān)注的比較的需要曼月。假設(shè)不同條件之間剪接異構(gòu)體（isoform）的使用沒有差別，差異表達(dá)分析研究同一基因在不同條件下的表達(dá)差異柔昼，而不是比較多個(gè)基因之間的表達(dá)或測定絕對(duì)表達(dá)量哑芹。換而言之，基因長度在我們關(guān)注的比較中始終不變岳锁，且任何觀測到的差異是來自于條件的變化而不是基因長度的變化绩衷。
在此處，使用edgeR中的cpm函數(shù)將原始計(jì)數(shù)轉(zhuǎn)換為CPM和log-CPM值激率。如果可以提供基因長度信息咳燕，可以使用edgeR中的rpkm函數(shù)計(jì)算RPKM值，就像計(jì)算CPM值那樣簡單乒躺。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)
對(duì)于一個(gè)基因招盲，CPM值為1相當(dāng)于在測序深度最低的樣品中（JMS9-P8c, 序列數(shù)量約2千萬）有20個(gè)計(jì)數(shù)，或者在測序深度最高的樣品中（JMS8-3嘉冒，序列數(shù)量約7.6千萬）有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ù)（以百萬計(jì)）的平均值礼饱，所以log-CPM值是根據(jù)CPM值通過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ù)變化帶來的巨大的偽差異性，這對(duì)于繪制探索性圖表很有用洒忧。在這個(gè)數(shù)據(jù)集中蝴韭，平均的樣本總序列數(shù)是4.55千萬，所以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-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函數(shù)時(shí)也會(huì)用到log-CPM值，但voom會(huì)默認(rèn)使用更小的預(yù)先計(jì)數(shù)重新計(jì)算自己的log-CPM值芝雪。
5.2刪除低表達(dá)基因
所有數(shù)據(jù)集中都混有表達(dá)的基因與不表達(dá)的基因减余。盡管我們想要檢測在一種條件中表達(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部分）抒抬。
在任何樣本中都沒有足夠多的序列片段的基因應(yīng)該從下游分析中過濾掉杨刨。這樣做的原因有好幾個(gè)。 從生物學(xué)的角度來看擦剑，在任何條件下的表達(dá)水平都不具有生物學(xué)意義的基因都不值得關(guān)注妖胀，因此最好忽略。 從統(tǒng)計(jì)學(xué)的角度來看惠勒，去除低表達(dá)計(jì)數(shù)基因使數(shù)據(jù)中的均值 - 方差關(guān)系可以得到更精確的估計(jì)赚抡，并且還減少了在觀察差異表達(dá)的下游分析中需要進(jìn)行的統(tǒng)計(jì)檢驗(yàn)的數(shù)量。
edgeR包中的filterByExpr函數(shù)提供了自動(dòng)過濾基因的方法纠屋，可保留盡可能多的有足夠表達(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)行過濾時(shí)使用CPM值而不是表達(dá)計(jì)數(shù)赁遗，以避免對(duì)總序列數(shù)大的樣本的偏向性。在這個(gè)數(shù)據(jù)集中族铆，總序列數(shù)的中位數(shù)是5.1千萬岩四，且10/51約等于0.2，所以filterByExpr函數(shù)保留在至少三個(gè)樣本中CPM值大于等于0.2的基因哥攘。此處剖煌，一個(gè)具有生物學(xué)意義的基因需要在至少三個(gè)樣本中表達(dá)，因?yàn)槿N細(xì)胞類型組內(nèi)各有三個(gè)重復(fù)献丑。所使用的閾值取決于測序深度和實(shí)驗(yàn)設(shè)計(jì)。如果樣本總表達(dá)計(jì)數(shù)量增大侠姑，那么可以選擇更低的CPM閾值创橄，因?yàn)楦蟮目偙磉_(dá)計(jì)數(shù)量提供了更好的分辨率來探究更多表達(dá)水平更低的基因。
使用這個(gè)標(biāo)準(zhǔn)莽红，基因的數(shù)量減少到了16624個(gè)妥畏，約為開始時(shí)數(shù)量的60%。過濾后的log-CPM值顯示出每個(gè)樣本的分布基本相同（下圖B部分）安吁。需要注意的是醉蚁，從整個(gè)DGEList對(duì)象中取子集時(shí)同時(shí)刪除了被過濾的基因的計(jì)數(shù)和其相關(guān)的基因信息。過濾后的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")
5.3歸一化基因表達(dá)分布
在樣品制備或測序過程中滥玷，不具備生物學(xué)意義的外部因素會(huì)影響單個(gè)樣品的表達(dá)氏身。比如說，在實(shí)驗(yàn)中第一批制備的樣品會(huì)總體上表達(dá)高于第二批制備的樣品惑畴。假設(shè)所有樣品表達(dá)值的范圍和分布都應(yīng)當(dāng)相似蛋欣，需要進(jìn)行歸一化來確保整個(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ù)被用作文庫大小的縮放系數(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 <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5
下圖顯示了沒有經(jīng)過歸一化的與經(jīng)過了歸一化的數(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")
5.4對(duì)樣本的無監(jiān)督聚類
在我們看來，用于檢查基因表達(dá)分析的最重要的探索性圖表之一便是MDS圖或其余類似的圖绢陌。這種圖表使用無監(jiān)督聚類方法展示出了樣品間的相似性和不相似性挨下，能讓我們在進(jìn)行正式的檢驗(yàn)之前對(duì)于能檢測到多少差異表達(dá)基因有個(gè)大致概念。理想情況下脐湾，樣本會(huì)在不同的實(shí)驗(yàn)組內(nèi)很好的聚類臭笆，且可以鑒別出遠(yuǎn)離所屬組的樣本，并追蹤誤差或額外方差的來源。如果存在技術(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è)因子获黔。如果在其中一些維度上樣本可按照某因子聚類，這說明該因子對(duì)于表達(dá)差異有影響在验，在線性模型中應(yīng)當(dāng)將其包括進(jìn)去玷氏。反之，沒有或者僅有微小影響的因子在下游分析時(shí)應(yīng)當(dāng)被剔除腋舌。
在這個(gè)數(shù)據(jù)集中盏触，可以看出樣本在維度1和2能很好地按照實(shí)驗(yàn)分組聚類，隨后在維度3按照測序道（樣品批次）分離（如下圖所示）块饺。請(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使用以測序泳道（批次）定義的色彩組合進(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)頁（如果設(shè)置launch=TRUE玉雾，將會(huì)在瀏覽器中打開）翔试，其左側(cè)面板含有一張MDS圖，而右側(cè)面板包含一張展示了各個(gè)維度所解釋的方差比例的柱形圖复旬。點(diǎn)擊柱形圖中的柱可切換MDS圖繪制時(shí)所使用的維度垦缅，且將鼠標(biāo)懸浮于單個(gè)點(diǎn)上可顯示相應(yīng)的樣本標(biāo)簽。也可切換配色方案驹碍，以突顯不同細(xì)胞類型或測序泳道（批次）壁涎。此數(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)
6差異表達(dá)分析
6.1創(chuàng)建設(shè)計(jì)矩陣和對(duì)比
在此研究中志秃，我們想知道哪些基因在我們研究的三組細(xì)胞之間以不同水平表達(dá)怔球。在我們的分析中，假設(shè)基礎(chǔ)數(shù)據(jù)是正態(tài)分布的浮还，為其擬合一個(gè)線性模型竟坛。在此之前，需要?jiǎng)?chuàng)建一個(gè)包含細(xì)胞類型以及測序泳道（批次）信息的設(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ì)矩陣。 比如說延刘，~0+group+lane去除了第一個(gè)因子group的截距漫试，但第二個(gè)因子lane的截距被保留。 此外也可以使用~group+lane碘赖，來自group和lane的截距均被保留驾荣。 此處的關(guān)鍵是理解如何解釋給定模型中估計(jì)得到的系數(shù)。 我們在此分析中選取第一種模型普泡，因?yàn)樵跊]有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)勢之一便是其適應(yīng)任意實(shí)驗(yàn)復(fù)雜程度的能力撼班。簡單的設(shè)計(jì)歧匈，比如此工作流程中關(guān)于細(xì)胞類型和批次的實(shí)驗(yàn)設(shè)計(jì)，直到更復(fù)雜的因子設(shè)計(jì)和含有交互作用項(xiàng)的模型砰嘁，都能夠被相對(duì)簡單地處理件炉。當(dāng)實(shí)驗(yàn)或技術(shù)效應(yīng)可被隨機(jī)效應(yīng)模型（random effect model）模擬時(shí)，limma中的另一種可能性是使用duplicateCorrelation函數(shù)來估計(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)分布來模擬計(jì)數(shù)的方法假設(shè)均值與方差間具有二次的關(guān)系。在limma中，假設(shè)log-CPM值符合正態(tài)分布秀撇，并使用由voom函數(shù)計(jì)算得到的精確權(quán)重來調(diào)整均值與方差的關(guān)系超棺，從而對(duì)log-CPM值進(jìn)行線性建模。
當(dāng)操作DGEList對(duì)象時(shí)呵燕，voom從x中自動(dòng)提取文庫大小和歸一化因子棠绘，以此將原始計(jì)數(shù)轉(zhuǎn)換為log-CPM值。在voom中再扭，對(duì)于log-CPM值額外的歸一化可以通過設(shè)定normalize.method參數(shù)來進(jìn)行弄唧。
下圖左側(cè)展示了這個(gè)數(shù)據(jù)集log-CPM值的均值-方差關(guān)系。通常而言霍衫，方差是測序?qū)嶒?yàn)中的技術(shù)差異和不同細(xì)胞類型的重復(fù)樣本之間的生物學(xué)差異的結(jié)合候引，而voom圖會(huì)顯示出一個(gè)在均值與方差之間遞減的趨勢。 生物學(xué)差異高的實(shí)驗(yàn)通常會(huì)有更平坦的趨勢敦跌，其方差值在高表達(dá)處穩(wěn)定澄干。 生物學(xué)差異低的實(shí)驗(yàn)更傾向于急劇下降的趨勢。
不僅如此柠傍，voom圖也提供了對(duì)于上游所進(jìn)行的過濾水平的可視化檢測麸俘。如果對(duì)于低表達(dá)基因的過濾不夠充分，在圖上表達(dá)低的一端惧笛，受到非常低的表達(dá)計(jì)數(shù)的影響从媚，可以觀察到方差水平的下降。如果觀察到了這種情況患整，應(yīng)當(dāng)回到最初的過濾步驟并提高用于該數(shù)據(jù)集的表達(dá)閾值拜效。
當(dāng)前面觀察的MDS圖中具有明顯的樣本水平的差異時(shí)者蠕，可以用voomWithQualityWeights函數(shù)來同時(shí)合并樣本水平的權(quán)重和voom(Liu et al. 2015)估算得到的豐度相關(guān)的權(quán)重赃梧。關(guā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")
值得注意的是鸠窗，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)行株婴，它們原先是為微陣列而設(shè)計(jì)的怎虫。這些函數(shù)不僅可以用于微陣列數(shù)據(jù)，也可以用于RNA-seq數(shù)據(jù)困介。它們會(huì)單獨(dú)為每個(gè)基因的表達(dá)值擬合一個(gè)模型大审。然后，通過利用所有基因的信息來進(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基因谅将，這與我們在MDS圖中觀察到的結(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)過經(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é)果提取并寫入到單個(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"))
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)過多重假設(shè)檢驗(yàn)校正的p值。列出前多少個(gè)基因的數(shù)量可由用戶指定特碳，如果設(shè)為n=Inf則會(huì)包括所有的基因诚亚。基因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頁面敏释，左側(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)
上方指令生成的均值-差異圖可以在線訪問（詳見http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MD-Plot.html）财喳。 Glimma提供的交互性使得單個(gè)圖形窗口內(nèi)能夠呈現(xiàn)出額外的信息。 Glimma是以R和Javascript實(shí)現(xiàn)的斩狱，使用R代碼生成數(shù)據(jù)耳高，并在之后使用Javascript庫D3（https://d3js.org）轉(zhuǎn)換為圖形，使用Bootstrap庫處理界面并生成互動(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")
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.wehi.edu.au/software/MSigDB/以RData對(duì)象格式獲取说庭。 此外然磷，對(duì)于人類和小鼠，來自MSigDB的其他有用的基因集也可從此網(wǎng)站獲取刊驴，比如標(biāo)志（hallmark）基因集姿搜。C2基因集的內(nèi)容收集自在線數(shù)據(jù)庫、出版物以及該領(lǐng)域?qū)＜依υ鳎鴺?biāo)志基因集的內(nèi)容來自MSigDB舅柜，從而獲得具有明確定義的生物狀態(tài)或過程。
load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123"))
idx <- ids2indices(Mm.c2,id=rownames(v))
cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1])
head(cam.BasalvsLP,5)
## NGenes Direction PValue FDR
## LIM_MAMMARY_STEM_CELL_UP 791 Up 1.77e-18 8.36e-15
## LIM_MAMMARY_STEM_CELL_DN 683 Down 4.03e-14 8.69e-11
## ROSTY_CERVICAL_CANCER_PROLIFERATION_CLUSTER 170 Up 5.52e-14 8.69e-11
## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 94 Down 2.74e-13 3.23e-10
## SOTIRIOU_BREAST_CANCER_GRADE_1_VS_3_UP 190 Up 5.16e-13 4.87e-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ù)通過比較假設(shè)檢驗(yàn)來評(píng)估一個(gè)給定基因集中的基因是否相對(duì)于不在集內(nèi)的基因而言在差異表達(dá)基因的排序中更靠前躲惰。 它使用limma的線性模型框架致份，并同時(shí)采用設(shè)計(jì)矩陣和對(duì)比矩陣（如果有的話），且在測試的過程中會(huì)使用來自voom的觀測水平權(quán)重础拨。 在通過基因間相關(guān)性（默認(rèn)設(shè)定為0.01氮块，但也可通過數(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)長im等人（2010）的成熟管腔基因集（上調(diào)及下調(diào)）繪制了條碼圖（barcodeplot）芳悲。需要注意的是立肘，由于我們的對(duì)比是將LP與ML相比而不是相反，這些基因集的方向在我們的數(shù)據(jù)集中是反過來的（如果將對(duì)比反過來名扛，基因集的方向?qū)?huì)與對(duì)比一致）谅年。
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP,
index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
limma也有其他的基因集檢驗(yàn)届良，比如mroast(Wu et al. 2010)的自包含檢驗(yàn)笆凌。雖然camera非常適合檢驗(yàn)基因集的大型數(shù)據(jù)庫并觀察其中哪些相對(duì)于其他的在排序上位次更高（如前文所示），自包含檢驗(yàn)更善于集中檢驗(yàn)一個(gè)或少個(gè)選中的集合是否本身差異表達(dá)士葫。換句話說乞而，camera更適用于搜尋具有意義的基因集，而mroast測試的是已經(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（可訪問https://bioconductor.org/packages/RNAseq123查看）內(nèi)包含此文章的英文和簡體中文版以及進(jìn)行整個(gè)分析流程所需要的代碼侦香。安裝此包即可管理以上提到的所有需要的包。對(duì)于RNA-seq數(shù)據(jù)分析實(shí)踐培訓(xùn)而言纽疟，此包也是非常有用的資源罐韩。
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
參考文獻(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.
