Analysis of microarray data: "eset" problem

Question

[orzech_mag]

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!

[Devon Ryan]

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).

[orzech_mag]

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")

Answer 1

[Devon Ryan]

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.

[orzech_mag]

It works now. Thank you very much!

[1234anjalianjali1234]

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.

[Devon Ryan]

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