Question: Analysis of microarray data: "eset" problem
0
gravatar for orzech_mag
4.6 years ago by
orzech_mag200
Poland/Łódź
orzech_mag200 wrote:

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 • 1.8k views
ADD COMMENTlink modified 4.6 years ago by Devon Ryan91k • written 4.6 years ago by orzech_mag200

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 REPLYlink written 4.6 years ago by Devon Ryan91k

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 REPLYlink written 4.6 years ago by orzech_mag200
0
gravatar for Devon Ryan
4.6 years ago by
Devon Ryan91k
Freiburg, Germany
Devon Ryan91k wrote:

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 COMMENTlink written 4.6 years ago by Devon Ryan91k

It works now. Thank you very much!

ADD REPLYlink written 4.6 years ago by orzech_mag200

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 REPLYlink written 3.2 years ago by 1234anjalianjali123430

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

ADD REPLYlink written 3.2 years ago by Devon Ryan91k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1824 users visited in the last hour