使用limma、Glimma和edgeR，RNA-seq數據分析易如反掌
使用limma、Glimma和edgeR，RNA-seq數據分析易如反掌
23 October 2018
Contents
1摘要
簡單且高效地分析RNA測序數據的能力正是Bioconductor的核心優點。 RNA-seq分析一般從基因水平的序列計數開始，涉及到數據預處理，探索性數據分析，差別表達檢驗以及通路分析，獲得的結果可用於指導進一步實驗和驗證研究。 在這篇工做流程文章中，咱們經過分析來自小鼠乳腺的RNA測序數據，示範瞭如何使用流行的edgeR包載入、整理、過濾和歸一化數據，而後用limma包的voom方法、線性模型和經驗貝葉斯調節（empirical Bayes moderation）來評估差別表達並進行基因集檢驗。經過使用Glimma包，此流程獲得了增進，實現告終果的互動探索，使用戶得以查看單個樣本與基因。 這三個軟件包提供的完整分析突出了研究人員可使用Bioconductor輕鬆地從RNA測序實驗的原始計數揭示生物學意義。ios
2背景介紹
RNA測序（RNA-seq）已成爲基因表達研究的主要技術。其中，基因組規模的多條件基因差別表達檢測是研究者最常探究的問題之一。對於RNA-seq數據，來自Bioconductor項目(Huber et al. 2015)的 edgeR (Robinson, McCarthy, and Smyth 2010)和limma包 (Ritchie et al. 2015)提供了一套用於處理此問題的完善的統計學方法。git
在這篇文章中，咱們描述了一個用於分析RNA-seq數據的edgeR - limma工做流程，使用基因水平計數做爲輸入，通過預處理和探索性數據處理，而後獲得差別表達（DE）基因和基因表達特徵（gene signatures）的列表。Glimma包(Su and Ritchie 2016)提供的交互式圖表能夠同時呈現總體樣本和單個基因水平的數據，使得咱們相對靜態的R圖表而言，能夠探索更多的細節。數據庫
此工做流程中所分析的實驗來自Sheridan等（2015）(Sheridan et al. 2015)，由三個細胞羣組成，即基底（basal）、管腔祖細胞（liminal progenitor, LP）和成熟管腔（mature luminal, ML）。細胞羣皆分選自雌性處女小鼠的乳腺，每種都設三個生物學重複。RNA樣品分三個批次使用Illumina HiSeq 2000進行測序，獲得長爲100鹼基對的單端序列片斷。 本文所描述的分析假設從RNA-seq實驗得到的序列片斷已經與適當的參考基因組比對，並已經在基因水平上對序列進行了統計計數。在本文條件下，使用Rsubread包提供的基於R的流程將序列片斷與小鼠參考基因組（mm20）比對（具體而言，先使用align函數(Liao, Smyth, and Shi 2013)，而後使用featureCounts (Liao, Smyth, and Shi 2014)進行基因水平的總結，利用其內置的mm10基於RefSeq的註釋），json
這些樣本的計數數據能夠從Gene Expression Omnibus (GEO)數據庫 http://www.ncbi.nlm.nih.gov/geo/使用GEO序列登記號GSE63310下載。更多關於實驗設計和樣品製備的信息也能夠在GEO使用該登記號查看。api
3初始配置
library(limma) library(Glimma) library(edgeR) library(Mus.musculus)
4數據整合
4.1讀入計數數據
爲開始此分析，從https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file在線下載文件GSE63310_RAW.tar，並從壓縮包中解壓出相關的文件。下方的代碼將完成此步驟，或者您也能夠手動進行這一步並繼續後續分析。瀏覽器
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file" utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb") utils::untar("GSE63310_RAW.tar", exdir = ".") files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt") for(i in paste(files, ".gz", sep="")) R.utils::gunzip(i, overwrite=TRUE)
每個文本文件均包含一個給定樣品的原始基因水平計數。須要注意的是，咱們的分析僅包含了此實驗中的basal、LP和ML樣品（請查看下方相關文件名）。markdown
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt") read.delim(files[1], nrow=5)
## EntrezID GeneLength Count ## 1 497097 3634 1 ## 2 100503874 3259 0 ## 3 100038431 1634 0 ## 4 19888 9747 0 ## 5 20671 3130 1
儘管這九個文本文件能夠分別讀入R而後合併爲一個計數矩陣，edgeR提供了更方便的途徑，使用readDGE函數便可一步完成。獲得的DGEList對象中包含一個計數矩陣，它的27179行分別對應惟一的Entrez基因標識（ID），九列分別對應此實驗中的每一個樣品。session
x <- readDGE(files, columns=c(1,3)) class(x)
## [1] "DGEList" ## attr(,"package") ## [1] "edgeR"
dim(x)
## [1] 27179 9
若是來自全部樣品的計數存儲在同一個文件中，數據能夠首先讀入R再使用DGEList函數轉換爲一個DGEList對象。
4.2組織樣品信息
爲進行下游分析，與實驗設計有關的樣品水平信息須要與計數矩陣的列關聯。這裏須要包括各類對錶達水平有影響的實驗變量，不管是生物變量仍是技術變量。例如，細胞類型（在這個實驗中是basal、LP和ML），基因型（野生型、敲除），表型（疾病狀態、性別、年齡），樣品處理（用藥、對照）和批次信息（若是樣品是在不一樣時間點進行收集和分析的，記錄進行實驗的時間）等。
咱們的DGEList對象中包含的samples數據框同時存儲了細胞類型（group）和批次（測序泳道 lane）信息，每種信息都包含三個不一樣的水平。須要注意的是，在x$samples中，程序會自動計算每一個樣品的文庫大小，歸一化係數會被設置爲1. 爲了簡單起見，咱們從咱們的DGEList對象x的列名中刪去了GEO樣品ID（GSM*）。
samplenames <- substring(colnames(x), 12, nchar(colnames(x))) samplenames
## [1] "10_6_5_11" "9_6_5_11" "purep53" "JMS8-2" "JMS8-3" "JMS8-4" "JMS8-5" ## [8] "JMS9-P7c" "JMS9-P8c"
colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", "Basal", "ML", "LP")) x$samples$group <- group lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2))) x$samples$lane <- lane x$samples
## files group lib.size norm.factors lane ## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004 ## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004 ## purep53 GSM1545538_purep53.txt Basal 57160817 1 L004 ## JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006 ## JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006 ## JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006 ## JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006 ## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008 ## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008
4.3組織基因註釋
咱們的DGEList對象中的第二個數據框名爲genes，用於存儲與計數矩陣的行相關聯的基因水平的信息。 爲檢索這些信息，咱們可使用包含特定物種信息的包，好比小鼠的Mus.musculus(Bioconductor Core Team 2016b)（或人類的Homo.sapiens (Bioconductor Core Team 2016a)）；或者也可使用biomaRt 包 (Durinck et al. 2005, 2009)，它經過接入Ensembl genome數據庫來進行基因註釋。
能夠檢索的信息類型包括基因符號（gene symbols）、基因名稱（gene names）、染色體名稱和位置（chromosome names and locations）、Entrez基因ID（Entrez gene IDs）、Refseq基因ID（Refseq gene IDs）和Ensembl基因ID（Ensembl gene IDs）等。biomaRt 主要使用Ensembl基因ID進行檢索，而因爲Mus.musculus包含多種不一樣來源的信息，它容許用戶從多種不一樣基因ID中選取檢索鍵。
咱們使用Mus.musculus包，利用咱們數據集中的Entrez基因ID來檢索相關的基因符號和染色體信息。
geneid <- rownames(x)
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"), keytype="ENTREZID") head(genes)
## ENTREZID SYMBOL TXCHROM ## 1 497097 Xkr4 chr1 ## 2 100503874 Gm19938 <NA> ## 3 100038431 Gm10568 <NA> ## 4 19888 Rp1 chr1 ## 5 20671 Sox17 chr1 ## 6 27395 Mrpl15 chr1
與任何基因ID同樣，Entrez基因ID可能不能一對一地匹配感興趣的基因信息。在處理以前，檢查重複的基因ID和弄清楚重複的來源很是重要。咱們的基因註釋中包含28個匹配到不一樣染色體的基因（好比基因Gm1987關聯於染色體chr4和chr4_JH584294_random，小RNA Mir5098關聯於chr2，chr5，chr8，chr11和chr17）。 爲了處理重複的基因ID，咱們能夠合併來自多重匹配基因的全部染色體信息，好比將基因Gm1987分配到chr4 and chr4_JH584294_random，或選取其中一條染色體來表明具備重讀註釋的基因。爲了簡單起見，咱們選擇後者，保留每一個基因ID第一次出現的信息。
genes <- genes[!duplicated(genes$ENTREZID),]
在此例子中，註釋與數據對象中的基因順序是相同的。若是因爲缺失和／或從新排列基因ID致使其順序不一致，能夠用match來正確排序基因。而後將基因註釋的數據框加入數據對象，數據即被整潔地整理入一個DGEList對象中，它包含原始計數數據和相關的樣品信息和基因註釋。
x$genes <- genes x
## An object of class "DGEList" ## $samples ## files group lib.size norm.factors lane ## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004 ## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004 ## purep53 GSM1545538_purep53.txt Basal 57160817 1 L004 ## JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006 ## JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006 ## JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006 ## JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006 ## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008 ## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008 ## ## $counts ## Samples ## Tags 10_6_5_11 9_6_5_11 purep53 JMS8-2 JMS8-3 JMS8-4 JMS8-5 JMS9-P7c JMS9-P8c ## 497097 1 2 342 526 3 3 535 2 0 ## 100503874 0 0 5 6 0 0 5 0 0 ## 100038431 0 0 0 0 0 0 1 0 0 ## 19888 0 1 0 0 17 2 0 1 0 ## 20671 1 1 76 40 33 14 98 18 8 ## 27174 more rows ... ## ## $genes ## ENTREZID SYMBOL TXCHROM ## 1 497097 Xkr4 chr1 ## 2 100503874 Gm19938 <NA> ## 3 100038431 Gm10568 <NA> ## 4 19888 Rp1 chr1 ## 5 20671 Sox17 chr1 ## 27174 more rows ...
5數據預處理
5.1原始數據尺度轉換
因爲更深的測序總會產生更多的序列，在差別表達相關的分析中，咱們不多使用原始的序列數。在實踐中，咱們一般將原始的序列數進行歸一化，來消除測序深度所致使的差別。一般被使用的方法有基於序列的CPM（counts per million）、log-CPM、FPKM（fragments per kilobase of transcript per million），和基於轉錄本數目的RPKM（reads per kilobase of transcript per million）。
儘管CPM和log-CPM轉換並不像RPKM和FPKM那樣考慮到基因長度區別的影響，但在咱們的分析中常常會用到它們。雖然也可使用RPKM和FPKM，但CPM和log-CPM只使用計數矩陣便可計算，且已足以知足咱們所感興趣的比較的須要。假設不一樣條件之間剪接異構體（isoform）的使用沒有差異，差別表達分析研究同一基因在不一樣條件下的表達差別，而不是比較多個基因之間的表達或測定絕對錶達量。換而言之，基因長度在咱們感興趣的比較中始終不變，且任何觀測到的差別是來自於條件的變化而不是基因長度的變化。
在此處，使用edgeR中的cpm函數將原始計數轉換爲CPM和log-CPM值。其中，在進行對數轉換前會給全部基因的計數加上2，以免對零取對數，且可減少低表達基因之間的差別。 若是能夠提供基因長度信息，可使用edgeR中的rpkm函數計算RPKM值，就像計算CPM值那樣簡單。
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)
5.2刪除低表達基因
全部數據集中都混有表達的基因與不表達的基因。儘管咱們想要檢測在一種條件中表達但再另外一種條件中不表達的基因，也有一些基因在全部樣品中都不表達。實際上，這個數據集中19%的基因在全部九個樣品中的計數都是零。
table(rowSums(x$counts==0)==9)
## ## FALSE TRUE ## 22026 5153
咱們應當忽略在任何條件都沒有表達到具備生物學意義的程度的基因，這樣能夠減少基因的子集而保留感興趣的基因，且能夠減少下游在進行差別表達測試時的計算量。 經過檢查log-CPM值，能夠看出每一個樣本中很大一部分基因都是不表達或低表達的（以下圖中A部分所示）。在此使用CPM值爲1做爲標準（至關於log-CPM值爲0），將表達量高於此閾值的基因看做表達，反之看做不表達。咱們只保留在至少一個組（或者至少在整個實驗的三個樣品中）表達的基因用於下游分析。
儘管任何合理的值都能用做表達的閾值，通常咱們在分析中使用CPM值爲1，由於它對於大多數數據集都能很好地分隔表達的基因與不表達的基因。在這裏，CPM值爲1意味着若是一個基因在測序深度最低的樣品中（JMS9-P8c, 序列數量約2千萬）有至少20個計數，或者在測序深度最高的樣品中（JMS8-3，序列數量約7.6千萬）有至少76個計數，那麼它是「表達的」。若是測序的序列片斷是在外顯子水平而不是基因水平統計的，且／或實驗的測序深度很低，可能須要考慮更低的CPM閾值。
keep.exprs <- rowSums(cpm>1)>=3 x <- x[keep.exprs,, keep.lib.sizes=FALSE] dim(x)
## [1] 14165 9
使用這個標準，基因的數量減小到了開始時數量的大約一半（14165個基因，下圖B部分）。須要注意的是，從整個DGEList對象中取子集不只會刪除基因計數，同時也刪除了相關的基因信息。下方給出的是繪圖所用代碼。
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.21), las=2, main="", xlab="") title(main="A. Raw data", xlab="Log-cpm") abline(v=0, 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, prior.count=2) plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.21), las=2, main="", xlab="") title(main="B. Filtered data", xlab="Log-cpm") abline(v=0, lty=3) for (i in 2:nsamples){ den <- density(lcpm[,i]) lines(den$x, den$y, col=col[i], lwd=2) } legend("topright", samplenames, text.col=col, bty="n")
Figure 1: 每一個樣本過濾前的原始數據（A）和過濾後（B）的數據的log-CPM值密度。豎直虛線標出了過濾步驟中所用閾值即log-CPM爲0（至關於CPM值爲1）。
5.3歸一化基因表達分佈
在樣品製備或測序過程當中，不具有生物學意義的外部因素會影響單個樣品的表達。好比說，在實驗中第一批製備的樣品會整體上表達高於第二批製備的樣品。假設全部樣品的表達值的範圍和分佈都應當類似，須要進行歸一化來確保整個實驗中每一個樣本的表達分佈都類似。
密度圖和箱線圖等展現每一個樣品基因表達量分佈的圖表能夠用於判斷是否有樣品和其餘樣品分佈有差別。在此數據集中，全部樣品的log-CPM分佈都很類似（上圖B部分）。
儘管如此，咱們依然須要使用edgeR中的calcNormFactors函數，用TMM(Robinson and Oshlack 2010)方法進行歸一化。此處計算獲得的歸一化係數被用做文庫大小的縮放係數。當咱們使用DGEList對象時，這些歸一化係數被自動存儲在x$samples$norm.factors。對此數據集而言，TMM歸一化的做用比較溫和，這體如今全部的縮放因子都相對接近1。
x <- calcNormFactors(x, method = "TMM") x$samples$norm.factors
## [1] 0.896 1.035 1.044 1.041 1.032 0.922 0.984 1.083 0.979
爲了更好地可視化表現出歸一化的影響，咱們複製了數據並進行了調整，使得第一個樣品的計數減小到了其原始值的5%，而第二個樣品增大到了5倍。
x2 <- x
x2$samples$norm.factors <- 1 x2$counts[,1] <- ceiling(x2$counts[,1]*0.05) x2$counts[,2] <- x2$counts[,2]*5
下圖顯示了沒有通過歸一化的與通過了歸一化的數據的樣本的表達分佈，其中歸一化前的分佈顯然不一樣，而歸一化後比較類似。此處，第一個樣品的TMM縮放係數0.05很是小，而第二個樣品的縮放係數6.13很大，它們都並不接近1。
par(mfrow=c(1,2)) lcpm <- cpm(x2, log=TRUE, prior.count=2) 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.0547 6.1306 1.2293 1.1705 1.2149 1.0562 1.1459 1.2613 1.1170
lcpm <- cpm(x2, log=TRUE, prior.count=2) boxplot(lcpm, las=2, col=col, main="") title(main="B. Example: Normalised data",ylab="Log-cpm")
Figure 2: 樣例數據：log-CPM值的箱線圖展現了未經歸一化的數據（A）及歸一化後的數據（B）中每一個樣本的表達分佈。數據集通過調整，樣本1和2中的表達計數分別被縮放到其原始值的5%和500%。
5.4對樣本的無監督聚類
在咱們看來，用於檢查基因表達分析的最重要的探索性圖表之一即是MDS圖或其他相似的圖。這種圖表使用無監督聚類方法展現出了樣品間的類似性和不類似性，能讓咱們在進行正式的檢驗以前對於能檢測到多少差別表達基因有個大體概念。理想狀況下，樣本會在不一樣的實驗組內很好的聚類，且能夠鑑別出遠離所屬組的樣本，並追蹤偏差或額外方差的來源。若是存在技術重複，它們應當互相很是接近。
這樣的圖能夠用limma中的plotMDS函數繪製。第一個維度表示可以最好地分離樣品且解釋最大比例的方差的引導性的倍數變化（leading-fold-change），然後續的維度的影響更小，並與以前的維度正交。當實驗設計涉及到多個因子時，建議在多個維度上檢查每一個因子。若是在其中一些維度上樣本可按照某因子聚類，這說明該因子對於表達差別有影響，在線性模型中應當將其包括進去。反之，沒有或者僅有微小影響的因子在下游分析時應當被剔除。
在這個數據集中，能夠看出樣本在維度1和2能很好地按照實驗分組聚類，隨後在維度3按照測序道（樣品批次）分離（以下圖所示）。請記住，第一維度解釋了數據中最大比例的方差，須要注意到，當咱們向高維度移動，維度上的取值範圍會變小。
儘管全部樣本都按組聚類，在維度1上最大的轉錄差別出如今basal和LP以及basal和ML之間。所以，預期在basal樣品與其餘之間的成對比較中可以獲得大量的DE基因，而在ML和LP之間的比較中獲得的DE基因數量略少。在其餘的數據集中，不按照實驗組聚類的樣本可能在下游分析中只表現出較小的或不表現出差別表達。
爲繪製MDS圖，咱們爲不一樣的感興趣的因子賦予不一樣的色彩組合。維度1和2使用以細胞類型定義的色彩組合進行檢查。
維度3和4使用以測序泳道（批次）定義的色彩組合進行檢查。
lcpm <- cpm(x, log=TRUE, prior.count=2) par(mfrow=c(1,2)) col.group <- group levels(col.group) <- brewer.pal(nlevels(col.group), "Set1") col.group <- as.character(col.group) col.lane <- lane levels(col.lane) <- brewer.pal(nlevels(col.lane), "Set2") col.lane <- as.character(col.lane) plotMDS(lcpm, labels=group, col=col.group) title(main="A. Sample groups") plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4)) title(main="B. Sequencing lanes")
Figure 3: log-CPM值在維度1和2的MDS圖，以樣品分組上色並標記（A）和維度3和4的MDS圖，以測序道上色並標記（B）。圖中的距離對應於最主要的倍數變化（fold change），默認狀況下也就是前500個在每對樣品之間差別最大的基因的平均（均方根）log2倍數變化。
做爲另外一種選擇，Glimma包也提供了便於探索多個維度的交互式MDS圖。其中的glMDSPlot函數可生成一個html網頁（若是設置launch=TRUE，將會在瀏覽器中打開），其左側面板含有一張MDS圖，而右側面板包含一張展現了各個維度所解釋的方差比例的柱形圖。點擊柱形圖中的柱可切換MDS圖繪製時所使用的維度，且將鼠標懸浮於單個點上可顯示相應的樣本標籤。也可切換配色方案，以突顯不一樣細胞類型或測序泳道（批次）。此數據集的交互式MDS圖能夠從http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MDS-Plot.html看到。
glMDSPlot(lcpm, labels=paste(group, lane, sep="_"), groups=x$samples[,c(2,5)], launch=FALSE)
6差別表達分析
6.1建立設計矩陣和對比
在此研究中，咱們想知道哪些基因在咱們研究的三組細胞之間以不一樣水平表達。在咱們的分析中，假設基礎數據是正態分佈的，爲其擬合一個線性模型。在此以前，須要建立一個包含細胞類型以及測序泳道（批次）信息的設計矩陣。
design <- model.matrix(~0+group+lane) colnames(design) <- gsub("group", "", colnames(design)) design
## Basal LP ML laneL006 laneL008 ## 1 0 1 0 0 0 ## 2 0 0 1 0 0 ## 3 1 0 0 0 0 ## 4 1 0 0 1 0 ## 5 0 0 1 1 0 ## 6 0 1 0 1 0 ## 7 1 0 0 1 0 ## 8 0 0 1 0 1 ## 9 0 1 0 0 1 ## attr(,"assign") ## [1] 1 1 1 2 2 ## attr(,"contrasts") ## attr(,"contrasts")$group ## [1] "contr.treatment" ## ## attr(,"contrasts")$lane ## [1] "contr.treatment"
對於一個給定的實驗，一般有幾種等價的方法能夠建立一個合適的設計矩陣。 好比說，~0+group+lane去除了第一個因子group的截距，但第二個因子lane的截距被保留。 此外也可使用~group+lane，來自group和lane的截距均被保留。 此處的關鍵是理解如何解釋給定模型中估計獲得的係數。 咱們在此分析中選取第一種模型，由於在沒有group的截距的狀況下能更直截了當地設定模型中的對比。用於細胞羣之間成對比較的對比能夠在limma中用makeContrasts函數設定。
contr.matrix <- makeContrasts( BasalvsLP = Basal-LP, BasalvsML = Basal - ML, LPvsML = LP - ML, levels = colnames(design)) contr.matrix
## Contrasts ## Levels BasalvsLP BasalvsML LPvsML ## Basal 1 1 0 ## LP -1 0 1 ## ML 0 -1 -1 ## laneL006 0 0 0 ## laneL008 0 0 0
limma的線性模型方法的核心優點之一即是其適應任意實驗複雜程度的能力。簡單的設計，好比此工做流程中關於細胞類型和批次的實驗設計，直到更復雜的因子設計和含有交互做用項的模型，都可以被相對簡單地處理。當實驗或技術效應可被隨機效應模型（random effect model）模擬時，limma中的另外一種可能性是使用duplicateCorrelation函數來估計交互做用，這須要在此函數以及lmFit的線性建模步驟均指定一個block參數。
6.2從表達計數數據中刪除異方差
據顯示對於RNA-seq計數數據而言，當使用原始計數或當其被轉換爲log-CPM值時，方差並不獨立於均值(Law et al. 2014)。使用負二項分佈來模擬計數的方法假設均值與方差間具備二次的關係。在limma中，假設log-CPM值符合正態分佈，並使用由voom函數計算獲得的精確權重來調整均值與方差的關係，從而對log-CPM值進行線性建模。
當操做DGEList對象時，voom從x中自動提取文庫大小和歸一化因子，以此將原始計數轉換爲log-CPM值。在voom中，對於log-CPM值額外的歸一化能夠經過設定normalize.method參數來進行。
下圖左側展現了這個數據集log-CPM值的均值-方差關係。一般而言，方差是測序實驗中的技術差別和不一樣細胞類型的重複樣本之間的生物學差別的結合，而voom圖會顯示出一個在均值與方差之間遞減的趨勢。 生物學差別高的實驗一般會有更平坦的趨勢，其方差值在高表達處穩定。 生物學差別低的實驗更傾向於急劇降低的趨勢。
不只如此，voom圖也提供了對於上游所進行的過濾水平的可視化檢測。若是對於低表達基因的過濾不夠充分，在圖上表達低的一端，受到很是低的表達計數的影響，能夠觀察到方差水平的降低。若是觀察到了這種狀況，應當回到最初的過濾步驟並提升用於該數據集的表達閾值。
當前面觀察的MDS圖中具備明顯的樣本水平的差別時，能夠用voomWithQualityWeights函數來同時合併樣本水平的權重和voom(Liu et al. 2015)估算獲得的丰度相關的權重。關於此種狀況的例子參見Liu等(2016) (Liu et al. 2016)。
par(mfrow=c(1,2)) v <- voom(x, design, plot=TRUE) v
## An object of class "EList" ## $genes ## ENTREZID SYMBOL TXCHROM ## 1 497097 Xkr4 chr1 ## 6 27395 Mrpl15 chr1 ## 7 18777 Lypla1 chr1 ## 9 21399 Tcea1 chr1 ## 10 58175 Rgs20 chr1 ## 14160 more rows ... ## ## $targets ## files group lib.size norm.factors lane ## 10_6_5_11 GSM1545535_10_6_5_11.txt LP 29409426 0.896 L004 ## 9_6_5_11 GSM1545536_9_6_5_11.txt ML 36528591 1.035 L004 ## purep53 GSM1545538_purep53.txt Basal 59598629 1.044 L004 ## JMS8-2 GSM1545539_JMS8-2.txt Basal 53382070 1.041 L006 ## JMS8-3 GSM1545540_JMS8-3.txt ML 78175314 1.032 L006 ## JMS8-4 GSM1545541_JMS8-4.txt LP 55762781 0.922 L006 ## JMS8-5 GSM1545542_JMS8-5.txt Basal 54115150 0.984 L006 ## JMS9-P7c GSM1545544_JMS9-P7c.txt ML 23043111 1.083 L008 ## JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19525423 0.979 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.87 2.52 3.30 -4.48 -3.99 3.31 -3.20 -5.29 ## 27395 3.88 4.40 4.52 4.57 4.32 3.79 3.92 4.35 4.13 ## 18777 4.71 5.56 5.40 5.17 5.63 5.08 5.08 5.76 5.15 ## 21399 4.78 4.74 5.37 5.13 4.85 4.94 5.16 5.04 4.99 ## 58175 3.94 3.29 -1.77 -1.88 2.99 3.36 -2.11 3.14 3.52 ## 14160 more rows ... ## ## $weights ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] ## [1,] 1.18 1.18 20.53 20.98 1.77 1.22 21.13 1.18 1.18 ## [2,] 20.88 26.56 31.60 29.66 32.56 26.75 29.79 21.90 17.15 ## [3,] 28.00 33.70 34.85 34.46 35.15 33.55 34.52 31.44 25.23 ## [4,] 27.67 29.60 34.90 34.43 34.84 33.16 34.49 26.14 24.50 ## [5,] 19.74 18.66 3.18 2.63 24.19 24.01 2.65 13.15 14.35 ## 14160 more rows ... ## ## $design ## Basal LP ML laneL006 laneL008 ## 1 0 1 0 0 0 ## 2 0 0 1 0 0 ## 3 1 0 0 0 0 ## 4 1 0 0 1 0 ## 5 0 0 1 1 0 ## 6 0 1 0 1 0 ## 7 1 0 0 1 0 ## 8 0 0 1 0 1 ## 9 0 1 0 0 1 ## attr(,"assign") ## [1] 1 1 1 2 2 ## attr(,"contrasts") ## attr(,"contrasts")$group ## [1] "contr.treatment" ## ## attr(,"contrasts")$lane ## [1] "contr.treatment"
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")
Figure 4: 圖中繪製了每一個基因的均值（x軸）和方差（y軸），顯示了在該數據上使用voom前它們之間的相關性（左），以及當運用voom的精確權重後這種趨勢是如何消除的（右）。左側的圖是使用voom函數繪製的，它爲進行log-CPM轉換後的數據擬合線性模型從而提取殘差方差。而後，對方差取平方根（或對標準差取平方根），並相對每一個基因的平均表達做圖。均值經過平均計數加上0
5再進行log2轉換計算獲得。右側的圖使用plotSA繪製了log2殘差標準差與log-CPM均值的關係。平均log2殘差標準差由水平藍線標出。在這兩幅圖中，每一個黑點表示一個基因，紅線爲對這些點的擬合。
值得注意的是，DGEList對象中存儲的另外一個數據框，即基因和樣本水平信息所存儲之處，保留在了voom建立的EList對象v中。v$genes數據框等同於x$genes，v$targets等同於x$samples，而v$E中所儲存的表達值相似於x$counts，儘管它進行了尺度轉換。此外，voom的EList對象中還有一個精確權重的矩陣v$weights，而設計矩陣存儲於v$design。
6.3擬合線性模型以進行比較
limma的線性建模使用lmFit和contrasts.fit函數進行，它們原先是爲微陣列而設計的。這些函數不只能夠用於微陣列數據，也能夠用於RNA-seq數據。它們會單獨爲每一個基因的表達值擬合一個模型。而後，經過利用全部基因的信息來進行經驗貝葉斯調整，這樣能夠得到更精確的基因水平的變異程度估計(Smyth 2004)。下一圖爲此模型的殘差關於平均表達值的圖。從圖中能夠看出，方差再也不與表達水平均值相關。
6.4檢查DE基因數量
爲快速查看差別表達水平，顯著上調或下調的基因能夠彙總到一個表格中。顯著性的判斷使用校訂p值閾值的默認值5%。在basal與LP的表達水平之間的比較中，發現了4127個在basal中相較於LP下調的基因和4298個在basal中相較於LP上調的基因，即共8425個DE基因。在basal和ML之間發現了一共8510個DE基因（4338個下調基因和4172個上調基因），而在LP和ML中發現了一共5340個DE基因（2895個下調基因和2445個上調基因）。在包括basal細胞類型的比較中皆找到了大量的DE基因，這與咱們在MDS圖中觀察到的結果相吻合。
summary(decideTests(efit))
## BasalvsLP BasalvsML LPvsML ## Down 4127 4338 2895 ## NotSig 5740 5655 8825 ## Up 4298 4172 2445
一些研究中不只僅須要使用校訂p值閾值，更爲嚴格定義的顯著性可能須要差別倍數的對數（log-FCs）也高於某個最小值。treat方法(McCarthy and Smyth 2009)能夠按照對最小log-FC值的要求，使用通過經驗貝葉斯調整的t統計值計算p值。當咱們的檢驗要求基因的log-FC顯著大於1（等同於在本來的尺度上不一樣細胞類型之間差兩倍）時，差別表達基因的數量獲得了降低，basal與LP相比只有3135個DE基因，basal與ML相比只有3270個DE基因，LP與ML相比只有385個DE基因。
tfit <- treat(vfit, lfc=1) dt <- decideTests(tfit) summary(dt)
## BasalvsLP BasalvsML LPvsML ## Down 1417 1512 203 ## NotSig 11030 10895 13780 ## Up 1718 1758 182
在多種比較中皆差別表達的基因能夠從decideTests的結果中提取，其中的0表明不差別表達的基因，1表明上調的基因，-1表明下調的基因。共有2409個基因在basal和LP以及basal和ML的比較中都差別表達，其中的20個於下方列出。write.fit函數可用於將三個比較的結果提取並寫入到單個輸出文件。
de.common <- which(dt[,1]!=0 & dt[,2]!=0) length(de.common)
## [1] 2409
head(tfit$genes$SYMBOL[de.common], n=20)
## [1] "Xkr4" "Rgs20" "Cpa6" "Sulf1" "Eya1" ## [6] "Msc" "Sbspon" "Pi15" "Crispld1" "Kcnq5" ## [11] "Ptpn18" "Arhgef4" "2010300C02Rik" "Aff3" "Npas2" ## [16] "Tbc1d8" "Creg2" "Il1r1" "Il18r1" "Il18rap"
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
Figure 5: 韋恩圖展現了僅basal和LP（左）、僅basal和ML（右）的對比的DE基因數量，還有兩種對比中共同的DE基因數量（中）。在任何對比中均不差別表達的基因數量標於右下。
write.fit(tfit, dt, file="results.txt")
6.5從上到下檢查單個DE基因
使用topTreat函數能夠列舉出使用treat獲得的結果中靠前的DE基因（對於eBayes的結果可使用topTable函數）。默認狀況下，topTreat將基因按照校訂p值從小到大排列，併爲每一個基因給出相關的基因信息、log-FC、平均log-CPM、校訂t值、原始及通過多重假設檢驗校訂的p值。列出前多少個基因的數量可由用戶指定，若是設爲n=Inf則會包括全部的基因。基因Cldn7和Rasef在basal與LP和basal於ML的比較中都位於DE基因的前幾名。
basal.vs.lp <- topTreat(tfit, coef=1, n=Inf) basal.vs.ml <- topTreat(tfit, coef=2, n=Inf) head(basal.vs.lp)
## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val ## 12759 12759 Clu chr14 -5.44 8.86 -33.4 3.99e-10 2.7e-06 ## 53624 53624 Cldn7 chr11 -5.51 6.30 -32.9 4.50e-10 2.7e-06 ## 242505 242505 Rasef chr4 -5.92 5.12 -31.8 6.06e-10 2.7e-06 ## 67451 67451 Pkp2 chr16 -5.72 4.42 -30.7 8.01e-10 2.7e-06 ## 228543 228543 Rhov chr2 -6.25 5.49 -29.5 1.11e-09 2.7e-06 ## 70350 70350 Basp1 chr15 -6.07 5.25 -28.6 1.38e-09 2.7e-06
head(basal.vs.ml)
## ENTREZID SYMBOL TXCHROM logFC AveExpr t P.Value adj.P.Val ## 242505 242505 Rasef chr4 -6.51 5.12 -35.5 2.57e-10 1.92e-06 ## 53624 53624 Cldn7 chr11 -5.47 6.30 -32.5 4.98e-10 1.92e-06 ## 12521 12521 Cd82 chr2 -4.67 7.07 -31.8 5.80e-10 1.92e-06 ## 71740 71740 Nectin4 chr1 -5.56 5.17 -31.3 6.76e-10 1.92e-06 ## 20661 20661 Sort1 chr3 -4.91 6.71 -31.2 6.76e-10 1.92e-06 ## 15375 15375 Foxa1 chr12 -5.75 5.63 -28.3 1.49e-09 2.28e-06
6.6差別表達結果的實用圖形表示
爲可視化地總結全部基因的結果，可以使用plotMD函數繪製均值-差別（MD）圖，其中展現了線性模型擬合所獲得的log-FC與log-CPM平均值間的關係，而差別表達的基因會被重點標出。
plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1], xlim=c(-8,13))
Glimma的glMDPlot函數提供了交互式的均值-差別圖，拓展了這種圖表的功能性。 此函數的輸出爲一個html頁面，左側面板爲結果的總結性圖表（與plotMD的輸出相似），右側面板包含各個樣本的log-CPM值，下方爲結果的表格。 這種交互式展現容許用戶使用基因名搜索特定基因，而這在R統計圖中是作不到的。 glMDPlot函數不只限於均值-差別圖，其默認版本容許讀入數據框，而用戶能夠選擇在左側面板繪圖所用的列。
glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1], side.main="ENTREZID", counts=lcpm, groups=group, launch=FALSE)
使用Glimma生成的均值-差別圖。左側面板顯示了總結性數據（log-FC與log-CPM值的關係），其中選中的基因在每一個樣本中的數值顯示於右側面板。下方爲結果的表格，其搜索框使用戶得以使用可行的註釋信息查找某個特定基因，如基因符號Clu。
上方指令生成的均值-差別圖能夠在線訪問（詳見http://bioinf.wehi.edu.au/folders/limmaWorkflow/glimma-plots/MD-Plot.html）。 Glimma提供的交互性使得單個圖形窗口內可以呈現出額外的信息。 Glimma是以R和Javascript實現的，使用R代碼生成數據，並在以後使用Javascript庫D3（https://d3js.org）轉換爲圖形，使用Bootstrap庫處理界面並生成互動性可搜索的表格的數據表。 這使得圖表能夠在任何現代的瀏覽器中查看，對於從Rmarkdown分析報告中將其做爲關聯文件而附加而言十分方便。
前文所展現的圖表中，一些展現了在任意一個條件下表達的全部基因（好比共同DE基因的韋恩圖或均值-差別圖），而另外一些展現單獨的基因（交互性均值-差別圖右邊面板中所展現的log-CPM值）。而熱圖使用戶得以查看一部分基因的表達。這對於查看單個組或樣本的表達頗有用，而不至於在關注於單個基因時失去對於研究總體的注意，也不會形成因爲上千個基因所取平均值而致使的失去分辨率。
使用gplots包的heatmap.2函數，咱們爲basal與LP的對照中前100個DE基因（按調整p值排序）繪製了一幅熱圖。熱圖中正確地將樣本按照細胞類型聚類，並從新排序了基因，造成了表達類似的塊狀。從熱圖中，咱們觀察到對於basal與LP之間的前100個DE基因，ML和LP樣本的表達很是類似。
library(gplots) basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100] i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes) mycol <- colorpanel(1000,"blue","white","red") heatmap.2(lcpm[i,], scale="row", labRow=v$genes$SYMBOL[i], labCol=group, col=mycol, trace="none", density.info="none", margin=c(8,6), lhei=c(2,10), dendrogram="column")
Figure 6: 在basal和LP的對比中前100個DE基因log-CPM值的熱圖。通過縮放調整後，每一個基因（每行）的表達均值爲0，而且標準差爲1。給定基因相對高表達的樣本被標記爲紅色，相對低表達的樣本被標記爲藍色。淺色和白色表明中等表達水平的基因。樣本和基因已經過分層聚類的方法從新排序。圖中顯示有樣本聚類的樹狀圖。
7使用camera的基因集檢驗
在這次分析的最後，咱們要進行一些基因集檢驗。爲此，咱們將camera方法(Wu and Smyth 2012)應用於Broad Institute的MSigDB c2中的(Subramanian et al. 2005)中適應小鼠的c2基因表達特徵，這可從http://bioinf.wehi.edu.au/software/MSigDB/以RData對象格式獲取。 此外，對於人類和小鼠，來自MSigDB的其餘有用的基因集也可今後網站獲取，好比表明性（hallmark）基因集。C2基因集的內容收集自在線數據庫、出版物以及該領域專家，而標誌基因集的內容來自MSigDB，從而對於生物狀態或過程的定義足夠明確。
load(system.file("extdata", "mouse_c2_v5p1.rda", package = "RNAseq123")) idx <- ids2indices(Mm.c2,id=rownames(v)) cam.BasalvsLP <- camera(v,idx,design,contrast=contr.matrix[,1]) head(cam.BasalvsLP,5)
## NGenes Direction PValue FDR ## LIM_MAMMARY_STEM_CELL_UP 739 Up 1.13e-18 5.36e-15 ## LIM_MAMMARY_STEM_CELL_DN 630 Down 1.57e-15 3.71e-12 ## ROSTY_CERVICAL_CANCER_PROLIFERATION_CLUSTER 163 Up 1.44e-13 2.26e-10 ## SOTIRIOU_BREAST_CANCER_GRADE_1_VS_3_UP 183 Up 2.18e-13 2.58e-10 ## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 87 Down 6.73e-13 6.36e-10
cam.BasalvsML <- camera(v,idx,design,contrast=contr.matrix[,2]) head(cam.BasalvsML,5)
## NGenes Direction PValue FDR ## LIM_MAMMARY_STEM_CELL_UP 739 Up 5.09e-23 2.40e-19 ## LIM_MAMMARY_STEM_CELL_DN 630 Down 5.13e-19 1.21e-15 ## LIM_MAMMARY_LUMINAL_MATURE_DN 166 Up 8.88e-16 1.40e-12 ## LIM_MAMMARY_LUMINAL_MATURE_UP 180 Down 6.29e-13 7.43e-10 ## ROSTY_CERVICAL_CANCER_PROLIFERATION_CLUSTER 163 Up 1.68e-12 1.59e-09
cam.LPvsML <- camera(v,idx,design,contrast=contr.matrix[,3]) head(cam.LPvsML,5)
## NGenes Direction PValue FDR ## LIM_MAMMARY_LUMINAL_MATURE_UP 180 Down 8.50e-14 3.40e-10 ## LIM_MAMMARY_LUMINAL_MATURE_DN 166 Up 1.44e-13 3.40e-10 ## LIM_MAMMARY_LUMINAL_PROGENITOR_UP 87 Up 3.84e-11 6.05e-08 ## REACTOME_RESPIRATORY_ELECTRON_TRANSPORT 91 Down 2.66e-08 3.14e-05 ## NABA_CORE_MATRISOME 222 Up 4.43e-08 4.19e-05
camera函數經過比較假設檢驗來評估一個給定基因集中的基因是否相對於不在集內的基於是言在差別表達基因的排序中更靠前。 它使用limma的線性模型框架，並同時採用設計矩陣和對比矩陣（若是有），且在測試的過程當中會適應來自voom的觀測水平權重。 在經過方差膨脹因子（variance inflation factor）調整獲得的基因集檢驗統計值的方差後，取決於基因間相關性（默認設定爲0.01，但也可經過數據估計）和基因集的規模，將會返回p值，且根據多重假設檢驗進行了校訂。
此實驗是等同於Lim等人(2010)(Lim et al. 2010)的RNA-seq，而他們使用Illumina微陣列分析了相同的分選細胞羣，所以當咱們看到該早期文獻中的基因表達特徵出如今每種對比的列表頂部時大可放心。在LP和ML的對比中，咱們爲Lim等人（2010）的成熟管腔基因集（上調及下調）繪製了條碼圖（barcodeplot）。須要注意的是，因爲咱們的對比是將LP與ML相比而不是相反，這些基因集的方向在咱們的數據集中是反過來的（若是將對比反過來，基因集的方向將會與對比一致）。
barcodeplot(efit$t[,3], index=idx$LIM_MAMMARY_LUMINAL_MATURE_UP, index2=idx$LIM_MAMMARY_LUMINAL_MATURE_DN, main="LPvsML")
Figure 7: LIM_MAMMARY_LUMINAL_MATURE_UP （紅色條形，圖表上方）和LIM_MAMMARY_LUMINAL_MATURE_DN（藍色條形，圖表下方）基因集在LP和ML的對比中的條碼圖，每一個基因集都有一條富集線展現了豎直條形在圖表每部分的相對富集程度。Lim等人的實驗(Lim et al
2010)很是相似於咱們的，用了相同的分選方式來獲取不一樣的細胞羣，只是他們使用的是微陣列而不是RNA-seq來測定基因表達。須要注意的是，上調基因集發生下調而下調基因集發生上調的逆相關性來自於對比的設定方式（LP相比於ML），若是將其對調，方向性將會吻合。
limma也有其餘的基因集檢驗，好比mroast(Wu et al. 2010)的自包含檢驗。雖然camera很是適合檢驗基因集的大型數據庫並觀察其中哪些相對於其餘的在排序上位次更高（如前文所示），自包含檢驗更善於集中檢驗一個或少個選中的集合是否自己差別表達。換句話說，camera更適用於搜尋具備意義的基因集，而mroast測試的是已經肯定有意義的基因集的顯著性。
8使用到的軟件和代碼
此RNA-seq工做流程使用了Bioconductor項目3.4版本中的多個軟件包，運行於R 3.3.0或更高版本。除了本文中重點提到的軟件（limma、Glimma以及edgeR），亦須要一些其餘軟件包，包括gplots和RColorBrewer還有基因註釋包Mus.musculus。 此文檔使用knitr編譯。全部用到的包的版本號以下所示。進行此分析所需代碼都可在Bioconductor工做流程包RNAseq123從http://www.bioconductor.org/help/workflows/RNAseq123/獲取。
sessionInfo()
## R Under development (unstable) (2018-10-16 r75448) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 18.04.1 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] 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.25.0 GenomicFeatures_1.35.0 ## [9] GenomicRanges_1.35.0 GenomeInfoDb_1.19.0 ## [11] AnnotationDbi_1.45.0 IRanges_2.17.0 ## [13] S4Vectors_0.21.0 Biobase_2.43.0 ## [15] BiocGenerics_0.29.1 edgeR_3.25.0 ## [17] Glimma_1.11.0 limma_3.39.1 ## [19] knitr_1.20 BiocStyle_2.11.0 ## ## loaded via a namespace (and not attached): ## [1] httr_1.3.1 bit64_0.9-7 jsonlite_1.5 ## [4] R.utils_2.7.0 gtools_3.8.1 assertthat_0.2.0 ## [7] BiocManager_1.30.3 highr_0.7 RBGL_1.59.0 ## [10] blob_1.1.1 GenomeInfoDbData_1.2.0 Rsamtools_1.35.0 ## [13] yaml_2.2.0 progress_1.2.0 RSQLite_2.1.1 ## [16] backports_1.1.2 lattice_0.20-35 digest_0.6.18 ## [19] XVector_0.23.0 htmltools_0.3.6 Matrix_1.2-15 ## [22] R.oo_1.22.0 XML_3.98-1.16 pkgconfig_2.0.2 ## [25] biomaRt_2.39.0 bookdown_0.7 zlibbioc_1.29.0 ## [28] gdata_2.18.0 BiocParallel_1.17.0 SummarizedExperiment_1.13.0 ## [31] magrittr_1.5 crayon_1.3.4 memoise_1.1.0 ## [34] evaluate_0.12 R.methodsS3_1.7.1 graph_1.61.0 ## [37] tools_3.6.0 prettyunits_1.0.2 hms_0.4.2 ## [40] matrixStats_0.54.0 stringr_1.3.1 locfit_1.5-9.1 ## [43] DelayedArray_0.9.0 Biostrings_2.51.0 compiler_3.6.0 ## [46] caTools_1.17.1.1 rlang_0.3.0.1 grid_3.6.0 ## [49] RCurl_1.95-4.11 bitops_1.0-6 rmarkdown_1.10 ## [52] DBI_1.0.0 R6_2.3.0 GenomicAlignments_1.19.0 ## [55] rtracklayer_1.43.0 bit_1.1-14 rprojroot_1.3-2 ## [58] KernSmooth_2.23-15 stringi_1.2.4 Rcpp_0.12.19 ## [61] xfun_0.4
參考文獻
Bioconductor Core Team. 2016a. Homo.sapiens: Annotation package for the Homo.sapiens object. https://bioconductor.org/packages/release/data/annotation/html/Homo.sapiens.html.
———. 2016b. Mus.musculus: Annotation package for the Mus.musculus object. https://bioconductor.org/packages/release/data/annotation/html/Mus.musculus.html.
Durinck, S., Y. Moreau, A. Kasprzyk, S. Davis, B. De Moor, A. Brazma, and W. Huber. 2005. 「BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis.」 Bioinformatics 21:3439–40.
Durinck, S., P. Spellman, E. Birney, and W. Huber. 2009. 「Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt.」 Nature Protocols 4:1184–91.
Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, et al. 2015. 「Orchestrating High-Throughput Genomic Analysis with Bioconductor.」 Nature Methods 12 (2):115–21. http://www.nature.com/nmeth/journal/v12/n2/full/nmeth.3252.html.
Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. 「Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-seq Read Counts.」 Genome Biology 15:R29.
Liao, Y., G. K. Smyth, and W. Shi. 2013. 「The Subread Aligner: Fast, Accurate and Scalable Read Mapping by Seed-and-Vote.」 Nucleic Acids Res 41 (10):e108.
———. 2014. 「featureCounts: an Efficient General-Purpose Program for Assigning Sequence Reads to Genomic Features.」 Bioinformatics 30 (7):923–30.
Lim, E., D. Wu, B. Pal, T. Bouras, M. L. Asselin-Labat, F. Vaillant, H. Yagita, G. J. Lindeman, G. K. Smyth, and J. E. Visvader. 2010. 「Transcriptome analyses of mouse and human mammary cell subpopulations reveal multiple conserved genes and pathways.」 Breast Cancer Research 12 (2):R21.
Liu, R., K. Chen, N. Jansz, M. E. Blewitt, and M. E. Ritchie. 2016. 「Transcriptional Profiling of the Epigenetic Regulator Smchd1.」 Genomics Data 7:144–7.
Liu, R., A. Z. Holik, S. Su, N. Jansz, K. Chen, H. S. Leong, M. E. Blewitt, M. L. Asselin-Labat, G. K. Smyth, and M. E. Ritchie. 2015. 「Why weight? Combining voom with estimates of sample quality improves power in RNA-seq analyses.」 Nucleic Acids Res 43:e97.
McCarthy, D. J., and G. K. Smyth. 2009. 「Testing significance relative to a fold-change threshold is a TREAT.」 Bioinformatics 25:765–71.
Ritchie, M. E., B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. 2015. 「limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies.」 Nucleic Acids Res 43 (7):e47.
Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. 「edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.」 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., and M. E. Ritchie. 2016. Glimma: Interactive HTML graphics for RNA-seq data. https://bioconductor.org/packages/devel/bioc/html/Glimma.html.
Subramanian, A., P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, et al. 2005. 「Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.」 Proc Natl Acad Sci U S A 102 (43):15545–50.
Wu, D., E. Lim, F. Vaillant, M. L. Asselin-Labat, J. E. Visvader, and G. K. Smyth. 2010. 「ROAST: rotation gene set tests for complex microarray experiments.」 Bioinformatics 26 (17):2176–82.
Wu, D., and G. K. Smyth. 2012. 「Camera: a competitive gene set test accounting for inter-gene correlation.」 Nucleic Acids Res 40 (17):e133.
- 1. 直通車上10分 易如反掌
- 2. 生物信息學入門 使用 RNAseq counts數據進行差異表達分析（DEG）——edgeR 算法 數據 代碼 結果解讀
- 3. 如何使用反射來分析PHP的數據結構
- 4. Python探索性數據分析，這樣才容易掌握！
- 5. cgmodel簡易數據分析
- 6. .Net反射機制分析和使用
- 7. RNAseq數據處理實驗記錄ydy_sir
- 8. 如何分析用戶反饋數據? | 下
- 9. 數據分析和數據分析師
- 10. 如何才能更好使用數據分析工具分析數據？
- 更多相關文章...
- • XSD 如何使用? - XML Schema 教程
- • 使用netwox構造IP數據包 - TCP/IP教程
- • Flink 數據傳輸及反壓詳解
- • TiDB 在摩拜單車在線數據業務的應用和實踐
-
每一个你不满意的现在，都有一个你没有努力的曾经。