GEO多芯片联合分析报错

#source("http://bioconductor.org/biocLite.R") #biocLite("limma")

logFoldChange=1 adjustP=0.05

library(limma) setwd("C:\Users\qiaog\Desktop\1\6.火山图") rt=read.table("normalize.txt",sep="\t",header=T,check.names=F) rt=as.matrix(rt) rownames(rt)=rt[,1] exp=rt[,2:ncol(rt)] dimnames=list(rownames(exp),colnames(exp)) rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)

#differential modType=c(rep("normal",10),rep("tumor",9),rep("normal",15),rep("tumor",15)) design <- model.matrix(~0+factor(modType)) colnames(design) <- c("con","treat") fit <- lmFit(rt,design) cont.matrix<-makeContrasts(treat-con,levels=design) fit2 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2)

allDiff=topTable(fit2,adjust='fdr',number=200000) write.table(allDiff,file="limmaTab.xls",sep="\t",quote=F,row.names=F)

#write table diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ] write.table(diffSig,file="diff.xls",sep="\t",quote=F,row.names=F) diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ] write.table(diffUp,file="up.xls",sep="\t",quote=F,row.names=F) diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ] write.table(diffDown,file="down.xls",sep="\t",quote=F,row.names=F)

#write expression level of diff gene **hmExp=rt[as.vector(diffSig[,1]),] 这里报错Error in rt[as.vector(diffSig[, 1]), ] : subscript out of bounds diffExp=rbind(id=colnames(hmExp),hmExp) write.table(diffExp,file="diffExp.txt",sep="\t",quote=F,col.names=F)

#volcano pdf(file="vol.pdf") xMax=max(-log10(allDiff$adj.P.Val)) yMax=max(abs(allDiff$logFC)) plot(-log10(allDiff$adj.P.Val), allDiff$logFC, xlab="-log10(adj.P.Val)",ylab="logFC", main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.8) diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange) points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.8) diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange)) points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="green",cex=0.8) abline(h=0,lty=2,lwd=3) dev.off()

你好,我是有问必答小助手。为了技术专家团更好地为您解答问题,烦请您补充下(1)问题背景详情,(2)您想解决的具体问题,(3)问题相关代码图片或者报错信息。便于技术专家团更好地理解问题,并给出解决方案。

您可以点击问题下方的【编辑】,进行补充修改问题。

img