[BioC] GO with R and Bioconductor

[Duke]

On 2/28/11 2:31 PM, Vincent Carey wrote:
> On Mon, Feb 28, 2011 at 2:07 PM, Duke<duke.lists at gmx.com>  wrote:
>> On 2/26/11 8:09 PM, Vincent Carey wrote:
>>> On Sat, Feb 26, 2011 at 6:40 PM, Duke<duke.lists at gmx.com>    wrote:
>>>> Hi Vincent,
>>>>
>>>> On 2/26/11 2:38 PM, Vincent Carey wrote:
>>>>> On Sat, Feb 26, 2011 at 12:00 PM,<duke.lists at gmx.com>      wrote:
>>>>>> Dear colleagues,
>>>>>>
>>>>>>   I used to download GO database from geneoncology.org and did some c++
>>>>>> coding to manipulate the data as I wished. Now I want to try my luck
>>>>>> with R
>>>>>> - Bioconductor. I have heard of tons of tools supporting GO such as
>>>>>> GO.db,
>>>>>> topGO, goseq, GOstats, biomart etc... and I have been reading their
>>>>>> description and examples, but honestly I am overhelmed and dont really
>>>>>> know
>>>>>> which package I should use to fulfill my task. So please advise me how
>>>>>> I can
>>>>>> do the following two simple tasks:
>>>>>>
>>>>>>   1. I have a list of genes (with gene names from UCSC such as Foxp3
>>>>>> etc...). How do I filter this list to get genes that have certain GO
>>>>>> term
>>>>>> such as transcription factor?
>>>>> since you said it was a simple task, consider the simple solution
>>>>> involving the "%annTo%" operator, which tells whether the symbols on
>>>>> the left have been annotated to the term on the right:
>>>>>
>>>>>> c("FOXP3", "BRCA2") %annTo% "mammary"
>>>>> FOXP3 BRCA2
>>>>> FALSE  TRUE
>>>>>> c("FOXP3", "BRCA2") %annTo% "transcription factor"
>>>>> FOXP3 BRCA2
>>>>>   TRUE FALSE
>>>>>
>>>>> you could use the named logical vectors generated in this way to
>>>>> perform the filtering you describe.  but see below.
>>>>>
>>>>>>   2. How do I know the capacity of the latest GO database on
>>>>>> bioconductor,
>>>>>> for example, how many genes available for mm9, and how many of them
>>>>>> have GO
>>>>>> term transcription factor?
>>>>> The "GO database" concerns the gene ontology, a structure of terms and
>>>>> relationships among them.  The association of GO terms to gene names
>>>>> for mouse is presented in various ways, but the most basic one is in
>>>>> the org.Mm.eg.db package.  With that, you could use
>>>>>
>>>>> org.Mm.eg()
>>>>>
>>>>> to find, among other statistics,
>>>>>
>>>>> org.Mm.egGO has 29984 mapped keys (of 63329 keys).  Your question
>>>>> concerning transcription factor mapping is not completely precise, and
>>>>> you might want to survey the family of GO terms to come up with a set
>>>>> of terms that meets your requirement.
>>>> You are right, I did not make it clear. What I wanted was exactly what
>>>> you
>>>> said: I need to get a list of genes that have a certain set of terms (for
>>>> example *transcription*) - which corresponds to a family of GO terms, not
>>>> just one GO term.
>>>>
>>>>>    Here's a
>>>>> demonstration of related queries:
>>>>>
>>>>>> get("GO:0003700", GOTERM)
>>>>> GOID: GO:0003700
>>>>> Term: sequence-specific DNA binding transcription factor activity
>>>>> Ontology: MF
>>>>> Definition: Interacting selectively and non-covalently with a specific
>>>>>      DNA sequence in order to modulate transcription. The transcription
>>>>>      factor may or may not also interact selectively with a protein or
>>>>>      macromolecular complex.
>>>>> Synonym: GO:0000130
>>>>> Secondary: GO:0000130
>>>>>> tfg = get("GO:0003700", revmap(org.Mm.egGO))
>>>>>> length(tfg)
>>>>> [1] 940
>>>> This works fine in case of one GO term (GO: 0003700). Is there a similar
>>>> function like getterm("transcription", revmap()) to get all the genes
>>>> that
>>>> their GO terms contain *transcription*?
>>> As far as I know there is no such function.  Solutions depend on the
>>> type of wild-card you are looking for.  If SQL LIKE operator is
>>> adequate
>>>
>>> termlikeToTags = function (liketok)
>>> {
>>>      require(GO.db)
>>>      gcon = GO_dbconn()
>>>      as.character(dbGetQuery(gcon, paste("select go_id from go_term
>>> where term like '%",
>>>          liketok, "%'", sep = ""))[[1]])
>>> }
>>>
>>> if you want to use full regular expressions
>>>
>>> reToTags = function (retok, ...) {
>>>      require(GO.db)
>>>      allt = dbGetQuery(GO_dbconn(), "select  go_id, term from go_term")
>>>      inds = grep(retok, as.character(allt[,"term"]), value=FALSE, ...)
>>> # will error if value set at call
>>>      if (length(inds)>0) return(as.character(allt[inds,"go_id"]))
>>>      stop("retok did not grep to any term")
>>> }
>>>
>>>> length(reToTags("transcription"))
>>> [1] 372
>>>> length(termlikeToTags("transcription"))
>>> [1] 372
>>>> length(reToTags("transcrip.*"))
>>> [1] 411
>>>> length(termlikeToTags("transcrip.*"))
>>> [1] 0
>>>
>> As far as I understand, these functions will list out all the GO Terms
>> containing "transcription" for example. In order to search for all of the
>> genes in the database that have these GO Terms, I suppose I will have to
>> loop
>>
>> tfg = get("GO:0003700", revmap(org.Mm.egGO))
>> length(tfg)
>>
>> for each of the term and then sum them up? Since your %annTo% function is to
>> check a certain gene list to see if any of them has the specified term,
>> would it be better to just run that function to the available genes in the
>> database? How do I list and use the gene list in the database?
> I'd suggest you read the vignettes at
> http://www.bioconductor.org/help/bioc-views/release/bioc/html/AnnotationDbi.html,
> particularly http://www.bioconductor.org/packages/2.7/bioc/vignettes/AnnotationDbi/inst/doc/AnnotationDbi.pdf
>
> The codes I provided are hints in the direction of solutions for the
> family of questions you pose.
>
>>>>> org.Mm.egGO is a mapping from mouse entrez gene ids to GO term tags.
>>>>> revmap reverses this mapping and takes a tag to the set of genes
>>>>> mapped to the tag by entrez.
>>>>>
>>>>> Now, to return to the first question -- it isn't simple and a lot of
>>>>> presuppositions have to be made explicit.  One of the most problematic
>>>>> is the commitment to use gene symbols.  If you don't read the docs
>>>>> about bioconductor annotation and R packages pertaining thereto, it's
>>>>> hard to make progress.
>>>> Can you elaborate a little bit more why using gene symbols is
>>>> problematic?
>>>> And if it is really a problem, what gene ID I should use to avoid that
>>>> problem?
>>> Entrez gene IDs are generally more specific.  If you haven't run into
>>> collisions and ambiguities with symbols by now, maybe you don't have
>>> to worry about it.
>>>
>> Well, our reference database is UCSC refFlat downloaded from UCSC genome
>> browser, and in that format there is only geneName (SYMBOL) and name
>> (refSEQ). There is no Entrez Gene ID in that file. We actually have a
>> problem of one same geneName (SYMBOL) has more than one isoforms (more than
>> one refSEQ), but in fact we dont have a good way of differentiate between
>> those (by mathematical methods). So we use either gene symbol or refSeq ID.
>>
> most of our annotation packages have mappings that employ refseq identifiers.
> the functions i provided can be readily altered to use refseq ids as
> keys, and if you
> introduce types for the identifier/term tokens, the functions can be
> implemented as methods
> that perform appropriately for inputs from different vocabularies.
> but you will have to
> learn how to use the mapping objects or the sqlite tables
> productively.  Package GSEABase is
> also a useful resource for working with collections of gene identifiers.
>

Thanks for both of your suggestion Vincent. I will definitely check them 
out. Your suggestions/hints do help me a lot :).

Bests,

D.