This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get the most expressed genes from one .CEL file in R.

Hi everybody!

LIMMA in R can give you a list of differentially expressed genes.

How can I simply get the probesets with highest signal intensity?

Can I get only the most expressed genes in an healty experiment, for example from one .CEL file?

Or the most expressed genes from a set of .CEL files of the same group (all of the control group, or all of the sample group).

Thanks! :)

limma r differential-expression
• 4,993 views • link
updated by Ram Ram • written by alevax • alevax 
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("GEOquery","affy","limma","gcrma"))
gse_number <- "GSE13887"
getGEOSuppFiles( gse_number )
COMPRESSED_CELS_DIRECTORY <- gse_number
untar( paste( gse_number , paste( gse_number , "RAW.tar" , sep="_") , sep="/" ), exdir=COMPRESSED_CELS_DIRECTORY)
cels <- list.files( COMPRESSED_CELS_DIRECTORY , pattern = "[gz]")
sapply( paste( COMPRESSED_CELS_DIRECTORY , cels, sep="/") , gunzip )
celData <- ReadAffy( celfile.path = gse_number )

gcrma.ExpressionSet <- gcrma(celData)

if(!require(panp)) { biocLite("panp") }
library(panp)
myGDS <- getGEO("GDS2697")
eset <- GDS2eSet(myGDS,do.log2=TRUE)
my_pa <- pa.calls(eset)

is(eset)
is(gcrma.ExpressionSet)

If you run the script you'll get an error while executing:

my_pa <- pa.calls(eset)

and not while executing

my_pa <- pa.calls(gcrma.ExpressionSet)

Why if they are both ExpressionSet?

> is(gcrma.ExpressionSet)
[1] "ExpressionSet"    "eSet"             "VersionedBiobase" "Versioned"
> is(eset)
[1] "ExpressionSet"    "eSet"             "VersionedBiobase" "Versioned"
• 0 views • link
updated by Ram Ram • written by alevax • alevax

The GDS does not contain enough information to make PA calls. The PA calls need match and mismatch probe data and the GDS contains only normalized probeset data, likely derived from match probes only.

• 0 views • link
updated by Ram Ram • written by Sean Davis Sean Davis

1 answer

The order() function in R is perhaps a useful place to start. Note, though, that microarrays measure relative expression between samples and not absolute expression, so simply ordering by intensity is not really a quantitative way to get the most expressed genes (though qualitatively it might do).

• 1 views • link
updated by Ram Ram • written by Sean Davis Sean Davis

Ok, I understand.

So, how can Iget the absolute expression of one microarray sample?

• 1 views • link
written by alevax • alevax

Unfortunately, microarrays (and RNA-seq, for that matter) do not measure absolute expression, so you cannot.

• 0 views • link
updated by Ram Ram • written by Sean Davis Sean Davis

Well, but if I take 2 .CEL files of the same group (e.g. 2 of the control group) and I perform the GCRMA obtaining the ExpessionSet object after the GCRMA, the gcrma.ExpressionSet, the pa.calls() works.

Thus, Should not I trust its result?

• 1 views • link
updated by Ram Ram • written by alevax • alevax

Thank you. They were interesting readings.

• 0 views • link
written by alevax • alevax

Log in to answer this question.

More posts like this