| 以下是我按照公众号【云生信】进行的一项试验，把再整理的代码贴上来分享一下。还有几个问题想请教各位大神： Q1：画出芯片直方图、箱式图、散点图、热图等的目的是什么？如何从图中得到相应数据？ Q2：我用GCBI在线实验室做出的结果和用此代码得到的结果并不完全相同，那到底该如何选择方法呢？ Q3：散点图中all/pm/mm究竟代表了什么呢（见下方标红部分）？ Q4：在这个代码中，比如散点图是比较两个数据的，每次都要修改代码才能画出下一个图，我试想用for循环语句，自动画出所需要的20张图，但是最终呈现的结果都是最后一张图（我已经用par（mfrow=c（））函数确定在一张画布上保留多图）。求如何用for循环实现？ library(affy) exprsSet.RMA<-rma(AffyData) ##用rma方法处理数据 exp.RMA<-exprs(exprsSet.RMA)## exprsSet.MAS5<-mas5(AffyData) ##用mas5方法处理数据 write.exprs(exprsSet.RMA,file="RMA-processed data.txt") ##结果保存 par(mfrow=c(1,2)) ##直方图为标准化前后数据分布的比较，可以观察出基因的以2为底对数转换后的信号值集中在什么区间 hist(log2(intensity(AffyData[,4])),breaks=100,col="blue") ##画原始芯片数据的直方图，因为实验数据变化范围大所以取log2，结果集中在6~8，大于10 的很少，说明样本仅少部分高表达，与实际吻合。 hist(log2(exp.RMA),breaks=100,col="green") ##画RMA方法处理后数据的芯片数据直方图 par(mfrow=c(1,2)) ##箱式图为RMA方法消除样本间差异前后的数据比较，红色图代表处理前数据分布，蓝色为处理后数据分布，可见RMA方法可以有效去除样本间差异。 boxplot(AffyData,col="red") ##画箱式图，比较数据分布情况 boxplot(data.frame(exp.RMA),col="blue") par(mfrow=c(2,2)) ##散点图评价两个样本之间的数据差异 plot(exp.RMA[,1:2],log="xy",pch=".",main="all") ##all,pm,mm分别代表什么?? plot(pm(AffyData)[,1:2],log="xy",pch=".",main="pm") plot(mm(AffyData)[,1:2],log="xy",pch=".",main="mm") plot(exp.RMA[,1:2],pch=".",main="exp.RMA") library(genefilter) ##函数rowSds在genenfliter中 rsd<-rowSds(exp.RMA) ##计算行的标准差 sel<-order(rsd,decreasing=TRUE)[1:20] ##选择表达差异高的top20，order返回行号 heatmap(exp.RMA[sel,],col=rainbow(256)) ##画热图 ##获得了标准化后的表达谱数据，接下来我们可以利用limma package的算法对表达谱数据进行差异基因提取。 ##需要设置的参数包括表达谱数据的文件目录，case组和control组的数据，以及输出文件的文件名，这里默认为case_control_limma.txt。那么对于我们的样例文件，我们将其放置于F：\盘下，case组（treated with**）为9个，control组（treated without**）为10个 library(limma) setwd("F:\\") exp<-read.table("RMA-processed data.txt",header=T,sep="\t") rownames(exp)<-exp[,1] exp<-exp[2:20] exp<-as.matrix(exp) ##构建余下19个样品的基因表达矩阵（除去某些样本） samps<-factor(c(rep("case",9),rep("control",10))) ##提取实验条件信息 design<-model.matrix(~0+samps) ##构建实验设计矩阵 colnames(design)<-c("case","control") fit<-lmFit(exp,design) ##线性模型拟合 cont.matrix<-makeContrasts(case-control,levels=design) ##构建对比模型，比较两个实验条件下表达数据 fit2<-contrasts.fit(fit,cont.matrix) ##根据对比模型进行差值计算 fit2<-eBayes(fit2) ##贝叶斯检验 final<-topTable(fit2,coef=1,number=dim(exp)[1],adjust.method="BH",sort.by="B",resort.by="M") ##生成所有基因的检验结果报表 write.table(final,paste("case_control","_limma.txt"),quote=FALSE,sep="\t") dif<-final[final[,"P.Value"]<0.01,] ##根据p.Value对结果进行筛选，得到全部差异表达基因 head(dif) ##显示结果的前六行 write.table(dif,paste("case_control","_limma-dif.txt"),quote=FALSE,sep="\t") source("http://Bioconductor.org/biocLite.R") affydb biocLite("hgu133plus2.db")##下载特定基因芯片注释包 library(annotate) ##加载注释工具包 affydb<-annPkgName(AffyData@annotation,type="db") ##获得基因芯片注释包名称 library(affydb,character.only=T) ##加载注释包"hgu95av2.db"，必须设定character.only final$Symbols<-getSYMBOL(rownames(final),affydb)##根据每个探针组的ID获取对应的基因Gene Symbol,并作为一个新的列，加到数据框dif最后 final$EntrezID<-getEG(rownames(final),affydb) ##根据每个探针组的ID获取对应的基因Entrez ID,同样加到数据框dif最后 write.table(final,paste("case_control","_limma.txt"),quote=FALSE,sep="\t") ##final和dif可直接在代码中手动转换。 谢谢各位大神！ |
联系客服
