Goals for this workshop

library(BiocAnno2016)
library(hugene20sttranscriptcluster.db)
library(EnsDb.Mmusculus.v79)
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(Homo.sapiens)
library(BSgenome)
library(BSgenome.Hsapiens.UCSC.hg19)
library(AnnotationHub)

What do we mean by annotation?

Map a known ID to other functional or positional information

Specific goal

We have data and statistics, and we want to add other useful information

The end result might be as simple as a data.frame or HTML table, or as complex as a RangedSummarizedExperiment

Data containers

ExpressionSet

load(system.file("data/eset.Rdata", package = "BiocAnno2016"))
eset

ExpressionSet (continued)

head(exprs(eset))
head(pData(phenoData(eset)))

ExpressionSet (continued)

head(pData(featureData(eset)))

BioC containers vs basic structures

Pros

Cons

Annotation sources

df <- data.frame("Package type" = c("ChipDb","OrgDb","TxDb/EnsDb","OrganismDb","BSgenome","Others","AnnotationHub","biomaRt"),
                 Example = c("hugene20sttranscriptcluster.db","org.Hs.eg.db","TxDb.Hsapiens.UCSC.hg19.knownGene; EnsDb.Hsapiens.v75",
                             "Homo.sapiens","BSgenome.Hsapiens.UCSC.hg19","GO.db; KEGG.db", "Online resource","Online resource"), check.names = FALSE)
knitr::kable(df)

Interacting with AnnoDb packages

The main function is select:

select(annopkg, keys, columns, keytype)

Where

Simple example

Say we have analyzed data from an Affymetrix Human Gene ST 2.0 array and want to know what the genes are. For purposes of this lab, we just select some IDs at random.

library(hugene20sttranscriptcluster.db)
set.seed(12345)
ids <- featureNames(eset)[sample(1:25000, 5)]
ids
select(hugene20sttranscriptcluster.db, ids, "SYMBOL")

Questions!

How do you know what the central keys are?

More questions!

What keytypes or columns are available for a given annotation package?

keytypes(hugene20sttranscriptcluster.db)
columns(hugene20sttranscriptcluster.db)

Another example

There is one issue with select however.

ids <- c('16737401','16657436' ,'16678303')
select(hugene20sttranscriptcluster.db, ids, c("SYMBOL","MAP"))

The mapIds function

An alternative to select is mapIds, which gives control of duplicates

mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID")

Choices for multiVals

Default is first, where we just choose the first of the duplicates. Other choices are list, CharacterList, filter, asNA or a user-specified function.

mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "list")

Choices for multiVals (continued)

mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "CharacterList")
mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "filter")
mapIds(hugene20sttranscriptcluster.db, ids, "SYMBOL", "PROBEID", multiVals = "asNA")

ChipDb/OrgDb questions

Using either the hugene20sttranscriptcluster.db or org.Hs.eg.db package,

TxDb packages

TxDb packages contain positional information; the contents can be inferred by the package name

TxDb.Species.Source.Build.Table

TxDb.Dmelanogaster.UCSC.dm3.ensGene TxDb.Athaliana.BioMart.plantsmart22

EnsDb packages

EnsDb packages are similar to TxDb packages, but based on Ensembl mappings

EnsDb.Hsapiens.v79
EnsDb.Mmusculus.v79
EnsDb.Rnorvegicus.v79

Transcript packages

As with ChipDb and OrgDb packages, select and mapIds can be used to make queries

select(TxDb.Hsapiens.UCSC.hg19.knownGene, c("1","10"),
       c("TXNAME","TXCHROM","TXSTART","TXEND"), "GENEID")
select(EnsDb.Hsapiens.v79, c("1", "10"),
       c("GENEID","GENENAME","SEQNAME","GENESEQSTART","GENESEQEND"), "ENTREZID")

But this is not how one normally uses them...

GRanges

The normal use case for transcript packages is to extract positional information into a GRanges or GRangesList object. An example is the genomic position of all genes:

gns <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
gns

GRangesList

Or the genomic position of all transcripts by gene:

txs <- transcriptsBy(TxDb.Hsapiens.UCSC.hg19.knownGene)
txs

Other accessors

Why *Ranges objects

The main rationale for *Ranges objects is to allow us to easily select and subset data based on genomic position information. This is really powerful!

GRanges and GRangesLists act like data.frames and lists, and can be subsetted using the [ function. As a really artificial example:

txs[txs %over% gns[1:2,]]

*Ranges use cases

SummarizedExperiment objects

SummarizedExperiment objects are like ExpressionSets, but the row-wise annotations are GRanges, so you can subset by genomic locations:

TxDb exercises

OrganismDb packages

OrganismDb packages are meta-packages that contain an OrgDb, a TxDb, and a GO.db package and allow cross-queries between those packages.

All previous accessors work; select, mapIds, transcripts, etc.

library(Homo.sapiens)
Homo.sapiens

OrganismDb packages

head(genes(Homo.sapiens, columns = c("ENTREZID","ALIAS","UNIPROT")),4)

OrganismDb exercises

BSgenome packages

BSgenome packages contain sequence information for a given species/build. There are many such packages - you can get a listing using available.genomes

library(BSgenome)
head(available.genomes())

BSgenome packages

We can load and inspect a BSgenome package

library(BSgenome.Hsapiens.UCSC.hg19)
Hsapiens

BSgenome packages

The main accessor is getSeq, and you can get data by sequence (e.g., entire chromosome or unplaced scaffold), or by passing in a GRanges object, to get just a region.

getSeq(Hsapiens, "chr1")
getSeq(Hsapiens, gns["5467",])

The Biostrings package contains most of the code for dealing with these *StringSet objects - please see the Biostrings vignettes and help pages for more information.

BSgenome exercises

AnnotationHub

AnnotationHub is a package that allows us to query and download many different annotation objects, without having to explicitly install them.

library(AnnotationHub)
hub <- AnnotationHub()
library(AnnotationHub)
hub <- AnnotationHub()
hub

Querying AnnotationHub

Finding the 'right' resource on AnnotationHub is like using Google - a well posed query is necessary to find what you are after. Useful queries are based on

names(mcols(hub))

AnnotationHub Data providers

unique(hub$dataprovider)

AnnotationHub Data classes

unique(hub$rdataclass)

AnnotationHub Species

head(unique(hub$species))
length(unique(hub$species))

AnnotationHub Data sources

unique(hub$sourcetype)

AnnotationHub query

qry <- query(hub, c("granges","homo sapiens","ensembl"))
qry

AnnotationHub query

qry$sourceurl

Selecting AnnotationHub resource

whatIwant <- qry[["AH50377"]]

We can use these data as they are, or convert to a TxDb format:

GRCh38TxDb <- makeTxDbFromGRanges(whatIwant)
GRCh38TxDb

AnnotationHub exercises

biomaRt

The biomaRt package allows queries to an Ensembl Biomart server. We can see the choices of servers that we can use:

library(biomaRt)
listMarts()

biomaRt data sets

And we can then check for the available data sets on a particular server.

mart <- useMart("ENSEMBL_MART_ENSEMBL")
head(listDatasets(mart))

biomaRt queries

After setting up a mart object pointing to the server and data set that we care about, we can make queries. We first set up the mart object.

mart <- useMart("ENSEMBL_MART_ENSEMBL","hsapiens_gene_ensembl")

Queries are of the form

getBM(attributes, filters, values, mart)

where

biomaRt attributes and filters

Both attributes and filters have rather inscrutable names, but a listing can be accessed using

atrib <- listAttributes(mart)
filts <- listFilters(mart)
head(atrib)
head(filts)

biomaRt query

A simple example query

afyids <- c("1000_at","1001_at","1002_f_at","1007_s_at")
getBM(c("affy_hg_u95av2", "hgnc_symbol"), c("affy_hg_u95av2"), afyids, mart)

biomaRt exercises



jmacdon/BiocAnno2016 documentation built on May 19, 2019, 12:49 p.m.