Lab 1.6: Annotation Resources
22 July 2019
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:
- Use
org.*packages for mapping between gene symbols. - Use
TxDb.*andensembldb(EnsDb.*) packages for working with gene models. - Use
AnnotationHubto easily obtain select consortium-level resources - Access
biomaRtand other internet-based resources for highly flexible annotation. - Use
VariantAnnotationandVariantFilteringfor annotating SNPs.
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.
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)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 informationThe idea is that there are
keytypes()that can be mapped to differentcolumns();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"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, andselect()(it’s often necessary to specify this function directly, to avoid a conflict with dplyr, asAnnotationDbi::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 2Some key - column mappings are 1:many, e.g., Entrez ID
"3812"maps to 44 Ensembl Ids. What doesmapIds()return when mapping Entrez ID"3812"to Ensembl ids? Use the additional argumentmultiVals = "CharacterList"to explore further. Compare results to those returned byselect().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 ... ENSG00000284381AnnotationDbi::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 ENSG00000284381It seems like it might often be useful to use the tidyverse on return values from
mapIds()andselect(); explore this usagelibrary(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 rowsAnnotationDbi::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
- 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.
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 97mart <- useMart("ENSEMBL_MART_ENSEMBL")Use
listDatasets()anduseDataset()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 ailMel1dataset <- useDataset("hsapiens_gene_ensembl", mart)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 Startfilters <- "ensembl_gene_id" # see `listFilters()`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_pageattrs <- c("ensembl_gene_id", "hgnc_symbol") # see `listAttributes()`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 foundtbl## # 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
Explore the KEGG web site https://www.genome.jp/kegg/ KEGG is a database of information on pathways.
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"Use
keggList()to query the pathway database for human pathways; present the result as a tibblehsa_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 rowsUse
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 rowshsa_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 rowsUpdate the
hsa_path_egtable to include information on gene symbol and Ensembl id from theorg.Hs.eg.dbpackage. Retrieve the relevant information usingmapIds(). 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 columnsUse
left_join()to append pathway descriptions to thehsa_kegg_annotable.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:
TxDb.*andEnsDb.*For transcript and other genome / coordinate annotation.- BSgenome For whole-genome representation. See
available.genomes()for pre-packaged genomes, and the vignette ‘How to forge a BSgenome data package’ in the - Homo.sapiens For integrating ‘TxDb’ and ’org.’ packages.
- ’SNPlocs.*’ For model organism SNP locations derived from dbSNP.
FaFile()(Rsamtools) for accessing indexed FASTA files.- [ensemblVEP][] Variant effect scores.
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
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.2txdb <- TxDb.Hsapiens.UCSC.hg38.knownGeneThe main purpose of this package is to provide genomic coordinates of genomic features such as
exons(), coding sequences (cds()),transcripts()andgenes(). 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 genomelibrary(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 genomeExtract all genes, and then keep only the ‘standard’ chromosomes 1:22, X, Y, and M. Use
table()ofseqnames()to determine how many genes are on each chromosome. Also do this in a dplyr way; note that theseqnames(gn)need to be coerced withas.factor().gn <- genes(txdb) length(gn)## [1] 26034std <- paste0("chr", c(1:22, "X", "Y", "M")) seqlevels(gn, pruning.mode = "coarse") <- std length(gn)## [1] 26012seqlevels(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 ## 2tibble(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 rowsexonsBy()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 1seqnames( 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>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)Note that the
TxDb.*packages also supportkeytypes(),columns(), andselect()for mapping between exon, cds, transcript, and gene identifiers.
Exercise 5: BSgenome.* packages
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.hg38Genomic 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.5728534Use the
org.*,TxDb.*, andBSgenome.*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 columnsbrca1_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.
Use AnnotationHub to discover and retrieve a current Ensembl annotation (‘EnsDb’) for Homo sapiens.
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 aGenenameFilteror a filter formula/expression to the function’sfilterparameter. This avoids to extract the coding region for all genes, which takes a long time.Visualize the transcripts in genomic coordinates using the Gviz package to construct a
GeneRegionTrack, and plotting it usingplotTracks().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] 29brca1cds[[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 genomecdswidth <- width(brca1cds) # width of each exon
all((sum(cdswidth) %% 3) == 0) # sum within cds, modulus 3## [1] FALSEThe 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] 8tx_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':
##
## Hsapiensgenome <- 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 ENST00000591849We 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 ENST00000489037Intron 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...ATTTTCTTTTTCTCCCCCCCTACCCTGCTAGExercise 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.
- Define our region of interest by creating a
GRangesinstance with appropriate genomic coordinates. Our region corresponds to 10Mb up- and down-stream of a particular gene. - Create a session for the UCSC genome browser
- 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.
- 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
- Easily access and import Roadmap Epigenomics files.
- ‘liftOver’ genomic range-based annotations from one coordinate system (e.g, hg38) to another (e.g., GRCh 38);
- Create TranscriptDb and BSgenome-style annotation resources ‘on the fly’ for a diverse set of organisms.
- Programmatically access the genomic coordiantes of clinically relevant variants cataloged in dbSNP.
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-35.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)