# Contents

# 1 Basic usage

Simply put, the `genbankr` package parses files in the NCBI’s GenBank (gb/gbk) format into R.

The primary workhorse is provided by the `readGenBank` function, which accepts either a GenBank file (via the `file` argument) or raw text in the GenBank format (via the `text` argument).

    suppressPackageStartupMessages(library(genbankr))
    smpfile = system.file("sample.gbk", package="genbankr")
    gb = readGenBank(smpfile)

    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column
    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column

    ## No exons read from genbank file. Assuming sections of CDS are full exons

    ## No transcript features (mRNA) found, using spans of CDSs

    ## Warning in fill_stack_df(feats[!typs %in% c("gene", "exon", "CDS",
    ## "variation", : Got unexpected multi-value field(s) [ rpt_type ]. The resulting
    ## column(s) will be of class CharacterList, rather than vector(s). Please contact
    ## the maintainer if multi-valuedness is expected/meaningful for the listed
    ## field(s).

    gb

    ## GenBank Annotations
    ## Human herpesvirus 5 strain VR1814, complete genome.
    ## Accession: GU179289
    ## 1 Sequence(s) with total length length: 235233
    ## 174 genes
    ## 170 transcripts
    ## 191 exons/cds elements
    ## 61 variations
    ## 24 other features

`readGenBank` generates a `GenBankFull` object, which contains the annotations and the origin sequence. genbankr provides methods for the AnnotationDbi API functions for retrieving information from the object.

We can retreive the genes via the function of the same name:

    genes(gb)

    ## GRanges object with 174 ranges and 5 metadata columns:
    ##         seqnames        ranges strand |        type        gene     loctype
    ##            <Rle>     <IRanges>  <Rle> | <character> <character> <character>
    ##     [1]   VR1814     1232-2164      + |        gene         RL1      normal
    ##     [2]   VR1814     2353-4863      - |        gene      RNA2.7      normal
    ##     [3]   VR1814     5252-5539      - |        gene        RL5A      normal
    ##     [4]   VR1814     5888-6223      - |        gene         RL6      normal
    ##     [5]   VR1814     6622-7677      - |        gene      RNA1.2      normal
    ##     ...      ...           ...    ... .         ...         ...         ...
    ##   [170]   VR1814 229242-229793      + |        gene        US32      normal
    ##   [171]   VR1814 229995-230168      + |        gene       US33A      normal
    ##   [172]   VR1814 230337-230828      + |        gene        US34      normal
    ##   [173]   VR1814 230822-231016      + |        gene       US34A      normal
    ##   [174]   VR1814 231986-234349      - |        gene        TRS1      normal
    ##                           note     gene_id
    ##                    <character> <character>
    ##     [1]                   <NA>         RL1
    ##     [2]                   <NA>      RNA2.7
    ##     [3]                   <NA>        RL5A
    ##     [4]                   <NA>         RL6
    ##     [5]                   <NA>      RNA1.2
    ##     ...                    ...         ...
    ##   [170]                   <NA>        US32
    ##   [171]                   <NA>       US33A
    ##   [172]                   <NA>        US34
    ##   [173]                   <NA>       US34A
    ##   [174] putative beta gene; ..        TRS1
    ##   -------
    ##   seqinfo: 1 sequence from GU179289.1 genome

We can also do the same for exons, cds elements, and transcripts (code run but output ommitted for brevity):

    cds(gb)
    exons(gb)
    transcripts(gb)

We can also access elements not in the standard TxDb API, such as variants and features which don’t fit into any of the other categories:

    variants(gb)
    otherFeatures(gb)

Furthermore, we can access header-level information via accessors as well:

    accession(gb)

    ## [1] "GU179289"

    vers(gb)

    ## accession.version         GenInfoID
    ##      "GU179289.1"       "270355759"

We can get the seqinfo for a `GenBankAnnot`/`GenBankFull` object:

    seqinfo(gb)

    ## Seqinfo object with 1 sequence from GU179289.1 genome:
    ##   seqnames seqlengths isCircular     genome
    ##   VR1814       235233      FALSE GU179289.1

Finally, we can ge the sequence itself from a `GenBankFull` object:

    getSeq(gb)

    ## DNAStringSet object of length 1:
    ##      width seq                                              names
    ## [1] 235233 CCATTCCGGGCCGCGTGCTGGGT...GCCGCCAGTGCGGGACAGGGCT VR1814

# 2 Low level parsing

While use and integration of the Bioconductor machinery is recommended, we also provide low-level parsing capabilities via the workhorse `parseGenBank` function. This function returns a list structure roughly corresponding to the top-level headings within the genbank format itself:

    pg = parseGenBank(smpfile)
    str(pg, max.level = 1)

    ## List of 10
    ##  $ LOCUS     : chr "LOCUS       GU179289              235233 bp    DNA     linear   VRL 09-MAY-2013"
    ##  $ FEATURES  :List of 429
    ##   ..- attr(*, "dim")= int 429
    ##   ..- attr(*, "dimnames")=List of 1
    ##  $ ORIGIN    :Formal class 'DNAStringSet' [package "Biostrings"] with 5 slots
    ##  $ ACCESSION : chr "GU179289"
    ##  $ COMMENT   : NULL
    ##  $ DEFINITION: chr "Human herpesvirus 5 strain VR1814, complete genome."
    ##  $ KEYWORDS  : chr NA
    ##  $ REFERENCE : NULL
    ##  $ SOURCE    :List of 3
    ##  $ VERSION   : Named chr [1:2] "GU179289.1" "270355759"
    ##   ..- attr(*, "names")= chr [1:2] "accession.version" "GenInfoID"

# 3 retaining the genomic sequence

If desired, `readGenBank` and `parseGenBank` can omit full sequence of the organism by specifying `ret.seq=FALSE`. In this case, `readGenBank` returns a `GenBankAnnot` object, rather than a `GenBankFull`

    gbf = readGenBank(smpfile, ret.seq = FALSE)

    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column
    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column

    ## No exons read from genbank file. Assuming sections of CDS are full exons

    ## No transcript features (mRNA) found, using spans of CDSs

    ## Warning in fill_stack_df(feats[!typs %in% c("gene", "exon", "CDS",
    ## "variation", : Got unexpected multi-value field(s) [ rpt_type ]. The resulting
    ## column(s) will be of class CharacterList, rather than vector(s). Please contact
    ## the maintainer if multi-valuedness is expected/meaningful for the listed
    ## field(s).

    gbf

    ## GenBank Annotations
    ## Human herpesvirus 5 strain VR1814, complete genome.
    ## Accession: GU179289
    ## 1 Sequence(s) with total length length: 235233
    ## 174 genes
    ## 170 transcripts
    ## 191 exons/cds elements
    ## 61 variations
    ## 24 other features

All of the accessor methods discussed above work for `GenBankFull` objects as they do for `GenBankAnnot` objects.

# 4 rtracklayer style import

We also provide a convenience method for using rtracklayer `import` style mechanics for reading GenBank files:

    gbkfile = GenBankFile(smpfile)

    ## Warning: This class is extending the deprecated RTLFile class from rtracklayer.
    ## Use BiocFile from BiocIO in place of RTLFile.

    gb2 = import(gbkfile)

    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column
    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column

    ## No exons read from genbank file. Assuming sections of CDS are full exons

    ## No transcript features (mRNA) found, using spans of CDSs

    ## Warning in fill_stack_df(feats[!typs %in% c("gene", "exon", "CDS",
    ## "variation", : Got unexpected multi-value field(s) [ rpt_type ]. The resulting
    ## column(s) will be of class CharacterList, rather than vector(s). Please contact
    ## the maintainer if multi-valuedness is expected/meaningful for the listed
    ## field(s).

# 5 Retrieving and parsing GenBank information by Versioned Accession

genbankr provides the `GBAccession` class and constructor for representing versioned Nuccore accession numbers.

    gba = GBAccession("U49845.1")
    gba

    ## GenBank Accession Number(s):  U49845.1

These accession objects can be passed directly to readGenBank:

    readGenBank(gba, partial=TRUE)

    ## Loading required namespace: rentrez

    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column
    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column

    ## genes not available for all CDS ranges, using internal grouping ids

    ## No exons read from genbank file. Assuming sections of CDS are full exons

    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column
    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column

    ## GenBank Annotations
    ## Saccharomyces cerevisiae TCP1-beta gene, partial cds; and Axl2p (AXL2) and Rev7p (REV7) genes, complete cds.
    ## Accession: U49845
    ## 1 Sequence(s) with total length length: 5028
    ## 2 genes
    ## 3 transcripts
    ## 3 exons/cds elements
    ## 0 variations
    ## 0 other features

# 6 Retrieving only origin sequence

`genbankr` also provides a fastpath for extracting only the sequence of an organism. We can call `getSeq` on a `GenBankFile` or `GBAccession` object

    getSeq(gbkfile)

    ## DNAStringSet object of length 1:
    ##      width seq                                              names
    ## [1] 235233 CCATTCCGGGCCGCGTGCTGGGT...GCCGCCAGTGCGGGACAGGGCT VR1814

Additionally, we can specify `ret.anno = FALSE` in `parseGenBank`

    parseGenBank(smpfile, ret.anno=FALSE)

    ## DNAStringSet object of length 1:
    ##      width seq                                              names
    ## [1] 235233 CCATTCCGGGCCGCGTGCTGGGT...GCCGCCAGTGCGGGACAGGGCT VR1814

# 7 Creating TxDb objects from genbank annotations

`genbankr` provides the `makeTxDbFromGenBank` function, which accepts a `GenBankRecord` or `GBAccession` object and returns a `TxDb` of the annotations.

    gbr = readGenBank(smpfile)

    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column
    ## Annotations don't have 'locus_tag' label, using 'gene' as gene_id column

    ## No exons read from genbank file. Assuming sections of CDS are full exons

    ## No transcript features (mRNA) found, using spans of CDSs

    ## Warning in fill_stack_df(feats[!typs %in% c("gene", "exon", "CDS",
    ## "variation", : Got unexpected multi-value field(s) [ rpt_type ]. The resulting
    ## column(s) will be of class CharacterList, rather than vector(s). Please contact
    ## the maintainer if multi-valuedness is expected/meaningful for the listed
    ## field(s).

    tx = makeTxDbFromGenBank(gbr)
    tx

    ## TxDb object:
    ## # Db type: TxDb
    ## # Supporting package: GenomicFeatures
    ## # accession.version: GU179289.1
    ## # GenInfoID: 270355759
    ## # genome: GU179289.1
    ## # source: NCBI (GenBank)
    ## # definition: Human herpesvirus 5 strain VR1814, complete genome.
    ## # Db content generated by: genbankr package from Bioconductor
    ## # genbankr version at creation time: 1.22.0
    ## # Nb of transcripts: 170
    ## # Db created by: GenomicFeatures package from Bioconductor
    ## # Creation time: 2021-10-26 21:36:24 -0400 (Tue, 26 Oct 2021)
    ## # GenomicFeatures version at creation time: 1.46.0
    ## # RSQLite version at creation time: 2.2.8
    ## # DBSCHEMAVERSION: 1.2

# 8 Details and caveats

Often times, GenBank files don’t contain exhaustive annotations. For example, files including CDS annotations often do not have separate transcript features. Furthermore, chromosomes are not always named, particularly in organisms that have only one. The details of how genbankr handles such cases are as follows:

In files where CDSs are annotated but individual exons are not, ‘approximate exons’ are defined as the individual contiguous elements within each CDS. Currently, no mixing of approximate and explicitly annotated exons is performed, even in cases where, e.g., exons are not annotated for some genes with CDS annotations.

In files where transcripts are not present, ‘approximate transcripts’ defined by the ranges spanned by groups of exons are used. Currently, we do not support generating approximate transcripts from CDSs in files that contain actual transcript annotations, even if those annotations do not cover all genes with CDS/exon annotations.

Features (gene, cds, variant, etc) are assumed to be contained within the most recent previous source feature (chromosome/physical piece of DNA). Chromosome name for source features (seqnames in the resulting GRanges}/VRanges is determined as follows:

1. The ‘chromosome’ attribute, as is (e.g., “chr1”);
2. the ‘strain’ attribute, combined with auto-generated count (e.g., “VR1814:1”);
3. the ‘organism’ attribute, combined with auto-generated count (e.g.“Human herpesvirus 5:1”)

Some GenBank files do not include origin sequence. In these cases, variation features are not supported, as there is no self-contained way to determine reference sequence and the features themselves typically contain only alt information (if that). In the case of files containing variation features but no origin sequence, those features are ignored with a warning.

Currently some information about from the header of a GenBank file, primarily reference and author based information, is not captured and returned. Please contact the maintainer if you have a direct use-case for this type of information.

# 9 Performance

We have taken pains to make the genbankr parser as effcient as easily possible. On our local machines, a 19MB genbank file takes 2-3 minutes to be parsed. That said, this package is not tested and likely is not suitable for processing extremely large genbank files. We suggest obtaining the annotations in a different format in such cases.