|
一、举例回顾 本节下载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分析以及火山图和热图就完啦，下一章是富集分析。 |
|
|
