差異表達分析通常作為根據基因表達矩陣進行生物資訊學分析的第一步，有助於我們觀察基因在不同樣本中的表達差異，從而確定要研究的基因和表型之間的聯絡。常用的基因表達資料來自基因晶片或高通量測序。雖然矩陣看起來差不多，但是由於服從不同的分佈，因此在進行差異表達的時候需要用不同的方法。對於一般的生命科學領域科研人員來說，瞭解晦澀的演算法並沒有太大價值。本文力求精簡，從資料——演算法——結果三個方面給出最簡單的示範。注意：文中程式碼僅適用於基因晶片的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處理的原因和判斷方法見：（還沒寫，過幾天補充）

如果資料不需要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將全部結果、滿足篩選條件（即差異表達倍數）的全部結果、上調結果、下調結果分別輸出。

這一步的篩選標準在程式碼剛開始時設定。