一、簡介

limma應(yīng)用于RNA-seq數(shù)據(jù)時，read counts被轉(zhuǎn)換為log2-counts-per-million（logCPM）?？梢杂袃煞N方式對均值-方差的關(guān)系（mean-variance relationship）進行建模：
（1）voom：精確權(quán)重（precision weights）；
（2）limma-trend：經(jīng)驗貝葉斯先驗趨勢（empirical Bayes prior trend）。
在這兩種情況下，RNA-seq數(shù)據(jù)可以被作為微陣列數(shù)據(jù)（microarray data）進行分析。這意味著limma包中的任何線性建?；蚧蚣瘻y試方法都可以應(yīng)用于RNA-seq數(shù)據(jù)。

limma使用線性模型（linear models）的方法來分析設(shè)計的微陣列實驗。該方法需要指定一個或兩個矩陣：
（1）第一種是設(shè)計矩陣（design matrix），它表明了每個陣列使用了哪些RNA樣本。
（2）第二種是對比矩陣（contrast matrix），它指定了希望在RNA樣本之間進行哪些比較。
對于非常簡單的實驗，可能不需要指定對比矩陣。

對于統(tǒng)計分析和評估差異表達，limma使用了一種經(jīng)驗貝葉斯方法（empirical Bayes method）來調(diào)整估計的倍數(shù)變化的標準誤差（moderate the standard errors of the estimated log-fold changes）。這導(dǎo)致了更穩(wěn)定的推理和更高的功效，特別是對于樣本量少的實驗。

基本過程：

• RNA-seq數(shù)據(jù)制作計數(shù)矩陣、樣本分組信息、基因注釋信息。
• 構(gòu)建DGEList對象，過濾filterByExpr()，TMM標準化calcNormFactors()。
• 準備設(shè)計矩陣model.matrix()和對比矩陣makeContrasts()。
• cpm()過濾和標準化后的數(shù)據(jù)轉(zhuǎn)換為logCPM。
• voom()轉(zhuǎn)換過濾和標準化后的數(shù)據(jù)。eBayes()和treat()不需要trend=T。
• lmFit()通過為每一個基因擬合線性模型來估計fold changes和standard errors。
• contrasts.fit()將原始模型重新定向到對比模型，得到新的系數(shù)和標準誤差。
• eBayes()使用trend=TRUE對標準誤差進行經(jīng)驗貝葉斯平滑，計算每個對比中每個基因的moderated t-statistic和log-odds。
• treat()使用trend=TRUE對log-fold-changes大于某一非零閾值進行檢驗。
• topTable()給出一個最有可能在給定對比下差異表達的基因列表。
• topTableF()給出一個最有可能在給定的一組對比中差異表達的基因列表。
• decideTests()展示所有對比結(jié)果。

二、RNA-seq數(shù)據(jù)的前處理

（1）制作計數(shù)矩陣

> head(countdata,3) 
          WTF1 WTF2 WTF3 WTM1 WTM2 WTM3 MF1 MF2 MF3 MF4 MF5 MMF1 MMF2
gene17427    0   46  210    0   26    0  40   0   0   0   0    0    0
gene17428   50  322  334   0  345  384  284  199  218  11   1    2    1
gene17429    0    0    0    6    0    0   0   0   0   0   5    0    0

（2）樣本分組信息

> Group
 [1] WTF WTF WTF WTM WTM WTM MF  MF  MF  MF  MF  MMF MMF
Levels: WTF WTM MF MMF

（3）基因注釋信息

> head(anno)
          chromosome ref_gene_name
gene17427          4       CR32011
gene17428          4       CR32010
gene17429          4       CR32009
gene17490          4       CR44028
gene17501          4       CR33797
gene17507          4       CR43361

（4）使用edgeR，將計數(shù)矩陣轉(zhuǎn)換為DGEList對象

> y <- DGEList(counts=countdata) 
> y 
An object of class "DGEList"
$counts
          WTF1 WTF2 WTF3 WTM1 WTM2 WTM3 MF1 MF2 MF3 MF4 MF5 MMF1 MMF2
gene17427    0   46  210    0   26    0  40   0   0   0   0    0    0
gene17428   50  322  334   0  345   384  284  199  218  11   1    2    1
gene17429    0    0    0    6    0    0   0   0   0   0   5    0    0
gene17490    4    0    3    0    0    0   0   3   0   0   0    0    0
gene17501    2    2    3    1    0    3   4   0   3   0   0    0    0
34367 more rows ...
$samples
     group lib.size norm.factors
WTF1     1 26485372            1
WTF2     1 21708533            1
WTF3     1 26038158            1
WTM1     1 14282144            1
WTM2     1 27198208            1
8 more rows ...

> y$samples$group <- Group #加入分組信息
> y$genes <- anno #加入基因注釋

（5）移除表達量非常低的基因

> #保留在至少在2個樣品中表達量大于min.CPM的基因
> keep <- rowSums(cpm(y)>1) >= 2
> yf <- y[keep,,keep.lib.sizes=FALSE] #重置了庫的大小lib.size
> y$samples$lib.size - yf$samples$lib.size
 [1]  26499  39607  80456  25102 125345 105390  63741  24067  46952
[10]  37984  44858  35573  34754
> nrow(y) - nrow(yf) #減少的基因數(shù)
[1] 11844

> #另一種過濾方法filterByExpr
> keep <- filterByExpr(y,group=Group)
> yf <- y[keep,,keep.lib.sizes=FALSE]
> y$samples$lib.size - yf$samples$lib.size
 [1] 11532 22890 58275 10496 88993 68069 35694 10642 24187 14905
[11] 22093 16030 19544
> nrow(y) - nrow(yf)
[1] 9521

用法：

filterByExpr(y, design = NULL, group = NULL, lib.size = NULL, ...)
filterByExpr(y, design = NULL, group = NULL, lib.size = NULL, min.count = 10, min.total.count = 15, large.n = 10, min.prop = 0.7, ...)

• y：DGEList對象
• design：設(shè)計矩陣，如果group不等于NULL時design會被忽略
• group：表示分組的向量或因子
• lib.size：默認為colSums(y)
• min.count：選取某基因表達量最小的樣本，這個樣品的最低計數(shù)要求。
• min.total.count：最小的總計數(shù)要求。
• large.n：每一個group中，樣本數(shù)多少會被認為是大樣本量。
• min.prop：最小分組中樣本數(shù)量的比例。

該方法保留在n個樣本中有大于k的CPM的基因。k由min.count 和lib.size決定，n由設(shè)計矩陣決定。n是最小的分組的樣本量。如果所有的分組樣本量都大于large.n，這就會稍微放松，但是n總是大于最小分組的min.prop（70%）。此外，每個保留的基因都需要在所有樣本中至少有min.total.count的計數(shù)。

（6）計算CPM和logCPM

> #cpmyf <- cpm(yf)
> lcpmyf <- cpm(yf,log=T) 

用法：

cpm(y,normalized.lib.sizes=T,log=F,prior.count=2)

• normalized.lib.sizes=T：使用標準化后的library sizes。effective library size = y$samples$lib.size * y$samples$norm.factors
• log=T：計算log2CPM
• prior.count=2：取對數(shù)前在每個觀測值上加一個數(shù)，避免對0取對數(shù)。這里使用prior count可以降低low counts取對數(shù)的方差。
• rpkm()：會尋找y$genes$Length或length

（7）箱線圖

> col.group <- Group
> levels(col.group) <- brewer.pal(nlevels(col.group),"Set1")
> col.group <- as.character(col.group)
> col.group
 [1] "#E41A1C" "#E41A1C" "#E41A1C" "#377EB8" "#377EB8" "#377EB8"
 [7] "#4DAF4A" "#4DAF4A" "#4DAF4A" "#4DAF4A" "#4DAF4A" "#984EA3"
[13] "#984EA3"
> boxplot(lcpmyf,las=2,col=col.group)

boxplot.JPG

（8）TMM標準化(trimmed mean of M-values)

> yf <- calcNormFactors(yf,method = "TMM") #計算歸一化系數(shù)
> #標準化后raw read counts值并沒改變，只有norm.factors變了
> yf$samples$norm.factors
 [1] 0.9027210 1.3505878 1.3102042 0.9202737 1.3636074 1.4155932
 [7] 0.9726943 1.0979219 0.9623963 0.9108200 0.7308098 0.7473150
[13] 0.6892828

> ##獲得標準化后的數(shù)據(jù)
> #cpmyf <- cpm(yf) #已經(jīng)考慮歸一化系數(shù)了，等于cpmyf/norm.factors
> #t(t(yf$counts)/yf$samples$lib.size/yf$samples$norm.factors)*1000000
> lcpmyf <- cpm(yf,log=T)
> boxplot(lcpmyf,las=2,col=col.group)

三、設(shè)計矩陣和對比矩陣

需要兩種矩陣：

• 設(shè)計矩陣（design matrix），表明了每個陣列使用了哪些RNA樣本。
• 對比矩陣（contrast matrix），指定了希望在RNA樣本之間進行哪些比較。對于非常簡單的實驗，可能不需要指定對比矩陣。

有兩種方法創(chuàng)建設(shè)計矩陣：

• create a design matrix which includes a coefficient for the mutant vs wild type difference, 將對照組作為截距，有一個系數(shù)表示實驗組與對照組的差異。
• create a design matrix which includes separate coefficients for wild type and mutant mice and then extract the difference as a contrast.截距為0，對照組和實驗組分別有自己的系數(shù)。

（1）對照組作為截距

> Group
 [1] WTF WTF WTF WTM WTM WTM MF  MF  MF  MF  MF  MMF MMF
Levels: WTF WTM MF MMF
> design <- model.matrix(~Group)
> colnames(design) <- c("WTF","WTMvsWTF","MFvsWTF","MMFvsWTF")
> design
   WTF WTMvsWTF MFvsWTF MMFvsWTF
1    1        0       0        0
2    1        0       0        0
3    1        0       0        0
4    1        1       0        0
5    1        1       0        0
6    1        1       0        0
7    1        0       1        0
8    1        0       1        0
9    1        0       1        0
10   1        0       1        0
11   1        0       1        0
12   1        0       0        1
13   1        0       0        1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Group
[1] "contr.treatment"

（2）截距為0，對照組和實驗組分別有自己的系數(shù)

> design <- model.matrix(~0+Group)
> colnames(design) <- c("WTF","WTM","MF","MMF")
> design
   WTF WTM MF MMF
1    1   0  0   0
2    1   0  0   0
3    1   0  0   0
4    0   1  0   0
5    0   1  0   0
6    0   1  0   0
7    0   0  1   0
8    0   0  1   0
9    0   0  1   0
10   0   0  1   0
11   0   0  1   0
12   0   0  0   1
13   0   0  0   1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Group
[1] "contr.treatment"

（3）對比矩陣

> cont.matrix <- makeContrasts(WTMvsWTF = WTM-WTF, MFvsWTF = MF-WTF, 
+                         MMFvsWTF = MMF-WTF, MMFvsMF = MMF-MF,
+                         levels = colnames(design))
> cont.matrix
      Contrasts
Levels WTMvsWTF MFvsWTF MMFvsWTF MMFvsMF
   WTF       -1      -1       -1       0
   WTM        1       0        0       0
   MF         0       1        0      -1
   MMF        0       0        1       1

用法：

makeContrasts(..., contrasts=NULL, levels)

• contrasts：指定對比的字符型向量。
• levels：字符向量或因子，給出所需對比的參數(shù)的名稱names，或設(shè)計矩陣等以參數(shù)名稱為列名的對象。參數(shù)通常是線性模型擬合的系數(shù)。
• 此函數(shù)的輸出通常用作contrasts.fit的輸入。

> makeContrasts(B-A,C-B,C-A,levels=c("A","B","C"))
      Contrasts
Levels B - A C - B C - A
     A    -1     0    -1
     B     1    -1     0
     C     0     1     1

（4）成對樣本的比較

成對樣本的比較出現(xiàn)在我們比較兩種處理條件，并且每個樣本給予一個處理自然地與一個給予另一種處理的特定的樣本配對的情況。與這種情況相關(guān)的經(jīng)典檢驗是配對t檢驗。Paired samples occur when we compare two treatments and each sample given one treatment is naturally paired with a particular sample given the other treatment. This is a special case of blocking with blocks of size two. The classical test associated with this situation is the paired t-test.

上述對于成對樣本的比較也可以用于存在batch effects或?qū)嶒瀼牟煌腷locks進行的情況。處理條件僅在每一個模塊內(nèi)進行比較。
The above approach used for paired samples can be applied in any situation where there are batch effects or where the experiment has been conducted in blocks. The treatments can be adjusted for differences between the blocks by using a model formula of the form:

> design <- model.matrix(~Block+Treatment) 

In this type of analysis, the treatments are compared only within each block.

四、差異表達分析

limma最初設(shè)計用于分析微陣列數(shù)據(jù)，但目前已擴展到RNA-seq數(shù)據(jù)、DNA甲基化等。limma應(yīng)用于RNA-seq數(shù)據(jù)，也是基于raw count的定量方式，但是它并不提供歸一化算法，在官方手冊中推薦采用edgeR的TMM歸一化算法。read counts被轉(zhuǎn)換為log2-counts-per-million（logCPM），并假設(shè)log-CPM值符合正態(tài)分布?？梢杂袃煞N方式調(diào)整均值-方差的關(guān)系（mean-variance relationship），從而對log-CPM值進行線性建模：

• voom：精確權(quán)重（precision weights）；
• limma-trend：經(jīng)驗貝葉斯先驗趨勢（empirical Bayes prior trend）。

在這兩種情況下，RNA-seq數(shù)據(jù)可以被作為微陣列數(shù)據(jù)（microarray data）進行分析。這意味著limma包中的任何線性建?；蚧蚣瘻y試方法都可以應(yīng)用于RNA-seq數(shù)據(jù)。

如果測序深度在RNA樣本之間差異不大（最大library size和最小library size相差3倍以內(nèi)），那么使用limma-trend是最簡單和最穩(wěn)健的方法。將計數(shù)數(shù)據(jù)轉(zhuǎn)化為logCPM，然后logCPM值可以在任何標準的limma pipeline中進行分析。在運行eBayes或treat時使用trend=TRUE。
當library sizes在樣本之間差距較大時，voom方法在理論上比limma-trend更有效。

（1）limma-trend

• 數(shù)據(jù)轉(zhuǎn)換為logCPM

> lcpmyf <- cpm(yf,log=T)

• 設(shè)計矩陣

> design <- model.matrix(~0+Group)
> colnames(design) <- gsub("Group","",colnames(design))
> rownames(design) <- colnames(yf)

• 對比矩陣

> cont.matrix <- makeContrasts(WTMvsWTF = WTM-WTF, MFvsWTF = MF-WTF, 
+                              MMFvsWTF = MMF-WTF, MMFvsMF = MMF-MF,
+                              levels = colnames(design))
> cont.matrix
      Contrasts
Levels WTMvsWTF MFvsWTF MMFvsWTF MMFvsMF
   WTF       -1      -1       -1       0
   WTM        1       0        0       0
   MF         0       1        0      -1
   MMF        0       0        1       1

• 擬合線性模型lmFit

> fit <- lmFit(lcpmyf, design) 

• 轉(zhuǎn)換為對比模型contrasts.fit

> fit2 <- contrasts.fit(fit, cont.matrix) 

• 經(jīng)驗貝葉斯平滑eBayes

> fit3 <- eBayes(fit2, trend=TRUE) 

• 提取某組比較的結(jié)果

> topTable(fit3, coef=ncol(design))
             logFC   AveExpr         t      P.Value    adj.P.Val        B
rna23260  9.733111 -2.401474  18.96688 6.770586e-10 1.525278e-05 7.506843
rna4289   8.671071 -2.564865  17.44008 1.691543e-09 1.905354e-05 7.248625
rna6574   7.687889 -2.716123  15.91767 4.555986e-09 2.608897e-05 6.933027

• 或者treat在基因排序中給予fold-changes更多權(quán)重

> fit4 <- treat(fit2, lfc=log2(1.2),trend=TRUE) 

• 提取有fold-changes閾值的某組比較的結(jié)果

> topTreat(fit4, coef=ncol(design))
            logFC     AveExpr        t      P.Value    adj.P.Val
rna23260 9.733111 -2.40147385 21.44573 2.000093e-10 4.505809e-06
rna4289  8.671071 -2.56486452 19.03239 7.024639e-10 7.912553e-06
rna34624 8.412221 -2.49050322 17.32005 1.908439e-09 1.433110e-05

• 提取所有比較結(jié)果

> dt <- decideTests(fit3)
> head(dt,3)
TestResults matrix
           Contrasts
            WTMvsWTF MFvsWTF MMFvsWTF MMFvsMF
  gene17427        0       0        0       0
  gene17428        0       0        0       0
  gene17507        0       0        0       0
> summary(dt)
       WTMvsWTF MFvsWTF MMFvsWTF MMFvsMF
Down         90      14       23      16
NotSig    20700   22477    22452   22481
Up         1738      37       53      31

• 保存為文件

> write.fit(fit3, dt, file="results.txt")

（2）voom

• voom轉(zhuǎn)換應(yīng)用于過濾后并標準化的DGEList對象

> v <- voom(yf, design, plot=TRUE) 
> v #EList對象

voom

• 在此之后，可以使用通常的limma差異表達分析工作流程

> fit <- lmFit(v, design) 
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit3 <- eBayes(fit2)

• 繪制log2殘差標準差與log-CPM均值的關(guān)系

> plotSA(fit3, main="Final model: Mean-variance trend")

plotSA

• 提取結(jié)果

> restop <- topTable(fit3, coef=ncol(design), n=Inf)

五、limma用法

（1）voom——Transform RNA-Seq Data Ready for Linear Modelling

voom(counts, design = NULL, lib.size = NULL, normalize.method = "none", block = NULL, correlation = NULL, weights = NULL, span = 0.5, plot = FALSE, save.plot = FALSE)

• counts：是一個包含raw counts的數(shù)值矩陣、或ExpressionSet、DGEList對象。counts必須是非負數(shù)，并且不含有NAs。
• design：設(shè)計矩陣。
• lib.size：每一個樣本庫的大小。默認為DGEList對象的normalized (effective) library sizes。

• plot：是否繪制mean-variance trend圖。

  
  
  voom用法

（2）lmFit

lmFit(object, design=NULL, ndups=1, spacing=1, block=NULL, correlation, weights=NULL, method="ls", ...)

• object：一種類似于矩陣的數(shù)據(jù)對象，其中包含log-expression values。
• design：設(shè)計矩陣。
• method：擬合方法，"ls"表示最小二乘，"robust"表示穩(wěn)健回歸。

（3）contrasts.fit

contrasts.fit(fit, contrasts=NULL, coefficients=NULL)

• fit：是lmFit的輸出。包含coefficients和stdev.unscaled。
• contrasts：數(shù)值矩陣，行對應(yīng)于擬合中的系數(shù)，列包含對比。如果只有一個對比，可能是一個向量?？梢允莔akeContrasts得到的對比矩陣。
• coefficients：向量，指示哪些系數(shù)將保留在修改后的擬合對象中。 指定對比的另一種方法。
• contrasts.fit之后會得到coefficients、stdev.unscaled、cov.coefficients，fit里面大多數(shù)組成成分會被保留，但是t、p.value、lods、F、F.p.value會被移除。

contrasts.fit

（4）eBayes

eBayes(fit, proportion = 0.01, stdev.coef.lim = c(0.1,4), trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1))
treat(fit, lfc = log2(1.2), trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1))

• 給定一個微陣列線性模型擬合，計算moderated t-statistics、moderated F-statistic和差異表達的logodds，通過經(jīng)驗Bayes將標準誤差調(diào)整到一個共同的值。
• fit：lmFit或contrasts.fit得到的擬合模型。
• proportion：0-1之間的數(shù)，差異表達基因的假定比例。
• trend：對于先驗方差是否允許一個intensity-trend。默認先驗方差是常數(shù)。
• robust：df.prior和var.prior的估計是否穩(wěn)健地反對離群樣本方差（be robustified against outlier sample variances）。
• lfc：被認為具有科學(xué)意義的最小log2倍數(shù)變化。

eBayes

（5）toptable

topTable(fit, coef=NULL, number=10, genelist=fitgenes, adjust.method="BH", sort.by="F", p.value=1, lfc=0)
topTreat(fit, coef=1, sort.by="p", resort.by=NULL, ...)

• fit：由lmFit和eBayes/treat產(chǎn)生。
• coef：列號或列名，指定線性模型的哪個系數(shù)或?qū)Ρ仁歉信d趣的。對于topTable，也可以是列下標的向量，在這種情況下，基因排序是按照這組對比的F-statistic。
• number：列出的基因的最大數(shù)量。
• sort.by：按照什么對基因進行排序。對于topTable可以是"logFC"/"M", "AveExpr"/"A"/"Amean", "t"/"T", "P"/"p", "B" or "none"。對于topTableF可以是"F" or "none"。對于topTreat，與topTable相同，除了"B"。
• p.value：校正后的p值的閾值，只有小于這個值的基因才會被列出。
• lfc：要求的最小absolute log2-fold-change。topTable和topTableF只列出log-fold-changes大于設(shè)置的lfc值的基因。在topTreat這個argument不可用，不會去除基因，但根據(jù)它們的對數(shù)變化超過lfc的證據(jù)對基因進行排序。
• confint：是否輸出logFC的95%置信區(qū)間。或者在0-1之間取值作為置信度。

toptable結(jié)果

（6）decideTests

decideTests(object, method = "separate", adjust.method = "BH", p.value = 0.05, lfc = 0, ...)

• object：一個p值的數(shù)值矩陣，或eBayes和treat產(chǎn)生的對象，從中可以提取p值和t統(tǒng)計量。
• method：指定如何在多個測試方案中組合基因和對比。可以是"separate", "global", "hierarchical" or "nestedF"。
• method="separate"將對p值的每一列分別應(yīng)用多重測試調(diào)整。設(shè)置method="separate"等價于對線性模型擬合中的每個系數(shù)分別使用topTable的方法。method="global"將把t統(tǒng)計量的整個矩陣看作一個不相關(guān)檢驗的單一向量。當需要在所有對比中有相同的t統(tǒng)計量閾值時，method="global"是有用的。
• adjust.method：可以是"none", "BH", "fdr" (equivalent to "BH"), "BY" and "holm"。
• p.value：0-1的數(shù)值，給出要求的family-wise error rate or false discovery rate。
• lfc：要求的最小absolute log2-fold-change。
• 結(jié)果用-1,0,1表示significantly negative, not significant, significantly positive。

（7）write.fit

write.fit(fit, results = NULL, file, digits = NULL, adjust = "none", method = "separate", F.adjust = "none", quote = FALSE, sep = "\t", row.names = TRUE, ...)

This function writes a tab-delimited text file containing for each gene (1) the average log2-intensity, (2) the coefficients or contrasts (log2-fold-changes), (3) moderated t-statistics, (4) t-statistic P-values, (5) F-statistic if available, (6) F-statistic P-values if available, (7) classification if available and (8) gene names and annotation.

> dt <- decideTests(fit3)
> write.fit(fit3, dt, file="resultsvoom.txt")
> restxt <- read.table("resultsvoom.txt")
> colnames(restxt)
 [1] "A"                   "Coef.WTMvsWTF"      
 [3] "Coef.MFvsWTF"        "Coef.MMFvsWTF"      
 [5] "Coef.MMFvsMF"        "t.WTMvsWTF"         
 [7] "t.MFvsWTF"           "t.MMFvsWTF"         
 [9] "t.MMFvsMF"           "p.value.WTMvsWTF"   
[11] "p.value.MFvsWTF"     "p.value.MMFvsWTF"   
[13] "p.value.MMFvsMF"     "F"                  
[15] "F.p.value"           "Res.WTMvsWTF"       
[17] "Res.MFvsWTF"         "Res.MMFvsWTF"       
[19] "Res.MMFvsMF"         "Genes.chromosome"   
[21] "Genes.ref_gene_name"