差異表達(dá)分析通常作為根據(jù)基因表達(dá)矩陣進(jìn)行生物信息學(xué)分析的第一步，有助于我們觀察基因在不同樣本中的表達(dá)差異，從而確定要研究的基因和表型之間的聯(lián)系。常用的基因表達(dá)數(shù)據(jù)來自基因芯片或高通量測序。雖然矩陣看起來差不多，但是由于服從不同的分布，因此在進(jìn)行差異表達(dá)的時候需要用不同的方法。對于一般的生命科學(xué)領(lǐng)域科研人員來說，了解晦澀的算法并沒有太大價值。本文力求精簡，從數(shù)據(jù)——算法——結(jié)果三個方面給出最簡單的示范。注意：文中代碼僅適用于基因芯片的counts數(shù)據(jù)！使用的是limma算法！

       基于TCGA的FPKM數(shù)據(jù)進(jìn)行差異表達(dá)的算法可以參考：（還沒寫，過幾天補(bǔ)充）

1.數(shù)據(jù)準(zhǔn)備

數(shù)據(jù)準(zhǔn)備包括表達(dá)矩陣和分組矩陣。

表達(dá)矩陣：

分組矩陣

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

2. 使用Limma包進(jìn)行差異分析

首先要安裝limma包和gplots包

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

讀取數(shù)據(jù)

#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)  
#讀取矩陣文件，這是輸入的數(shù)據(jù)路徑，改成自己的文件名#
dim(rawexprSet)
exprSet=log2(rawexprSet)
par(mfrow=c(1,2))
boxplot(data.frame(exprSet),col="blue") ## 畫箱式圖，比較數(shù)據(jù)分布情況
exprSet[1:5,1:5]
group <- read.csv("datTraits.csv",header=TRUE,row.names=1,check.names = FALSE)
group <- group[,1] #定義比較組，按照癌癥和正常樣品數(shù)目修改#
design <- model.matrix(~0+factor(group))#把group設(shè)置成一個model matrix#
colnames(design)=levels(factor(group))
rownames(design)=colnames(exprSet)

這里需要注意，從GEO下載的表達(dá)矩陣中，并非所有的數(shù)據(jù)都是已經(jīng)log處理，對于沒有l(wèi)og處理的數(shù)據(jù)需要自己log.

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

GEO芯片數(shù)據(jù)差異表達(dá)分析時需要log2處理的原因

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

GEO芯片數(shù)據(jù)差異表達(dá)分析時是否需要log2以及標(biāo)準(zhǔn)化的問題

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

如果數(shù)據(jù)不需要log處理，只要將圖中所示的代碼前面加上#，即注釋掉

注釋后：

右下角的箱線圖表明數(shù)據(jù)還是比較整齊的，可以進(jìn)行下一步分析

計算步驟

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) 

 輸出結(jié)果：

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)#輸出有顯著差異表達(dá)的到diffSig這個文件#
write.csv(diffSig, "diffSig.csv")
diffUp = diff[(diff$P.Value < padj & (diff$logFC>foldChange)),]#foldchange>0是上調(diào)，foldchange<0是下調(diào)#
#write.table(diffUp, file="up.xls",sep="\t",quote=F)#39-42把上調(diào)和下調(diào)分別輸入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將全部結(jié)果、滿足篩選條件（即差異表達(dá)倍數(shù)）的全部結(jié)果、上調(diào)結(jié)果、下調(diào)結(jié)果分別輸出。

這一步的篩選標(biāo)準(zhǔn)在代碼剛開始時設(shè)置。

GEO芯片數(shù)據(jù)差異表達(dá)分析時需要log2處理的原因

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

GEO芯片數(shù)據(jù)差異表達(dá)分析時是否需要log2以及標(biāo)準(zhǔn)化的問題

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

差異表達(dá)矩陣制作教程

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

差異表達(dá)的熱圖繪制詳見

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

使用edgeR對RNAseq數(shù)據(jù)進(jìn)行差異表達(dá)分析教程

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

差異表達(dá)分析（DEG）時 row.names'里不能有重復(fù)的名字 的解決方案

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

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

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

富集分析與蛋白質(zhì)互作用網(wǎng)絡(luò)（PPI）的可視化 Cystocape入門指南

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

進(jìn)階版Venn plot：Upset plot入門實(shí)戰(zhàn)代碼詳解——UpSetR包介紹

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

使用R語言ggplot2包繪制pathway富集分析氣泡圖（Bubble圖）:數(shù)據(jù)結(jié)構(gòu)及代碼

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