[BioC] problems with Affymetrix probes id using exonmap and xmapcore

Gomez Moruno, Antonio agomezm at idibell.cat
Wed May 4 11:17:36 CEST 2011


Hi Tim,

sorry i am very busy and i had no time to answer you, by the way, that's my code

library(gplots) 
#Load libraries
library(exonmap)
library(plier)
library(limma)
library(hugene10stprobeset.db)
library(hugene10sttranscriptcluster.db)


raw.data<-read.exon(covdesc="covdesc")
dat <- ReadAffy(cdfname="hugene10stv1cdf") #alternative using hugene10stv1cdf in order to avoid problems with the probesetidd
####DON'T USE IT!!!raw.data at cdfName<-"exon.pmcdf"# which CDF file defines the array by setting the name of the .cdf
#checking consistency of the data

summary(exprs(raw.data))

#RMA normalization and generate probeset summaries

x.rma<-rma(raw.data)

#PLIER normalization and generate probeset summaries
x.pli<-justPlier(raw.data, usemm=F, normalize=T, norm.type="pmonly")

#calculating fold changes and p-values

pc.rma <-pc(x.rma,"group",c("a","b"))## for rma normalization
pc.pli <-pc(x.pli,"group",c("a","b"))## for plier normalization

#find all probesets with an absolute fold change greater than 2 and unadjusted p-value less than 10e-4

keep.data<-(abs(fc(pc.rma))>2) &tt(pc.rma)<1e-4### only 2

keep.data.rma<-(abs(fc(pc.rma))>2) &tt(pc.rma)<1e-2### >260 Use this!!!!!!!!!!! For rma
keep.data.pli<-(abs(fc(pc.pli))>2) &tt(pc.pli)<1e-2### >106 Use this!!!!!!!!!!! for plier

#fishes thenames of these probesets out form raw.data objects

sigs.rma<-featureNames(x.rma)[keep.data.rma]### Caution!!! must be determined if this command exists in xmapcore or has been overwrited/substituted!!!!

sigs.pli<-featureNames(x.pli)[keep.data.pli]### Caution!!! must be determined if this command exists in xmapcore or has been overwrited/substituted!!!!

####ALTERNATIVELY CALCULATING FOLD CHANGE USING LIMMA

design<-as.matrix(cbind(DKO=c(1,1,1,1,1,1,0,0,0,0,0,0),WT=c(0,0,0,0,0,0,1,1,1,1,1,1)))

fit.rma<-lmFit(x.rma,design)## For rma
fit.plier<-lmFit(x.pli,design)# For plier

cont.matrix<-makeContrasts(DKOvsWT=DKO-WT,levels=design)

fit2.rma<-contrasts.fit(fit.rma,cont.matrix)
fit2.pli<-contrasts.fit(fit.plier,cont.matrix)

fit3.rma<-eBayes(fit2.rma)
fit3.pli<-eBayes(fit2.pli)

result.fold.limma.rma<-decideTests(fit3.rma, lfc=2)
result.fold.limma.pli<-decideTests(fit3.pli, lfc=2)

keep.limma.rma<-result.fold.limma.rma@".Data" !=0
keep.limma.pli<-result.fold.limma.pli@".Data" !=0


limma.sigs.rma<-rownames(result.fold.limma.rma@".Data")[keep.limma.rma]
limma.sigs.pli<-rownames(result.fold.limma.pli@".Data")[keep.limma.pli]

#####
##### Exporting the data in order to be used with xmapcore script

v<-as.data.frame(limma.sigs.pli)#### coercing the vector into a data frame
sink("significant_probesets_limma_plier.csv")
write.table(v, quote = FALSE, sep = ",")# exporting the data as csv file
sink()


detach(exonmap)

#####MAPPING TO ANNOTATION
##### Which transcripts are differentially expressed???
library(xmapcore)
library(exon.pmcdf)
Sys.setenv(R_XMAP_CONF_DIR="~/.exonmap")#configuring XMAP_DIR
xmap.connect()##connecting to the database

probeset.to.trancript(sigs) #### HERE ARE THE NULL RESULTS!!!!!!!!!!!


I guess that the problem is that i'm using the results from a Human Gene 1.0 ST array and i don't know if exonmap is able to handle this data in order to find exons. I checked the probesets that i obtained using featureNames in HuGene-1_0-st-v1.na30.1.hg19.probeset files and i found that the code obtained is    transcript_cluster_id and no from probesets. Is xmapcore able to handle those transcripts ids or may i have to use another package as xps or affygui in order to translate this ids to probesets or gene ids previous to the use of xmapcore?

Thanks in advance


________________________________________
De: Tim Yates [tyates at picr.man.ac.uk]
Enviat el: diumenge, 1 / maig / 2011 12:22
Per a: Gomez Moruno, Antonio; bioconductor at r-project.org
Tema: Re: [BioC] problems with Affymetrix probes id using exonmap and xmapcore

Hi again

Sorry for the late responses, it's a big public holiday here in the UK...

Can you give an example of the code which is returning NULL?

Cheers,

Tim


On 28/04/2011 16:05, "Gomez Moruno, Antonio" <agomezm at idibell.cat> wrote:

> Hi Tim,
>
> the problem is that i need to process the data previously with exonmap and
> then try to annotate this data with xmapcore. Could it be possible? From your
> reply i assume that i can normalize the data and work with probesets using
> only xmapcore?? But, i read the doc concerning xmapcore and  i didn't find
> nothing about it.
>
> I tried to use another Affymetrix related package, oligo, and i obtained this
> info about my data
>
> ExpressionSet (storageMode: lockedEnvironment)
> assayData: 257430 features, 12 samples
>   element names: exprs
> protocolData
>   rowNames: IG100120JMI035-3gst_01.CEL IG100120JMI036-3gst_01.CEL ...
>     IG100120JMI046-3gst_01.CEL (12 total)
>   varLabels: exprs dates
>   varMetadata: labelDescription channel
> phenoData
>   rowNames: IG100120JMI035-3gst_01.CEL IG100120JMI036-3gst_01.CEL ...
>     IG100120JMI046-3gst_01.CEL (12 total)
>   varLabels: index
>   varMetadata: labelDescription channel
> featureData
>   featureNames: 7892501 7892502 ... 8180418 (257430 total)
>   fvarLabels: probesetid seqname ... probesettype (39 total)
>   fvarMetadata: labelDescription
> experimentData: use 'experimentData(object)'
> Annotation: pd.hugene.1.0.st.v1
>
> So, the numbers that i got here are the same that i got in exonmpa, how can i
> handle this on xmapcore??
>
> Thanks
>
> ________________________________________
> De: Tim Yates [TYates at picr.man.ac.uk]
> Enviat el: dijous, 28 / abril / 2011 15:20
> Per a: Gomez Moruno, Antonio
> Tema: Re: [BioC] problems with Affymetrix probes id using exonmap and xmapcore
>
> Hiya,
>
> You seem to have both the deprecated exonmap package, and the xmapcore package
> loaded at the same time.
>
> This will cause misbehavior at best, and crashes at worst.
>
> Can you try with just xmapcore loaded, and if that still doesn't work, can you
> post a failing example of your code?
>
> Cheers,
>
> Tim
>
>
>
> ----- Reply message -----
> From: "Gomez Moruno, Antonio" <agomezm at idibell.cat>
> Date: Thu, Apr 28, 2011 08:58
> Subject: [BioC] problems with Affymetrix probes id using exonmap and xmapcore
> To: "bioconductor at r-project.org" <bioconductor at r-project.org>
>
> I got several .CEL files to do RMA using exonmap library, when i got my list
> significant differentially expressed probesets and i try to use xmapcore to
> map to annotation using "exonic", all i got is a NULL vector. I used xmapcore
> examples and all worked perfectly, but my list of Probesets have no
> correspondance in xmapcore. Anyone has any idea?
>
> This is my session info
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
>  [1] hugene10stv1cdf_2.8.0                hugene10sttranscriptcluster.db_7.0.1
>  [3] hugene10stprobeset.db_7.0.1          org.Hs.eg.db_2.5.0
>  [5] RSQLite_0.9-4                        AnnotationDbi_1.14.1
>  [7] limma_3.8.1                          plier_1.22.0
>  [9] exon.pmcdf_1.1                       xmapcore_1.6.0
> [11] IRanges_1.10.0                       RMySQL_0.7-5
> [13] DBI_0.2-5                            RColorBrewer_1.0-2
> [15] genefilter_1.34.0                    affy_1.30.0
> [17] Biobase_2.12.1                      exonmap_2.0.03
>
> loaded via a namespace (and not attached):
> [1] affyio_1.20.0         annotate_1.30.0       digest_0.4.2
> [4] preprocessCore_1.14.0 splines_2.13.0
> [7] survival_2.36-8       tools_2.13.0          xtable_1.5-6
>
> this is my raw data
>
> AffyBatch object
> size of arrays=1050x1050 features (21 kb)
> cdf=exon.pmcdf (1411189 affyids)
> number of samples=12
> number of genes=1411189
> annotation=hugene10stv1
> notes=
>
>
>
>
> Thanks in advance
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
> --------------------------------------------------------
> This email is confidential and intended solely for the...{{dropped:29}}



More information about the Bioconductor mailing list