[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