生物信息学入门使用 GEO基因芯片数据进行差异表达分析（DEG）

imtravelinghah 2022-07-16 发布于广西

展开全文

差异表达分析通常作为根据基因表达矩阵进行生物信息学分析的第一步，有助于我们观察基因在不同样本中的表达差异，从而确定要研究的基因和表型之间的联系。常用的基因表达数据来自基因芯片或高通量测序。虽然矩阵看起来差不多，但是由于服从不同的分布，因此在进行差异表达的时候需要用不同的方法。对于一般的生命科学领域科研人员来说，了解晦涩的算法并没有太大价值。本文力求精简，从数据——算法——结果三个方面给出最简单的示范。注意：文中代码仅适用于基因芯片的counts数据！使用的是limma算法！

基于TCGA的FPKM数据进行差异表达的算法可以参考：（还没写，过几天补充）

1.数据准备

数据准备包括表达矩阵和分组矩阵。

表达矩阵：

分组矩阵

第一列为样本名称，第二列为组名称，注意每一列都要有列名

2. 使用Limma包进行差异分析

首先要安装limma包和gplots包

source("http:///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")