Contents

Original Authors: Martin Morgan, Sonali Arora
Presenting Authors: Martin Morgan, Lori Shepherd
Date: 22 July, 2019
Back: Monday labs

Objective: Learn about Bioconductor resources for gene and genome annotation.

Lessons learned:

1 Gene annotation

1.1 Data packages

Organism-level (‘org’) packages contain mappings between a central identifier (e.g., Entrez gene ids) and other identifiers (e.g. GenBank or Uniprot accession number, RefSeq id, etc.). The name of an org package is always of the form org.<Sp>.<id>.db (e.g. org.Sc.sgd.db) where <Sp> is a 2-letter abbreviation of the organism (e.g. Sc for Saccharomyces cerevisiae) and <id> is an abbreviation (in lower-case) describing the type of central identifier (e.g. sgd for gene identifiers assigned by the Saccharomyces Genome Database, or eg for Entrez gene ids). The “How to use the ‘.db’ annotation packages” vignette in the AnnotationDbi package (org packages are only one type of “.db” annotation packages) is a key reference. The ‘.db’ and most other Bioconductor annotation packages are updated every 6 months.

Annotation packages usually contain an object named after the package itself. These objects are collectively called AnnotationDb objects, with more specific classes named OrgDb, ChipDb or TranscriptDb objects. Methods that can be applied to these objects include cols(), keys(), keytypes() and select(). Common operations for retrieving annotations are summarized in the table.

Category Function Description
Discover columns() List the kinds of columns that can be returned
keytypes() List columns that can be used as keys
keys() List values that can be expected for a given keytype
select() Retrieve annotations matching keys, keytype and columns
Manipulate setdiff(), union(), intersect() Operations on sets
duplicated(), unique() Mark or remove duplicates
%in%, match() Find matches
any(), all() Are any TRUE? Are all?
merge() Combine two different based on shared keys
GRanges* transcripts(), exons(), cds() Features (transcripts, exons, coding sequence) as GRanges.
transcriptsBy() , exonsBy() Features group by gene, transcript, etc., as GRangesList.
cdsBy()

1.2 Internet resources

A short summary of select Bioconductor packages enabling web-based queries is in following Table.

Package Description
AnnotationHub Ensembl, Encode, dbSNP, UCSC data objects
biomaRt Ensembl and other annotations
PSICQUIC Protein interactions
uniprot.ws Protein annotations
KEGGREST KEGG pathways
SRAdb Sequencing experiments.
rtracklayer genome tracks.
GEOquery Array and other data
ArrayExpress Array and other data

1.3 Exercises

Exercise 1: This exercise illustrates basic use of the `select’ interface to annotation packages.

  1. Install and attach the org.Hs.eg.db annotation package; it contains ‘symbol mapping’ information for Homo sapiens, based on NCBI ‘Entrez’ identifiers.

    library(org.Hs.eg.db)
  2. Take a quick look at a summary of data in this package

    org.Hs.eg.db
    ## OrgDb object:
    ## | DBSCHEMAVERSION: 2.1
    ## | Db type: OrgDb
    ## | Supporting package: AnnotationDbi
    ## | DBSCHEMA: HUMAN_DB
    ## | ORGANISM: Homo sapiens
    ## | SPECIES: Human
    ## | EGSOURCEDATE: 2019-Apr26
    ## | EGSOURCENAME: Entrez Gene
    ## | EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
    ## | CENTRALID: EG
    ## | TAXID: 9606
    ## | GOSOURCENAME: Gene Ontology
    ## | GOSOURCEURL: ftp://ftp.geneontology.org/pub/go/godatabase/archive/latest-lite/
    ## | GOSOURCEDATE: 2019-Apr24
    ## | GOEGSOURCEDATE: 2019-Apr26
    ## | GOEGSOURCENAME: Entrez Gene
    ## | GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
    ## | KEGGSOURCENAME: KEGG GENOME
    ## | KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
    ## | KEGGSOURCEDATE: 2011-Mar15
    ## | GPSOURCENAME: UCSC Genome Bioinformatics (Homo sapiens)
    ## | GPSOURCEURL: 
    ## | GPSOURCEDATE: 2018-Dec3
    ## | ENSOURCEDATE: 2019-Apr08
    ## | ENSOURCENAME: Ensembl
    ## | ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
    ## | UPSOURCENAME: Uniprot
    ## | UPSOURCEURL: http://www.UniProt.org/
    ## | UPSOURCEDATE: Fri Apr 26 20:12:58 2019
    ## 
    ## Please see: help('select') for usage information
  3. The idea is that there are keytypes() that can be mapped to different columns(); keys() can be used to see available keys. Explore the package to see what sorts of information is available, e.g.,

    keytypes(org.Hs.eg.db)
    ##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
    ##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
    ##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
    ## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
    ## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
    ## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
    ## [25] "UNIGENE"      "UNIPROT"
    columns(org.Hs.eg.db)
    ##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
    ##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
    ##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
    ## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
    ## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
    ## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
    ## [25] "UNIGENE"      "UNIPROT"
    head(keys(org.Hs.eg.db, "SYMBOL"))
    ## [1] "A1BG"  "A2M"   "A2MP1" "NAT1"  "NAT2"  "NATP"
  4. There are two basic ways of extracting data from an org.* package – mapIds() to create a 1:1 mapping between key and a single column, and select() (it’s often necessary to specify this function directly, to avoid a conflict with dplyr, as AnnotationDbi::select()). Explore these functions, e.g.,

    set.seed(123)
    egid <- sample(keys(org.Hs.eg.db), 6)
    mapIds(org.Hs.eg.db, egid, "SYMBOL", "ENTREZID")
    ## 'select()' returned 1:1 mapping between keys and columns
    ##      106480363      110120927           3654      110485085      100129570 
    ##    "RNU7-186P" "LOC110120927"        "IRAK1" "LOC110485085"      "PEBP1P3" 
    ##      100129013 
    ##      "ACTR3P2"
    AnnotationDbi::select(
        org.Hs.eg.db, egid, c("SYMBOL", "ENSEMBL", "GENENAME"), "ENTREZID"
    )
    ## 'select()' returned 1:1 mapping between keys and columns
    ##    ENTREZID       SYMBOL         ENSEMBL
    ## 1 106480363    RNU7-186P            <NA>
    ## 2 110120927 LOC110120927            <NA>
    ## 3      3654        IRAK1 ENSG00000184216
    ## 4 110485085 LOC110485085            <NA>
    ## 5 100129570      PEBP1P3            <NA>
    ## 6 100129013      ACTR3P2            <NA>
    ##                                                  GENENAME
    ## 1                    RNA, U7 small nuclear 186 pseudogene
    ## 2                                    VISTA enhancer hs721
    ## 3              interleukin 1 receptor associated kinase 1
    ## 4                 BRCA1P1 intergenic recombination region
    ## 5 phosphatidylethanolamine binding protein 1 pseudogene 3
    ## 6                                      ACTR3 pseudogene 2
  5. Some key - column mappings are 1:many, e.g., Entrez ID "3812" maps to 44 Ensembl Ids. What does mapIds() return when mapping Entrez ID "3812" to Ensembl ids? Use the additional argument multiVals = "CharacterList" to explore further. Compare results to those returned by select().

    egid <- "3812"
    mapIds(org.Hs.eg.db, egid, "ENSEMBL", "ENTREZID")
    ## 'select()' returned 1:many mapping between keys and columns
    ##              3812 
    ## "ENSG00000240403"
    mapIds(
        org.Hs.eg.db, egid, "ENSEMBL", "ENTREZID",
        multiVals = "CharacterList"
    )
    ## 'select()' returned 1:many mapping between keys and columns
    ## CharacterList of length 1
    ## [["3812"]] ENSG00000240403 ENSG00000278850 ... ENSG00000284381
    AnnotationDbi::select(
        org.Hs.eg.db, egid, c("SYMBOL", "ENSEMBL"),
        multiVals = "CharacterList"
    )
    ## 'select()' returned 1:many mapping between keys and columns
    ##    ENTREZID  SYMBOL         ENSEMBL
    ## 1      3812 KIR3DL2 ENSG00000240403
    ## 2      3812 KIR3DL2 ENSG00000278850
    ## 3      3812 KIR3DL2 ENSG00000276424
    ## 4      3812 KIR3DL2 ENSG00000275566
    ## 5      3812 KIR3DL2 ENSG00000276004
    ## 6      3812 KIR3DL2 ENSG00000278809
    ## 7      3812 KIR3DL2 ENSG00000278726
    ## 8      3812 KIR3DL2 ENSG00000278403
    ## 9      3812 KIR3DL2 ENSG00000275629
    ## 10     3812 KIR3DL2 ENSG00000276882
    ## 11     3812 KIR3DL2 ENSG00000278710
    ## 12     3812 KIR3DL2 ENSG00000275083
    ## 13     3812 KIR3DL2 ENSG00000277181
    ## 14     3812 KIR3DL2 ENSG00000278656
    ## 15     3812 KIR3DL2 ENSG00000274722
    ## 16     3812 KIR3DL2 ENSG00000275416
    ## 17     3812 KIR3DL2 ENSG00000278474
    ## 18     3812 KIR3DL2 ENSG00000278707
    ## 19     3812 KIR3DL2 ENSG00000275838
    ## 20     3812 KIR3DL2 ENSG00000276739
    ## 21     3812 KIR3DL2 ENSG00000277709
    ## 22     3812 KIR3DL2 ENSG00000273911
    ## 23     3812 KIR3DL2 ENSG00000278361
    ## 24     3812 KIR3DL2 ENSG00000277982
    ## 25     3812 KIR3DL2 ENSG00000275626
    ## 26     3812 KIR3DL2 ENSG00000276357
    ## 27     3812 KIR3DL2 ENSG00000275511
    ## 28     3812 KIR3DL2 ENSG00000273735
    ## 29     3812 KIR3DL2 ENSG00000275262
    ## 30     3812 KIR3DL2 ENSG00000278442
    ## 31     3812 KIR3DL2 ENSG00000278758
    ## 32     3812 KIR3DL2 ENSG00000284384
    ## 33     3812 KIR3DL2 ENSG00000283975
    ## 34     3812 KIR3DL2 ENSG00000284213
    ## 35     3812 KIR3DL2 ENSG00000284295
    ## 36     3812 KIR3DL2 ENSG00000284063
    ## 37     3812 KIR3DL2 ENSG00000284466
    ## 38     3812 KIR3DL2 ENSG00000284101
    ## 39     3812 KIR3DL2 ENSG00000284046
    ## 40     3812 KIR3DL2 ENSG00000284053
    ## 41     3812 KIR3DL2 ENSG00000283951
    ## 42     3812 KIR3DL2 ENSG00000284192
    ## 43     3812 KIR3DL2 ENSG00000284528
    ## 44     3812 KIR3DL2 ENSG00000284381
  6. It seems like it might often be useful to use the tidyverse on return values from mapIds() and select(); explore this usage

    library(tidyverse)
    egid <- keys(org.Hs.eg.db)    # all ENTREZIDs
    mapIds(org.Hs.eg.db, egid, "SYMBOL", "ENTREZID") %>% 
       as_tibble() %>% 
       rownames_to_column("ENTREZID")
    ## Warning: Calling `as_tibble()` on a vector is discouraged, because the behavior is likely to change in the future. Use `tibble::enframe(name = NULL)` instead.
    ## This warning is displayed once per session.
    ## # A tibble: 61,521 x 2
    ##    ENTREZID value   
    ##    <chr>    <chr>   
    ##  1 1        A1BG    
    ##  2 2        A2M     
    ##  3 3        A2MP1   
    ##  4 4        NAT1    
    ##  5 5        NAT2    
    ##  6 6        NATP    
    ##  7 7        SERPINA3
    ##  8 8        AADAC   
    ##  9 9        AAMP    
    ## 10 10       AANAT   
    ## # … with 61,511 more rows
    AnnotationDbi::select(
        org.Hs.eg.db, egid, c("SYMBOL", "GO", "GENENAME"), "ENTREZID"
    ) %>% as_tibble()
    ## # A tibble: 309,097 x 6
    ##    ENTREZID SYMBOL GO         EVIDENCE ONTOLOGY GENENAME              
    ##    <chr>    <chr>  <chr>      <chr>    <chr>    <chr>                 
    ##  1 1        A1BG   GO:0002576 TAS      BP       alpha-1-B glycoprotein
    ##  2 1        A1BG   GO:0003674 ND       MF       alpha-1-B glycoprotein
    ##  3 1        A1BG   GO:0005576 HDA      CC       alpha-1-B glycoprotein
    ##  4 1        A1BG   GO:0005576 IDA      CC       alpha-1-B glycoprotein
    ##  5 1        A1BG   GO:0005576 TAS      CC       alpha-1-B glycoprotein
    ##  6 1        A1BG   GO:0005615 HDA      CC       alpha-1-B glycoprotein
    ##  7 1        A1BG   GO:0008150 ND       BP       alpha-1-B glycoprotein
    ##  8 1        A1BG   GO:0031093 TAS      CC       alpha-1-B glycoprotein
    ##  9 1        A1BG   GO:0034774 TAS      CC       alpha-1-B glycoprotein
    ## 10 1        A1BG   GO:0043312 TAS      BP       alpha-1-B glycoprotein
    ## # … with 309,087 more rows

Exercise 2: biomaRt.

Internet access required for this exercise

  1. Explore the Biomart web site https://www.ensembl.org/biomart for retrieving all kinds of genomic annotations.

Start by choosing a database (e.g., ‘Ensembl Genes 92’), dataset (e.g., ‘Human genes (GRCh38.p12)’), filter (e.g., ‘GENE’ / ‘Input external reference’ / ‘Gene stable id’ and enter ‘ENSG00000000003’), attributes (default is ok), then press ‘Results’ to map from Ensembl identifier to transcript identifier.

  1. Install (if necessary) and load the biomaRt package. Use listMarts() to see availble databases, useMart() to select the mart you’re interested in.

    library(biomaRt)
    head(listMarts())
    ##                biomart               version
    ## 1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 97
    ## 2   ENSEMBL_MART_MOUSE      Mouse strains 97
    ## 3     ENSEMBL_MART_SNP  Ensembl Variation 97
    ## 4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 97
    mart <- useMart("ENSEMBL_MART_ENSEMBL")
  2. Use listDatasets() and useDataset() to select the Homo sapiens gene dataset.

    head(listDatasets(mart))
    ##                        dataset                           description
    ## 1 abrachyrhynchus_gene_ensembl Pink-footed goose genes (ASM259213v1)
    ## 2     acalliptera_gene_ensembl      Eastern happy genes (fAstCal1.2)
    ## 3   acarolinensis_gene_ensembl        Anole lizard genes (AnoCar2.0)
    ## 4    acitrinellus_gene_ensembl        Midas cichlid genes (Midas_v5)
    ## 5        ahaastii_gene_ensembl    Great spotted kiwi genes (aptHaa1)
    ## 6    amelanoleuca_gene_ensembl                 Panda genes (ailMel1)
    ##       version
    ## 1 ASM259213v1
    ## 2  fAstCal1.2
    ## 3   AnoCar2.0
    ## 4    Midas_v5
    ## 5     aptHaa1
    ## 6     ailMel1
    dataset <- useDataset("hsapiens_gene_ensembl", mart)
  3. Use listFilters() to see available filters. The filter is the type of data that you are querying with. Choose one.

    head(listFilters(dataset))
    ##              name              description
    ## 1 chromosome_name Chromosome/scaffold name
    ## 2           start                    Start
    ## 3             end                      End
    ## 4      band_start               Band Start
    ## 5        band_end                 Band End
    ## 6    marker_start             Marker Start
    filters <- "ensembl_gene_id"                    # see `listFilters()`
  4. Use listAttrbutes() to see available attributes. Attributes represent the information you’d like to retrieve. Choose some!

    head(listAttributes(dataset))
    ##                            name                  description         page
    ## 1               ensembl_gene_id               Gene stable ID feature_page
    ## 2       ensembl_gene_id_version       Gene stable ID version feature_page
    ## 3         ensembl_transcript_id         Transcript stable ID feature_page
    ## 4 ensembl_transcript_id_version Transcript stable ID version feature_page
    ## 5            ensembl_peptide_id            Protein stable ID feature_page
    ## 6    ensembl_peptide_id_version    Protein stable ID version feature_page
    attrs <- c("ensembl_gene_id", "hgnc_symbol")    # see `listAttributes()`
  5. Create a character vector of Ensembl gene ids, compose and execute the query, transforming the result to a tibble.

    ids <- c(
        "ENSG00000000003", "ENSG00000000005", "ENSG00000000419", 
        "ENSG00000000457", "ENSG00000000460", "ENSG00000000938"
    )
    tbl <- getBM(attrs, filters, ids, dataset) %>% as_tibble()
    ## Cache found
    tbl
    ## # A tibble: 6 x 2
    ##   ensembl_gene_id hgnc_symbol
    ##   <chr>           <chr>      
    ## 1 ENSG00000000003 TSPAN6     
    ## 2 ENSG00000000005 TNMD       
    ## 3 ENSG00000000419 DPM1       
    ## 4 ENSG00000000457 SCYL3      
    ## 5 ENSG00000000460 C1orf112   
    ## 6 ENSG00000000938 FGR

Exercise 3: KEGGREST

Internet access required for this exercise

  1. Explore the KEGG web site https://www.genome.jp/kegg/ KEGG is a database of information on pathways.

  2. Load the KEGGREST package and discover available databases

    library(KEGGREST)
    KEGGREST::listDatabases()
    ##  [1] "pathway"  "brite"    "module"   "ko"       "genome"   "vg"      
    ##  [7] "ag"       "compound" "glycan"   "reaction" "rclass"   "enzyme"  
    ## [13] "disease"  "drug"     "dgroup"   "environ"  "genes"    "ligand"  
    ## [19] "kegg"
  3. Use keggList() to query the pathway database for human pathways; present the result as a tibble

    hsa_pathways <- keggList("pathway", "hsa") %>% 
        tibble(pathway = names(.), description = .)
    hsa_pathways
    ## # A tibble: 333 x 2
    ##    pathway       description                                                 
    ##    <chr>         <chr>                                                       
    ##  1 path:hsa00010 Glycolysis / Gluconeogenesis - Homo sapiens (human)         
    ##  2 path:hsa00020 Citrate cycle (TCA cycle) - Homo sapiens (human)            
    ##  3 path:hsa00030 Pentose phosphate pathway - Homo sapiens (human)            
    ##  4 path:hsa00040 Pentose and glucuronate interconversions - Homo sapiens (hu…
    ##  5 path:hsa00051 Fructose and mannose metabolism - Homo sapiens (human)      
    ##  6 path:hsa00052 Galactose metabolism - Homo sapiens (human)                 
    ##  7 path:hsa00053 Ascorbate and aldarate metabolism - Homo sapiens (human)    
    ##  8 path:hsa00061 Fatty acid biosynthesis - Homo sapiens (human)              
    ##  9 path:hsa00062 Fatty acid elongation - Homo sapiens (human)                
    ## 10 path:hsa00071 Fatty acid degradation - Homo sapiens (human)               
    ## # … with 323 more rows
  4. Use keggLink() to recover the genes in each pathway.

    hsa_path_eg  <- keggLink("pathway", "hsa") %>% 
        tibble(pathway = ., egid = sub("hsa:", "", names(.)))
    hsa_path_eg
    ## # A tibble: 30,967 x 2
    ##    pathway       egid  
    ##    <chr>         <chr> 
    ##  1 path:hsa00010 10327 
    ##  2 path:hsa00010 124   
    ##  3 path:hsa00010 125   
    ##  4 path:hsa00010 126   
    ##  5 path:hsa00010 127   
    ##  6 path:hsa00010 128   
    ##  7 path:hsa00010 130   
    ##  8 path:hsa00010 130589
    ##  9 path:hsa00010 131   
    ## 10 path:hsa00010 160287
    ## # … with 30,957 more rows
    hsa_path_eg %>% group_by(pathway) %>% summarize(genes = list(egid))
    ## # A tibble: 333 x 2
    ##    pathway       genes     
    ##    <chr>         <list>    
    ##  1 path:hsa00010 <chr [68]>
    ##  2 path:hsa00020 <chr [30]>
    ##  3 path:hsa00030 <chr [30]>
    ##  4 path:hsa00040 <chr [34]>
    ##  5 path:hsa00051 <chr [33]>
    ##  6 path:hsa00052 <chr [31]>
    ##  7 path:hsa00053 <chr [27]>
    ##  8 path:hsa00061 <chr [18]>
    ##  9 path:hsa00062 <chr [27]>
    ## 10 path:hsa00071 <chr [44]>
    ## # … with 323 more rows
  5. Update the hsa_path_eg table to include information on gene symbol and Ensembl id from the org.Hs.eg.db package. Retrieve the relevant information using mapIds(). How would you deal with entrez gene ids that map to multiple Ensembl ids?

    hsa_kegg_anno <- hsa_path_eg %>%
        mutate(
            symbol = mapIds(org.Hs.eg.db, egid, "SYMBOL", "ENTREZID"),
            ensembl = mapIds(org.Hs.eg.db, egid, "ENSEMBL", "ENTREZID")
        )
    ## 'select()' returned 1:1 mapping between keys and columns
    ## 'select()' returned 1:many mapping between keys and columns
  6. Use left_join() to append pathway descriptions to the hsa_kegg_anno table.

    left_join(hsa_kegg_anno, hsa_pathways)
    ## Joining, by = "pathway"
    ## # A tibble: 30,967 x 5
    ##    pathway     egid   symbol  ensembl     description                        
    ##    <chr>       <chr>  <chr>   <chr>       <chr>                              
    ##  1 path:hsa00… 10327  AKR1A1  ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  2 path:hsa00… 124    ADH1A   ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  3 path:hsa00… 125    ADH1B   ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  4 path:hsa00… 126    ADH1C   ENSG000002… Glycolysis / Gluconeogenesis - Hom…
    ##  5 path:hsa00… 127    ADH4    ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  6 path:hsa00… 128    ADH5    ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  7 path:hsa00… 130    ADH6    ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  8 path:hsa00… 130589 GALM    ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ##  9 path:hsa00… 131    ADH7    ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ## 10 path:hsa00… 160287 LDHAL6A ENSG000001… Glycolysis / Gluconeogenesis - Hom…
    ## # … with 30,957 more rows

2 Genome annotation

There are a diversity of packages and classes available for representing large genomes. Several include:

2.1 Transcript annotation packages

Genome-centric packages are very useful for annotations involving genomic coordinates. It is straight-forward, for instance, to discover the coordinates of coding sequences in regions of interest, and from these retrieve corresponding DNA or protein coding sequences. Other examples of the types of operations that are easy to perform with genome-centric annotations include defining regions of interest for counting aligned reads in RNA-seq experiments and retrieving DNA sequences underlying regions of interest in ChIP-seq analysis, e.g., for motif characterization.

2.2 rtracklayer

The rtracklayer package allows us to query the UCSC genome browser, as well as providing import() and export() functions for common annotation file formats like GFF, GTF, and BED. The exercise below illustrates some of the functionality of rtracklayer.

2.3 Exercises

Exercise 4: TxDb.* packages

  1. Install and attach the TxDb.Hsapiens.UCSC.hg38.knownGene package. This contains the gene models for Homo sapiens based on the ‘hg38’ build of the human genome, using gene annotations in the UCSC ‘knownGene’ annotation track; TxDb’s for more recent builds and for different annotation tracks are available. Take a look at a summary of the package, and create an alias for easy typing

    library(TxDb.Hsapiens.UCSC.hg38.knownGene)
    TxDb.Hsapiens.UCSC.hg38.knownGene
    ## TxDb object:
    ## # Db type: TxDb
    ## # Supporting package: GenomicFeatures
    ## # Data source: UCSC
    ## # Genome: hg38
    ## # Organism: Homo sapiens
    ## # Taxonomy ID: 9606
    ## # UCSC Table: knownGene
    ## # UCSC Track: GENCODE v29
    ## # Resource URL: http://genome.ucsc.edu/
    ## # Type of Gene ID: Entrez Gene ID
    ## # Full dataset: yes
    ## # miRBase build ID: NA
    ## # transcript_nrow: 226811
    ## # exon_nrow: 647025
    ## # cds_nrow: 301804
    ## # Db created by: GenomicFeatures package from Bioconductor
    ## # Creation time: 2019-04-25 15:16:59 +0000 (Thu, 25 Apr 2019)
    ## # GenomicFeatures version at creation time: 1.35.11
    ## # RSQLite version at creation time: 2.1.1
    ## # DBSCHEMAVERSION: 1.2
    txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
  2. The main purpose of this package is to provide genomic coordinates of genomic features such as exons(), coding sequences (cds()), transcripts() and genes(). Explore, for example,

    ex <- exons(txdb)
    ex
    ## GRanges object with 647025 ranges and 1 metadata column:
    ##                    seqnames        ranges strand |   exon_id
    ##                       <Rle>     <IRanges>  <Rle> | <integer>
    ##        [1]             chr1   11869-12227      + |         1
    ##        [2]             chr1   12010-12057      + |         2
    ##        [3]             chr1   12179-12227      + |         3
    ##        [4]             chr1   12613-12697      + |         4
    ##        [5]             chr1   12613-12721      + |         5
    ##        ...              ...           ...    ... .       ...
    ##   [647021] chrUn_GL000220v1 155997-156149      + |    647021
    ##   [647022] chrUn_KI270442v1 380608-380726      + |    647022
    ##   [647023] chrUn_KI270442v1 217250-217401      - |    647023
    ##   [647024] chrUn_KI270744v1   51009-51114      - |    647024
    ##   [647025] chrUn_KI270750v1 148668-148843      + |    647025
    ##   -------
    ##   seqinfo: 595 sequences (1 circular) from hg38 genome
    library(ggplot2)
    qplot(log10(width(ex)))
    ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

    ex[ which.max(width(ex)) ]
    ## GRanges object with 1 range and 1 metadata column:
    ##       seqnames              ranges strand |   exon_id
    ##          <Rle>           <IRanges>  <Rle> | <integer>
    ##   [1]    chr12 102197585-102402596      + |    350923
    ##   -------
    ##   seqinfo: 595 sequences (1 circular) from hg38 genome
  3. Extract all genes, and then keep only the ‘standard’ chromosomes 1:22, X, Y, and M. Use table() of seqnames() to determine how many genes are on each chromosome. Also do this in a dplyr way; note that the seqnames(gn) need to be coerced with as.factor().

    gn <- genes(txdb)
    length(gn)
    ## [1] 26034
    std <- paste0("chr", c(1:22, "X", "Y", "M"))
    seqlevels(gn, pruning.mode = "coarse") <- std
    length(gn)
    ## [1] 26012
    seqlevels(gn)
    ##  [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9" 
    ## [10] "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18"
    ## [19] "chr19" "chr20" "chr21" "chr22" "chrX"  "chrY"  "chrM"
    table( seqnames(gn) )
    ## 
    ##  chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9 chr10 chr11 chr12 
    ##  2625  1711  1497  1043  1206  1146  1231   922  1047  1039  1577  1315 
    ## chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22  chrX  chrY 
    ##   592   869   892  1036  1407   421  1655   750   347   582  1038    62 
    ##  chrM 
    ##     2
    tibble(chr = as.factor(seqnames(gn))) %>% 
        group_by(chr) %>% 
        summarize(n = n())
    ## # A tibble: 25 x 2
    ##    chr       n
    ##    <fct> <int>
    ##  1 chr1   2625
    ##  2 chr2   1711
    ##  3 chr3   1497
    ##  4 chr4   1043
    ##  5 chr5   1206
    ##  6 chr6   1146
    ##  7 chr7   1231
    ##  8 chr8    922
    ##  9 chr9   1047
    ## 10 chr10  1039
    ## # … with 15 more rows
  4. exonsBy() groups exons by gene or transcript; extract exons grouped by gene. (Challenging!) can you identify genes with exons on different chromosomes? Are there any of these genes on the standard chromosomes?

    exByGn <- exonsBy(txdb, "gene")
    ##
    trans <- lengths(unique(seqnames(exByGn)))
    table( trans )
    ## trans
    ##     1     2     3     4     5     6     7     8     9    10    13    16 
    ## 26175   629   111    27    34    39    79    72     2    17     1     1 
    ##    22    24    27    28    30    33 
    ##     1     1     1     1     3     1
    seqnames( exByGn[ trans > 1 ] )
    ## RleList of length 1020
    ## $`10000`
    ## factor-Rle of length 56 with 2 runs
    ##   Lengths:                  30                  26
    ##   Values :                chr1 chr1_KI270763v1_alt
    ## Levels(595): chr1 chr2 chr3 ... chrUn_KI270756v1 chrUn_KI270757v1
    ## 
    ## $`100037417`
    ## factor-Rle of length 6 with 2 runs
    ##   Lengths:                    3                    3
    ##   Values :                chr22 chr22_KI270879v1_alt
    ## Levels(595): chr1 chr2 chr3 ... chrUn_KI270756v1 chrUn_KI270757v1
    ## 
    ## $`100049076`
    ## factor-Rle of length 20 with 2 runs
    ##   Lengths:                  11                   9
    ##   Values :                chr5 chr5_GL339449v2_alt
    ## Levels(595): chr1 chr2 chr3 ... chrUn_KI270756v1 chrUn_KI270757v1
    ## 
    ## $`100128260`
    ## factor-Rle of length 4 with 2 runs
    ##   Lengths:    2    2
    ##   Values : chrX chrY
    ## Levels(595): chr1 chr2 chr3 ... chrUn_KI270756v1 chrUn_KI270757v1
    ## 
    ## $`100128292`
    ## factor-Rle of length 2 with 2 runs
    ##   Lengths:                    1                    1
    ##   Values :                chr10 chr10_KI270825v1_alt
    ## Levels(595): chr1 chr2 chr3 ... chrUn_KI270756v1 chrUn_KI270757v1
    ## 
    ## ...
    ## <1015 more elements>
    ##
    std <- paste0("chr", c(1:22, "X", "Y", "M"))
    unames <- unique(seqnames(exByGn[ trans > 1 ]))
    transstd <- all(unames %in% std)
    unames[transstd]
    ## FactorList of length 40
    ## [["100128260"]] chrX chrY
    ## [["100132062"]] chr1 chr5
    ## [["100302278"]] chr1 chr15
    ## [["100359394"]] chrX chrY
    ## [["100422831"]] chr9 chr19
    ## [["100422872"]] chr2 chr5
    ## [["100500894"]] chrX chrY
    ## [["100616345"]] chr1 chr16
    ## [["101929127"]] chr2 chr9
    ## [["101954271"]] chr1 chr2 chr3 chr7 chr10 chr15 chr19
    ## ...
    ## <30 more elements>
  5. The previous exercise indicated that gene "22947" has exons on both chromosomes 4 and 10. Find out more about this gene using the org.Hs.eg.db package and by searching for the gene symbol on the NCBI web site.

    egid <- "22947"
    AnnotationDbi::select(
        org.Hs.eg.db, egid, c("SYMBOL", "GENENAME"), "ENTREZID"
    )
    ## 'select()' returned 1:1 mapping between keys and columns
    ##   ENTREZID SYMBOL                              GENENAME
    ## 1    22947 DUX4L1 double homeobox 4 like 1 (pseudogene)
    url <- paste0("https://www.ncbi.nlm.nih.gov/gene/", egid)
    browseURL(url)
  6. Note that the TxDb.* packages also support keytypes(), columns(), and select() for mapping between exon, cds, transcript, and gene identifiers.

Exercise 5: BSgenome.* packages

  1. Install (if necessary) and load the BSgenome.Hsapiens.UCSC.hg38 package, containing the entire sequence of the hg38 build of Homo sapiens. Check out it’s contents, and create a simple alias.

    library(BSgenome.Hsapiens.UCSC.hg38)
    BSgenome.Hsapiens.UCSC.hg38
    ## Human genome:
    ## # organism: Homo sapiens (Human)
    ## # provider: UCSC
    ## # provider version: hg38
    ## # release date: Dec. 2013
    ## # release name: Genome Reference Consortium GRCh38
    ## # 455 sequences:
    ## #   chr1                    chr2                    chr3                   
    ## #   chr4                    chr5                    chr6                   
    ## #   chr7                    chr8                    chr9                   
    ## #   chr10                   chr11                   chr12                  
    ## #   chr13                   chr14                   chr15                  
    ## #   ...                     ...                     ...                    
    ## #   chrUn_KI270744v1        chrUn_KI270745v1        chrUn_KI270746v1       
    ## #   chrUn_KI270747v1        chrUn_KI270748v1        chrUn_KI270749v1       
    ## #   chrUn_KI270750v1        chrUn_KI270751v1        chrUn_KI270752v1       
    ## #   chrUn_KI270753v1        chrUn_KI270754v1        chrUn_KI270755v1       
    ## #   chrUn_KI270756v1        chrUn_KI270757v1                               
    ## # (use 'seqnames()' to see all the sequence names, use the '$' or '[['
    ## # operator to access a given sequence)
    hg38 <- BSgenome.Hsapiens.UCSC.hg38
  2. Genomic sequence can be retrieved by chromosome, e.g., hg38[["chr1"]], or by genomic range, e.g., getSeq(hg38, GRanges("chr1:1000000-2000000")). Retrieve your favorite chunk(s) of DNA and calculate GC content.

    dna <- getSeq(hg38, GRanges("chr1:1000000-2000000"))
    letterFrequency(dna, "GC", as.prob=TRUE)
    ##            G|C
    ## [1,] 0.5728534
  3. Use the org.*, TxDb.*, and BSgenome.* packages to retrieve the BRCA1 exon DNA sequence.

    brca1_egid <- mapIds(org.Hs.eg.db, "BRCA1", "ENTREZID", "SYMBOL")
    ## 'select()' returned 1:1 mapping between keys and columns
    brca1_exons <- exonsBy(txdb, "gene")[[brca1_egid]]
    getSeq(hg38, brca1_exons)
    ##   A DNAStringSet instance of length 80
    ##      width seq
    ##  [1]  1508 CAATTGGGCAGATGTGTGAGGCACCTGTGGTG...CTGCAAATAAACTTGGTAGCAAACACTTCCA
    ##  [2]   998 CAATTGGGCAGATGTGTGAGGCACCTGTGGTG...AGACTGTGGCTCAAAAAAAAAAAAAAAAAAA
    ##  [3]   717 CAATTGGGCAGATGTGTGAGGCACCTGTGGTG...AAGGAAACTTGAAACCTGGGCATGGTGGCTC
    ##  [4]   240 CAATTGGGCAGATGTGTGAGGCACCTGTGGTG...ACTCTTCAGTCCTTCTACTGTCCTGGCTACT
    ##  [5]   174 CAATTGGGCAGATGTGTGAGGCACCTGTGGTG...ACAGAGCCACAGGACCCCAAGAATGAGCTTA
    ##  ...   ... ...
    ## [76]    94 GAGCTCGCTGAGACTTCCTGGACGGGGGACAG...GCGCTCAGGAGGCCTTCACCCTCTGCTCTGG
    ## [77]   126 ACAGATAAATTAAAACTGCGACTGCGCGGCGT...GCGCTCAGGAGGCCTTCACCCTCTGCTCTGG
    ## [78]   174 TTAGCGGTAGCCCCTTGGTTTCCGTGGCAACG...GCGCTCAGGAGGCCTTCACCCTCTGCTCTGG
    ## [79]   175 CTTAGCGGTAGCCCCTTGGTTTCCGTGGCAAC...GCGCTCAGGAGGCCTTCACCCTCTGCTCTGG
    ## [80]   120 AAAGCGTGGGAATTACAGATAAATTAAAACTG...CCGCGTTGGGGTGAGACCCTCACTTCATCCG

Exercise 6

This exercise uses annotation resources to go from a gene symbol ‘BRCA1’ through to the genomic coordinates of each transcript associated with the gene, and finally to the DNA sequences of the transcripts. This can be achieved using an EnsDb package along with a BSgenome package, or with a combination of TxDb, Homo.sapiens and BSgenome packages. We will focus here on the former approach.

  1. Use AnnotationHub to discover and retrieve a current Ensembl annotation (‘EnsDb’) for Homo sapiens.

  2. Use the cdsBy() function to retrieve the genomic coordinates of all coding sequences for the gene ‘BRCA1’ from the EnsDb.Hsapiens.v86 package. To retrieve only data for the specified gene, submit either a GenenameFilter or a filter formula/expression to the function’s filter parameter. This avoids to extract the coding region for all genes, which takes a long time.

  3. Visualize the transcripts in genomic coordinates using the Gviz package to construct a GeneRegionTrack, and plotting it using plotTracks().

  4. Use the Bsgenome.Hsapiens.UCSC.hg38 package and extractTranscriptSeqs() function to extract the DNA sequence of each transcript.

Solution

Retrieve the coding sequences grouped by transcript for the gene of interest and verify that each coding sequence is a multiple of 3.

library(EnsDb.Hsapiens.v86)
edb <- EnsDb.Hsapiens.v86

brca1cds <- cdsBy(edb, by = "tx", filter = ~ genename == "BRCA1")

class(brca1cds)
## [1] "CompressedGRangesList"
## attr(,"package")
## [1] "GenomicRanges"
length(brca1cds)
## [1] 29
brca1cds[[1]]                           # exons in cds
## GRanges object with 21 ranges and 3 metadata columns:
##        seqnames            ranges strand |   gene_name         exon_id
##           <Rle>         <IRanges>  <Rle> | <character>     <character>
##    [1]       17 43124017-43124096      - |       BRCA1 ENSE00003559512
##    [2]       17 43115726-43115779      - |       BRCA1 ENSE00003510592
##    [3]       17 43106456-43106533      - |       BRCA1 ENSE00003541068
##    [4]       17 43104868-43104956      - |       BRCA1 ENSE00003531836
##    [5]       17 43104122-43104261      - |       BRCA1 ENSE00003513709
##    ...      ...               ...    ... .         ...             ...
##   [17]       17 43057052-43057135      - |       BRCA1 ENSE00003458468
##   [18]       17 43051063-43051117      - |       BRCA1 ENSE00003477922
##   [19]       17 43049121-43049194      - |       BRCA1 ENSE00003628864
##   [20]       17 43047643-43047703      - |       BRCA1 ENSE00003687053
##   [21]       17 43045678-43045802      - |       BRCA1 ENSE00001814242
##        exon_rank
##        <integer>
##    [1]         2
##    [2]         3
##    [3]         4
##    [4]         5
##    [5]         6
##    ...       ...
##   [17]        18
##   [18]        19
##   [19]        20
##   [20]        21
##   [21]        22
##   -------
##   seqinfo: 1 sequence from GRCh38 genome
cdswidth <- width(brca1cds)             # width of each exon
all((sum(cdswidth) %% 3) == 0)          # sum within cds, modulus 3
## [1] FALSE

The CDS for some transcripts is not of the expected length, how come? Get the transcript ID of the first transcript that does have a CDS of the wrong size and look this transcript up in the Ensembl genome browser (http://www.ensembl.org).

tx_cds_fail <- names(brca1cds)[(sum(cdswidth) %% 3) != 0]

length(tx_cds_fail)
## [1] 8
tx_cds_fail[1]
## [1] "ENST00000412061"

In the description of the transcript it says CDS 5’ incomplete. Thus, in addition to known protein coding transcripts, Ensembl provides annotations for transcripts known to be targeted for nonsense mediated mRNA decay or that have incomplete CDS. Such transcripts would however not be listed in e.g. the TxDb.Hsapiens.UCSC.hg38.knownGene package.

Next we visualize the BRCA1 transcripts using Gviz (this package has an excellent vignette, vignette("Gviz"))

library(Gviz)

## Use the function from the ensembldb package to extract the data in the
## format suitable for Gviz
grt <- getGeneRegionTrackForGviz(edb, filter = ~genename == "BRCA1")
plotTracks(list(GenomeAxisTrack(), GeneRegionTrack(grt)))

Extract the coding sequences of each transcript. EnsDb databases provide annotations from Ensembl and use hence Ensembl style chromosome names (such as “Y”) while the BSgenome package is based on UCSC annotations that use a naming style that prepends a “chr” to each chromosome name (e.g. “chrY”). Change thus the seqlevelsStyle from the default UCSC chromosome naming to Ensembl naming style.

library(BSgenome.Hsapiens.UCSC.hg19)
## 
## Attaching package: 'BSgenome.Hsapiens.UCSC.hg19'
## The following object is masked from 'package:BSgenome.Hsapiens.UCSC.hg38':
## 
##     Hsapiens
genome <- BSgenome.Hsapiens.UCSC.hg19

## Change the seqlevelsStyle from UCSC to Ensembl
seqlevelsStyle(genome) <- "Ensembl"
tx_seq <- extractTranscriptSeqs(genome, brca1cds)
tx_seq
##   A DNAStringSet instance of length 29
##      width seq                                           names               
##  [1]  4704 ATGAATGTAGAAAAGGCTGAA...ATCCCCCACAGCCACTACTGA ENST00000309486
##  [2]  4875 ATGGATTTATCTGCTCTTCGC...ATCCCCCACAGCCACTACTGA ENST00000346315
##  [3]  2043 ATGGATTTATCTGCTCTTCGC...ATCCCCCACAGCCACTACTGA ENST00000351666
##  [4]  2166 ATGGATTTATCTGCTCTTCGC...ATCCCCCACAGCCACTACTGA ENST00000352993
##  [5]  4797 ATGGATTTATCTGCTCTTCGC...ATCCCCCACAGCCACTACTGA ENST00000354071
##  ...   ... ...
## [25]  1419 ATGGATTTATCTGCTCTTCGC...AAGAAGGCAAGCCTCCCCAAC ENST00000494123
## [26]   531 ATGAATGTAGAAAAGGCTGAA...AAGAAGGCAAGCCTCCCCAAC ENST00000497488
## [27]   522 ATGGATGCTGAGTTTGTGTGT...ATCCCCCACAGCCACTACTGA ENST00000586385
## [28]  1065 ATGCACAGTTGCTCTGGGAGT...ATCCCCCACAGCCACTACTGA ENST00000591534
## [29]   291 ATGAGTGACAGCAAGAAAACC...ATCCCCCACAGCCACTACTGA ENST00000591849

We can also inspect the CDS sequence for the transcripts with incomplete CDS. Many of them do not start with a start codon hence indicating that the CDS is incomplete on their 5’ end.

tx_seq[tx_cds_fail]
##   A DNAStringSet instance of length 8
##     width seq                                            names               
## [1]  1312 GTTTGGATTCTGCAAAAAAGGC...TGAAGAGATAAAGAAAAAAAA ENST00000412061
## [2]   958 TTCAGCTTGACACAGGTTTGGA...TCACTCCAAATCAGTAGAGAG ENST00000473961
## [3]   667 ATGGATTTATCTGCTCTTCGCG...GTTTGGATTCTGCAAAAAAGG ENST00000476777
## [4]  1867 ATGGATTTATCTGCTCTTCGCG...ATAGTTGTTCTAGCAGTGAAG ENST00000477152
## [5]  1870 ATGGATTTATCTGCTCTTCGCG...AGTCTATTAAAGAAAGAAAAA ENST00000478531
## [6]  1495 GAGCTATTGAAAATCATTTGTG...AGTCTATTAAAGAAAGAAAAA ENST00000484087
## [7]   800 GAGCTATTGAAAATCATTTGTG...TAAAGAACCAGGAGTGGAAAG ENST00000487825
## [8]   296 ATGGATTTATCTGCTCTTCGCG...AAAAGATGAAGTTTCTATCAT ENST00000489037

Intron coordinates can be identified by first calculating the range of the genome (from the start of the first exon to the end of the last exon) covered by each transcript, and then taking the (algebraic) set difference between this and the genomic coordinates covered by each exon

introns <- psetdiff(unlist(range(brca1cds)), brca1cds)

Retrieve the intronic sequences with getSeq() (these are not assembled, the way that extractTranscriptSeqs() assembles exon sequences into mature transcripts); note that introns start and end with the appropriate acceptor and donor site sequences. Unfortunately, UCSC and Ensembl do also use different names for the genome assembly. Change the genome name for the introns object to matche the one from the genome object.

unique(genome(genome))
## [1] "hg19"
genome(introns)
##       17 
## "GRCh37"
## Change the genome name on introns to match the one from the
## BSgenome package
genome(introns) <- c(`17` = unique(genome(genome)))

seq <- getSeq(genome, introns)
names(seq)
##  [1] "ENST00000309486" "ENST00000346315" "ENST00000351666" "ENST00000352993"
##  [5] "ENST00000354071" "ENST00000357654" "ENST00000412061" "ENST00000461221"
##  [9] "ENST00000461574" "ENST00000461798" "ENST00000468300" "ENST00000470026"
## [13] "ENST00000471181" "ENST00000473961" "ENST00000476777" "ENST00000477152"
## [17] "ENST00000478531" "ENST00000484087" "ENST00000487825" "ENST00000489037"
## [21] "ENST00000491747" "ENST00000492859" "ENST00000493795" "ENST00000493919"
## [25] "ENST00000494123" "ENST00000497488" "ENST00000586385" "ENST00000591534"
## [29] "ENST00000591849"
seq[["ENST00000352993"]]                     # 20 introns
##   A DNAStringSet instance of length 20
##      width seq
##  [1]  1840 GTAAGGTGCCTGCATGTACCTGTGCTATATGG...ACTAATCTCTGCTTGTGTTCTCTGTCTCCAG
##  [2]  1417 GTAAGTATTGGGTGCCCTGTCAGAGAGGGAGG...TTTGAATGCTCTTTCCTTCCTGGGGATCCAG
##  [3]  1868 GTAAGAGCCTGGGAGAACCCCAGAGTTCCAGC...CAGTGATTTTACATCTAAATGTCCATTTTAG
##  [4]  5934 GTAAAGCTCCCTCCCTCAAGTTGACAAAAATC...CCTGTCCCTCTCTCTTCCTCTCTTCTTCCAG
##  [5]  6197 GTAAGTACTTGATGTTACAAACTAACCAGAGA...ATCCTGATGGGTTGTGTTTGGTTTCTTTCAG
##  ...   ... ...
## [16]  4241 GTAAAACCATTTGTTTTCTTCTTCTTCTTCTT...TTGCTTGACTGTTCTTTACCATACTGTTTAG
## [17]   606 GTAAGTGTTGAATATCCCAAGAATGACACTCA...CAAACATAATGTTTTCCCTTGTATTTTACAG
## [18]  1499 GTATATAATTTGGTAATGATGCTAGGTTGGAA...TGAGTGTGTTTCTCAAACAATTTAATTTCAG
## [19]  9192 GTAAGTTTGAATGTGTTATGTGGCTCCATTAT...AAATTGTTCTTTCTTTCTTTATAATTTATAG
## [20]  8237 GTAAGTCAGCACAAGAGTGTATTAATTTGGGA...ATTTTCTTTTTCTCCCCCCCTACCCTGCTAG

Exercise 7

Internet access required for this exercise

Here we use rtracklayer to retrieve estrogen receptor binding sites identified across cell lines in the ENCODE project. We focus on binding sites in the vicinity of a particularly interesting region.

  1. Define our region of interest by creating a GRanges instance with appropriate genomic coordinates. Our region corresponds to 10Mb up- and down-stream of a particular gene.
  2. Create a session for the UCSC genome browser
  3. Query the UCSC genome browser for ENCODE estrogen receptor ERalphaa transcription marks; identifying the appropriate track, table, and transcription factor requires biological knowledge and detective work.
  4. Visualize the location of the binding sites and their scores; annotate the mid-point of the region of interest.

Solution

Define the region of interest

library(GenomicRanges)
roi <- GRanges("chr10", IRanges(92106877, 112106876, names="ENSG00000099194"))

Create a session

library(rtracklayer) 
session <- browserSession()

Query the UCSC for a particular track, table, and transcription factor, in our region of interest

trackName <- "wgEncodeRegTfbsClusteredV2"
tableName <- "wgEncodeRegTfbsClusteredV2"
trFactor <- "ERalpha_a"
ucscTable <- getTable(ucscTableQuery(session, track=trackName,
    range=roi, table=tableName, name=trFactor))

Visualize the result

plot(score ~ chromStart, ucscTable, pch="+")
abline(v=start(roi) + (end(roi) - start(roi) + 1) / 2, col="blue")

3 AnnotationHub

AnnotationHub is a data base of large-scale whole-genome resources, e.g., regulatory elements from the Roadmap Epigenomics project, Ensembl GTF and FASTA files for model and other organisms, and the NHLBI grasp2db data base of GWAS results. There are many interesting ways in which these resources can be used. Examples include

Unfortunately, AnnotationHub makes extensive use of internet resources and so we will not pursue it in this course; see the vignettes that come with the pacakge, for instance AnnotationHub HOW-TOs.

4 Annotating variants

Bioconductor provides facilities for reading VCF files. These work very well with the annotation resources described above, so for instance it is straight-forward to identify variants in coding or other regions of interest.

To develop a sense of the capabilities available, work through the VariantAnnotation vignette ‘Introduction to Variant Annotation’, and the VariantFiltering vignette.

5 End matter

5.1 Session Info

sessionInfo()
## R version 3.6.1 Patched (2019-07-15 r76837)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.6 LTS
## 
## Matrix products: default
## BLAS:   /home/lori/bin/R-3-6-branch/lib/libRblas.so
## LAPACK: /home/lori/bin/R-3-6-branch/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] BSgenome.Hsapiens.UCSC.hg19_1.4.0      
##  [2] KEGGREST_1.25.0                        
##  [3] forcats_0.4.0                          
##  [4] stringr_1.4.0                          
##  [5] dplyr_0.8.3                            
##  [6] purrr_0.3.2                            
##  [7] readr_1.3.1                            
##  [8] tidyr_0.8.3                            
##  [9] tibble_2.1.3                           
## [10] ggplot2_3.2.0                          
## [11] tidyverse_1.2.1                        
## [12] AnnotationHub_2.17.5                   
## [13] BiocFileCache_1.9.1                    
## [14] dbplyr_1.4.2                           
## [15] Gviz_1.29.0                            
## [16] biomaRt_2.41.7                         
## [17] BSgenome.Hsapiens.UCSC.hg38_1.4.1      
## [18] BSgenome_1.53.0                        
## [19] rtracklayer_1.45.1                     
## [20] Biostrings_2.53.2                      
## [21] XVector_0.25.0                         
## [22] EnsDb.Hsapiens.v75_2.99.0              
## [23] ensembldb_2.9.2                        
## [24] AnnotationFilter_1.9.0                 
## [25] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.6
## [26] GenomicFeatures_1.37.4                 
## [27] GenomicRanges_1.37.14                  
## [28] GenomeInfoDb_1.21.1                    
## [29] org.Hs.eg.db_3.8.2                     
## [30] AnnotationDbi_1.47.0                   
## [31] IRanges_2.19.10                        
## [32] S4Vectors_0.23.17                      
## [33] Biobase_2.45.0                         
## [34] BiocGenerics_0.31.5                    
## [35] BiocStyle_2.13.2                       
## [36] praise_1.0.0                           
## [37] BiocManager_1.30.4                     
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1                  backports_1.1.4              
##   [3] Hmisc_4.2-0                   lazyeval_0.2.2               
##   [5] splines_3.6.1                 BiocParallel_1.19.0          
##   [7] digest_0.6.20                 htmltools_0.3.6              
##   [9] fansi_0.4.0                   magrittr_1.5                 
##  [11] checkmate_1.9.4               memoise_1.1.0                
##  [13] cluster_2.1.0                 modelr_0.1.4                 
##  [15] matrixStats_0.54.0            askpass_1.1                  
##  [17] prettyunits_1.0.2             colorspace_1.4-1             
##  [19] blob_1.2.0                    rvest_0.3.4                  
##  [21] rappdirs_0.3.1                haven_2.1.1                  
##  [23] xfun_0.8                      crayon_1.3.4                 
##  [25] RCurl_1.95-4.12               jsonlite_1.6                 
##  [27] zeallot_0.1.0                 survival_2.44-1.1            
##  [29] VariantAnnotation_1.31.3      glue_1.3.1                   
##  [31] gtable_0.3.0                  zlibbioc_1.31.0              
##  [33] DelayedArray_0.11.4           scales_1.0.0                 
##  [35] DBI_1.0.0                     Rcpp_1.0.1                   
##  [37] xtable_1.8-4                  progress_1.2.2               
##  [39] htmlTable_1.13.1              foreign_0.8-71               
##  [41] bit_1.1-14                    Formula_1.2-3                
##  [43] htmlwidgets_1.3               httr_1.4.0                   
##  [45] RColorBrewer_1.1-2            acepack_1.4.1                
##  [47] pkgconfig_2.0.2               XML_3.98-1.20                
##  [49] nnet_7.3-12                   utf8_1.1.4                   
##  [51] tidyselect_0.2.5              labeling_0.3                 
##  [53] rlang_0.4.0                   later_0.8.0                  
##  [55] munsell_0.5.0                 cellranger_1.1.0             
##  [57] tools_3.6.1                   cli_1.1.0                    
##  [59] generics_0.0.2                RSQLite_2.1.1                
##  [61] broom_0.5.2                   evaluate_0.14                
##  [63] yaml_2.2.0                    knitr_1.23                   
##  [65] bit64_0.9-7                   nlme_3.1-140                 
##  [67] mime_0.7                      xml2_1.2.0                   
##  [69] compiler_3.6.1                rstudioapi_0.10              
##  [71] curl_3.3                      png_0.1-7                    
##  [73] interactiveDisplayBase_1.23.0 stringi_1.4.3                
##  [75] lattice_0.20-38               ProtGenerics_1.17.2          
##  [77] Matrix_1.2-17                 vctrs_0.2.0                  
##  [79] pillar_1.4.2                  data.table_1.12.2            
##  [81] bitops_1.0-6                  httpuv_1.5.1                 
##  [83] R6_2.4.0                      latticeExtra_0.6-28          
##  [85] bookdown_0.12                 promises_1.0.1               
##  [87] gridExtra_2.3                 codetools_0.2-16             
##  [89] dichromat_2.0-0               assertthat_0.2.1             
##  [91] SummarizedExperiment_1.15.5   openssl_1.4.1                
##  [93] withr_2.1.2                   GenomicAlignments_1.21.4     
##  [95] Rsamtools_2.1.3               GenomeInfoDbData_1.2.1       
##  [97] hms_0.5.0                     rpart_4.1-15                 
##  [99] rmarkdown_1.14                biovizBase_1.33.0            
## [101] shiny_1.3.2                   lubridate_1.7.4              
## [103] base64enc_0.1-3

5.2 Acknowledgements

Research reported in this tutorial was supported by the National Human Genome Research Institute and the National Cancer Institute of the National Institutes of Health under award numbers U41HG004059 and U24CA180996.

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 633974)