This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Analysis of microarray data: "eset" problem

Hello,

I was performing analysis of microarray data according to the tutorial from Center for Research Informatics. There is a step-by-step description and codes, but I found an error and I cannot solve my problem on my own. So, I will be greatful for any help.

The codes look like:

setwd('~/Desktop/test')

source("http://www.bioconductor.org/biocLite.R")
biocLite("affy")
biocLite("affycoretools")
biocLite("GEOquery")

library(affy)
library(affycoretools)
library(GEOquery)

mydata<-ReadAffy()
pData(mydata)<-read.table("phenod.txt", header=T, row.names=1, sep="\t")

eset<-rma(mydata)
eset

write.exprs(eset, file="Expression_values.xls")

biocLite("limma")
library("limma")

Group<-factor(pData(mydata)[,1], levels=levels(pData(mydata)[,1]))
design<-model.matrix(~Group)

fit<-lmFit(eset,design)

fit<-eBayes(fit)

tab<-topTable(fit, coef=2, adjust="fdr", n=50)

# annotation? ------> here is the problem: After below code I get the error:

Error in rows[i] : only 0's may be mixed with negative subscripts
eset2 <- eset[tab[,1]]

So, what to do in such situation?

Best regards!

software-error r
• 3,977 views • link
updated by Ram Ram • written by orzech_mag orzech_mag

What's the output of quantile(tab[,1])? I don't think topTable produces a data.frame whose first column is guaranteed to be a proper row index number (especially with n=50).

• 0 views • link
written by Devon Ryan Devon Ryan

The output is:

> quantile(tab[,1])
         0%         25%         50%         75%        100% 
-0.85955683 -0.53124254 -0.01184481  0.24215869  1.80077207

But, please, look at the subsequent codes which are using the eset2<-eset[tab[,1]] command:

library(hgu133a.db)
library(annotate)
library(R2HTML)

ls("package:hgu133a.db")

ID<-featureNames(eset2)

Symbol<-getSYMBOL(ID, "hgu133a.db")
Name<-as.character(lookUp(ID, "hgu133a.db", "GENENAME"))
Ensembl<-as.character(lookUp(ID, "hgu133a.db", "ENSEMBL"))

Ensembl<- ifelse(Ensembl=="NA", NA,
                 paste("<a href='http://useast.ensembl.org/Homo_sapiens/Gene/Summary?g=",
                       Ensembl, "'>",Ensembl, "</a>", sep=""))
tmp<-data.frame(ID=ID, Symbol=Symbol, Name=Name, Ensembl=Ensembl, stringsAsFactors=F)
tmp[tmp=="NA"]<-NA #poprawia błąd komendy stringsAsFactors
HTML(tmp, "out.html", append=F)
write.table(tmp, file="target.txt", row.names=F, sep="\t")
• 0 views • link
updated by Ram Ram • written by orzech_mag orzech_mag

1 answer

So the problem is that you're using non-row indices as row indices. That's not going to work well. What you want to do is something along the lines of:

tab<-topTable(fit, coef=2, adjust="fdr", n=inf, sort.by="none")
eset2 <- eset[which(tab$padj<0.1)] #Using which() in case there are NAs

Something along those lines should work better for you. If you really do just want to the top 50 or so then you need to match the gene names/IDs/probe annotations/whatever from the columns of tab to the info held in eset. That will allow you to get row numbers and to then subset eset.

• 0 views • link
written by Devon Ryan Devon Ryan

It works now. Thank you very much!

• 0 views • link
written by orzech_mag orzech_mag

I want to read phenodata. but I don't know where to get this phenod.txt file... it would be great if you help me.

• 0 views • link
written by 1234anjalianjali1234 • 1234anjalianjali1234

It's not something you get from somewhere, you need to create it.

• 0 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this