劉小澤寫于19.12.29 + 2020.1.1
這一篇不扯別的，只列我認(rèn)為比較重點(diǎn)的limma包的知識(shí)，所以代碼部分可能比較多。但也會(huì)用不同的知識(shí)點(diǎn)來(lái)區(qū)分

前言

limma的全稱是：Linear Models for Microarray Data

需要閱讀limma的官方說(shuō)明：https://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf

尤其是第8章 Linear Models Overview

常用知識(shí)點(diǎn)

知識(shí)點(diǎn)一 （文字解釋）

進(jìn)行差異分析時(shí)常用limma。雖然它是針對(duì)芯片數(shù)據(jù)開發(fā)的，但也有l(wèi)imma-voom可以分析轉(zhuǎn)錄組數(shù)據(jù)，可以作為金標(biāo)準(zhǔn)。

它需要的輸入文件有：

• 表達(dá)矩陣 （exprSet）（這個(gè)容易獲得），芯片數(shù)據(jù)可以通過(guò)exprs() ，常規(guī)的轉(zhuǎn)錄組可以通過(guò)read.csv()等導(dǎo)入
• 分組矩陣 (design) ：就是將表達(dá)矩陣的列（各個(gè)樣本）分成幾組（例如最簡(jiǎn)單的case - control，或者一些時(shí)間序列的樣本day0, day1, day2 ...）【通過(guò)model.matrix()得到】
• 比較矩陣（contrast）：意思就是如何指定函數(shù)去進(jìn)行組間比較【通過(guò)makeContrasts()得到】

它的主要流程有：

• lmFit()：線性擬合模型構(gòu)建【需要兩個(gè)東西：exprSet和design】 ，得到的結(jié)果再和contrast一起導(dǎo)入contrasts.fit()函數(shù)
• eBayes()：利用上一步contrasts.fit()的結(jié)果
• topTable()：利用上一步eBayes()的結(jié)果，并最終導(dǎo)出差異分析結(jié)果

知識(shí)點(diǎn)二（代碼演示）

搭配上面??的解釋來(lái)看，效果更好

分開展示 =》 構(gòu)建三個(gè)輸入文件

# 輸入文件一：例如我現(xiàn)在已經(jīng)有了一個(gè)表達(dá)矩陣eset

# 輸入文件二：分組矩陣【假設(shè)9個(gè)樣本分成了3組】
design <- model.matrix(~ 0+factor(c(1,1,1,2,2,3,3,3)))
colnames(design) <- c("group1", "group2", "group3")
rownames(design) <- colnames(eset)

# 輸入文件三：比較矩陣【如果要進(jìn)行三組之間的兩兩比較】
contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)

• 注意一個(gè)問(wèn)題：樣本數(shù)量和分組數(shù)量不一樣，因?yàn)榧词褂?00個(gè)樣本，最后還可能依然分成2組（處理和對(duì)照組，我們只關(guān)心分組）
• 比較矩陣的組之間用-連接

分開展示 =》 進(jìn)行三個(gè)主要流程

# 第一步：lmFit
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrast.matrix)

# 第二步：eBayes
fit2 <- eBayes(fit2)

# 第三步：topTable【最后例如要挑出第一個(gè)比較組：group2-group1的差異分析結(jié)果】
topTable(fit2, coef=1, adjust="BH")

• 使用coef參數(shù)，這里設(shè)為1，也就是表示??根據(jù)上面makeContrasts的第一個(gè)（group2-group1）來(lái)提取結(jié)果

• adjust="BH"表示對(duì)p值的校正方法，包括了："none", "BH", "BY" and "holm"。

  那么為啥要對(duì)P值進(jìn)行校正呢？
  p值是針對(duì)單次檢驗(yàn)，假設(shè)采用的p值為小于0.05，我們通常認(rèn)為這個(gè)基因在兩組樣本中的表達(dá)是有顯著差異的，但是仍舊有5%的概率表示這個(gè)基因并不是差異基因。

  但是，當(dāng)兩組樣本中有20000個(gè)基因采用同樣的檢驗(yàn)方式進(jìn)行統(tǒng)計(jì)檢驗(yàn)時(shí)，就會(huì)遇到一個(gè)問(wèn)題：?jiǎn)未畏稿e(cuò)的概率為0.05， 如果進(jìn)行20000次檢驗(yàn)，那么就會(huì)有0.05*20000=1000 個(gè)基因在組間的差異被錯(cuò)誤估計(jì)

最后整合展示代碼

# 還是假設(shè)有了表達(dá)矩陣eset
design <- model.matrix(~ 0+factor(c(1,1,1,2,2,3,3,3)))
colnames(design) <- c("group1", "group2", "group3")
fit <- lmFit(eset, design)
contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
# 【例如要挑出第一個(gè)比較組：group2-group1的差異分析結(jié)果】
output <- topTable(fit2, coef=1, adjust="BH")

最后的最后，別忘了去掉那些NA值

limma_DEG = na.omit(output) 
# 然后可以選擇保存
# write.csv(limma_DEG,"limma_results.csv",quote = F)

知識(shí)點(diǎn)三 一定要使用比較矩陣嗎？

答案是不一定
看這里：https://github.com/bioconductor-china/basic/blob/master/makeContrasts.md

然后可以再結(jié)合說(shuō)明書Chapter9.2 (第42頁(yè))

Group <- factor(targets$Target, levels=c("WT","Mu"))

• 如果使用：

  design <- model.matrix(~Group)
  colnames(design) <- c("WT","MUvsWT")
  

  那么它已經(jīng)把要比較的組放在了第一列，然后其余的列都與第一列進(jìn)行比較，而結(jié)果使用coef進(jìn)行指定提取

• 或者可以用

  design <- model.matrix(~0+Group)
  colnames(design) <- c("WT","MU")
  

  它沒有設(shè)置默認(rèn)如何比較，僅僅是做了一個(gè)分組，后續(xù)還需要使用makeContrasts來(lái)定義

其實(shí)差異分析，一個(gè)使用難點(diǎn)就是：分組
limma包關(guān)于如何針對(duì)特殊情況進(jìn)行分組描述了很大的篇幅

例如 如果包含多個(gè)組多個(gè)處理

參考：limma說(shuō)明書 P43

• 實(shí)驗(yàn)設(shè)計(jì)背景

  如果包含多個(gè)組多個(gè)處理

  有6個(gè)樣本，分成3組（1、2、3），并且每組包括兩個(gè)處理方式：一個(gè)對(duì)照（C），一個(gè)處理（T）

• 我們自己來(lái)構(gòu)建這樣一個(gè)數(shù)據(jù)

  targets <- data.frame(FileName=paste0("File",1:6),
                        SibShip=c(1,1,2,2,3,3),
                        Treatment=rep(c("C","T"),3))
  

• 分組矩陣這樣設(shè)計(jì)

  library(limma)
  
  > (SibShip <- factor(targets$SibShip))
  [1] 1 1 2 2 3 3
  Levels: 1 2 3
  
  > (Treat <- factor(targets$Treatment, levels=c("C","T")))
  [1] C T C T C T
  Levels: C T
  # 注意這種設(shè)計(jì)方式：
  design <- model.matrix(~SibShip+Treat)
  

• 最后還是三步走

  fit <- lmFit(eset, design)
  fit <- eBayes(fit)
  topTable(fit, coef="TreatT")
  

上面分組矩陣的設(shè)計(jì)規(guī)律就是：

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

將每個(gè)組作為一個(gè)block，其中只比較組內(nèi)的處理和對(duì)照（Treatment）

例如 時(shí)間順序的樣本

參考P47 的 Chapter 9.6

比方說(shuō)，有兩種基因型（wt、mu），各自測(cè)了3個(gè)時(shí)間點(diǎn)（0、6、24h）

時(shí)間順序的樣本

可以這樣操作：

lev <- c("wt.0hr","wt.6hr","wt.24hr","mu.0hr","mu.6hr","mu.24hr")

f <- factor(Target, levels=lev)

design <- model.matrix(~0+f)
colnames(design) <- lev

fit <- lmFit(eset, design)

如果要探索野生型中從0到6、從6到24h等待后發(fā)生了什么變化

cont.wt <- makeContrasts(
      "wt.6hr-wt.0hr",
      "wt.24hr-wt.6hr", levels=design)

fit2 <- contrasts.fit(fit, cont.wt)
fit2 <- eBayes(fit2)
topTableF(fit2, adjust="BH")
# 那么對(duì)于mu型也是一樣的

如果要探索從0-6、從6到24h處理后，mu相對(duì)wt發(fā)生了什么變化

cont.dif <- makeContrasts(
     Dif6hr =(mu.6hr-mu.0hr)-(wt.6hr-wt.0hr),
     Dif24hr=(mu.24hr-mu.6hr)-(wt.24hr-wt.6hr),
 levels=design)

fit3 <- contrasts.fit(fit, cont.dif)
fit3 <- eBayes(fit3)
topTableF(fit3, adjust="BH")

當(dāng)然，還有很多很多的例子，都在說(shuō)明書有體現(xiàn)。
這里只是列出來(lái)一個(gè)思路：凡是想用一個(gè)包，它的幫助文檔就是最好的答案

歡迎關(guān)注我們的公眾號(hào)～_～　　
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個(gè)不拽術(shù)語(yǔ)、通俗易懂的生信知識(shí)平臺(tái)。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!