Re: Questions about annotations
68 views
Skip to first unread message
Yu Chuan
Nov 2, 2009, 8:36:08 PM11/2/09
to aroma.affymetrix
Hi all,
Follow up the below discussion....I am also trying to get the
chromosome locations for the genes on a HuEx-1_0-st-v2 dataset.
I used the core CDF file HuEx-1_0-st-v2,coreR3,A20071112,EP.cdf made
by Elizabeth for pre-processing using aroma.affymetrix. At the end I
got a matrix of transcript-level expression values with the first 5
columns/rows as
> trFit[1:5,1:5]
unitName groupName unit group cell
1 2315251 2315252 1 1 1
2 2315373 2315374 2 1 3
3 2315554 2315586 3 1 7
4 2315633 2315638 4 1 22
5 2315674 2315690 5 1 30
Based on what EP suggested below, is it correct that I should have
used the Ensembl cdf that Mark has created during the pre-processing
step if I want to use biomaRt package to get the gene annotations? Or,
is there any way that I can get the chromosome locations for each
unitName using biomaRt package directly based on what I already have
here? It seems that biomaRt package can also take Affy HuEx 1_0 st v2
IDs for searching Ensembl database, however, I don't think the
unitName above are Affy HuEx 1_0 st v2 IDs, right? I am still confused
about what's the most straightforward approach to get chromosome
locations using what I already have, and wonder if anyone has any
suggestion?
Thanks!
Yu Chuan
Hi Shuying,
So for your main question -- how to get a gene name for a transcript
cluster id -- you will probably need to use the Affymetrix annotation
for this, regardless of any flaws it might have. See their support
pages
for files with names like HuEx-1_0-st-v2.na24.hg18.transcript.csv.
While
I'm sure there are many solutions, I don't know of any other
*straight-forward* mapping. However, I am not familiar with what the
Genome Browser is providing these days.
If you are only working with the 25% of the array that make up the
core
probesets (so I assume you just want the well-annotated genome and
don't
want to search the extended and full probesets) you could instead try
the Ensembl cdf that Mark has created -- he has matched the Affy
probesets to the Ensembl gene/transcript annotation. This has the
advantage that you can take the unit name which is the Ensembl id
(e.g.
ENSG00000136527) and then use biomaRt package in R to get the gene
name.
Furthermore, Steffen Durinck has created a nice package that (among
other things) allow you to plot data from the exon array against the
transcript information in Ensembl. It will be released on bioconductor
soon. Moreover, genes with overlap are dealt with slightly differently
than Affymetrix, so I think there is less combining of separate genes
into one 'gene'. So there are a lot of advantages to using this cdf.
However, it will not give you a way to extend your analysis to the
other
75%(!) of the chip that are not included.
Another, more specific comment regarding your code. The *much* better
way to get the unit/group names to match the output of the
extractMatrix
command is to use extractDataFrame() command which has the (non-
default)
option to return the unit(gene) and group(exon) names.
Best,
Elizabeth
- Hide quoted text -
- Show quoted text -
Shuying wrote:
> Hi all,
> As a beginner, I copied some code (see the following)
> and did the preprocessing for my exon arrays. But I
> do not know how to intepret those normalized files.
> Here are my specific questions:
> As you see in the following code, I generated two
> files, one is to merge all the exon probe sets and
> get an overall gene expression level for a transcript
> cluster ID, that is the file "plmCore.BrCL.txt" and it
> has 18708 rows; another is to NOT merge those exon
> probe sets, file name is "plmCoreBrCL.NoMerge.txt",
> and I got 284258 rows. My goal is to map those exon
> probe sets and transcript cluster ID to official gene
> symbols (e.g. genes MYOD1, PR, ER and so on) for my
> further analysis. How can I do that? I know how to do
> this for affymetrix U133 data, but for the exon array
> data I do not know how to extract the gene ID (or
> transcript ID) information from the affy CDF file
> and/or annotation file. I searched UCSC genome browser
> web and found one file (named "affyHumanExon.txt") in
> their annotation directory, but I am not sure how will
> that help. I searched in this aroma.affymetrix group
> list, did not find exact answer. Any feedback is
> helpful.
> Thanks! and here is my code:
> library(aroma.affymetrix)
> verbose<-Arguments$getVerbose(-8)
> timestampOn(verbose)
> chipType <- "HuEx-1_0-st-v2"
> cdf <-
> AffymetrixCdfFile$fromChipType("HuEx-1_0-st-v2")
> print(cdf)
> # cdf <- AffymetrixCdfFile$fromChipType(chipType,
> tags="core")
> # print(cdf)
> cdfCore<-AffymetrixCdfFile$fromChipType(paste(chipType,"coreR3,A20071112,EP",sep=","))
> BC.exon <- AffymetrixCelSet$fromName(name=name,
> chipType="HuEx-1_0-st-v2" )
> setCdf(BC.exon, cdfCore)
> bc<-RmaBackgroundCorrection(BC.exon)
> system.time(process.bc<-process(bc, verbose=verbose) )
> # How come it works in my laptop, but not bilbo, is it
> becasue the version of R there is not correct?
> #quantile normalization:
> BcExon.quan <- QuantileNormalization(process.bc,
> typesToUpdate="pm",tags=c("*"))
> process.BcExon.quan <- process(BcExon.quan,
> verbose=verbose,force=F)
> # takes about 30 to 40 minutes
> #RMA model fitting:
> #----
> #Fit to core transcript clusters
> #----
> setCdf(BC.exon, cdfCore)
> plmCore.BrCL <- ExonRmaPlm(process.BcExon.quan ,
> mergeGroups=TRUE,tags=c("*","core"))
> fit(plmCore.BrCL, verbose=verbose)
> # About 50 minutes in my laptop
> #----
> #Fit to core probesets
> #----
> setCdf(BC.exon, cdfCore)
> plmCoreBrCL.NoMerge <- ExonRmaPlm(process.BcExon.quan,
> mergeGroups=FALSE,tags=c("*","core"))
> fit(plmCoreBrCL.NoMerge, verbose=verbose)
> # take about 70 minutes
> #----------------- output
> celNames <- paste(getNames(BC.exon),".CEL",sep="");
> header <- t(c("Probeset", celNames));
> ces<-getChipEffectSet(plmCoreBrCL.NoMerge);
> pbsetSumm<-extractMatrix(ces,returnUgcMap=TRUE,units=NULL);
> # pbsetSumm<-log2(pbsetSumm);
> pbsetSumm<-round(log2(pbsetSumm), 3) ; # Mar 31,
> notes, to save space, here I only save 3 digits
> pbsetNames<-readCdfGroupNames(getPathname(getCdf(ces)),unit=NULL,truncate=F)
> rownames(pbsetSumm)<-unlist(pbsetNames)
> write.table(header, file="plmCoreBrCL.NoMerge.txt",
> sep="\t", quote=FALSE, row.names=FALSE,
> col.names=FALSE);
> write.table(pbsetSumm, file="plmCoreBrCL.NoMerge.txt",
> sep="\t", quote=FALSE, append=TRUE, row.names=TRUE,
> col.names=FALSE);
> ces<-getChipEffectSet(plmCore.BrCL);
> pbsetSumm<-extractMatrix(ces,returnUgcMap=TRUE,units=NULL);
> # pbsetSumm<-log2(pbsetSumm);
> pbsetSumm<-round(log2(pbsetSumm), 3) ; # Mar 31,
> notes, to save space, here I only save 3 digits
> pbsetNames<-readCdfGroupNames(getPathname(getCdf(ces)),unit=NULL,truncate=F)
> rownames(pbsetSumm)<-names(pbsetNames)
> write.table(header, file="plmCore.BrCL.txt", sep="\t",
> quote=FALSE, row.names=FALSE, col.names=FALSE);
> write.table(pbsetSumm, file="plmCore.BrCL.txt",
> sep="\t", quote=FALSE, append=TRUE, row.names=TRUE,
> col.names=FALSE);
> save.image()
> ____________________________________________________________________________________
> You rock. That's why Blockbuster's offering you one month of Blockbuster Total Access, No Cost.
> http://tc.deals.yahoo.com/tc/blockbuster/text5.com
Follow up the below discussion....I am also trying to get the
chromosome locations for the genes on a HuEx-1_0-st-v2 dataset.
I used the core CDF file HuEx-1_0-st-v2,coreR3,A20071112,EP.cdf made
by Elizabeth for pre-processing using aroma.affymetrix. At the end I
got a matrix of transcript-level expression values with the first 5
columns/rows as
> trFit[1:5,1:5]
unitName groupName unit group cell
1 2315251 2315252 1 1 1
2 2315373 2315374 2 1 3
3 2315554 2315586 3 1 7
4 2315633 2315638 4 1 22
5 2315674 2315690 5 1 30
Based on what EP suggested below, is it correct that I should have
used the Ensembl cdf that Mark has created during the pre-processing
step if I want to use biomaRt package to get the gene annotations? Or,
is there any way that I can get the chromosome locations for each
unitName using biomaRt package directly based on what I already have
here? It seems that biomaRt package can also take Affy HuEx 1_0 st v2
IDs for searching Ensembl database, however, I don't think the
unitName above are Affy HuEx 1_0 st v2 IDs, right? I am still confused
about what's the most straightforward approach to get chromosome
locations using what I already have, and wonder if anyone has any
suggestion?
Thanks!
Yu Chuan
Hi Shuying,
So for your main question -- how to get a gene name for a transcript
cluster id -- you will probably need to use the Affymetrix annotation
for this, regardless of any flaws it might have. See their support
pages
for files with names like HuEx-1_0-st-v2.na24.hg18.transcript.csv.
While
I'm sure there are many solutions, I don't know of any other
*straight-forward* mapping. However, I am not familiar with what the
Genome Browser is providing these days.
If you are only working with the 25% of the array that make up the
core
probesets (so I assume you just want the well-annotated genome and
don't
want to search the extended and full probesets) you could instead try
the Ensembl cdf that Mark has created -- he has matched the Affy
probesets to the Ensembl gene/transcript annotation. This has the
advantage that you can take the unit name which is the Ensembl id
(e.g.
ENSG00000136527) and then use biomaRt package in R to get the gene
name.
Furthermore, Steffen Durinck has created a nice package that (among
other things) allow you to plot data from the exon array against the
transcript information in Ensembl. It will be released on bioconductor
soon. Moreover, genes with overlap are dealt with slightly differently
than Affymetrix, so I think there is less combining of separate genes
into one 'gene'. So there are a lot of advantages to using this cdf.
However, it will not give you a way to extend your analysis to the
other
75%(!) of the chip that are not included.
Another, more specific comment regarding your code. The *much* better
way to get the unit/group names to match the output of the
extractMatrix
command is to use extractDataFrame() command which has the (non-
default)
option to return the unit(gene) and group(exon) names.
Best,
Elizabeth
- Hide quoted text -
- Show quoted text -
Shuying wrote:
> Hi all,
> As a beginner, I copied some code (see the following)
> and did the preprocessing for my exon arrays. But I
> do not know how to intepret those normalized files.
> Here are my specific questions:
> As you see in the following code, I generated two
> files, one is to merge all the exon probe sets and
> get an overall gene expression level for a transcript
> cluster ID, that is the file "plmCore.BrCL.txt" and it
> has 18708 rows; another is to NOT merge those exon
> probe sets, file name is "plmCoreBrCL.NoMerge.txt",
> and I got 284258 rows. My goal is to map those exon
> probe sets and transcript cluster ID to official gene
> symbols (e.g. genes MYOD1, PR, ER and so on) for my
> further analysis. How can I do that? I know how to do
> this for affymetrix U133 data, but for the exon array
> data I do not know how to extract the gene ID (or
> transcript ID) information from the affy CDF file
> and/or annotation file. I searched UCSC genome browser
> web and found one file (named "affyHumanExon.txt") in
> their annotation directory, but I am not sure how will
> that help. I searched in this aroma.affymetrix group
> list, did not find exact answer. Any feedback is
> helpful.
> Thanks! and here is my code:
> library(aroma.affymetrix)
> verbose<-Arguments$getVerbose(-8)
> timestampOn(verbose)
> chipType <- "HuEx-1_0-st-v2"
> cdf <-
> AffymetrixCdfFile$fromChipType("HuEx-1_0-st-v2")
> print(cdf)
> # cdf <- AffymetrixCdfFile$fromChipType(chipType,
> tags="core")
> # print(cdf)
> cdfCore<-AffymetrixCdfFile$fromChipType(paste(chipType,"coreR3,A20071112,EP",sep=","))
> BC.exon <- AffymetrixCelSet$fromName(name=name,
> chipType="HuEx-1_0-st-v2" )
> setCdf(BC.exon, cdfCore)
> bc<-RmaBackgroundCorrection(BC.exon)
> system.time(process.bc<-process(bc, verbose=verbose) )
> # How come it works in my laptop, but not bilbo, is it
> becasue the version of R there is not correct?
> #quantile normalization:
> BcExon.quan <- QuantileNormalization(process.bc,
> typesToUpdate="pm",tags=c("*"))
> process.BcExon.quan <- process(BcExon.quan,
> verbose=verbose,force=F)
> # takes about 30 to 40 minutes
> #RMA model fitting:
> #----
> #Fit to core transcript clusters
> #----
> setCdf(BC.exon, cdfCore)
> plmCore.BrCL <- ExonRmaPlm(process.BcExon.quan ,
> mergeGroups=TRUE,tags=c("*","core"))
> fit(plmCore.BrCL, verbose=verbose)
> # About 50 minutes in my laptop
> #----
> #Fit to core probesets
> #----
> setCdf(BC.exon, cdfCore)
> plmCoreBrCL.NoMerge <- ExonRmaPlm(process.BcExon.quan,
> mergeGroups=FALSE,tags=c("*","core"))
> fit(plmCoreBrCL.NoMerge, verbose=verbose)
> # take about 70 minutes
> #----------------- output
> celNames <- paste(getNames(BC.exon),".CEL",sep="");
> header <- t(c("Probeset", celNames));
> ces<-getChipEffectSet(plmCoreBrCL.NoMerge);
> pbsetSumm<-extractMatrix(ces,returnUgcMap=TRUE,units=NULL);
> # pbsetSumm<-log2(pbsetSumm);
> pbsetSumm<-round(log2(pbsetSumm), 3) ; # Mar 31,
> notes, to save space, here I only save 3 digits
> pbsetNames<-readCdfGroupNames(getPathname(getCdf(ces)),unit=NULL,truncate=F)
> rownames(pbsetSumm)<-unlist(pbsetNames)
> write.table(header, file="plmCoreBrCL.NoMerge.txt",
> sep="\t", quote=FALSE, row.names=FALSE,
> col.names=FALSE);
> write.table(pbsetSumm, file="plmCoreBrCL.NoMerge.txt",
> sep="\t", quote=FALSE, append=TRUE, row.names=TRUE,
> col.names=FALSE);
> ces<-getChipEffectSet(plmCore.BrCL);
> pbsetSumm<-extractMatrix(ces,returnUgcMap=TRUE,units=NULL);
> # pbsetSumm<-log2(pbsetSumm);
> pbsetSumm<-round(log2(pbsetSumm), 3) ; # Mar 31,
> notes, to save space, here I only save 3 digits
> pbsetNames<-readCdfGroupNames(getPathname(getCdf(ces)),unit=NULL,truncate=F)
> rownames(pbsetSumm)<-names(pbsetNames)
> write.table(header, file="plmCore.BrCL.txt", sep="\t",
> quote=FALSE, row.names=FALSE, col.names=FALSE);
> write.table(pbsetSumm, file="plmCore.BrCL.txt",
> sep="\t", quote=FALSE, append=TRUE, row.names=TRUE,
> col.names=FALSE);
> save.image()
> ____________________________________________________________________________________
> You rock. That's why Blockbuster's offering you one month of Blockbuster Total Access, No Cost.
> http://tc.deals.yahoo.com/tc/blockbuster/text5.com
Mark Robinson
Nov 6, 2009, 7:14:59 AM11/6/09
to aroma-af...@googlegroups.com
Hi Yu Chuan.
Some comments below.
On 3-Nov-09, at 12:36 PM, Yu Chuan wrote:
> Hi all,
>
> Follow up the below discussion....I am also trying to get the
> chromosome locations for the genes on a HuEx-1_0-st-v2 dataset.
Have a look at the section "Current NetAffx Annotation Files" at:
http://www.affymetrix.com/products_services/arrays/specific/exon.affx#1_4
Those CSV files go a long way to getting you what you need.
> I used the core CDF file HuEx-1_0-st-v2,coreR3,A20071112,EP.cdf made
> by Elizabeth for pre-processing using aroma.affymetrix. At the end I
> got a matrix of transcript-level expression values with the first 5
> columns/rows as
>
>> trFit[1:5,1:5]
> unitName groupName unit group cell
> 1 2315251 2315252 1 1 1
> 2 2315373 2315374 2 1 3
> 3 2315554 2315586 3 1 7
> 4 2315633 2315638 4 1 22
> 5 2315674 2315690 5 1 30
>
> Based on what EP suggested below, is it correct that I should have
> used the Ensembl cdf that Mark has created during the pre-processing
> step if I want to use biomaRt package to get the gene annotations?
It doesn't matter what CDF you use, you can still use biomaRt/
GenomeGraphs to pull down annotation information. You just need
chromosome coordinates, which you can get from the CSV file from
Affy. Or, see below.
> Or,
> is there any way that I can get the chromosome locations for each
> unitName using biomaRt package directly based on what I already have
> here? It seems that biomaRt package can also take Affy HuEx 1_0 st v2
> IDs for searching Ensembl database, however, I don't think the
> unitName above are Affy HuEx 1_0 st v2 IDs, right?
You'll have to look up the biomaRt details to do this exactly, but
here is a quick example (using groupNames i.e. probeset ids, not
transcript cluster ids):
> mart<-useDataset("hsapiens_gene_ensembl",mart)Checking
attributes ... ok
Checking filters ... ok
>
getBM
(attributes
=
c
("affy_huex_1_0_st_v2
","hgnc_symbol","strand","transcript_start","chromosome_name"),
filters="affy_huex_1_0_st_v2",values=c("2315252","2315374"), mart=mart)
affy_huex_1_0_st_v2 hgnc_symbol strand transcript_start
chromosome_name
1 2315374 FAM87B 1
752751 1
2 2315374 FAM87A -1
325934 8
Cheers,
Mark
> --
> When reporting problems on aroma.affymetrix, make sure 1) to run the
> latest version of the package, 2) to report the output of
> sessionInfo() and traceback(), and 3) to post a complete code example.
>
>
> You received this message because you are subscribed to the Google
> Groups "aroma.affymetrix" group.
> To post to this group, send email to aroma-af...@googlegroups.com
> To unsubscribe from this group, send email to aroma-affymetr...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/aroma-affymetrix?hl=en
------------------------------
Mark Robinson, PhD (Melb)
Epigenetics Laboratory, Garvan
Bioinformatics Division, WEHI
e: m.rob...@garvan.org.au
e: mrob...@wehi.edu.au
p: +61 (0)3 9345 2628
f: +61 (0)3 9347 0852
------------------------------
Some comments below.
On 3-Nov-09, at 12:36 PM, Yu Chuan wrote:
> Hi all,
>
> Follow up the below discussion....I am also trying to get the
> chromosome locations for the genes on a HuEx-1_0-st-v2 dataset.
http://www.affymetrix.com/products_services/arrays/specific/exon.affx#1_4
Those CSV files go a long way to getting you what you need.
> I used the core CDF file HuEx-1_0-st-v2,coreR3,A20071112,EP.cdf made
> by Elizabeth for pre-processing using aroma.affymetrix. At the end I
> got a matrix of transcript-level expression values with the first 5
> columns/rows as
>
>> trFit[1:5,1:5]
> unitName groupName unit group cell
> 1 2315251 2315252 1 1 1
> 2 2315373 2315374 2 1 3
> 3 2315554 2315586 3 1 7
> 4 2315633 2315638 4 1 22
> 5 2315674 2315690 5 1 30
>
> Based on what EP suggested below, is it correct that I should have
> used the Ensembl cdf that Mark has created during the pre-processing
> step if I want to use biomaRt package to get the gene annotations?
GenomeGraphs to pull down annotation information. You just need
chromosome coordinates, which you can get from the CSV file from
Affy. Or, see below.
> Or,
> is there any way that I can get the chromosome locations for each
> unitName using biomaRt package directly based on what I already have
> here? It seems that biomaRt package can also take Affy HuEx 1_0 st v2
> IDs for searching Ensembl database, however, I don't think the
> unitName above are Affy HuEx 1_0 st v2 IDs, right?
here is a quick example (using groupNames i.e. probeset ids, not
transcript cluster ids):
> mart<-useDataset("hsapiens_gene_ensembl",mart)Checking
attributes ... ok
Checking filters ... ok
>
getBM
(attributes
=
c
("affy_huex_1_0_st_v2
","hgnc_symbol","strand","transcript_start","chromosome_name"),
filters="affy_huex_1_0_st_v2",values=c("2315252","2315374"), mart=mart)
affy_huex_1_0_st_v2 hgnc_symbol strand transcript_start
chromosome_name
1 2315374 FAM87B 1
752751 1
2 2315374 FAM87A -1
325934 8
Cheers,
Mark
> When reporting problems on aroma.affymetrix, make sure 1) to run the
> latest version of the package, 2) to report the output of
> sessionInfo() and traceback(), and 3) to post a complete code example.
>
>
> You received this message because you are subscribed to the Google
> Groups "aroma.affymetrix" group.
> To post to this group, send email to aroma-af...@googlegroups.com
> To unsubscribe from this group, send email to aroma-affymetr...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/aroma-affymetrix?hl=en
------------------------------
Mark Robinson, PhD (Melb)
Epigenetics Laboratory, Garvan
Bioinformatics Division, WEHI
e: m.rob...@garvan.org.au
e: mrob...@wehi.edu.au
p: +61 (0)3 9345 2628
f: +61 (0)3 9347 0852
------------------------------
Yu Chuan
Nov 8, 2009, 1:09:46 PM11/8/09
to aroma.affymetrix
Hi Mark,
Thanks a lot!! Very useful. I think I know what to do now.
Best,
Yu Chuan
On Nov 6, 4:14 am, Mark Robinson <mrobin...@wehi.EDU.AU> wrote:
> Hi Yu Chuan.
>
> Some comments below.
>
> On 3-Nov-09, at 12:36 PM, Yu Chuan wrote:
>
> > Hi all,
>
> > Follow up the below discussion....I am also trying to get the
> > chromosome locations for the genes on a HuEx-1_0-st-v2 dataset.
>
> Have a look at the section "Current NetAffx Annotation Files" at:http://www.affymetrix.com/products_services/arrays/specific/exon.affx...
> > For more options, visit this group athttp://groups.google.com/group/aroma-affymetrix?hl=en
> e: mrobin...@wehi.edu.au
Thanks a lot!! Very useful. I think I know what to do now.
Best,
Yu Chuan
On Nov 6, 4:14 am, Mark Robinson <mrobin...@wehi.EDU.AU> wrote:
> Hi Yu Chuan.
>
> Some comments below.
>
> On 3-Nov-09, at 12:36 PM, Yu Chuan wrote:
>
> > Hi all,
>
> > Follow up the below discussion....I am also trying to get the
> > chromosome locations for the genes on a HuEx-1_0-st-v2 dataset.
>
>
> ------------------------------
> Mark Robinson, PhD (Melb)
> Epigenetics Laboratory, Garvan
> Bioinformatics Division, WEHI
> e: m.robin...@garvan.org.au
> ------------------------------
> Mark Robinson, PhD (Melb)
> Epigenetics Laboratory, Garvan
> Bioinformatics Division, WEHI
> e: mrobin...@wehi.edu.au
Reply all
Reply to author
Forward
0 new messages