差异表达分析通常作为根据基因表达矩阵进行生物信息学分析的第一步,有助于我们观察基因在不同样本中的表达差异,从而确定要研究的基因和表型之间的联系。常用的基因表达数据来自基因芯片或高通量测序。虽然矩阵看起来差不多,但是由于服从不同的分布,因此在进行差异表达的时候需要用不同的方法。对于一般的生命科学领域科研人员来说,了解晦涩的算法并没有太大价值。本文力求精简,从数据——算法——结果三个方面给出最简单的示范。注意:文中代码仅适用于基因芯片的counts数据!使用的是limma算法!
基于TCGA的FPKM数据进行差异表达的算法可以参考:(还没写,过几天补充)
1.数据准备
数据准备包括表达矩阵和分组矩阵。
表达矩阵:
分组矩阵
第一列为样本名称,第二列为组名称,注意每一列都要有列名
2. 使用Limma包进行差异分析
首先要安装limma包和gplots包
source("/")
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=("",header=TRUE,=1, = FALSE)
#读取矩阵文件,这是输入的数据路径,改成自己的文件名#
dim(rawexprSet)
exprSet=log2(rawexprSet)
par(mfrow=c(1,2))
boxplot((exprSet),col="blue") ## 画箱式图,比较数据分布情况
exprSet[1:5,1:5]
group <- ("",header=TRUE,=1, = FALSE)
group <- group[,1] #定义比较组,按照癌症和正常样品数目修改#
design <- (~0+factor(group))#把group设置成一个model matrix#
colnames(design)=levels(factor(group))
rownames(design)=colnames(exprSet)
这里需要注意,从GEO下载的表达矩阵中,并非所有的数据都是已经log处理,对于没有log处理的数据需要自己log.
log处理的原因和判断方法见:
GEO芯片数据差异表达分析时需要log2处理的原因
/tuanzide5233/article/details/88542805
GEO芯片数据差异表达分析时是否需要log2以及标准化的问题
/tuanzide5233/article/details/88542558
如果数据不需要log处理,只要将图中所示的代码前面加上#,即注释掉
注释后:
右下角的箱线图表明数据还是比较整齐的,可以进行下一步分析
计算步骤
fit <- lmFit(exprSet,design)
<-makeContrasts(paste0(unique(group),collapse = "-"),levels = design)
fit2=(fit,)
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
tempOutput = topTable(fit2,coef=1,n=Inf,adjust="BH")
nrDEG = (tempOutput)
输出结果:
allDiff <- nrDEG
diff=allDiff
(diff, "")
diffSig = diff[(diff$ < padj & (diff$logFC>foldChange | diff$logFC<(-foldChange))),]#筛选有显著差异的#
#(diffSig, file="",sep="\t",quote=F)#输出有显著差异表达的到diffSig这个文件#
(diffSig, "")
diffUp = diff[(diff$ < padj & (diff$logFC>foldChange)),]#foldchange>0是上调,foldchange<0是下调#
#(diffUp, file="",sep="\t",quote=F)#39-42把上调和下调分别输入up和down两个文件#
(diffUp, "")
diffDown = diff[(diff$ < padj & (diff$logFC<(-foldChange))),]
#(diffDown, file="",sep="\t",quote=F)
(diffDown, "")
这里可以看到按照padj将全部结果、满足筛选条件(即差异表达倍数)的全部结果、上调结果、下调结果分别输出。
这一步的筛选标准在代码刚开始时设置。
GEO芯片数据差异表达分析时需要log2处理的原因
/tuanzide5233/article/details/88542805
GEO芯片数据差异表达分析时是否需要log2以及标准化的问题
/tuanzide5233/article/details/88542558
差异表达矩阵制作教程
/tuanzide5233/article/details/83659768
差异表达的热图绘制详见
/tuanzide5233/article/details/83659501
使用edgeR对RNAseq数据进行差异表达分析教程
/tuanzide5233/article/details/88785486
差异表达分析(DEG)时 '里不能有重复的名字 的解决方案
/tuanzide5233/article/details/86568155
生存分析系列教程(一)使用生信人工具盒进行生存分析
/tuanzide5233/article/details/83685403
富集分析与蛋白质互作用网络(PPI)的可视化 Cystocape入门指南
/tuanzide5233/article/details/88048439
进阶版Venn plot:Upset plot入门实战代码详解——UpSetR包介绍
/tuanzide5233/article/details/83109527
使用R语言ggplot2包绘制pathway富集分析气泡图(Bubble图):数据结构及代码
/tuanzide5233/article/details/82141817