阅读 16 SEO

13.GEO数据集的R语言差异分析和代码解析--2.差异分析

一、举例回顾

本节下载GSE1009数据集,使用limma包进行差异分析举例。

GSE1009  

样本量:共6个样本,其中后3为糖尿病肾病(DN)肾小球样本,前3个为正常肾小球样本。

使用芯片:Affymetrix Human Genome U95 Version 2 Array。

平台:GPL8300。

二、需要准备的文件:

上一节中保存的GSE1009.Rdata和GSE1009expressionmetrix_GSE.csv(请把这两个文件放在一个文件夹中。)

三、差异分析代码解析:

>setwd("D:\\Rfile")

>rm(list = ls())

>options(stringsAsFactors=F)

#老规矩,先设置工作目录。


>load(file = 'GSE1009.Rdata')

#加载上一节保存的R文件,里面包含表达矩阵和分组信息,不记得了可以看看上一节代码。


>table(group_list)

#计数每组样本的个数。


>library(limma)

#调用limma包,注意调用包之前,请先安装。

#limma包是biocondutor中的包,搜索安装命令的方式见上一节,这里就不重复了。最新的安装命令如下:


>logFoldChange=1

>adjustP=0.05

#设置差异基因DEGs的筛选阈值。

#关于筛选阈值:一般有三个筛选条件,log2FC、P、FDR(调整的P),大家可以根据自己的需要设置。


>rt=read.csv("GSE1009expressionmetrix_GSE.csv")

#读取整理好的表达矩阵数据,读入R中的矩阵的数据结构为列表,将整个列表赋值给变量rt.


>rt=as.matrix(rt)

#将列表rt转换为矩阵。

#这时表达矩阵如下:


注意行名

>rownames(rt)=rt[,1]

#将第一列数据(基因名)作为行名。


再看看行名,已经变了

>exp=rt[,2:ncol(rt)]

#选择第2列到最后一列,即表达矩阵,并赋值给exp。这样就得到行名为基因名,列名为样本名的矩阵。

>dimnames=list(rownames(exp),colnames(exp))

#提取出exp的行名和列名。这里的dimnames就是矩阵创建matrix里面的参数,在数据结构矩阵那一节有讲。

>rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)

#合并整个矩阵,得到一个行名为基因名,列名为样本名的表达矩阵。

>modType=c(rep("Control",3),rep("DN",3))

#构建分组模型,根据GSE1009样本的排列,前3个为正常样本,后3个为糖肾样本,构建一个分组的向量模型。

>design <- model.matrix(~0+factor(modType))

>colnames(design) <- c("Control","DN")

#创建线性模型。

#注意这里,在构建模型时,有一些数据集的处理组在前,对照组在后,且因为R里面默认的字符存储顺序为字母表顺序,如果处理组的首写字母的顺序在对照组的后面,构建模型时前三个又是处理组,那design出来的结果会上下调就完全相反。所以如果是这种情况的话,需要定义因子顺序。如:“modType=c(rep("DN",40),rep("Control",40))

design <- model.matrix(~0+factor(modType,levels = c("DN","Control"),ordered = F)).”



>fit <- lmFit(rt,design)

>cont.matrix<-makeContrasts(DN-Control,levels=design)

>fit2 <- contrasts.fit(fit, cont.matrix)

>fit2 <- eBayes(fit2)

>allDiff=topTable(fit2,adjust='fdr',number=200000)

#贝叶斯检验。

>diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]

#根据筛选条件筛选出差异基因。

>diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ]

#筛选上调的基因。

>diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ]

#筛选下调的基因。

##如果想把这些基因的结果输出,可以用write.table函数write.table(diffUp,file="up.xls",sep="\t",quote=F,row.names=F),要基因名字就把参数row.names设置为T.



##绘制热图

>group<-modType

#设置热图的分组。

>group<-as.data.frame(group)

>rownames(group)<-colnames(rt)

#按照分组,将rt的列名改为DN或control

>heat<-rt[rownames(rt) %in% c(head(rownames(subset(diffSig,diffSig$logFC>0)),20),head(rownames(subset(diffSig,diffSig$logFC<0)),20)),]

>library(pheatmap)

>x <- t(scale(t(heat)))

>pheatmap(x,annotation_col=group)

#绘制前20个基因的热图,自己调整基因数、图形颜色等等。


这里GSE1009的热图就绘制成功,实际处理时大家可根据自己的需求调整图形参数。

##绘制火山图

>png(file="火山图.png",width = 600,height = 600)

#设置图片格式

>yMax=max(-log10(allDiff$adj.P.Val))

#定义y轴最大值

>xMax=max(abs(allDiff$logFC))

#定义x轴最大值

>plot(allDiff$logFC,-log10(allDiff$adj.P.Val), xlab="logFC",ylab="-log10(adj.P.Val)",main="Volcano", xlim=c(-xMax,xMax),ylim=c(0,yMax),yaxs="i",pch=19, cex=1.2)

#绘图

>diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange)

>points( diffSub$logFC,-log10(diffSub$adj.P.Val), pch=19, col="red",cex=1.2)

#修改上调基因的图形参数

>diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange))

>points(diffSub$logFC,-log10(diffSub$adj.P.Val),  pch=19, col="green",cex=1.2)

#修改下调基因的图形参数

>abline(h=-log10(adjustP),lty=2,lwd=2)

>abline(v=c(-1,1),lty=2,lwd=2)

#拼接图形

>dev.off()

#返回终端。





所有代码如下:





整个分析中的变量如下:


可以看到,共有66个DEGs(diffsig),22个上调(diffup),44个上调(diffDown)

至此,DEGs分析以及火山图和热图就完啦,下一章是富集分析。

作者:NiKasu

原文链接:https://www.jianshu.com/p/44c879cff1cb

文章分类
后端
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 gxwowoo@163.com 举报,一经查实,本站将立刻删除。
相关推荐