Get the most expressed genes from one .CEL file in R.

Question

[alevax]

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! :)

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

[Sean Davis]

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.

[zx8754]

http://stackoverflow.com/questions/25724486/get-the-most-expressed-genes-from-one-cel-file-in-r

Answer 1

[Sean Davis]

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

[alevax]

Ok, I understand.

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

[Sean Davis]

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

[alevax]

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?

[Sean Davis]

See these papers for some background on expressed/unexpressed gene calling:

• http://www.ncbi.nlm.nih.gov/pubmed/21177656
• http://www.ncbi.nlm.nih.gov/pubmed/24271388

[alevax]

Thank you. They were interesting readings.