生物資訊學入門 使用 GEO基因晶片資料進行差異表達分析（DEG）——Limma 算法 資料 代碼 結果解讀

       差異表達分析通常作為根據基因表達矩陣進行生物資訊學分析的第一步，有助于我們觀察基因在不同樣本中的表達差異，進而确定要研究的基因和表型之間的聯系。常用的基因表達資料來自基因晶片或高通量測序。雖然矩陣看起來差不多，但是由于服從不同的分布，是以在進行差異表達的時候需要用不同的方法。對于一般的生命科學領域科研人員來說，了解晦澀的算法并沒有太大價值。本文力求精簡，從資料——算法——結果三個方面給出最簡單的示範。注意：文中代碼僅适用于基因晶片的counts資料！使用的是limma算法！

       基于TCGA的FPKM資料進行差異表達的算法可以參考：（還沒寫，過幾天補充）

1.資料準備

資料準備包括表達矩陣和分組矩陣。

表達矩陣：

分組矩陣

第一列為樣本名稱，第二列為組名稱，注意每一列都要有列名

2. 使用Limma包進行差異分析

首先要安裝limma包和gplots包

source("http://bioconductor.org/biocLite.R")
biocLite("Limma")
biocLite("gplots")
           

讀取資料

#DGE for microarray by limma
library('gplots')
library('limma')
setwd("C:/Users/lenovo/DEG")
foldChange=0.5 #fold change=1意思是差異是兩倍
padj=0.01#padj=0.05意思是矯正後P值小于0.05
rawexprSet=read.csv("express-counts2.csv",header=TRUE,row.names=1,check.names = FALSE)  
#讀取矩陣檔案，這是輸入的資料路徑，改成自己的檔案名#
dim(rawexprSet)
exprSet=log2(rawexprSet)
par(mfrow=c(1,2))
boxplot(data.frame(exprSet),col="blue") ## 畫箱式圖，比較資料分布情況
exprSet[1:5,1:5]
group <- read.csv("datTraits.csv",header=TRUE,row.names=1,check.names = FALSE)
group <- group[,1] #定義比較組，按照癌症和正常樣品數目修改#
design <- model.matrix(~0+factor(group))#把group設定成一個model matrix#
colnames(design)=levels(factor(group))
rownames(design)=colnames(exprSet)
           

這裡需要注意，從GEO下載下傳的表達矩陣中，并非所有的資料都是已經log處理，對于沒有log處理的資料需要自己log.

log處理的原因和判斷方法見：

GEO晶片資料差異表達分析時需要log2處理的原因

https://blog.csdn.net/tuanzide5233/article/details/88542805

GEO晶片資料差異表達分析時是否需要log2以及标準化的問題

https://blog.csdn.net/tuanzide5233/article/details/88542558

如果資料不需要log處理，隻要将圖中所示的代碼前面加上#，即注釋掉

注釋後：

右下角的箱線圖表明資料還是比較整齊的，可以進行下一步分析

計算步驟

fit <- lmFit(exprSet,design)
cont.matrix<-makeContrasts(paste0(unique(group),collapse = "-"),levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
tempOutput = topTable(fit2,coef=1,n=Inf,adjust="BH") 
nrDEG = na.omit(tempOutput) 
           

 輸出結果：

allDiff <- nrDEG
diff=allDiff
write.csv(diff, "limmaOut.csv")
diffSig = diff[(diff$P.Value < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]#篩選有顯著差異的#
#write.table(diffSig, file="diffSig.xls",sep="\t",quote=F)#輸出有顯著差異表達的到diffSig這個檔案#
write.csv(diffSig, "diffSig.csv")
diffUp = diff[(diff$P.Value < padj & (diff$logFC>foldChange)),]#foldchange>0是上調，foldchange<0是下調#
#write.table(diffUp, file="up.xls",sep="\t",quote=F)#39-42把上調和下調分别輸入up和down兩個檔案#
write.csv(diffUp, "diffUp.csv")
diffDown = diff[(diff$P.Value < padj & (diff$logFC<(-foldChange))),]
#write.table(diffDown, file="down.xls",sep="\t",quote=F)
write.csv(diffDown, "diffDown.csv")
           

這裡可以看到按照padj将全部結果、滿足篩選條件（即差異表達倍數）的全部結果、上調結果、下調結果分别輸出。

這一步的篩選标準在代碼剛開始時設定。

GEO晶片資料差異表達分析時需要log2處理的原因

https://blog.csdn.net/tuanzide5233/article/details/88542805

GEO晶片資料差異表達分析時是否需要log2以及标準化的問題

https://blog.csdn.net/tuanzide5233/article/details/88542558

差異表達矩陣制作教程

https://blog.csdn.net/tuanzide5233/article/details/83659768

差異表達的熱圖繪制詳見

https://blog.csdn.net/tuanzide5233/article/details/83659501

使用edgeR對RNAseq資料進行差異表達分析教程

https://blog.csdn.net/tuanzide5233/article/details/88785486

差異表達分析（DEG）時 row.names'裡不能有重複的名字 的解決方案

https://blog.csdn.net/tuanzide5233/article/details/86568155

生存分析系列教程（一）使用生信人工具盒進行生存分析

https://blog.csdn.net/tuanzide5233/article/details/83685403

富集分析與蛋白質互作用網絡（PPI）的可視化 Cystocape入門指南

https://blog.csdn.net/tuanzide5233/article/details/88048439

進階版Venn plot：Upset plot入門實戰代碼詳解——UpSetR包介紹

https://blog.csdn.net/tuanzide5233/article/details/83109527

使用R語言ggplot2包繪制pathway富集分析氣泡圖（Bubble圖）:資料結構及代碼

https://blog.csdn.net/tuanzide5233/article/details/82141817