Analysis of microarray data: "eset" problem
1
0
Entering edit mode
7.0 years ago
orzech_mag ▴ 230

Hallo,

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!

 

R software error • 2.7k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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")
ADD REPLY
0
Entering edit mode
7.0 years ago

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.

ADD COMMENT
0
Entering edit mode

It works now. Thank you very much!

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 2294 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6