[BioC] problem generating Chinese hamster txdb from GFF file

Marc Carlson mcarlson at fhcrc.org
Wed Aug 21 01:53:47 CEST 2013


Hi Guido,

Sorry for the delay.  Your reply seems to have snuck past my radar.

Below is the code I ended up using to process this file.  I looked more 
closely at the file you are trying to parse and I noticed that there is 
something very strange about it.  Many different features are labeled as 
being of type "exon", but only SOME of them actually have the exon 
ranking information.  :(

But what is really odd here is that the features which don't have an 
exon ranking don't seem to be the "same" as the ones that do. I would 
recommend that you get an explanation for that from the person who 
generated this file.

Anyhow I also made another quick update to GenomicFeatures to allow you 
to better make use of the other IDs that are in there. So with version 
1.13.33 you should be able to do this (for example), and get gene ids 
that may be closer to what you want.  I aim to do a little more in this 
department, but I have some other things I have to do 1st.

library(GenomicFeatures)

txdb <- makeTranscriptDbFromGFF(file="ref_CriGri_1.0_scaffolds.gff3",
                                 format="gff3",
                                 gffGeneIdAttributeName="Name",
dataSource="ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/GFF",
                                 species="Cricetulus griseus")


   Marc



On 08/16/2013 01:59 AM, Hooiveld, Guido wrote:
> Hi Marc,
> Thanks for helpful comments!
> Keys, keytypes, columns, select,.... still difficult for me...  Moreover, it is the first time I am working with these type of files, so I still have to get used to all concepts and do's & don'ts.
> Anyway, I tried to regenerate the txdb utilizing the info on exon ranks. Unfortunately, this doesn't work in my hands... Am I missing something obviously, or is something going seriously wrong? Note that the error thrown is slightly different between the production and dev versions of the packages (see below).
>
> I will also check the functions in rtracklayer.
>
> Regards,
> Guido
>
>
> Dev:
>> txdb <- makeTranscriptDbFromGFF(file="ref_CriGri_1.0_scaffolds.gff3", format="gff3", exonRankAttributeName="exon_number",
> + dataSource="ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/GFF", species="Cricetulus griseus")
> extracting transcript information
> Extracting gene IDs
> extracting transcript information
> Processing splicing information for gff3 file.
> Error in if (any(splicings$exon_rank == 0) && !any(splicings$exon_rank <  :
>    missing value where TRUE/FALSE needed
>> sessionInfo()
> R Under development (unstable) (2013-08-10 r63532)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C
> [5] LC_TIME=English_United States.1252
>
> attached base packages:
> [1] parallel  stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] GenomicFeatures_1.13.29 AnnotationDbi_1.23.18   Biobase_2.21.6          rtracklayer_1.21.9      GenomicRanges_1.13.36   XVector_0.1.0           IRanges_1.19.24
> [8] BiocGenerics_0.7.4      BiocInstaller_1.11.4
>
> loaded via a namespace (and not attached):
>   [1] biomaRt_2.17.2     Biostrings_2.29.15 bitops_1.0-5       BSgenome_1.29.1    DBI_0.2-7          RCurl_1.95-4.1     Rsamtools_1.13.29  RSQLite_0.11.4     stats4_3.1.0
> [10] tools_3.1.0        XML_3.98-1.1       zlibbioc_1.7.0
>
> Production:
>> txdb <- makeTranscriptDbFromGFF(file="ref_CriGri_1.0_scaffolds.gff3", format="gff3", exonRankAttributeName="exon_number",
> + dataSource="ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/GFF", species="Cricetulus griseus")
> extracting transcript information
> Extracting gene IDs
> extracting transcript information
> Processing splicing information for gff3 file.
> Prepare the 'metadata' data frame ... metadata: OK
> Now generating chrominfo from available sequence names. No chromosome length information is available.
> Error in .normargSplicings(splicings, transcripts_tx_id) :
>    'splicings$exon_rank' must be an integer vector with no NAs
> In addition: Warning message:
> In matchCircularity(chroms, circ_seqs) :
>    None of the strings in your circ_seqs argument match your seqnames.
>> sessionInfo()
> R version 3.0.1 Patched (2013-06-05 r62877)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C
> [5] LC_TIME=English_United States.1252
>
> attached base packages:
> [1] parallel  stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] TxDb.Hsapiens.UCSC.hg19.knownGene_2.9.2 GenomicFeatures_1.12.3                  AnnotationDbi_1.22.6                    Biobase_2.20.1
> [5] GenomicRanges_1.12.4                    IRanges_1.18.2                          BiocGenerics_0.6.0
>
> loaded via a namespace (and not attached):
>   [1] biomaRt_2.16.0     Biostrings_2.28.0  bitops_1.0-5       BSgenome_1.28.0    DBI_0.2-7          RCurl_1.95-4.1     Rsamtools_1.12.3   RSQLite_0.11.4
>   [9] rtracklayer_1.20.4 stats4_3.0.1       tools_3.0.1        XML_3.98-1.1       zlibbioc_1.6.0
>
> -----Original Message-----
> From: bioconductor-bounces at r-project.org [mailto:bioconductor-bounces at r-project.org] On Behalf Of Marc Carlson
> Sent: Thursday, August 15, 2013 22:44
> To: bioconductor at r-project.org
> Subject: Re: [BioC] problem generating Chinese hamster txdb from GFF file
>
> Hi Guido,
>
> The problem is that your keys are not valid.  The keys for GENEID can be retrieved with the keys method like this:
>
> k <- keys(txdb, keytype="GENEID")
>
> And if we look at them we can see that they don't look anything like the keys you are trying to extract data for:
>
> head(k)
> [1] "gene10"    "gene1000"  "gene10000" "gene10001" "gene10002" "gene10003"
>
>
> So if you use keys like that (valid keys), then everything works fine.
> For example:
>
> select(txdb, keys = head(k), columns="TXNAME", keytype="GENEID")
>
>        GENEID  TXNAME
> 1    gene10    rna8
> 2    gene10    rna7
> 3  gene1000  rna917
> 4 gene10000 rna9527
> 5 gene10001 rna9528
> 6 gene10002 rna9530
> 7 gene10002 rna9529
> 8 gene10003 rna9531
>
>
> Also I noticed that you are getting the warning about deducing the exon ranks.  Are you sure you want to do that with a mammal?  This file contains a field called "exon_number" that looks like something you should probably try to use.
>
> And finally, since you seem to want to use a real gene ID, you will have to extract those from the file separately.  This is because whoever created this file decided that they wanted to encode the actual gene symbols here in a peculiar way.  So you can do that step separately by using the rtracklayer import() method like this:
>
> library(rtracklayer)
> result <- import("ref_CriGri_1.0_scaffolds.gff3", format="gff3")
>
> If you look at what is in result, you will see (for example) that gene0 is cross referenced to GeneID:100754456 (for example).
>
>
> Hope this clarifies things,
>
>
>     Marc
>
>
>
>
> On 08/15/2013 10:20 AM, Hooiveld, Guido wrote:
>> Hello,
>> I would like to make a transcript database for Chinese hamster (Cricetulus griseus) using the function makeTranscriptDbFromGFF. Since it is the first time I use this function, I do this per the code in the vignette. However, I face some problems and would appreciate any pointer.
>>
>> I started by downloading the GFF (ref_CriGri_1.0_scaffolds.gff3) from NCBI.
>> [ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/GFF/]
>>
>> Creating the txdb seems to go fine:
>>
>>> txdb <- makeTranscriptDbFromGFF(file="ref_CriGri_1.0_scaffolds.gff3", format="gff3", dataSource="ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/GFF/", species="Cricetulus griseus")
>> extracting transcript information
>> Extracting gene IDs
>> extracting transcript information
>> Processing splicing information for gff3 file.
>> Deducing exon rank from relative coordinates provided
>> Prepare the 'metadata' data frame ... metadata: OK
>> Now generating chrominfo from available sequence names. No chromosome length information is available.
>> Warning messages:
>> 1: In .deduceExonRankings(exs, format = "gff") :
>>     Infering Exon Rankings.  If this is not what you expected, then please be sure that you have provided a valid attribute for exonRankAttributeName
>> 2: In matchCircularity(chroms, circ_seqs) :
>>     None of the strings in your circ_seqs argument match your seqnames.
>>
>> # Characterization txdb
>>> txdb
>> TranscriptDb object:
>> | Db type: TranscriptDb
>> | Supporting package: GenomicFeatures
>> | Data source: ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/GFF/
>> | Organism: Cricetulus griseus
>> | miRBase build ID: NA
>> | transcript_nrow: 21612
>> | exon_nrow: 212163
>> | cds_nrow: 183531
>> | Db created by: GenomicFeatures package from Bioconductor
>> | Creation time: 2013-08-15 18:26:33 +0200 (Thu, 15 Aug 2013)
>> | GenomicFeatures version at creation time: 1.13.29
>> | RSQLite version at creation time: 0.11.4
>> | DBSCHEMAVERSION: 1.0
>>> columns(txdb)
>> [1] "CDSID"      "CDSNAME"    "CDSCHROM"   "CDSSTRAND"  "CDSSTART"   "CDSEND"     "EXONID"
>>    [8] "EXONNAME"   "EXONCHROM"  "EXONSTRAND" "EXONSTART"  "EXONEND"    "GENEID"     "TXID"
>> [15] "EXONRANK"   "TXNAME"     "TXCHROM"    "TXSTRAND"   "TXSTART"    "TXEND"
>>> keytypes(txdb)
>> [1] "GENEID"   "TXID"     "TXNAME"   "EXONID"   "EXONNAME" "CDSID"    "CDSNAME"
>> As a test, I would like to retrieve the transcript names for 3 Entrez Gene IDs. Note: I selected these 3 IDs randomly; they are in the top5 -listed IDs in the GFF.
>> Unfortunately, this doesn't seem to work:
>>> keys <- c("100754456", "100759368", "100760238")
>>> select(txdb, keys = keys, columns="TXNAME", keytype="GENEID")
>> [1] GENEID TXNAME
>> <0 rows> (or 0-length row.names)
>> Note: the human example from the vignette worked fine.
>>
>> I think this is somehow due to the fact that during the parsing of the GFF the keytypes (?) are not properly recognized/linked. This is also reflected by the plain use of 'rna' and 'gene' as names/IDs:
>>> GR <- transcripts(txdb, columns=c("tx_id", "tx_name", "gene_id", "cds_id","cds_name"))
>>> head(GR)
>> GRanges with 6 ranges and 5 metadata columns:
>>               seqnames             ranges strand |     tx_id     tx_name         gene_id        cds_id
>>                  <Rle>          <IRanges>  <Rle> | <integer> <character> <CharacterList> <IntegerList>
>>     [1] NW_003613580.1 [ 736932,  738552]      + |         1        rna0           gene2           1,2
>>     [2] NW_003613580.1 [2073869, 2085513]      + |         2        rna4           gene6     3,4,5,...
>>     [3] NW_003613580.1 [2143936, 2241368]      + |         3        rna6           gene8    8,9,10,...
>>     [4] NW_003613580.1 [2390123, 2449470]      + |         4        rna8          gene10  22,23,24,...
>>     [5] NW_003613580.1 [2390123, 2555467]      + |         5        rna7          gene10  22,23,24,...
>>     [6] NW_003613580.1 [3070833, 3147709]      + |         6        rna9          gene12  43,44,45,...
>>                cds_name
>>         <CharacterList>
>>     [1]              NA
>>     [2]              NA
>>     [3]              NA
>>     [4]              NA
>>     [5]              NA
>>     [6]              NA
>>     ---
>>     seqlengths:
>>      NW_003613580.1 NW_003613581.1 NW_003613582.1 ... NW_003698075.1 NW_003704321.1 NW_003717424.1
>>                  NA             NA             NA ...             NA             NA             NA
>>> tail(GR)
>> GRanges with 6 ranges and 5 metadata columns:
>>               seqnames    ranges strand |     tx_id     tx_name         gene_id        cds_id
>>                  <Rle> <IRanges>  <Rle> | <integer> <character> <CharacterList> <IntegerList>
>>     [1] NW_003694533.1 [29, 205]      - |     21607    rna22504       gene25148        183527
>>     [2] NW_003697941.1 [ 4, 135]      - |     21608    rna22505       gene25149        183528
>>     [3] NW_003698068.1 [ 1, 267]      + |     21609    rna22506       gene25150        183529
>>     [4] NW_003698075.1 [ 1, 267]      - |     21610    rna22507       gene25151            NA
>>     [5] NW_003704321.1 [24, 203]      - |     21611    rna22508       gene25152        183530
>>     [6] NW_003717424.1 [10, 189]      + |     21612    rna22509       gene25155        183531
>>                cds_name
>>         <CharacterList>
>>     [1]              NA
>>     [2]              NA
>>     [3]              NA
>>     [4]              NA
>>     [5]              NA
>>     [6]              NA
>>     ---
>>     seqlengths:
>>      NW_003613580.1 NW_003613581.1 NW_003613582.1 ... NW_003698075.1 NW_003704321.1 NW_003717424.1
>>                  NA             NA             NA ...             NA             NA             NA
>>
>>
>> Any suggestions on how to properly generate a txdb from this GFF would be appreciated.
>>
>> Thanks,
>> Guido
>>
>>> sessionInfo()
>> R Under development (unstable) (2013-08-10 r63532)
>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>>
>> locale:
>> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252
>> [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
>> [5] LC_TIME=English_United States.1252
>>
>> attached base packages:
>> [1] parallel  stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] GenomicFeatures_1.13.29 AnnotationDbi_1.23.18   Biobase_2.21.6          GenomicRanges_1.13.36
>> [5] XVector_0.1.0           IRanges_1.19.21         BiocGenerics_0.7.3
>>
>> loaded via a namespace (and not attached):
>> [1] biomaRt_2.17.2     Biostrings_2.29.14 bitops_1.0-5       BSgenome_1.29.1    DBI_0.2-7
>>    [6] RCurl_1.95-4.1     Rsamtools_1.13.28  RSQLite_0.11.4     rtracklayer_1.21.9 stats4_3.1.0
>> [11] tools_3.1.0        XML_3.98-1.1       zlibbioc_1.7.0
>>
>> ---------------------------------------------------------
>> Guido Hooiveld, PhD
>> Nutrition, Metabolism & Genomics Group
>> Division of Human Nutrition
>> Wageningen University
>> Biotechnion, Bomenweg 2
>> NL-6703 HD Wageningen
>> the Netherlands
>> tel: (+)31 317 485788
>> fax: (+)31 317 483342
>> email:      guido.hooiveld at wur.nl
>> internet:   http://nutrigene.4t.com
>> http://scholar.google.com/citations?user=qFHaMnoAAAAJ
>> http://www.researcherid.com/rid/F-4912-2010
>>
>>
>> 	[[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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
> _______________________________________________
> 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
>
>



More information about the Bioconductor mailing list