---
title: "RNA-seq Analysis 1"
author: "Pietà Schofield"
subtitle: Introduction to Bioconductor
output:
  html_document:
    code_folding: show
    fig_caption: yes
    highlight: textmate
    theme: spacelab
    toc: yes
    toc_depth: 3
    toc_float:
      collapsed: yes
      smooth_scroll: no
  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)
}
```

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

# Rmarkdown

[Donald Knuth](https://en.wikipedia.org/wiki/Donald_Knuth) coined the phrase [literate programming](https://en.wikipedia.org/wiki/Literate_programming) `r Cite(bib,c("knuth1984","knuth1992"))` and he made a distinction

- **Literate proramming** embeds the code within the natural language description of the logic behind the code 

$$
versus
$$

- **Documentation generation** structured comments embedded in the code are extracted to produce documentation 

*"Let us change our traditional attitude to the construction of programs: Instead of imagining that our main task is to instruct a computer what to do, let us concentrate rather on explaining to human beings what we want a computer to do."* `r TextCite(bib,"knuth1984")`

This can be particularly useful way of looking at things when you are using a computer to perform data analysis. 

## RStudio 

RStudio is designed to make literate programming simple. It has a few options but the main two are provided by the packages `sweave` and `knitr`, which enable the writing of documents in a markup language such as [HTML](https://en.wikipedia.org/wiki/HTML) , [LaTeX](https://en.wikipedia.org/wiki/LaTeX) or [markdown](https://en.wikipedia.org/wiki/Markdown), with embedded code chunks, that can then be compiled into a final document in a range of formats such as portable document format (PDF), Microsoft Word (docx) or HTML. During the compilation process code chunks are executed and output such as figures are generated and incorporated into the final document.

In this course we will focus of [Rmarkdown](http://rmarkdown.rstudio.com) there is extensive documentation available as [http://rmarkdown.rstudio.com](http://rmarkdown.rstudio.com). You can start a new `rmarkdown` document from with-in RStudio from the **File - New File - R Markdown...** menu option that opens the dialogue box below

![R markdown dialogue box in RStudio](http://www.compbio.dundee.ac.uk/user/pschofield/Projects/teaching_BS32010/figures/rmarkdowndialogue.jpg)

Select Document option on the left hand side and the PDF option on the right hand side, fill in a title of your choice and your name and a  to create a demo rmarkdown document, that will open in the RStudio editor panel. Without making any changes to this file you can compile it into a PDF by clicking the **Knit PDF** button on the editor panel button bar. You will be prompted to save the file, so you should chose and name and location for the file.

![R markdown knit PDF button](http://www.compbio.dundee.ac.uk/user/pschofield/Projects/teaching_BS32010/figures/knitbutton.jpg)

If you examine the file in the editor you will see it is structured in a particular ways. 

### YAML header

In typical recursive geek speak [YAML](https://en.wikipedia.org/wiki/YAML) is an acronym that stands for **Y**AML **A**in't **M**arkup **L**anguage. It is a header block of parameters used through out the document compilation process. It uses a standard format for many programming and markup languages, though the permissable option vary. It is separated from the rest of the document being surrounded by three dashes `---`

For example the YAML block for this document is 

```
---
title: 'RNA-seq Analysis 1: Introduction to Bioconductor'
author: "Pietà Schofield"
output: 
  html_document:
    fig_caption: yes
    toc: yes
    toc_depth: 2
    toc_float:
      collapsed: no
      smooth_scroll: no
    code_folding: show
    css: workshop.css
---
```

### Code Chunks

Sections of R (and other language) code can be embedded in the document surrounded by the three back tick delimites as in the example below. 

```{r chuckName, fig.caption="Some Random Stuff", verbatim=TRUE, eval=FALSE}
#
# Plot anything R can plot
#
plot(1:10+rnorm(10),1:10+rnorm(10),pch="x", 
     xlab="expected", ylab="measured",main="Demo Plot")
#
# One option is to just send it to the default device and knitr
# captures it and put it in a temporary place
#
abline(0,1,col="red")
# 
# alternatively write it to a file and link to the file in the markdown
#
```

The first line holds paramters or **chunk options** between the parentheses that control how the chunk will be processed during compilation and how the output from the chuck if it is executed will be incorporated and formated in the document. For example whether the chunk is executed is controlled by the **eval** option and whether the chunk is displayed in the output document is controlled by the **echo** option

### Markdown Text

The rest of the document is in markdown and certain codes will produce various formatting styles in the output document. For example various numbers of # characters at the start of a line control heading levels * or _ will control italic and bold text. 


### Graphics

The code chunk above produces a graph in the output. It is possible to change the code in the chunk to change the figure. It is also possible to set chunk options  so that only the graph is displayed and not the code, or only the code and not the graph.

```{r plotgraph, fig.cap="Some Random Stuff", echo=FALSE}
#
# Plot anything R can plot
#
plot(1:10+rnorm(10),1:10+rnorm(10),pch="x", 
     xlab="expected", ylab="measured",main="Demo Plot")
#
# One option is to just send it to the default device and knitr
# captures it and put it in a temporary place
#
abline(0,1,col="red")
# 
# alternatively write it to a file and link to the file in the markdown
#
```

## Exercise 1.1

Create an rmarkdown document that will compile to PDF. 

- Alter the document to hide the code that displays the summary of the car dataset.
- Alter the document to add a table of contents.
- Add a code chunk to the document to use the `head()` function to display the first 10 items in the pressure data frame.
- Add a final section heading and code chunk using the `sessionInfo()` function to include the R configuration information in the document.

#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](http://www.compbio.dundee.ac.uk/user/pschofield/Projects/teaching_BS32010/figures/biocHP.png)

## Installation

The initial installation of bioconductor

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

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

Subsequently you can install packages with the BiocInstaller package

```{r bioc5, eval=FALSE}
require(BiocInstaller)
biocLite("NameOfPackage")
```

## 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. 


# 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/)

## 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. 

# Bibliography

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

# Session Info

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



