使用limma包進行差異基因分析時，做最多的是兩分類的，例如control組和disease組，但也會碰到按照序列進行的分組。這時，如果逐一使用兩兩比較求差異基因則略顯復雜。其實開發limma包的大神們已經替我們考慮到。我自己當下limma包的PDF，仔細研讀并將代碼運行后分享給有需要的同道，相互學習。

library(limma)##加載包 load("GSE37761exp_groupfile-00.RData")#加載匹配過的基因表達矩陣 eset=exp ##將基因表達矩陣賦值給eset targets<-read.csv("group_file.csv",row.names = 1)##讀入樣本數據，包括兩列， #第一列GSM號，第二列為樣本分組#該數據集中實際樣本分組中存在24h_1、24h_5、24h_21等，然而這不符合R的命名原則 targets$Target=gsub("_",".",targets$Target)##將"_"替換成“.”，也可以不替換 ##該數據集中實際樣本分組中存在24h_1、24h_5、24h_21等，然而這不符合R的命名原則，所以在沒個分類前加一個“F”，具體自己定 targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep="")) colnames(targets)=c("FileName","Target")#更改列名，為了和limma包中的一致 lev<-unique(targets$Target)##使用unique（）函數進行去重 f <- factor(targets$Target, levels=lev) design <- model.matrix(~0+f) #樣本矩陣 colnames(design) <- lev #更改列名為levels名 ###兩兩之間的比較，求差異基因使用topTable函數 cont.wt <- makeContrasts("F12h.1-Fcontrol",#舉例levels=design) fit <- lmFit(eset, design) fit2 <- contrasts.fit(fit, cont.wt) fit2 <- eBayes(fit2) tT=topTable(fit2, adjust="BH",sort.by="logFC",n=Inf) tT = subset(tT, select=c("adj.P.Val","P.Value","logFC")) colnames(tT)=c("FDR","P.Value","logFC") ### ##后面可以設置fdr及logFC的閾值進行篩選，進而得到差異基因列表，此處略。##兩組以上的比較，尋找差異基因，使用topTableF函數 cont.wt <- makeContrasts("F12h.1-Fcontrol",##舉例"F24h.1-F12h.1","F5d.1-F24h.1","F10d.1-F5d.1",levels=design) fit <- lmFit(eset, design) fit2 <- contrasts.fit(fit, cont.wt) fit2 <- eBayes(fit2) tT=topTableF(fit2, adjust="BH",sort.by="F",n=Inf) ##tT = subset(tT, select=c("adj.P.Val","P.Value","logFC")) 這一條和下一條代碼需要根據自己需要進行更改，topTableF函數得到的結果中包括F值，而logFC則是沒兩組的比較，均得出了 ##colnames(tT)=c("FDR","P.Value","logFC")

因為不支持PDF上傳，最后附上limma包usersguide.pdf的下載地址。http://www.bioconductor.org/packages/3.2/bioc/vignettes/limma/inst/doc/usersguide.pdf

---

學習過程的就是分享的過程，分享的過程也是交流的過程，交流的過程就是進步的過程。

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的python 分析两组数据的差异_R语言limma包差异基因分析（两组或两组以上）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。