Overview

• Introduction to Bioconductor
• GenomicRanges
• Reading Sequence alignments
• Annotations
• Exercises

Bioconductor

Bioconductor (BioC) is an open source, open development software project to provide tools for the analysis and comprehension of high-throughput genomics data

• Started in 2001
• Gained popularity through microarray analysis packages
• Colloborative effort by developers across the world
• Distributed as R packages
• Two releases every year

Objectives

• Provide access to powerful statistical and graphical methods for analysing genomics data.
• Facilitate retrieval and integration of annotation data (GenBank, GO, Entrez Gene, PubMed).
• Allow the rapid development of extensible, interoperable, and scalable software.
• Promote high-quality documentation and reproducible research.
• Provide training in computational and statistical methods.

www.bioconductor.org

BioC webpage

Bioconductor Release 3.14 - it works with R version 4.1.0

• Software (2083)

  • Provides implementation of analysis methods

• AnnotationData (904)

  • mapping between microarray probe, gene, pathway, gene ontology, homology and other annotations
  • Representations of GO, KEGG and other annotations, and can easily access NCBI, Biomart, UCSC and other sources

• ExperimentData (408)

  • code, data and documentation for specific experiments or projects

• Workflows (29)

  • Common analysis work flows for genomics data

Installating BioC Packages

Installing & Using a package for R version < 3.5.0:

source("https://bioconductor.org/biocLite.R")
BiocInstaller::biocLite(c("GenomicFeatures", "AnnotationDbi"))

Installing & Using a package for R version after 3.5.0:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("GenomicFeatures", "AnnotationDbi"))

Finding Help

• ?function

library("limma")
?lmFit

• Browse Vignettes browseVignettes(package=“limma”)
• Bioconductor Course Materials: http://www.bioconductor.org/help/course-materials/
• Bioconductor Support Forum - https://support.bioconductor.org/

Prints version information about R and all loaded packages

sessionInfo()

GenomicRanges

• Genomic Ranges provides data structure for efficiently storing genomic coordinates
  • Collection of genes coordinates
  • Transcription factor binding sites (ChIP-Seq peaks)
  • Collection of aligned sequencing reads
• Builds on top of Interval Ranges (IRanges) package and lays foundation for sequencing analysis.
• IRanges are collection of integer interval and GenomicRanges extends IRanges by including chromosome and strand.
• Provides collection of functions for accessing and manipulating Genomic coordinates
• Use cases: Identifying TF binding overlap, counting sequencing reads overlap with a gene
• Main classes: GRanges and GRangesList

Run Length Encoding (Rle)

• Run length encoding is a data compression technique
• Efficiently encoding the redundant information

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GenomicRanges")

## Warning: package(s) not installed when version(s) same as current; use `force = TRUE` to
##   re-install: 'GenomicRanges'

# Orginial vector
# chr1, chr2, chr2, chr2, chr1, chr1, chr3, chr3

library(GenomicRanges)
chr <- Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 2))
chr

## character-Rle of length 8 with 4 runs
##   Lengths:      1      3      2      2
##   Values : "chr1" "chr2" "chr1" "chr3"

The above Rle can be interpreted as a run of length 1 of chr1, followed by run length of 3 of chr2, followed by run length of 2 of chr1 and followed by run length of 2 of chr3

Constructing GRanges object from data frame

mm9genes <- read.table("./data/mm9Genes.txt",sep="\t",header=T)
head(mm9genes)

##     chr     start       end strand                ens        Symbol
## 1  chr9 105729415 105731415      - ENSMUSG00000043719        Col6a6
## 2 chr18  43636450  43638450      - ENSMUSG00000043424        Gm9781
## 3  chr7  70356455  70358455      - ENSMUSG00000030525        Chrna7
## 4  chrX 100818093 100820093      - ENSMUSG00000086370 B230206F22Rik
## 5 chr12  82881157  82883157      - ENSMUSG00000042724        Map3k9
## 6  chr7 152124774 152126774      - ENSMUSG00000070348         Ccnd1

mm9genes.GR <- GRanges(seqnames=mm9genes$chr,
                       ranges=IRanges(start=mm9genes$start,end=mm9genes$end),
                       strand=mm9genes$strand,
                       ENSID=mm9genes$ens,
                       Symbol=mm9genes$Symbol)

Another option: makeGRangesFromDataFrame()

Converting GRanges object to data frame: as.data.frame(GRanges)

Constructing GRanges object from data frame

head(mm9genes.GR)

## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames              ranges strand |              ENSID        Symbol
##          <Rle>           <IRanges>  <Rle> |        <character>   <character>
##   [1]     chr9 105729415-105731415      - | ENSMUSG00000043719        Col6a6
##   [2]    chr18   43636450-43638450      - | ENSMUSG00000043424        Gm9781
##   [3]     chr7   70356455-70358455      - | ENSMUSG00000030525        Chrna7
##   [4]     chrX 100818093-100820093      - | ENSMUSG00000086370 B230206F22Rik
##   [5]    chr12   82881157-82883157      - | ENSMUSG00000042724        Map3k9
##   [6]     chr7 152124774-152126774      - | ENSMUSG00000070348         Ccnd1
##   -------
##   seqinfo: 21 sequences from an unspecified genome; no seqlengths

GenomicRangesList

• To represent hierarchical structured data, ex: Exons in a transcript
• List-like data structure
• Each element of the list is GRanges instance

gr2 <- GRanges(seqnames = Rle(c("chr1", "chr3","chr2", "chr1", "chr3"), c(1, 2,1, 2, 4)),
                ranges = IRanges(start=55:64, end = 94:103, names = letters[11:20]),
                strand = Rle(c("+", "-", "+", "-"), c(1, 4, 3, 2)),
                score = 1:10, GC = runif(10,0,1))

GRL <- GRangesList("Peak1" = gr1, "Peak2" = gr2)

Operations on GenomicRanges

Operations on GenomicRanges

Operations on GenomicRanges

Operations on GenomicRanges

head(ranges(gr1))

## IRanges object with 6 ranges and 0 metadata columns:
##         start       end     width
##     <integer> <integer> <integer>
##   a        11        50        40
##   b        12        51        40
##   c        13        52        40
##   d        14        53        40
##   e        15        54        40
##   f        16        55        40

start(gr1)

##  [1] 11 12 13 14 15 16 17 18 19 20

width(gr1)

##  [1] 40 40 40 40 40 40 40 40 40 40

Finding overlapping regions

• One of the common tasks in Sequencing data analysis
• Ex: Identifying transcription factor binding sites overlap with promoters
• findOverlaps function finds intervals overlap between two GRanges object.
• Usage: function(query,subject)

gr1_overlaps <- findOverlaps(gr1,gr2,ignore.strand=F)
gr1_overlaps

## Hits object with 5 hits and 0 metadata columns:
##       queryHits subjectHits
##       <integer>   <integer>
##   [1]         6           1
##   [2]         9           2
##   [3]         9           3
##   [4]        10           2
##   [5]        10           3
##   -------
##   queryLength: 10 / subjectLength: 10

Output of findOverlaps is a ‘Hits’ object indicating which of the query and subject intervals overlap.

Finding overlapping regions

Convert the ‘Hits’ object to 2 column matrix using as.matrix(). Values in the first column are indices of the query and values in second column are indices of the subject.

gr1_overlaps.m <- as.matrix(gr1_overlaps)
gr1[gr1_overlaps.m[,"queryHits"], ]

## GRanges object with 5 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score         GC
##        <Rle> <IRanges>  <Rle> | <integer>  <numeric>
##   f     chr1     16-55      + |         6 0.00232794
##   i     chr3     19-58      - |         9 0.60866204
##   i     chr3     19-58      - |         9 0.60866204
##   j     chr3     20-59      - |        10 0.75465976
##   j     chr3     20-59      - |        10 0.75465976
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

Finding overlapping regions

subsetByOverlaps extracts the query intervals overlap with subject intervals

subsetByOverlaps(gr1,gr2,ignore.strand=F)

## GRanges object with 3 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score         GC
##        <Rle> <IRanges>  <Rle> | <integer>  <numeric>
##   f     chr1     16-55      + |         6 0.00232794
##   i     chr3     19-58      - |         9 0.60866204
##   j     chr3     20-59      - |        10 0.75465976
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

Finding overlapping regions - %in% and %over%

Alternate ways of find overlapping regions

gr1[gr1 %in% gr2]

## GRanges object with 0 ranges and 2 metadata columns:
##    seqnames    ranges strand |     score        GC
##       <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

gr1[gr1 %over% gr2]

## GRanges object with 3 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score         GC
##        <Rle> <IRanges>  <Rle> | <integer>  <numeric>
##   f     chr1     16-55      + |         6 0.00232794
##   i     chr3     19-58      - |         9 0.60866204
##   j     chr3     20-59      - |        10 0.75465976
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths

Note the difference between %in% and %over%

Finding overlapping regions

Other interesting functions: nearest() and distanceToNearest()

distanceToNearest(gr1,gr2,ignore.strand=F)

## Hits object with 8 hits and 1 metadata column:
##       queryHits subjectHits |  distance
##       <integer>   <integer> | <integer>
##   [1]         1           5 |         8
##   [2]         4           4 |         4
##   [3]         5           5 |         4
##   [4]         6           1 |         0
##   [5]         7           7 |         4
##   [6]         8           7 |         3
##   [7]         9           3 |         0
##   [8]        10           3 |         0
##   -------
##   queryLength: 10 / subjectLength: 10

precede, follow

Find nearest range in gr2 that precede or follow each range in gr1

precede(gr1,gr2)

##  [1] NA NA NA NA NA  6  7  7 NA NA

follow(gr1,gr2)

##  [1]  5 NA NA  4  5 NA NA NA  9  9

precede, follow returns the index of previous/next ranges. Overlapping ranges are excluded

Counting overlapping regions

countOverlaps() tabulates number of subject intervals overlap with each interval in query, ex: counting number of sequencing reads overlap with genes in RNA-Seq

countOverlaps(gr1,gr2,ignore.strand=T) # note the strand!

## a b c d e f g h i j 
## 0 0 0 0 0 1 1 2 2 2

Computing Coverage

coverage calculates how many ranges overlap with individual positions in the genome. coverage function returns the coverage as Rle instance.

coverage(gr1)

## RleList of length 3
## $chr1
## integer-Rle of length 55 with 6 runs
##   Lengths: 10  4  1 35  4  1
##   Values :  0  1  2  3  2  1
## 
## $chr2
## integer-Rle of length 53 with 6 runs
##   Lengths: 11  1  1 38  1  1
##   Values :  0  1  2  3  2  1
## 
## $chr3
## integer-Rle of length 59 with 8 runs
##   Lengths: 16  1  1  1 37  1  1  1
##   Values :  0  1  2  3  4  3  2  1

Time for Exercises!

• Exercises Part1
  
  

• Exercises Part1 Solutions

File formats in NGS

• FASTQ (Fasta with quality)
• SAM/BAM - Sequence Alignment/Map format
• BED
• Wiggle/bedgraph

FASTQ

@ERR590398.1 HWI-ST1146:148:C2FVTACXX:8:1101:1172:2059
NAAAATGCATATTCCTAGCATACTTCCCAAACATACTGAATTATAATCTC
+ERR590398.1 HWI-ST1146:148:C2FVTACXX:8:1101:1172:2059
A1BDFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJIJJ

BAM/SAM

BAM/SAM consist two sections: Header and Alignment

Header section:

Meta data (reference genome, aligner), starts with “@”

Reading Sequence alignments (BAM/SAM)

Methods for reading BAM/SAM

• readAligned from ShortRead package – Accept multiple formats – BAM, export
  • Reads all files in a directory
  • Reads base call qualities, chromosome, position, and strand
• scanBam from Rsamtools package
  • scanBam reads BAM files into list structure
  • Options to select what fields and which records to import using ScanBamParam
• readGAlignments from GenomicAlignments package

Reading Sequence alignments (BAM/SAM)

library("Rsamtools")
BamFile <- system.file("extdata", "ex1.bam", package="Rsamtools")
which <- IRangesList(seq1=IRanges(1000, 2000),seq2=IRanges(c(100, 1000), c(1000, 2000)))
what <- c("rname", "strand", "pos", "qwidth", "seq")
param <- ScanBamParam(which=which, what=what)
bamReads <-  scanBam(BamFile, param=param)

length(bamReads) # Each element of the list corresponds to a range specified by the which argument in ScanBamParam

## [1] 3

names(bamReads[[1]]) # elements specified by the what and tag arguments to ScanBamParam

## [1] "rname"  "strand" "pos"    "qwidth" "seq"

Reading Sequence alignments (BAM/SAM)

ScanBamParam:

# Constructor
ScanBamParam(flag = scanBamFlag(), what = character(0), which)

# Constructor helpers
scanBamFlag(isPaired = NA, isProperPair = NA, isUnmappedQuery = NA, 
    hasUnmappedMate = NA, isMinusStrand = NA, isMateMinusStrand = NA,
    isFirstMateRead = NA, isSecondMateRead = NA, isSecondaryAlignment = NA,
    isNotPassingQualityControls = NA, isDuplicate = NA)

Reading Sequence alignments (BAM/SAM)

• GenomicRanges package defines the GAlignments class – a specialised class for storing set of genomic alignments (ex: sequencing data)
• Only BAM support now – future version may include other formats
• The readGAlignments function takes an additional argument, param allowing the user to customise which genomic regions and which fields to read from BAM -param can be constructed using ScanBamParam function

library(GenomicAlignments)
SampleAlign <- readGAlignments(BamFile)

Reading Sequence alignments (BAM/SAM)

We can also customise which regions to read

region <- IRangesList(seq1=IRanges(1000, 2000),seq2=IRanges(1000, 2000))
param1 <- ScanBamParam(what=c("rname", "pos", "cigar","qwidth"),which=region)
SampleAlign1 <- readGAlignments(BamFile,param=param1)

Time for Exercises!

• Exercises Part2
  
  

• Exercises Part2 Solutions

AnnotationDbi Accessor Functions

• columns What kind of annotation available in AnnotationDb object.
• keytypes Displays which type of identifiers can be passed in to select function.
• keys returns keys (index) for the database contained in the AnnotationDb object. Used along with keytypes in select function to retrieve interested annotation
• select will retrieve the annotation data as a data.frame based on the supplied keys, keytypes and columns.

Note: Package name = Annotation object

We will explore how to retrieve annotation from gene-centric organism level annotation package (org.Hs.eg.db)

TranscriptDB (TxDb) packages

• TxDb packages provide access genomic coordinates to various transcript-related features from UCSC and Biomart data sources.
• TxDb objects contains relationship between mRNA transcripts, exons, CDS and their associated identifiers
• TxDb packages follows specific naming scheme, ex: TxDb.Mmusculus.UCSC.mm9.knownGene

We will explore TxDb package for Mouse mm9 genome from UCSC. We will first install the package and load in to our current working space.

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("TxDb.Mmusculus.UCSC.mm9.knownGene")

TxDb.Mmusculus.UCSC.mm9.knownGene

By default, the annotation object will have same name as package name. Create an alias for convenience.

library("TxDb.Mmusculus.UCSC.mm9.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene

Since TxDb are inherited from AnnotationDb object, we can use columns, keys, select, and keytypes functions.

keys <- c("100009600", "100009609", "100009614")
select(txdb,keys=keys,columns=c("GENEID","TXNAME"),keytype="GENEID")

## 'select()' returned 1:1 mapping between keys and columns

##      GENEID     TXNAME
## 1 100009600 uc009veu.1
## 2 100009609 uc012fog.1
## 3 100009614 uc011xhj.1

biomaRt - Example

library("biomaRt")
marts <- listMarts() # List available marts
marts[1,]

##                biomart           version
## 1 ENSEMBL_MART_ENSEMBL Ensembl Genes 106

ensembl <- useMart("ensembl") # select ensembl
ens_datasets <- listDatasets(ensembl) # list datasets 
ens_human <- useDataset("hsapiens_gene_ensembl",mart=ensembl) # select human dataset
ens_human_Attr <- listAttributes(ens_human) # list available annotation
ens_human_filters <- listFilters(ens_human) # list availabel filters

Extract genomic coordinates, ensembl gene identifier and gene symbol for genes in chromsome X

chrXGenes <- getBM(attributes = c("chromosome_name","start_position","end_position","ensembl_gene_id","strand","external_gene_name"), filter="chromosome_name",values="X", mart=ens_human)

Using rtracklayer for creating bigwig format files

library(rtracklayer)

myRleList <- RleList(chr1=Rle(c(0,0,0,0,0,1,1,1,0,0,0,0,0)),chr2=Rle(c(0,0,0,0,0,1,1,1,2,2,2,2,2)))

 export.bw(myRleList,con="chr1_Myc.bw")

Time for Exercises!

• Exercises Part3
  
  

• Exercises Part3 Solutions