1 摘要
簡單且高效地分析RNA測序數(shù)據(jù)的能力正是Bioconductor的核心優(yōu)勢之一。在獲得RNA-seq基因表達矩陣后竭缝，通常需要對數(shù)據(jù)進行預處理膜蛔、探索性數(shù)據(jù)分析呜呐、差異表達檢驗以及通路分析颠通，以得到可以幫助進一步實驗和驗證研究的結果顿锰。在本工作流程中谨垃，我們將通過分析來自小鼠乳腺的RNA測序數(shù)據(jù)，演示如何使用流行的edgeR包載入硼控、整理刘陶、過濾和歸一化數(shù)據(jù)，然后用limma包的voom方法牢撼、線性模型和經(jīng)驗貝葉斯調節(jié)來評估差異表達并進行基因集檢驗匙隔。通過Glimma包，本流程進一步實現(xiàn)了結果的互動探索熏版，便于用戶查看特定樣本與基因的分析結果纷责。通過使用這三個Bioconductor包，研究者可以輕松地運行完整的RNA-seq數(shù)據(jù)分析流程撼短，從原始計數(shù)（raw counts）中挖掘出其中蘊含的生物學意義再膳。
2 背景介紹
RNA測序（RNA-seq）是用于研究基因表達的重要技術。其中曲横，在基因組規(guī)模下檢測多條件之間基因的差異表達是研究者最常探究的問題之一喂柒。對于RNA-seq數(shù)據(jù)，來自Bioconductor項目(Huber et al. 2015)的edgeR?(Robinson, McCarthy, and Smyth 2010)和limma包(Ritchie et al. 2015)提供了一套用于處理此問題的完善的統(tǒng)計學方法。
在這篇文章中胳喷，我們描述了一個用于分析RNA-seq數(shù)據(jù)的edgeR?-?limma工作流程，使用基因水平的計數(shù)（gene-level counts）作為輸入夭织，經(jīng)過預處理和探索性數(shù)據(jù)分析吭露，然后得到差異表達（DE）基因和基因表達特征（gene signatures）的列表。Glimma包(Su et al. 2017)提供的交互式圖表可以同時呈現(xiàn)整體樣本層面與單個基因層面的數(shù)據(jù)尊惰，相對靜態(tài)的R圖表而言讲竿，更便于我們探索更多的細節(jié)。
此工作流程中我們分析的數(shù)據(jù)來自Sheridan等人的實驗（2015）(Sheridan et al. 2015)弄屡，它包含三個細胞群题禀，即基底（basal）、管腔祖細胞（luminal progenitor, LP）和成熟管腔（mature luminal, ML）膀捷。細胞群皆分選自雌性處女小鼠的乳腺迈嘹，每種都設三個生物學重復。RNA樣品分三個批次使用Illumina HiSeq 2000進行測序全庸，得到長為100堿基對的單端序列片段秀仲。
本文所述的分析流程假設從RNA-seq實驗獲得的序列片段已經(jīng)與適當?shù)膮⒖蓟蚪M比對，并已經(jīng)在基因水平上對序列進行了統(tǒng)計計數(shù)壶笼。在本文條件下神僵，使用Rsubread包提供的基于R的流程將序列片段與小鼠參考基因組（mm10）比對（具體而言，先使用align函數(shù)(Liao, Smyth, and Shi 2013)進行比對覆劈，然后使用featureCounts?(Liao, Smyth, and Shi 2014)函數(shù)保礼，利用其內置的基于RefSeq的mm10注釋進行基因水平的總結）。
這些樣本的計數(shù)數(shù)據(jù)可以從Gene Expression Omnibus (GEO)數(shù)據(jù)庫http://www.ncbi.nlm.nih.gov/geo/使用GEO序列登記號;. 下載责语。更多關于實驗設計和樣品制備的信息也可以在GEO使用該登記號查看炮障。
3 初始配置
> if (!requireNamespace("BiocManager")) {
+ ????install.packages("BiocManager")
+ }
> BiocManager::install("edgeR")以此類推Mus.musculus等包?
library(limma)
library(Glimma)
library(edgeR)
library(Mus.musculus)
4 數(shù)據(jù)整合
4.1讀入計數(shù)數(shù)據(jù)
為開始此分析，從https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file在線下載文件GSE63310_RAW.tar鹦筹，并從壓縮包中解壓出相關的文件铝阐。下方的代碼將完成此步驟，或者您也可以手動進行這一步并繼續(xù)后續(xù)分析铐拐。
url <-?"https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir =?".")
files <-?c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
??"GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
??"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in?paste(files, ".gz", sep=""))
??R.utils::gunzip(i, overwrite=TRUE)#gunzip解壓縮
每一個文本文件均包含一個給定樣品的原始基因水平計數(shù)徘键。需要注意的是，我們的分析僅包含了此實驗中的basal遍蟋、LP和ML樣品（請查看下方相關文件名）吹害。
每一個文本文件均為對應樣品的原始基因水平計數(shù)矩陣。需要注意我們的這次分析僅包含了此實驗中的basal虚青、LP和ML樣品（可見下方所示文件名）它呀。
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
???"GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
???"GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
???"GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
???"GSM1545545_JMS9-P8c.txt")
read.delim(files[1], nrow=5)
## ???EntrezID GeneLength Count
## 1 ???497097 ??????3634 ????1
## 2 100503874 ??????3259 ????0
## 3 100038431 ??????1634 ????0
## 4 ????19888 ??????9747 ????0
## 5 ????20671 ??????3130 ????1
相比于分別讀入這九個文本文件然后合并為一個計數(shù)矩陣，edgeR提供了更方便的途徑，使用readDGE函數(shù)即可一步完成纵穿。得到的DGEList對象中包含一個計數(shù)矩陣下隧，它的27179行分別對應每個基因不重復的Entrez基因ID，九列分別對應此實驗中的每個樣品谓媒。
x <- readDGE(files, columns=c(1,3))
class(x)
## [1] "DGEList"
## attr(,"package")
## [1] "edgeR"
dim(x)
dim()檢索或設置對象的尺寸淆院。
## [1] 27179 ????9
如果數(shù)據(jù)不是每個樣品一個文件的形式，而是一個包含所有樣品的計數(shù)的文件句惯，則可以先將文件讀入R土辩，再使用DGEList函數(shù)轉換為一個DGEList對象。
4.2組織樣品信息
為進行下游分析抢野，需要將有關實驗設計的樣品信息與計數(shù)矩陣的列關聯(lián)起來拷淘。這里需要包括各種對表達水平有影響的實驗變量，無論是生物變量還是技術變量指孤。例如启涯，細胞類型（在這個實驗中是basal、LP和ML）恃轩、基因型（野生型逝嚎、敲除）、表型（疾病狀態(tài)详恼、性別补君、年齡）、樣品處理（用藥昧互、對照）和批次信息（如果樣品是在不同時間點進行收集和分析的挽铁，需要記錄進行實驗的時間）等。
我們的DGEList對象中包含的samples數(shù)據(jù)框同時存儲了細胞類型（group）和批次（測序泳道lane）信息敞掘，每種信息都包含三個不同的水平叽掘。在x$samples中，程序會自動計算每個樣品的文庫大芯裂恪（即樣品的總序列計數(shù)）更扁，歸一化系數(shù)會被預先設置為1。為了方便閱讀赫冬，我們從DGEList對象x的列名中刪去了GEO樣品ID（GSM*）浓镜。
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))
samplenames
## [1] "10_6_5_11" "9_6_5_11" ?"purep53" ??"JMS8-2" ???"JMS8-3" ???"JMS8-4" ???"JMS8-5" ??
## [8] "JMS9-P7c" ?"JMS9-P8c"
colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
?????????????????????"Basal", "ML", "LP"))
x$samples$group <- group
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$samples
## ?????????????????????????????files group lib.size norm.factors lane
## 10_6_5_11 GSM1545535_10_6_5_11.txt ???LP 32863052 ???????????1 L004
## 9_6_5_11 ??GSM1545536_9_6_5_11.txt ???ML 35335491 ???????????1 L004
## purep53 ????GSM1545538_purep53.txt Basal 57160817 ???????????1 L004
## JMS8-2 ??????GSM1545539_JMS8-2.txt Basal 51368625 ???????????1 L006
## JMS8-3 ??????GSM1545540_JMS8-3.txt ???ML 75795034 ???????????1 L006
## JMS8-4 ??????GSM1545541_JMS8-4.txt ???LP 60517657 ???????????1 L006
## JMS8-5 ??????GSM1545542_JMS8-5.txt Basal 55086324 ???????????1 L006
## JMS9-P7c ??GSM1545544_JMS9-P7c.txt ???ML 21311068 ???????????1 L008
## JMS9-P8c ??GSM1545545_JMS9-P8c.txt ???LP 19958838 ???????????1 L008
4.3組織基因注釋
我們的DGEList對象中的第二個數(shù)據(jù)框名為genes，用于存儲與計數(shù)矩陣的行相關聯(lián)的基因信息劲厌。為檢索這些信息膛薛，我們可以使用特定物種的注釋包，比如小鼠的Mus.musculus?(Bioconductor Core Team 2016b)（或人類的Homo.sapiens?(Bioconductor Core Team 2016a)）补鼻；或者也可以使用biomaRt?包?(Durinck et al. 2005, 2009)哄啄，它通過接入Ensembl genome數(shù)據(jù)庫來進行基因注釋雅任。
可以檢索的信息類型包括基因符號（gene symbols）、基因名稱（gene names）咨跌、染色體名稱和位置沪么、Entrez基因ID、Refseq基因ID和Ensembl基因ID等锌半。biomaRt主要通過Ensembl基因ID進行檢索成玫，而Mus.musculus包含來自不同來源的信息，允許用戶從不同基因ID中選擇某一種作為檢索鍵拳喻。
我們使用Mus.musculus包，利用我們數(shù)據(jù)集中的Entrez基因ID來檢索相關的基因符號和染色體信息猪腕。
geneid <- rownames(x)
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),
????????????????keytype="ENTREZID")
head(genes)
## ???ENTREZID ?SYMBOL TXCHROM
## 1 ???497097 ???Xkr4 ???chr1
## 2 100503874 Gm19938 ???<NA>
## 3 100038431 Gm10568 ???<NA>
## 4 ????19888 ????Rp1 ???chr1
## 5 ????20671 ??Sox17 ???chr1
## 6 ????27395 ?Mrpl15 ???chr1
與任何基因ID一樣冗澈，Entrez基因ID可能不能一對一地匹配我們想獲得的基因信息。在處理之前陋葡，檢查重復的基因ID和弄清楚重復的來源非常重要亚亲。我們的基因注釋中包含28個能匹配到多個不同染色體的基因（比如基因Gm1987關聯(lián)于染色體chr4和chr4_JH584294_random，小RNA Mir5098關聯(lián)于chr2腐缤，chr5捌归，chr8，chr11和chr17）岭粤。為了處理重復的基因ID惜索，我們可以合并來自多重匹配基因的所有染色體信息，比如將基因Gm1987分配到chr4 and chr4_JH584294_random剃浇，或選取其中一條染色體來代表具有重復注釋的基因巾兆。為了簡單起見，我們選擇后者虎囚，保留每個基因ID第一次出現(xiàn)的信息角塑。
genes <- genes[!duplicated(genes$ENTREZID),]去除重復注釋的基因
在此例子中，注釋與數(shù)據(jù)對象中的基因順序是相同的淘讥。如果由于缺失和／或重新排列基因ID導致其順序不一致圃伶，我們可以用match函數(shù)來正確排序基因。然后蒲列，我們將基因注釋的數(shù)據(jù)框添加到DGEList對象窒朋，數(shù)據(jù)的整合就完成了，此時的數(shù)據(jù)對象中含有原始計數(shù)數(shù)據(jù)以及相關的樣品信息和基因注釋蝗岖。
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ù)預處理
5.1原始數(shù)據(jù)尺度轉換
由于更深的測序總會產(chǎn)生更多的序列片段炼邀，在差異表達及相關的分析中，我們很少直接使用序列數(shù)剪侮。在實際操作時拭宁，我們通常將原始的序列數(shù)進行歸一化洛退，來消除測序深度所導致的差異。通常被使用的方法有基于序列的CPM（counts per million）杰标、log-CPM兵怯、FPKM（fragments per kilobase of transcript per million），和基于轉錄本數(shù)目的RPKM（reads per kilobase of transcript per million）腔剂。
我們在分析中通常使用CPM和log-CPM轉換媒区。雖然RPKM和FPKM可以校正基因長度區(qū)別的影響，但CPM和log-CPM只使用計數(shù)矩陣即可計算掸犬，且已足以滿足我們所關注的比較的需要袜漩。假設不同條件之間剪接異構體（isoform）的表達比例沒有變化，差異表達分析關注的是同一基因在不同條件之間表達水平的相對差異湾碎，而不是比較多個基因之間的差異或測定絕對表達量宙攻。換而言之，基因長度在我們進行比較的不同組之間是始終不變的介褥，且任何觀測到的差異都來自于不同組的條件的變化而不是基因長度的變化座掘。
我們使用edgeR中的cpm函數(shù)將原始計數(shù)轉換為CPM和log-CPM值。如果可以提供基因長度信息柔滔，RPKM值的計算也和CPM值的計算一樣簡單溢陪，只需使用edgeR中的rpkm函數(shù)。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)
對于一個基因睛廊，CPM值為1相當于在本實驗測序深度最低的樣品中（JMS9-P8c, 文庫大小約2千萬）有20個計數(shù)形真，或者在測序深度最高的樣品中（JMS8-3，文庫大小約7.6千萬）有76個計數(shù)超全。
log-CPM值將被用于探索性圖表中没酣。當設置log=TRUE時，cpm函數(shù)會給CPM值加上一個彌補值并進行l(wèi)og2轉換卵迂。默認的彌補值是2/L裕便，其中2是“預先計數(shù)”，而L是樣本文庫大屑洹（以百萬計）的平均值偿衰，所以log-CPM值是根
CPM值通過log2(CPM + 2/L)計算得到的。這樣的計算方式可以確保任意兩個具有相同CPM值的序列片段計數(shù)的log-CPM值也相同改览。彌補值的使用可以避免對零取對數(shù)下翎，并能使所有樣本間的對數(shù)倍數(shù)變化（log-fold-change）向0推移而減小低表達基因間微小計數(shù)變化帶來的巨大的偽差異性，這對于繪制探索性圖表很有幫助宝当。在這個數(shù)據(jù)集中视事，平均的樣本文庫大小是4.55千萬，所以L約等于45.5庆揩，且每個樣本中的最小log-CPM值為log2(2/45.5) = -4.51俐东。換而言之跌穗，在加上了預先計數(shù)彌補值后，此數(shù)據(jù)集中的零表達計數(shù)對應的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ù)時也會用到log-CPM值蚌吸，但voom會默認使用更小的預先計數(shù)重新計算自己的log-CPM值。
5.2刪除低表達基因
所有數(shù)據(jù)集中都混有表達的基因與不表達的基因砌庄。我們想要檢測在一種條件中表達但在另一種條件中不表達的基因羹唠，但也有一些基因在所有樣品中都不表達。實際上娄昆，這個數(shù)據(jù)集中19%的基因在所有九個樣品中的計數(shù)都是零佩微。
table(rowSums(x$counts==0)==9)
##
## FALSE ?TRUE
## 22026 ?5153
Figure 1:?每個樣本過濾前的原始數(shù)據(jù)（A）和過濾后（B）的數(shù)據(jù)的log-CPM值密度。豎直虛線標出了過濾步驟中所用閾值（相當于CPM值為約0.2）萌焰。
log-CPM值的分布圖表顯示每個樣本中很大一部分基因都是不表達或者表達程度相當?shù)偷牟该校鼈兊膌og-CPM值非常小甚至是負的（圖1A）。
在任何樣本中都沒有足夠多的序列片段的基因應該從下游分析中過濾掉杆怕。這樣做的原因有好幾個。從生物學的角度來看壳贪，在任何條件下的表達水平都不具有生物學意義的基因都不值得關注陵珍，因此最好忽略。從統(tǒng)計學的角度來看违施，去除低表達計數(shù)基因使數(shù)據(jù)中的均值- 方差關系可以得到更精確的估計互纯，并且還減少了下游的差異表達分析中需要進行的統(tǒng)計檢驗的數(shù)量。
edgeR包中的filterByExpr函數(shù)提供了自動過濾基因的方法磕蒲，可保留盡可能多的有足夠表達計數(shù)的基因留潦。
keep.exprs <- filterByExpr(x, group=group)
x <- x[keep.exprs,, keep.lib.sizes=FALSE]#keep.lib.sizes=FALSE參數(shù)在過濾結束后重新計算樣本庫size，一般推薦是這樣做的辣往。盡管做不做對后續(xù)影響都不大...
dim(x)
## [1] 16624 ????9
此函數(shù)默認選取最小的組內的樣本數(shù)量為最小樣本數(shù)兔院，保留至少在這個數(shù)量的樣本中有10個或更多計數(shù)的基因。實際進行過濾時站削，使用的是CPM值而不是表達計數(shù)坊萝，以避免對總序列數(shù)大的樣本的偏向性。在這個數(shù)據(jù)集中许起，總序列數(shù)的中位數(shù)是5.1千萬十偶，且10/51約等于0.2，所以filterByExpr函數(shù)保留在至少三個樣本中CPM值大于等于0.2的基因园细。在我們的此次實驗中惦积，一個具有生物學意義的基因需要在至少三個樣本中表達，因為三種細胞類型組內各有三個重復猛频。過濾的閾值取決于測序深度和實驗設計狮崩。如果樣本總表達計數(shù)量增大蛛勉，那么可以選擇更低的CPM閾值，因為更大的總表達計數(shù)量提供了更好的分辨率來探究更多表達水平更低的基因厉亏。
使用這個標準董习，基因的數(shù)量減少到了16624個，約為開始時數(shù)量的60%爱只。過濾后的log-CPM值顯示出每個樣本的分布基本相同（下圖B部分）皿淋。需要注意的是，從整個DGEList對象中取子集時同時刪除了被過濾的基因的計數(shù)和其相關的基因信息恬试。留下的基因相對應的基因信息和計數(shù)在過濾后的DGEList對象中被保留窝趣。
下方給出的是繪圖所用代碼。
lcpm.cutoff <- log2(10/M + 2/L)
library(RColorBrewer)
nsamples <- ncol(x)???#nsamples是9
col <- brewer.pal(nsamples, "Paired")#顏色名為?"Paired"训柴，支持12種顏色區(qū)分哑舒，種類為qual型（區(qū)分色，幾種區(qū)分度很高的顏色組合）幻馁，colorblind為TRUE（對色盲友好）
par(mfrow=c(1,2))?# par(mfrow=c(1,2))實現(xiàn)一頁多圖的功能洗鸵，其中，通過設定函數(shù)par()的各個參數(shù)來調整我們的圖形仗嗦，mfrow=c(2,2) 就是畫4幅圖膘滨，mfrow=c(3,5),是畫15幅圖，
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)#加lcpm.cutoff這條分界線稀拐，設置樣式為虛線
for?(i in?2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)#繪制從2到9的其余樣本的散點圖
}
legend("topright", samplenames, text.col=col, bty="n")?#"topright"定位圖例火邓，圖例名稱為?samplenames，text.col=col設置圖例顏色德撬，bty="n"圖例框不劃出
lcpm <- cpm(x, log=TRUE)此時x已被削減至16624個
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歸一化基因表達分布
在樣品制備或測序過程中铲咨，不具備生物學意義的外部因素會影響單個樣品的表達。比如說蜓洪，在實驗中第一批制備的樣品會總體上表達高于第二批制備的樣品纤勒。差異表達分析假設所有樣品表達值的范圍和分布都應當相似。我們需要進行歸一化來確保整個實驗中每個樣本的表達分布都相似隆檀。
密度圖和箱線圖等展示每個樣品基因表達量分布的圖表可以用于判斷是否有樣品和其他樣品分布有差異踊东。在此數(shù)據(jù)集中，所有樣品的log-CPM分布都很相似（上圖B部分）刚操。
盡管如此闸翅，我們依然需要使用edgeR中的calcNormFactors函數(shù)，用TMM(Robinson and Oshlack 2010)方法進行歸一化菊霜。此處計算得到的歸一化系數(shù)被用作文庫大小的縮放系數(shù)坚冀。當我們使用DGEList對象時，這些歸一化系數(shù)被自動存儲在x$samples$norm.factors鉴逞。對此數(shù)據(jù)集而言记某，TMM歸一化的作用比較溫和司训，這體現(xiàn)在所有的縮放因子都相對接近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
在這里液南，為了更好地展示出歸一化的效果壳猜，我們復制了數(shù)據(jù)并進行了人工調整，使得第一個樣品的計數(shù)減少到了其原始值的5%滑凉，而第二個樣品增大到了5倍统扳。要注意在實際的數(shù)據(jù)分析流程中，不應當進行這樣的操作畅姊。
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ù)的表達分布咒钟，其中歸一化前不同樣本的分布明顯不同，而歸一化后比較相似若未。此處朱嘴，經(jīng)過我們人工調整的第一個樣品的TMM縮放系數(shù)0.06非常小，而第二個樣品的縮放系數(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對樣本的無監(jiān)督聚類
在我們看來，用于檢查基因表達分析的最重要的探索性圖表之一便是MDS圖甚淡，或類似的圖大诸。這種圖表使用無監(jiān)督聚類方法展示出了樣品間的相似性和不相似性捅厂，能讓我們在進行正式的檢驗之前對于能檢測到多少差異表達基因有個大致概念贯卦。理想情況下，樣本會在各個實驗組內很好的聚類焙贷，且我們可以鑒別出遠離所屬組的樣本撵割，并追蹤誤差或額外方差的來源。如果存在的話辙芍，技術重復應當互相非常接近啡彬。
這樣的圖可以用limma中的plotMDS函數(shù)繪制。第一個維度表示能夠最好地分離樣品且解釋最大比例的方差的領先倍數(shù)變化（leading-fold-change）故硅，而后續(xù)的維度的影響更小庶灿，并與之前的維度正交。當實驗設計涉及到多個因子時吃衅，建議在多個維度上檢查每個因子往踢。如果在其中一些維度上樣本可按照某因子聚類，這說明該因子對于表達差異有影響徘层，在線性模型中應當將其包括進去峻呕。反之利职，沒有或者僅有微小影響的因子在下游分析時應當被剔除。
在這個數(shù)據(jù)集中瘦癌，可以看出樣本在維度1和2能很好地按照實驗分組聚類猪贪，隨后在維度3按照測序泳道（樣品批次）分離（如下圖所示）。由于第一維度解釋了數(shù)據(jù)中最大比例的方差讯私，我們會發(fā)現(xiàn)當關注更高維度時热押，維度上的取值范圍會變小。
盡管所有樣本都按組聚類妄帘，在維度1上最大的轉錄差異出現(xiàn)在basal和LP以及basal和ML之間楞黄。因此，預期在basal樣品與其他之間的成對比較中能夠得到大量的DE基因抡驼，而在ML和LP之間的比較中得到的DE基因數(shù)量略少鬼廓。在其他的數(shù)據(jù)集中扫步，不按照實驗組聚類的樣本可能在下游分析中只表現(xiàn)出較小的或不表現(xiàn)出差異表達坐儿。
為繪制MDS圖，我們?yōu)椴煌囊蜃釉O立不同的配色啃匿。維度1和2以細胞類型上色馏锡，而維度3和4以測序泳道（批次）上色雷蹂。
lcpm <- cpm(x, log=TRUE)
par(mfrow=c(1,2))
col.group <- group# col.group為 LP ?ML ??Basal Basal ML ??LP ?Basal ML ???LP ??Levels: Basal LP ML
levels(col.group) <- ?brewer.pal(nlevels(col.group), "Set1")??
#nlevels(col.group)是3，此時col.group
變成#377EB8 #4DAF4A #E41A1C #E41A1C #4DAF4A #377EB8 #E41A1C #4DAF4A #377EB8
Levels: #E41A1C #377EB8 #4DAF4A
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")
作為另一種繪制MDS圖的方式党巾，Glimma包提供了便于探索多個維度的交互式MDS圖萎庭。其中的glMDSPlot函數(shù)可生成一個html網(wǎng)頁（如果設置launch=TRUE參數(shù)，將會在生成后直接在瀏覽器中打開）齿拂，其左側面板含有一張MDS圖驳规，而右側面板包含一張展示了各個維度所解釋的方差比例的柱形圖。點擊柱形圖中的柱可切換MDS圖繪制時所使用的維度署海，且將鼠標懸浮于單個點上可顯示相應的樣本標簽吗购。也可切換配色方案，以突顯不同細胞類型或測序泳道（批次）砸狞。此數(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 差異表達分析
6.1創(chuàng)建設計矩陣和對比
在此研究中，我們想知道哪些基因在我們研究的三組細胞之間以不同水平表達刀森。我們的分析中所用到的線性模型假設數(shù)據(jù)是正態(tài)分布的踱启。首先，我們要創(chuàng)建一個包含細胞類型以及測序泳道（批次）信息的設計矩陣。
design <- model.matrix(~0+group+lane)
colnames(design) <- gsub("group", "", colnames(design))
design
## ??Basal LP ML laneL006 laneL008
## 1 ????0 ?1 ?0 ???????0 ???????0
## 2 ????0 ?0 ?1 ???????0 ???????0
## 3 ????1 ?0 ?0 ???????0 ???????0
## 4 ????1 ?0 ?0 ???????1 ???????0
## 5 ????0 ?0 ?1 ???????1 ???????0
## 6 ????0 ?1 ?0 ???????1 ???????0
## 7 ????1 ?0 ?0 ???????1 ???????0
## 8 ????0 ?0 ?1 ???????0 ???????1
## 9 ????0 ?1 ?0 ???????0 ???????1
## attr(,"assign")
## [1] 1 1 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
##
## attr(,"contrasts")$lane
## [1] "contr.treatment"
對于一個給定的實驗禽捆，通常有多種等價的方法都能用來創(chuàng)建合適的設計矩陣笙什。比如說，~0+group+lane去除了第一個因子group的截距胚想，但第二個因子lane的截距被保留琐凭。此外也可以使用~group+lane，來自group和lane的截距均被保留浊服。
> design2 <- model.matrix(~group+lane)
> colnames(design2) <- gsub("group", "", colnames(design))
> design2
??Basal LP ML laneL006 laneL008
1 ????1 ?1 ?0 ???????0 ???????0
2 ????1 ?0 ?1 ???????0 ???????0
3 ????1 ?0 ?0 ???????0 ???????0
4 ????1 ?0 ?0 ???????1 ???????0
5 ????1 ?0 ?1 ???????1 ???????0
6 ????1 ?1 ?0 ???????1 ???????0
7 ????1 ?0 ?0 ???????1 ???????0
8 ????1 ?0 ?1 ???????0 ???????1
9 ????1 ?1 ?0 ???????0 ???????1
attr(,"assign")
[1] 0 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
attr(,"contrasts")$lane
[if !supportLists][1]?[endif]"contr.treatment"
理解如何解釋模型中估計的系數(shù)是創(chuàng)建合適的設計矩陣的關鍵统屈。我們在此分析中選取第一種模型，因為在沒有group的截距的情況下能更直截了當?shù)卦O定模型中的對比牙躺。用于細胞群之間成對比較的對比可以在limma中用makeContrasts函數(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)勢之一便是其適應任意實驗復雜程度的能力。簡單的實驗設計孽拷，比如此流程中關于細胞類型和批次的實驗設計吨掌，直到更復雜的析因設計和含有交互作用項的模型，都能夠被較簡單地處理脓恕。當實驗或技術效應可被隨機效應模型（random effect model）模擬時膜宋，可以使用limma中的duplicateCorrelation函數(shù)來估計交互作用，這需要在此函數(shù)以及l(fā)mFit的線性建模步驟均指定一個block參數(shù)炼幔。
6.2從表達計數(shù)數(shù)據(jù)中刪除異方差
對于RNA-seq計數(shù)數(shù)據(jù)而言秋茫，原始計數(shù)或log-CPM值的方差并不獨立于均值(Law et al. 2014)。有些差異表達分析方法使用負二項分布模型乃秀，假設均值與方差間具有二次的關系肛著。而在limma中，假設log-CPM值符合正態(tài)分布跺讯，因此我們在對RNA-seq的log-CPM值進行線性建模時枢贿，需要使用voom函數(shù)計算每個基因的權重從而調整均值與方差的關系，否則分析得到的結果可能是錯誤的抬吟。
當操作DGEList對象時萨咕，voom從x中自動提取文庫大小和歸一化因子统抬，以此將原始計數(shù)轉換為log-CPM值火本。在voom中，對于log-CPM值額外的歸一化可以通過設定normalize.method參數(shù)來進行聪建。
下圖左側展示了這個數(shù)據(jù)集log-CPM值的均值-方差關系钙畔。通常而言，方差是測序實驗操作中的技術差異和來自不同細胞類群的重復樣本之間的生物學差異的結合金麸，而voom圖會顯示出一個在均值與方差之間遞減的趨勢。生物學差異高的實驗通常會有更平坦的趨勢桨醋，其方差值在高表達處穩(wěn)定。生物學差異低的實驗更傾向于急劇下降的趨勢现斋。
不僅如此喜最，voom圖也提供了對于上游所進行的過濾水平的可視化檢測。如果對于低表達基因的過濾不夠充分庄蹋，在圖上表達低的一端瞬内，受到非常低的表達計數(shù)的影響，會出現(xiàn)方差的下降限书。如果觀察到了這種情況虫蝶，應當回到最初的過濾步驟并提高用于該數(shù)據(jù)集的表達閾值。
當先前繪制的MDS圖中發(fā)現(xiàn)組內重復樣本的聚集與分離程度出現(xiàn)明顯的差異時倦西，可以用voomWithQualityWeights函數(shù)(Liu et al. 2015)來代替voom能真，在計算基因權重值以外還能計算每個樣本的權重值。關于使用此種方式的例子參見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對象中的其他數(shù)據(jù)框，即基因和樣本信息普办，也保留在了voom創(chuàng)建的EList對象v中工扎。v$genes數(shù)據(jù)框等同于x$genes，v$targets等同于x$samples衔蹲，而v$E中所儲存的表達值類似于進行了log-CPM轉換后的x$counts肢娘。此外，voom的EList對象中還有一個權重值的矩陣v$weights，而設計矩陣存儲于v$design橱健。
6.3擬合線性模型以進行比較
limma的線性建模使用lmFit和contrasts.fit函數(shù)進行而钞，它們原先是為微陣列而設計的。這些函數(shù)不僅可以用于微陣列數(shù)據(jù)拘荡，也可以用于使用voom計算了基因權重后的RNA-seq數(shù)據(jù)笨忌。每個基因的表達值都會單獨擬合一個模型。然后通過借用全體基因的信息來進行經(jīng)驗貝葉斯調整（empirical Bayes moderation）俱病，這樣可以更精確地估算各基因的差異性(Smyth 2004)官疲。圖4為此模型的殘差關于平均表達值的圖。從圖中可以看出亮隙，方差不再與表達水平均值相關途凫。
6.4檢查DE基因數(shù)量
為快速查看表達差異水平，顯著上調或下調的基因可以匯總到一個表格中溢吻。顯著性的判斷使用默認的校正p值閾值维费，即5%。在basal與LP的表達水平之間的比較中促王，發(fā)現(xiàn)了4648個在basal中相較于LP下調的基因犀盟，和4863個在basal中相較于LP上調的基因，即共9511個差異表達基因蝇狼。在basal和ML之間發(fā)現(xiàn)了一共9598個差異表達基因（4927個下調基因和4671個上調基因）阅畴，而在LP和ML中發(fā)現(xiàn)了一共5652個差異表達基因（3135個下調基因和2517個上調基因）。在basal類群所參與的比較中皆找到了大量的差異表達基因迅耘，這與我們在MDS圖中觀察到的結果相吻合贱枣。
summary(decideTests(efit))
## ???????BasalvsLP BasalvsML LPvsML
## Down ???????4648 ?????4927 ??3135
## NotSig ?????7113 ?????7026 ?10972
## Up ?????????4863 ?????4671 ??2517
在某些時候，不僅僅需要用到校正p值閾值颤专，還需要差異倍數(shù)的對數(shù)（log-FCs）也高于某個最小值來更為嚴格地定義顯著性纽哥。treat方法(McCarthy and Smyth 2009)可以按照對最小log-FC值的要求，使用經(jīng)過經(jīng)驗貝葉斯調整的t統(tǒng)計值計算p值栖秕。當我們要求差異表達基因的log-FC顯著大于1（等同于不同細胞類群之間表達量差兩倍）并對此-進行檢驗時春塌，差異表達基因的數(shù)量得到了下降，basal與LP相比只有3684個差異表達基因簇捍，basal與ML相比只有3834個差異表達基因只壳，LP與ML相比只有414個差異表達基因。
tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)
## ???????BasalvsLP BasalvsML LPvsML
## Down ???????1632 ?????1777 ???224
## NotSig ????12976 ????12790 ?16210
## Up ?????????2016 ?????2057 ???190
在多個對比中皆差異表達的基因可以從decideTests的結果中提取垦写，其中的0代表不差異表達的基因吕世，1代表上調的基因彰触，-1代表下調的基因梯投。共有2784個基因在basal和LP以及basal和ML的比較中都差異表達，其中的20個于下方列出。write.fit函數(shù)可用于將三個比較的結果一起提取并寫入一個輸出文件分蓖。
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" ????????"Cracdl" ???????"Aff3" ?????????"Npas2"
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
write.fit(tfit, dt, file="results.txt")
6.5從頭到尾檢查單個DE基因
使用topTreat函數(shù)可以列舉出使用treat得到的結果中靠前的DE基因（對于eBayes的結果則應使用topTable函數(shù)）。默認情況下蒸甜，topTreat將基因按照經(jīng)過多重假設檢驗校正的p值從小到大排列棠耕，并為每個基因給出相關的基因信息、log-FC柠新、平均log-CPM窍荧、校正t值、原始及校正p值恨憎。列出前多少個基因的數(shù)量可由用戶設定蕊退，如果設為n=Inf則會包括所有的基因°究遥基因Cldn7和Rasef在basal與LP和basal于ML的比較中都位于DE基因的前幾名瓤荔。
basal.vs.lp <- topTreat(tfit, coef=1, n=Inf)
basal.vs.ml <- topTreat(tfit, coef=2, n=Inf)
head(basal.vs.lp)
## ???????ENTREZID SYMBOL TXCHROM logFC AveExpr ????t ?P.Value adj.P.Val
## 12759 ????12759 ???Clu ??chr14 -5.46 ???8.86 -33.6 1.72e-10 ?1.71e-06
## 53624 ????53624 ?Cldn7 ??chr11 -5.53 ???6.30 -32.0 2.58e-10 ?1.71e-06
## 242505 ??242505 ?Rasef ???chr4 -5.94 ???5.12 -31.3 3.08e-10 ?1.71e-06
## 67451 ????67451 ??Pkp2 ??chr16 -5.74 ???4.42 -29.9 4.58e-10 ?1.74e-06
## 228543 ??228543 ??Rhov ???chr2 -6.26 ???5.49 -29.1 5.78e-10 ?1.74e-06
## 70350 ????70350 ?Basp1 ??chr15 -6.08 ???5.25 -28.3 7.27e-10 ?1.74e-06
head(basal.vs.ml)
## ???????ENTREZID ?SYMBOL TXCHROM logFC AveExpr ????t ?P.Value adj.P.Val
## 242505 ??242505 ??Rasef ???chr4 -6.53 ???5.12 -35.1 1.23e-10 ?1.24e-06
## 53624 ????53624 ??Cldn7 ??chr11 -5.50 ???6.30 -31.7 2.77e-10 ?1.24e-06
## 12521 ????12521 ???Cd82 ???chr2 -4.69 ???7.07 -31.4 2.91e-10 ?1.24e-06
## 20661 ????20661 ??Sort1 ???chr3 -4.93 ???6.70 -30.7 3.56e-10 ?1.24e-06
## 71740 ????71740 Nectin4 ???chr1 -5.58 ???5.16 -30.6 3.72e-10 ?1.24e-06
## 12759 ????12759 ????Clu ??chr14 -4.69 ???8.86 -28.0 7.69e-10 ?1.48e-06
6.6實用的差異表達結果可視化
為可視化地總結所有基因的結果，可使用plotMD函數(shù)繪制均值-差異關系（MD）圖钥组，其中展示了線性模型擬合所得到的每個基因log-FC與log-CPM平均值間的關系输硝，而差異表達的基因會被重點標出。
plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1],
???????xlim=c(-8,13))
Glimma的glMDPlot函數(shù)提供了交互式的均值-差異圖程梦，拓展了這種圖表的功能腔丧。此函數(shù)的輸出為一個html頁面，左側面板為結果的總結性圖表（與plotMD的輸出類似）作烟，右側面板展示了各個樣本的log-CPM值愉粤，而下方面板為結果的匯總表。這使得用戶可以使用提供的注釋中的信息（比如基因符號）搜索特定基因拿撩，而這在R統(tǒng)計圖中是做不到的衣厘。
glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1],
?????????side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)文件生成在當前目錄
上方指令生成的均值-差異圖可以在線訪問（詳見http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MD-Plot.html）。Glimma提供的交互性使得單個圖形窗口內能夠呈現(xiàn)出額外的信息伦吠。Glimma是以R和Javascript實現(xiàn)的妆兑，使用R代碼生成數(shù)據(jù)魂拦，并在之后使用Javascript庫D3（https://d3js.org）轉換為圖形，使用Bootstrap庫處理界面并生成互動性可搜索的表格的數(shù)據(jù)表搁嗓。這使得圖表可以在任意當代瀏覽器中直接查看而無需后段服務器時刻運行芯勘，對于將其作為關聯(lián)文件附加在Rmarkdown分析報告中而言非常方便。
前文所展示的圖表中腺逛，一些僅展示了在任意一個條件下表達的所有基因而缺少單個基因的具體信息（比如不同對比中共同差異表達基因的韋恩圖或均值-差異圖）荷愕，而另一些僅展示單個基因（交互式均值-差異圖右邊面板中所展示的log-CPM值）。而熱圖使用戶得以查看某一小組基因的表達情況棍矛，既便于查看單個組或樣本的表達安疗，又不至于在關注于單個基因時失去對于整體的注意，也不會造成由于對上千個基因取平均值而導致的分辨率丟失够委。
使用gplots包的heatmap.2函數(shù)茂契，我們?yōu)閎asal與LP的對照中前100個差異表達基因（按校正p值排序）繪制了一幅熱圖。此熱圖中正確地將樣本按照細胞類型聚類慨绳，并重新排序了基因掉冶，表達相似的基因形成了塊狀。從熱圖中脐雪，我們發(fā)現(xiàn)basal與LP之間的前100個差異表達基因在ML和LP樣本中的表達非常相似厌小。
library(gplots)
basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100]
i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes)
mycol <- colorpanel(1000,"blue","white","red")
heatmap.2(lcpm[i,], scale="row", labRow=v$genes$SYMBOL[i], labCol=group, col=mycol, trace="none", density.info="none", margin=c(8,6), lhei=c(2,10), dendrogram="column")
7使用camera的基因集檢驗
在此次分析的最后，我們要進行一些基因集檢驗埋泵。為此幔欧，我們將camera方法(Wu and Smyth 2012)應用于Broad Institute的MSigDB c2中的(Subramanian et al. 2005)中適應小鼠的c2基因表達特征，這可從http://bioinf.wehi.edu.au/software/MSigDB/以RData對象格式獲取丽声。 此外礁蔗，對于人類和小鼠，來自MSigDB的其他有用的基因集也可從此網(wǎng)站獲取雁社，比如標志（hallmark）基因集浴井。C2基因集的內容收集自在線數(shù)據(jù)庫、出版物以及該領域專家霉撵，而標志基因集的內容來自MSigDB磺浙，從而獲得具有明確定義的生物狀態(tài)或過程洪囤。
load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123"))出錯，改用load("C:/Users/葛鵬/Documents/R/mouse_c2_v5p1.rdata")
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ù)通過比較假設檢驗來評估一個給定基因集中的基因是否相對于不在集內的基因而言在差異表達基因的排序中更靠前屠缭。它使用limma的線性模型框架箍鼓，并同時采用設計矩陣和對比矩陣（如果有的話）崭参，且在檢驗的過程中會運用到來自voom的權重值呵曹。在通過基因間相關性（默認設定為0.01，但也可通過數(shù)據(jù)估計）和基因集的規(guī)模得到方差膨脹因子（variance inflation factor）何暮，并使用它調整基因集檢驗統(tǒng)計值的方差后奄喂，將會返回根據(jù)多重假設檢驗進行了校正的p值。
Lim等人(2010)(Lim et al. 2010)使用Illumina微陣列分析了與此實驗相同的分選細胞群海洼，而我們的實驗是與他們的數(shù)據(jù)集等價的RNA-seq跨新，因此我們預期來自該早期文獻的基因表達特征將會出現(xiàn)在每種對比的列表頂端，而結果正符合我們的預期坏逢。在LP和ML的對比中域帐，我們?yōu)長im等人（2010）的成熟管腔基因集（上調及下調）繪制了條碼圖（barcodeplot）。由于我們的對比是將LP與ML相比而不是相反是整，這些基因集的方向在我們的數(shù)據(jù)集中是反過來的（如果將對比反過來肖揣，基因集的方向將會與對比一致）。
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP, index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
與我們的非常相似彤断，用了相同的方式分選獲取細胞群，只是他們使用了微陣列而不是RNA-seq來測定基因表達易迹。需要注意上調基因集發(fā)生下調而下調基因集發(fā)生上調的逆相關性來自于對比的設定方式（LP相比于ML）宰衙，如果將其對調，方向性將會吻合睹欲。
limma還提供了另外的基因集檢驗方法菩浙，比如mroast(Wu et al. 2010)的自包含檢驗。camera適用于檢驗包含許多基因集的大型數(shù)據(jù)庫中哪些基因集相對于其他基因集整體變化更為顯著（如前文所示）句伶，自包含檢驗則更善于集中檢驗一個或少個選中的基因集自身是否差異表達劲蜻。換句話說，camera更適用于找出具有意義的基因集考余，而mroast測試的是已經(jīng)確定有意義的基因集的顯著性先嬉。
參考文獻
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): 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.
