---
title: "RNA-seq Analysis 1"
author: "Pietà Schofield"
subtitle: Introduction to Bioconductor Packages
output:
  html_document:
    code_folding: show
    fig_caption: yes
    highlight: textmate
    theme: spacelab
    toc: yes
    toc_depth: 3
    toc_float:
      collapsed: no
      smooth_scroll: yes
  pdf_document:
    toc: yes
    toc_depth: '3'
---


```{r setup,echo=FALSE,include=FALSE}
require(RefManageR)
BibOptions(style="html",longnamesfirst=F)
if(!file.exists(".pieta.bib")){
 download.file("http://www.compbio.dundee.ac.uk/user/pschofield/biblio.bib",".pieta.bib")
}
bib <- ReadBib(".pieta.bib",check=FALSE)
require(knitr)
opts_chunk$set( message=F, warning=F, comment=NA)
#
# This little bit of trickery is rather complex knitr-foo thanks to Ramnath Vaidyanathan
#
# http://ramnathv.github.io/posts/verbatim-chunks-knitr/index.html
#
# It enables code chunks to be included with the surrounding brackets and options
# for illustrative purposes
#
knitr::opts_chunk$set(tidy = F)
knit_hooks$set(source = function(x, options){
  if (!is.null(options$verbatim) && options$verbatim){
    opts = gsub(",\\s*verbatim\\s*=\\s*TRUE\\s*", "", options$params.src)
    bef = sprintf('\n\n    ```{r %s}\n', opts, "\n")
    stringr::str_c(
      bef, 
      knitr:::indent_block(paste(x, collapse = '\n'), "    "), 
      "\n    ```\n"
    )
  } else {
    stringr::str_c("\n\n```", tolower(options$engine), "\n", 
      paste(x, collapse = '\n'), "\n```\n\n"
    )
  }
})
#
# a little bit of R from Yihui Xieto put in verbatim inline R
#
# http://stackoverflow.com/questions/20409172/how-to-display-verbatim-inline-r-code-with-backticks-using-rmarkdown
#
rinline <- function(code) {
  sprintf('``` `r %s` ```', code)
}

urlStub <- "http://www.compbio.dundee.ac.uk/user/pschofield/Projects/teaching_pg/"
```

- [Return to the  main page](http://www.compbio.dundee.ac.uk/user/pschofield/Projects/teaching_pg)
- [Full Rmarkdown script for this workshop](http://www.compbio.dundee.ac.uk/user/pschofield/Projects/teaching_pg/code/biocClasses.Rmd)


#Bioconductor

[Bioconductor](https://www.bioconductor.org) 
`r Cite(bib,"huber2015",.opts=list(max.names=3))` is a curated repository of open source packages for the R statistical environment specifically aimed at biological information processing. There are three main strands of packages

- Software packages that contain libraries of functions and data object classes for the processing and statistical analysis of biological data
- Annotation packages that contain data files that hold curated biological information databases for example lists of genes in a species, or the DNA sequence of an organism.
- Experimental Data packages that contain the results files in various states of processing and analysis form biological experiment.

The project has a website which is well worth exploring to discover the range of available packages.

![Bioconductor Home Page](`r paste0(urlStub,"/figures/biocHP.png")`)

## Installation

If you are installing R RStudio and bioconductor on a computer the initial installation of bioconductor requires the installation of a script from the bioconductor website. 

```{r bioc3}
source("http://bioconductor.org/biocLite.R")
```

Once this script has been installed it can be run to install the basic bioconductor installation package

```{r bioc4, eval=FALSE}
biocLite("BiocInstaller")
```

Subsequently you can install packages with the BiocInstaller packages without the sourcing of the script

```{r bioc5, eval=FALSE}
# with either
BiocInstaller::biocLite("NameOfPackage")

# or load the installer package first

require(BiocInstaller)
biocLite("NameOfPackage")
```

**However if you want a bioconductor package on the RStudio Server which is not installed it is probably better to request SLS Computing install the package centrally for you and everyone else can also then access it. The same is probably true of CRAN packages. This has the added advantage that you are less likely to get incompatibilities between packages which can arise if you install a package that is incompatible with a version of another package installed globally**

## Core Data Classes

One of the most useful things Bioconductor provides is a collection of common data classes that have storage structures and processing methods that as used extensively throughout the library. 

### Biostrings

The [Biostrings]() package 
`r Cite(bib,"pages2016",.opts=list(max.names=3))` provides several classes for the storage of sequences of single characters. It provides a general virtual class [XString]() for storing big string sequences of any character. It provides derived classes [BString], [DNAString], [RNAString] and [AAString] for general strings, DNA sequences, RNA sequences and Amino Acid sequences such as peptides and proteins. I also provides a class for collections or sets of such objects derived from the [XStringSet] class. It also provides classes for holding multiple sequence alignments and masked and quality scaled sequences and sets of aligned sequences.

Also provided are several standard constants such as `DNA_ALPHABET` with the full set of [IUPAC]() letters for DNA, and `DNA_BASES` that contains just the unambiguous base alphebet.

```{r dna}
require(Biostrings)
DNA_ALPHABET
DNA_BASES
```

So for example we can generate a vector of 10 bases chosen at random from the `DNA_BASES` alphabet  using the `sample()` function

```{r randomDNA}
# sample the alphabet
vcDNA <- sample(DNA_BASES,10,replace=T)
str(vcDNA)
```

The constructor for an object of the DNAString class wants a single string of chararters not a vector of single characters so we can use the `paste()` function to turn the vector of bases to a single string

```{r paste}
stDNA <- paste(vcDNA,collapse="")
str(stDNA)
```

We can then use the `DNAString()` constructor function to make a DNAString object.

```{r makeDNAString}
DNA <- DNAString(stDNA)
str(DNA)
```

The `str()` function here reveals the structure of the DNA objects you can see it is of class `DNAString` and you can see the name of the package that defined this class. It also has 5 `slots`, there are the objcts that appear in the indented list starting with an `@` sign.

You can retrieve the content of these slots with the syntax _objectname_ @ _slotname_ .

```{r slots}
DNA@shared
```

In the case of XString objects you can also retrieve the sequence as a character string using the casting function `as.character()`

```{r casrting}
seq <- as.character(DNA)
seq
```

As well as the `str()` function used above we can use the `is()` function to see what type of thing R knows an object to be.

```{r showTypes}
is(DNA)
```

Biostrings provides methods for reading and writing common file formats containing sequences searching such as FASTA and FASTQ, for searching with pattern and position weight matrices, and for pairwise sequence alignments.

```{r pattermatch}
bigSeq <- DNAString("AGCTGGTCGATGACGTTAGGATCAGTACGATTGAGGCAGGGTCAGTAATTGGATTGAGGATTGACCCCGATG")
smallSeq <- DNAString("TCAGTA")
pos <- matchPattern(smallSeq,bigSeq)
pos
```

### Execise 1.2

- What kind of object is the pos object above (the result of the `patternMatch()` call)
    - use the `str()` function to view its structure.
    - recover content of the ranges slot
    - use the `is()` function discover the full list of object types and classes the ranges slot 
    belongs to

### IRanges

```{r echo=FALSE}
require(GenomicAlignments)
plotRanges <- function(ir){
  bins <- disjointBins(IRanges(start(ir), end(ir) + 1))
  dat <- cbind(as.data.frame(ir), bin = bins, seqid = seq(length(bins)))
  ggplot2::ggplot(dat) + 
    ggplot2::geom_rect(ggplot2::aes(xmin = start, xmax = end, ymin = bin/4, ymax = bin/4 + 0.25,
                                    fill = as.factor(seqid)) )+
   ggplot2::scale_fill_discrete(guide=FALSE)
}
```

The IRanges package `r Cite(bib,"lawrence2013",.opts=list(max.names=3))` provides a way of storing, manipulating and calculating with interval ranges. Interval ranges represent regions of a sequence. 

You can construct an _IRange_ object with the `IRanges()` constructor, by specifying either of two lists, 

- start and end positions 
- start positions and widths

```{r ranges, fig.height=3}
library(IRanges)   
ir <- IRanges(start = c(1, 8, 14, 15, 19, 34, 40),
              width = c(12, 6, 6, 15, 6, 2, 7))
plotRanges(ir)
```

The IRanges package provides functions for manipulating and analysing IRanges objects, for example there is a function `overLapsAny()` that will tell you which ranges in an IRanges object overlap

```{r overlaps}
overLaps <- overlapsAny(ir)
overLaps
```

and exactly which the overlap with

```{r findOver}
overLaps <- findOverlaps(ir)
overLaps
```

We can split the ranges in to groups (or bins) of non-overlapping (or disjoint) ranges

```{r disjoint}
bins <- disjointBins(ir)
bins
```

```{r niceplot, echo=FALSE, fig.height=3}
dat <- cbind(as.data.frame(ir), bin = bins)
ggplot2::ggplot(dat) + 
  ggplot2::geom_rect(ggplot2::aes(xmin = start, xmax = end, ymin = bin/4, ymax = bin/4 + 0.25,
                                    fill = as.factor(bin)) ) +
  ggplot2::scale_fill_discrete(name="bin")
```

### GenomicRanges

The GenomicRanges package `r Cite(bib,"lawrence2013",.opts=list(max.names=3))` makes use of a compression technique called [Run Length Encoding(RLE)]() to store extra information with interval ranges to hold data such as genomic locations 

#### Run Length Encoding

```{r}
# generate a repetative sequence
mySeq <- unlist(sapply(rnorm(10,5,1),
                  function(n){
                    rep(sample(DNA_BASES,1),n)
                  }))
# what does it look like
mySeq
# how big is it
paste0("Size of mySeq ",format(object.size(mySeq),units = "Kb"))
# turn it into an RLE object
myRLE <- Rle(mySeq)
# what does that look like
myRLE
# how big is that
paste0("Size of myRLE ",format(object.size(myRLE),units = "Kb"))
```

Hmmmm! not much advantage "I thought you said it was a compression technique", well if you a looking at longer sequences with more repeats

```{r}
# generate a much bigger repetative sequence
mySeq <- unlist(sapply(rnorm(100,10,2),
                  function(n){
                    rep(sample(DNA_BASES,1),n)
                  }))
paste0("Size of mySeq ",format(object.size(mySeq),units = "Kb"))
# and now the RLE
myRLE <- Rle(mySeq)
paste0("Size of myRLE ",format(object.size(myRLE),units = "Kb"))
```

GenomicRanges (also called GRanges) are excelent for genomic annotations. The minimal requirements for a GRanges object are an RLE object with sequences names and a set of ranges associated with the seqences in the form of an IRanges object. The GRanges object also contains a slot for strand (sense on the DNA), however this defaults to `*` which represents unspecified if not specifically given.

```{r grangesConstructor}
grObj <- GRanges(seqnames=Rle(c(rep("seq1",4),rep("seq2",2),"seq3")),
                 ranges = IRanges(start = c(110, 200 , 250, 350, 100, 300, 40),
                                  width = c(50, 40, 60, 150, 60, 200, 70)))
grObj
```

GRanges objects can store other information along with this basic data, typically they are used for information in GTF and BED annotation files. The package `rtracklayer` has many functions for importing and exporting data from these file formats

```{r rtracklayer}
require(rtracklayer)
track <- import(system.file("tests", "v1.gff", package = "rtracklayer"))
export(track, "my.gff3", "gff3")
track
```

R has a few packages that are designed for visualiation of genomic range data. Two of the more flexible are the package `Gviz` and the package `ggbio`. 

For example we can plot use Gviz to plot some mocked up annotation 

```{r gvisexample, fig.height=3}
require(Gviz)
# if we mock up some data
grObj <- GRanges(seqnames=Rle(rep("7",4)),
                 ranges=IRanges(start=c(2000000, 2070000, 2100000, 2160000),
                                end=c(2050000, 2130000, 2150000, 2170000)),
                 strand=c("-", "+", "-", "-"),
                 group=c("Group1","Group2","Group1", "Group3"))
# make and annotation track
annTrack <- AnnotationTrack(grObj, genome="hg19", feature="test",
                            id=paste("annTrack item", 1:4),
                            name="annotation track",
                            stacking="squish")
# make and axis track
axis <- GenomeAxisTrack()
# and an ideogram (chromosome map track)
ideo <- IdeogramTrack(genome="hg19",chromosome=7)
## Now plot the tracks
res <- plotTracks(list(ideo, axis, annTrack))
```

### GenomicAlignments

There are several packages in bioconductor that are specifically designed for the handling of short read data, in the form of FASTQ raw files and alignment files such as BAM and SAM files. There are differences in the underlying data structures they use and the functions the provide.

The `ShortRead` package `r Cite(bib,"morgan2009",.opts=list(max.names=3))` is particularly suited for processing reads in FASTQ files and provides functions for sampling and QC tasks, it imports the reads into `Biostring` objects. It also contains some functions for single ended alignments. However, the packages `GenomicAlignments` `r Cite(bib,"lawrence2013",.opts=list(max.names=3))` provides much more extensive support for paired end and gapped alignments in BAM files. It uses the package `Rsamtools` to import and export alignment data files.

When alignments are read in the reads are held in GRanges objects with all the functions of the GenomicRanges and IRanges packages available for processing.

```{r genomicAlignments}
alignFile <- system.file(package="Gviz", "extdata", "gapped.bam")
gaData <- readGAlignmentPairs(alignFile)
gaData
```

`Gvis` also has functions to make plotting of alignment data relatively painless.

```{r genomicAlignmentPlot, fig.height=3}
atrack <- AlignmentsTrack(alignFile,isPaired=T)
plotTracks(atrack,from=2960000,to=3160000,chromosome="chr12")
```

## Data Import Export

Bioinformatics is settling on a fairly steady set of file formats, FASTA and FASTQ for sequence data for example and SAM/BAM for NGS short read alignment data. While many of packages mentioned in this workshop possess functions for importing and exporting data that the package data structures are designed to handle they often don't possess functions for other files. There are also some standard formats for processed data for example read depths, peak locations and genetic annotations. This type of data often appears as track in a genome browser.

### rtracklayer

The package [rtracklayer](https://bioconductor.org/packages/release/bioc/html/rtracklayer.html) 
`r Cite(bib, "lawrence2009")` has two useful generic functions `import()` and `export()` that will handle several track formats including FASTA, GTF, BED, WIG. 

## Annotations

There are many packages in bioconductor that contain curated annotation data for numerous and diverse organisms. See the list on the [BiocViews AnnotationData](https://www.bioconductor.org/packages/release/BiocViews.html#___AnnotationData) page

## Annotation Packages

There are annotation datasets that are available as Bioconductor packages many of them use the `RSQLite` package to interface to `sqlite` database files that hold the information locally. Bioconductor provides some standard packages for a range of annotation databases. Examples of these are:

- [BSGenome](https://bioconductor.org/packages/release/BiocViews.html#___BSgenome) packages `r Cite(bib,"pages2016", .opts=list(max.names=3))` store the full genomic sequence in Biostrings format for various organisms. 
- `AnnotationDbi` packages `r Cite(bib,"pages2016a", .opts=list(max.names=3))` that store several type of information
    - [OrgDb](https://bioconductor.org/packages/release/BiocViews.html#___OrgDb) packages  store genome wide annotation data in a  format for many organisms. 
    - [ChipDb](https://bioconductor.org/packages/release/BiocViews.html#___ChipDb`) and `probe` annotations that hold information on probe-sets on various micro array platforms
    - [InparanoidDb](https://bioconductor.org/packages/release/BiocViews.html#___InparanoidDb) packages that hold homology information on proteins
- [Txdb](https://bioconductor.org/packages/release/BiocViews.html#___TxDb) uses the [GenomicFeatures](https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html) package `r Cite(bib,"lawrence2013", .opts=list(max.names=3))` to manage database that hold transcript model annotations.
- [MeSHDb](https://bioconductor.org/packages/release/BiocViews.html#___MeSHDb) `r Cite(bib,"tsuyuzaki2015")` provide links between genes and the Medical Subject Headings curated database for many organisms


## Online Packages

Bioconductor also provides some library packages that are designed to ease searching various on-line databases for annotation data. These can be used to programmatically retrieve data that is often also avaliable to be downloaded from web sites such as [ENSEMBL](http://www.ensembl.org/info/data/ftp/index.html) and [UCSC](http://hgdownload.soe.ucsc.edu/downloads.html)

- [BiomaRt](https://bioconductor.org/packages/release/bioc/html/biomaRt.html) 
Interface to BioMart databases  `r Cite(bib,c("durinck2009","durinck2005"))` (e.g. Ensembl, COSMIC ,Wormbase and Gramene)
- [AnnotationHub](https://bioconductor.org/packages/release/bioc/html/AnnotationHub.html) The home page for the AnnotationHub package states
"The AnnotationHub web resource provides a central location where genomic files (e.g., VCF, bed, wig) and other resources from standard locations (e.g., UCSC, Ensembl) can be discovered"
- [GEOQuery](http://bioconductor.org/packages/release/bioc/html/GEOquery.html) 
`r Cite(bib,"davis2007")` This package provides access to the NCBI Gene Expression Omnibus (GEO) is a public repository of micro array data
- [SRAdb](http://bioconductor.org/packages/release/bioc/html/SRAdb.html) 
`r Cite(bib, "zhu2013")` provides similar access to GEOQuery except to the NCBI Small Reads Archive, the largest public repository of sequencing data from the next generation of sequencing platforms.
- [STRINGdb](https://www.bioconductor.org/packages/release/bioc/html/STRINGdb.html)
`r Cite(bib, "franceschini2013")` provides an interface to the STRING protein-protein interactions database. It also gives access to KEGG and GO annotations.
- [pathview](http://bioconductor.org/packages/release/bioc/html/pathview.html) 
`r Cite(bib, "luo2013")` is a tool set for pathway based data integration and visualization and will download pathway information from KEGG using the [KEGGgraph](http://bioconductor.org/packages/release/bioc/html/KEGGgraph.html) 
`r Cite(bib, "zhang2015")` package


###Biomart Example

Provided you are not working in bacteria the [biomaRt] package provides a useful interface onto the [ENSEMBL](http://www.ensembl.org/index.html) databases. [Biomart](http://www.biomart.org) provides a standardised data access protocols and ENSEMBL have implemented it on their genome servers for

- [Metazoa](http://metazoa.ensembl.org/index.html)
- [Protist](http://protists.ensembl.org/index.html)
- [Fungi](http://fungi.ensembl.org/index.html)
- [Plants](http://plants.ensembl.org/index.html)

as well as their popular genomes server, but not [bacteria](http://bacteria.ensembl.org/index.html)

###AnnotationHub Example

We know there is some _Saccharomyces cerevisiae_ data coming our way so lets look to see what is available through AnnotationHub and lets look at the data provider ENSEMBL.

```{r ah, eval=F}
require(AnnotationHub)
ah <- AnnotationHub()
sc <- query(ah,c("ENSEMBL","Saccharomyces cerevisiae"))
sc
```

This displays a list of the data items available. A fuller list in the form of a data frame is available with the `mocls()` function, this displays more information including the URL of the on-line data source.

```{r ahmcols, eval=FALSE}
mcols(sc)
```

Item `AH49366` is the full genomic sequence for _Saccharomyces cerevisiae_ release 81 

```{r ahGenome, eval=FALSE}
scFaFile <- sc[["AH49366"]]
scFaFile
```

The sequences in the FaFile object can be retrieved with the `getSeq()` command

```{r getSeq, eval=FALSE}
getSeq(scFaFile)
```

The last item in the list of available _S. cerevisiae_ annotation is the GTF file containing the genetic features and we can retrieve this into a `GenomicRanges` object

```{r ahGTF, eval=FALSE}
genes <- sc[["AH50407"]]
head(genes)
```

It is possible to retrieve the sequences for a gene of interest using these two annotations

```{r getGeneSeq, eval=FALSE}
# get the GenomicRanges item that has the gene name of interest and is of type gene 
ser3 <- genes[which(genes$gene_name=="SER3" & genes$type=="gene"),]
# pass this as a parameter to the getSeq function
seqSer3 <- getSeq(scFaFile,ser3)
names(seqSer3) <- ser3$gene_id
seqSer3
```

## Exercise 2.1

- Which chromosome is the SER3 gene on.
- In the example above what other types genetic features are listed for the SER3 gene?
- How many different bio-types are there in the genetic features data and what are they?
- Find the name, length and sequence for the longest snoRNA.

# Learning Bioconductor

Bioconductor provides many resources for assisting in learning how to do various analyses. It has a very active [support forum](https://support.bioconductor.org) where many bioconductor uses and many of the package authors will answer questions regarding the packages. 

## Package Vignettes

All Bioconductor packages are supposed to come with a user manual which is a concise list of the functions and data structures available in the package listed alphabetically and with all the function parameters documented. They are also expected to have examples for all the functions and at least one vignette document that is meant to be a more recipe, How-to, document. You can list the names of the available vignettes with the command `vignette()`

```{r vignette}
vignette(package="Biostrings")
```

You can then view the vignette with the same command only this time passing the name of the vignette you want.

```{r openVignette, eval=FALSE}
vignette("BiostringsQuickOverview")
```

## Workflows

One of the ways this is done is via [workflows](https://www.bioconductor.org/help/workflows/)

For example the [RNA-seq Workflow](http://www.bioconductor.org/help/workflows/rnaseqGene/) 
introduces several alternative packages for performing RNA-seq. We have concentrated on the
[edgeR]() package whereas the workflow concentrates on the [DESeq2]()

## Tutorials

The support forum has a specific section of links to [tutorial](https://support.bioconductor.org/t/Tutorials/)

## CRAN

Not all R packages for handling biological data are in bioconductor, some are in the [Comprehensive R Archive Network (CRAN)](https://cran.r-project.org/web/packages) and can be installed using the `install.packages()` function. Several of the [CRAN Task Views](https://cran.r-project.org/web/views/) are of interest to biological data analysis most specifically the [Phylogenetic Task View](https://cran.r-project.org/web/views/Phylogenetics.html) that list available packages on CRAN for phylogenetic studies. 

These are installed through the base packages `install.packages()` function. 

```{r cran, eval=FALSE}
install.packages("packagename")
```

- [Analyses of Phylogenetics and Evolution (ape) ](https://www.bioconductor.org/packages/release/bioc/html/ape.html) `r Cite(bib,"paradis2004")` is a very useful package that has several functions for importing data formats and even accessing online information. For example 

    - the `read.GenBank()` function can be used to retrieve sequence information by GenBank accession number 
- [genoPlotR](https://cran.r-project.org/web/packages/genoPlotR/index.html) `r Cite(bib, "guy2010")` also has a functions for reading feature and sequence annotation data 
    - the `read_dna_seg_from_genbank()` function will read the features from a GenBank file

There are also numerous packages available on GitHub and the [devtools]() package has a useful 
funtion `install.git_hub()` for installing packages from GitHub

# Bibliography

```{r bibliography, echo=FALSE, results="asis"}
PrintBibliography(bib)
```

# Session Info

```{r session, echo=FALSE}
sessionInfo()
```



