劉小澤寫于19.12.29 + 2020.1.1
這一篇不扯別的，只列我認(rèn)為比較重點(diǎn)的limma包的知識熬芜，所以代碼部分可能比較多社裆。但也會用不同的知識點(diǎn)來區(qū)分
前言
limma的全稱是：Linear Models for Microarray Data
需要閱讀limma的官方說明：https://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf
尤其是第8章 Linear Models Overview
常用知識點(diǎn)
知識點(diǎn)一 （文字解釋）
進(jìn)行差異分析時(shí)常用limma效览。雖然它是針對芯片數(shù)據(jù)開發(fā)的，但也有l(wèi)imma-voom可以分析轉(zhuǎn)錄組數(shù)據(jù)刹孔，可以作為金標(biāo)準(zhǔn)啡省。
它需要的輸入文件有：
-
表達(dá)矩陣 （exprSet）（這個(gè)容易獲得），芯片數(shù)據(jù)可以通過
exprs()髓霞，常規(guī)的轉(zhuǎn)錄組可以通過read.csv()等導(dǎo)入 -
分組矩陣 (design) ：就是將表達(dá)矩陣的列（各個(gè)樣本）分成幾組（例如最簡單的
case - control卦睹，或者一些時(shí)間序列的樣本day0, day1, day2 ...）【通過model.matrix()得到】 -
比較矩陣（contrast）：意思就是如何指定函數(shù)去進(jìn)行組間比較【通過
makeContrasts()得到】
它的主要流程有：
-
lmFit()：線性擬合模型構(gòu)建【需要兩個(gè)東西：
exprSet和design】 ，得到的結(jié)果再和contrast一起導(dǎo)入contrasts.fit()函數(shù) -
eBayes()：利用上一步
contrasts.fit()的結(jié)果 -
topTable()：利用上一步
eBayes()的結(jié)果酸茴，并最終導(dǎo)出差異分析結(jié)果
知識點(diǎn)二（代碼演示）
搭配上面??的解釋來看分预，效果更好
分開展示 =》 構(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è)問題：樣本數(shù)量和分組數(shù)量不一樣，因?yàn)榧词褂?00個(gè)樣本薪捍，最后還可能依然分成2組（處理和對照組笼痹，我們只關(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）來提取結(jié)果-
adjust="BH"表示對p值的校正方法凳干，包括了："none","BH","BY"and"holm"晴裹。那么為啥要對P值進(jìn)行校正呢？
p值是針對單次檢驗(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í)，就會遇到一個(gè)問題：單次犯錯的概率為0.05预厌， 如果進(jìn)行20000次檢驗(yàn)阿迈，那么就會有0.05*20000=1000 個(gè)基因在組間的差異被錯誤估計(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)
知識點(diǎn)三 一定要使用比較矩陣嗎轧叽？
答案是不一定
看這里：https://github.com/bioconductor-china/basic/blob/master/makeContrasts.md
然后可以再結(jié)合說明書Chapter9.2 (第42頁)
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來定義
其實(shí)差異分析，一個(gè)使用難點(diǎn)就是：分組
limma包關(guān)于如何針對特殊情況進(jìn)行分組描述了很大的篇幅
例如 如果包含多個(gè)組多個(gè)處理
參考：limma說明書 P43
-
實(shí)驗(yàn)設(shè)計(jì)背景
如果包含多個(gè)組多個(gè)處理有6個(gè)樣本网严，分成3組（1识樱、2、3）震束，并且每組包括兩個(gè)處理方式：一個(gè)對照（C）牺荠，一個(gè)處理（T）
-
我們自己來構(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)的處理和對照（Treatment）
例如 時(shí)間順序的樣本
參考P47 的 Chapter 9.6
比方說驴一，有兩種基因型（wt、mu）灶壶，各自測了3個(gè)時(shí)間點(diǎn)（0肝断、6、24h）
可以這樣操作：
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")
# 那么對于mu型也是一樣的
如果要探索從0-6胸懈、從6到24h處理后，mu相對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)然恰响，還有很多很多的例子趣钱，都在說明書有體現(xiàn)。
這里只是列出來一個(gè)思路：凡是想用一個(gè)包胚宦，它的幫助文檔就是最好的答案
歡迎關(guān)注我們的公眾號～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩首有，打造生信星球燕垃，想讓它成為一個(gè)不拽術(shù)語、通俗易懂的生信知識平臺井联。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
