劉小澤寫于18.11.18
好資料值得花時間學習，即便是英文版鳖谈。limma的學習資料有很多岁疼，我每次更新的也是不同的內(nèi)容，希望能多涉及一些關于limma的知識缆娃，充實大腦
最近經(jīng)常從花花那里聽說.jpg捷绒，被感染了哈哈??
配置
之前果子更新了一篇推送，講的就是bioconductor的安裝方式發(fā)生了改變贯要。我們之前直接source+install.package安裝暖侨，現(xiàn)在用到了biocmanager。剛開始轉換還遇到了一點問題崇渗，可能是網(wǎng)絡原因吧字逗，總提示安不上biocmanager。但是有問題不能怕宅广，相信“一切軟件問題都不是問題”葫掉，只要電腦還健在，就可以解決
# 先安裝BiocManager跟狱，它位于CRAN
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
# 然后添加BioC的國內(nèi)源俭厚，可以選清華或者中科大
if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
# 這就是最新的安裝方式了
if(!require("limma")) BiocManager::install("limma",update = F,ask = F)
# 再安裝一個數(shù)據(jù)包（關于白血病的）、
if(!require("leukemiasEset")) BiocManager::install("leukemiasEset",update = F,ask = F)
前言
LIMMA的意思原來是“l(fā)inear models for microarray data”驶臊，這也就解釋了為什么其中有構建線性模型的步驟【線性模型包括了線性回歸挪挤、多重線性回歸叼丑、方差分析等】另外limma是為芯片而生。但不僅止于此扛门，后來衍生的limma-voom也依然成了今天轉錄組差異分析的金標準鸠信。另外DNA甲基化也可以用limma去分析（當然我現(xiàn)在還沒有涉及這部分）
使用limma需要有一個代表某種測量的數(shù)值矩陣，比如基因表達矩陣（一般就是行為feature论寨，列為sample）或者組蛋白修飾星立、Hi-C contact matrix。這里feature可以代表許多含義葬凳，在樣本信息中（比如GEO得到的樣本信息數(shù)據(jù)框）一個feature表示一個樣本贞铣，在表達矩陣中一個feature就等于一個基因，另外還可以表示基因組區(qū)間（比如不同的啟動子promoters或者基因重疊區(qū)域 genomic bins）沮明。不論是limma還是deseq2或者是edgeR，都采用了經(jīng)驗貝葉斯方法獲取feature信息窍奋。
我們這里只討論基因表達數(shù)據(jù)荐健，包括芯片和轉錄組數(shù)據(jù)，這種數(shù)據(jù)中一般樣本數(shù)量在100以內(nèi)琳袄，最多一般也不超過1000江场，并且樣本是有分組的，最常見的就是對照組control和處理組case窖逗。我們分析的目標也就是找到不同組之間差異的features（genes）址否。
但是并非每次都是兩組比較，有時涉及到更多信息碎紊，比如年齡佑附、性別；然后還有更復雜的仗考，就是時間序列信息音同，所有的樣本是按照時間點劃分的，每個時間點為一個組秃嗜，包括了幾個樣本权均。因此，一個實驗的design矩陣 就是組間樣本是如何分布的锅锨。樣本一般是相互獨立的叽赊，但是有時也是成對出現(xiàn)，比如同一個病人的正常與癌變組織就是這樣一對樣本必搞。最常見的design就是兩個未成對的組之間的比較必指。
兩組比較 A two group comparison
獲取數(shù)據(jù)
# 使用的是一個表達數(shù)據(jù)集ExpressionSet=》leukemiasEset
library(leukemiasEset)
data(leukemiasEset) #使用data來加載這個數(shù)據(jù)集
table(leukemiasEset$LeukemiaType)
得到結果如下：這是4組不同的白血病和1組對照的樣本信息，其中NoL意思是不是白血补嘶（Not leukemia）即正常對照組
##
## ALL AML CLL CML NoL
## 12 12 12 12 12
如果想知道ALL型白血病與對照組之間的差異基因取劫，可以先取出相關的子集ALL和NoL
myData <- leukemiasEset[, leukemiasEset$LeukemiaType %in% c("ALL", "NoL")]
myData$LeukemiaType <- factor(myData$LeukemiaType)
# 這里你會看到匆笤，雖然myData只包含了ALL和NoL兩種type，但是顯示的因子還是5種谱邪，因此需要重新賦值一下因子類型
myData$LeukemiaType <- factor(myData$LeukemiaType)
線性模型
現(xiàn)在做一個標準的limma model fit
design <- model.matrix(~ myData$LeukemiaType)
fit <- lmFit(myData, design)
fit <- eBayes(fit)
topTable(fit)
首先炮捧，為了比較感興趣組的樣本（這里的ALL和NoL），我們用model.matrix構建了design矩陣（叫design可以理解為是從實驗設計角度出發(fā)）惦银，這里的design矩陣有兩組【多組的需要稍微變一下構建方法】咆课；
其次，擬合線性模型（line model fitting）扯俱，然后利用經(jīng)驗貝葉斯方法【這是limma的核心所在】 得到基因間強度關系（gene strength）书蚪，因為這個實驗矩陣（design）只有兩組，因此只有一種比較方法：即比較他們之間的基因差別迅栅；
然后殊校，使用topTable函數(shù)列出差異表達最顯著的基因【如果實驗設計比較復雜，還需要告訴topTable比較哪些interest特征】
得到的結果部分如下：
logFC AveExpr t P.Value adj.P.Val
ENSG00000163751 4.089587 5.819472 22.51729 9.894742e-18 1.733025e-13
ENSG00000104043 4.519488 5.762115 21.98550 1.718248e-17 1.733025e-13
ENSG00000008394 5.267835 7.482490 20.08250 1.374231e-16 9.240332e-13
- 其中一個重要的部分是logFC读存，也就是case與control組之間的log Fold change值为流，這里有個注意事項：你需要知道樣本組比較順序（是ALL-NoL還是NoL-ALL）。如果是前者让簿，那么logFC為正就說明相對NoL敬察，基因在ALL是高表達的。因此判斷哪個樣本為reference是很重要的尔当。那么如何判斷莲祸？
myData$LeukemiaType
## [1] ALL ALL ALL ALL ALL ALL ALL ALL ALL ALL ALL ALL NoL NoL NoL NoL NoL
## [18] NoL NoL NoL NoL NoL NoL NoL
## Levels: ALL NoL
? 就看上面這個代碼的結果中的Levels部分哪個因子首先出現(xiàn)，它就作為reference
? 這里的結果是ALL先出現(xiàn)椭迎，因此ALL組是reference锐帜，也就是NoL-ALL的模式，這樣的話畜号，如果得到的logFC是正值抹估，那么就說明在ALL組中基因表達下調(diào)
? 如果想要改變reference組，可以用relevel()函數(shù)
- 另外結果還包括t檢驗數(shù)據(jù)弄兜、P值（P.Value 對于多組比較未校正）药蜻、校正的P值（adj.P.Val 默認使用Benjamini-Horchberg方法對多組進行P值校正）
歡迎關注我們的公眾號～_～
我們是兩個農(nóng)轉生信的小碩，打造生信星球替饿，想讓它成為一個不拽術語语泽、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com
