Lecture 3

Set up

using <- function(...) {
    libs <- unlist(list(...))
    need <- libs[!unlist(lapply(libs, require, character.only = T))]
    if(length(need) > 0){ 
        install.packages(need, repos = "https://cloud.r-project.org")
        need <- need[!unlist(lapply(need, require, character.only = T))]
        if (length(need) > 0) {
          if (!requireNamespace("BiocManager", quietly = T))
            install.packages("BiocManager")
          BiocManager::install(need)
        }
    }
}
using("ggplot2", "dplyr", "ggseqlogo", "GenomicRanges",
      "BSgenome.Hsapiens.UCSC.hg19", "biomaRt")

Gene annotations

As an R implementation of the BioMart project, biomaRt facilitates access to a variety of resources ranging from Ensembl Genomes to the COSMIC cancer database. The chunk below retrieves the Ensembl gene set, subsets autosomal genes with valid HGNC symbols, converts the chromosome names to UCSC scheme, and finally derives a GRanges object from it.

# the useMart function connects to a specified BioMart database and dataset within this database.
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", host = "grch37.ensembl.org",
                 path = "/biomart/martservice", dataset = "hsapiens_gene_ensembl")
# seqinfo extracts the sequence information stored in a DNAStringSet object
si <- Seqinfo(genome = "hg19")[paste0("chr", 1:22)]
# Given a set of filters and corresponding values, getBM retrieves the user specified attributes from the BioMart database
genes <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "start_position",
                              "end_position", "strand", "ensembl_gene_id"),
        filters = "chromosome_name", values = as.character(1:22), mart = mart) %>%
  arrange(chromosome_name, start_position) %>%
  filter(hgnc_symbol != "") %>% # remove any element without a gene symbol
  mutate(strand = ifelse(strand == 1, '+', '-'),
         chromosome_name = paste0('chr', chromosome_name)) %>%
  makeGRangesFromDataFrame(keep.extra.columns = T,
                           start.field = "start_position",
                           end.field = "end_position",
                           seqinfo = si) # make granges object
# The GRanges class is a container for the genomic locations and their associated annotations.
genes

## GRanges object with 24827 ranges and 2 metadata columns:
##           seqnames            ranges strand | hgnc_symbol ensembl_gene_id
##              <Rle>         <IRanges>  <Rle> | <character>     <character>
##       [1]     chr1       11869-14412      + |     DDX11L1 ENSG00000223972
##       [2]     chr1       14363-29806      - |      WASH7P ENSG00000227232
##       [3]     chr1       29554-31109      + |  MIR1302-11 ENSG00000243485
##       [4]     chr1       29554-31109      + |  MIR1302-10 ENSG00000243485
##       [5]     chr1       29554-31109      + |   MIR1302-9 ENSG00000243485
##       ...      ...               ...    ... .         ...             ...
##   [24823]    chr22 50528308-50600119      + |     MOV10L1 ENSG00000073146
##   [24824]    chr22 50609160-50618723      + |       PANX2 ENSG00000073150
##   [24825]    chr22 50624344-50638027      + |       TRABD ENSG00000170638
##   [24826]    chr22 50656118-50683421      - |     TUBGCP6 ENSG00000128159
##   [24827]    chr22 50683612-50689834      - |      HDAC10 ENSG00000100429
##   -------
##   seqinfo: 22 sequences from hg19 genome

Sequence analysis

set.seed(42)
genes.sub <- genes[strand(genes) == '+'] %>%
  sample(1000) # sample 1000 genes from positive strand
# get exons from bioMart database
exons.sub <- getBM(attributes = c("ensembl_exon_id", "exon_chrom_start", "exon_chrom_end",
                              "chromosome_name", "strand"), filters = "ensembl_gene_id",
               values = genes.sub$ensembl_gene_id, mart = mart) %>%
  arrange(chromosome_name, exon_chrom_start) %>%
  mutate(strand = ifelse(strand == 1, '+', '-'),
         chromosome_name = paste0('chr', chromosome_name)) %>%
  makeGRangesFromDataFrame(keep.extra.columns = T,
                           start.field = "exon_chrom_start",
                           end.field = "exon_chrom_end",
                           seqinfo = si)
# resize exons
resize(exons.sub, fix = 'start', width = 10) %>% # resize returns an object of the same type and length as x containing intervals that have been resized to width width based on the strand(x) values
  resize(fix = 'end', width = 20) %>%
  getSeq(BSgenome.Hsapiens.UCSC.hg19, .) %>%
  consensusMatrix() %>%
  .[c('A', 'C', 'G', 'T'),] %>%
  ggseqlogo() # ggseqlogo is a shortcut for generating sequence logos

Additional exercise

Plot the sequence of splice donor sites using the same exons.sub object

---