%\VignetteIndexEntry{Rsubread Vignette}
%\VignetteDepends{}
%\VignetteKeywords{Read mapping}
%\VignettePackage{Rsubread}
\documentclass[12pt]{article}

\usepackage{hyperref}

\newcommand{\R}[1]{{\texttt{#1}}}
\newcommand{\C}[1]{{\texttt{#1}}}

\textwidth=6.2in
\textheight=8.5in
\oddsidemargin=0.2in
\evensidemargin=0.2in
\headheight=0in
\headsep=0in

\begin{document}
\title{Rsubread package: high-performance read alignment, quantification and mutation discovery}
\author{Wei Shi}
\date{14 September 2015}
\maketitle

\section{Introduction}
This vignette provides a brief description to the Rsubread package.
For more details, please refer to the Users Guide which can brought up in your R session via the following commands: \\

\noindent\R{> library(Rsubread)}\\
\R{> RsubreadUsersGuide()}\\

The Rsubread package provides facilities for processing the read data generated by next-gen sequencing technologies.
These facilities include quality assessment, read alignment, read summarization, exon-exon junction detection, absolute expression calling and SNP discovery.
They can be used to analyze data generated from all major sequencing platforms including Illumina GA/HiSeq, Roche 454, ABI SOLiD and Ion Torrent.

The Subread aligner (\R{align} function) is a highly efficient and accurate aligner for mapping genomic DNA and RNA sequencing reads.
It adopts a novel mapping paradigm called ``seed-and-vote".
Under this paradigm, a number of 16mers (called seeds or subreads) are extracted from each read and they were mapped to the reference genome to vote for the mapping location of the read.
Read mapping performed under this paradigm has been found to be more efficient and accurate than that carried out under the conventional ``seed-and-extend" paradigm (Liao et al. 2013).
This package also includes a program for detecting exon-exon junctions, \R{subjunc}, that makes use of the ``seed-and-vote" paradigm as well.

An important step in processing next-gen sequencing data is to assign mapped reads to genomic features such as genes, exons, and promoters.
This package includes a general-purpose read summarization function \R{featureCounts} that takes mapped reads as input and assigns them to genomic features.
In-built annotations are provided for users convenience.

Different from microarray technologies, the next-gen sequencing technologies do not provide Present/Absent calls for genomic features such as genes.
We have developed an algorithm to use the background noise measured from the RNA-seq data to call absolutely expressed genes.
The function \R{detectionCall} returns a detection p value for each gene from the read mapping results.

We have also developed a new SNP calling algorithm which is being implemented in function \R{exactSNPs}.
Our results showed that it compared favorably to competing methods, but was an order of magnitude faster.

This package also includes some other useful functions such as quality assessment (\R{qualityScores}, \R{atgcContent}), duplicate read removal (\R{removeDupReads}) and mapping percentage calculation (\R{propmapped}).

\section{Read alignment}
An index needs to be built first and then alignments can be carried out.
Building the index is an one-off operation.
The generated index can be re-used in subsequent read alignments.\\

{\noindent\bf Step 1: Index building \\}

The Rsubread package includes a dummy reference sequence that was generated by concatenating 900 100bp reads that were taken from a pilot dataset generated from the SEquencing Quality Control (SEQC) project.
We further extracted 100 reads from the same dataset and combine them with the 900 reads to make a read dataset for mapping.
Below is the command for building an index for the reference sequence:

\begin{scriptsize}
<<>>=
library(Rsubread)
ref <- system.file("extdata","reference.fa",package="Rsubread")
buildindex(basename="reference_index",reference=ref)
@
\end{scriptsize}

The generated index files were saved to the current working directory.
Rsubread creates a hash table for indexing the reference genome. 
Keys in the hash table are the 16bp sequences and hash values are their corresponding chromosomal locations. 
Color space index can be built by setting the \R{colorsapce} argument to \R{TRUE}.\\

A unique feature of Rsubread is that it allows users to control the computer memory usage in read mapping process.
Users can do this by tuning the amount of memory (in MB) to be used in read mapping.\\

{\noindent\bf Step 2: read mapping\\}

After the index was successfully built, we map the read dataset (including 1,000 reads) to the reference sequence:

\begin{scriptsize}
<<>>=
reads <- system.file("extdata","reads.txt.gz",package="Rsubread")
align(index="reference_index",readfile1=reads,output_file="alignResults.BAM",phredOffset=64)
@
\end{scriptsize}

Map paired-end reads:

\begin{scriptsize}
<<>>=
reads1 <- system.file("extdata","reads1.txt.gz",package="Rsubread")
reads2 <- system.file("extdata","reads2.txt.gz",package="Rsubread")
align(index="reference_index",readfile1=reads1,readfile2=reads2,
output_file="alignResultsPE.BAM",phredOffset=64)
@
\end{scriptsize}

\section{Counting mapped reads for genomic features}

The \R{featureCounts} function is a general-purpose read summarization function that assigns mapped reads (RNA-seq or gDNA-seq reads) to genomic features such as genes, exons, promoters, gene bodies and genomic bins.

This function takes as input a set of files that contain read mapping results and an annotation file that includes genomic features.
It automatically detects the format of input read files (supported formats include SAM and BAM).
For paired end reads, it automatically re-orders the reads if reads from the same pair were found not to be next to each other in the input.

In-built NCBI RefSeq gene annotations for genomes mm9, mm10 and hg19 are provided for convenience.
These annotations include chromosomal coordinates of exons of each gene.
When these annotations are used for summarization, only reads overlapping with exons will be counted by \R{featureCounts}.
Users can use {\texttt{getInBuiltAnnotation}} function to retrieve these annotations.

Below gives the example code of assigning reads and fragments, generated in the last section, to two artificial genes.
Assign single end reads to genes:

\begin{scriptsize}
<<>>=
ann <- data.frame(
GeneID=c("gene1","gene1","gene2","gene2"),
Chr="chr_dummy",
Start=c(100,1000,3000,5000),
End=c(500,1800,4000,5500),
Strand=c("+","+","-","-"),
stringsAsFactors=FALSE)
ann
fc_SE <- featureCounts("alignResults.BAM",annot.ext=ann)
fc_SE
@
\end{scriptsize}

Assign fragments (read pairs) to the two genes:

\begin{scriptsize}
<<>>=
fc_PE <- featureCounts("alignResultsPE.BAM",annot.ext=ann,isPairedEnd=TRUE)
fc_PE
@
\end{scriptsize}

\section{Finding exon junctions}

The RNA-seq technology provides a unique opportunity to identify the alternative splicing events that occur during the gene transcription process.
The \R{subjunc} function can be used to detect exon-exon junctions.
It first extracts a number of subreads (16mers) from each read, maps them to the reference genome and identifies the two best mapping locations for each read (representing potential locations of exons spanned by the read).
Then, it builds a junction table including all putative junctions.
Finally, it carries out a verification step to remove false positives in junction detection by realigning all the reads. 
The donor (`GT') and receptor sites(`AG'), are required to be present when calling exon-exon junctions. 
Output of this function includes the discovered exon-exon junctions and also read mapping results.

\section{Base quality scores}

Quality scores give the probabilities of read bases being incorrectly called, which is useful for examining the quality of sequencing data.
The \R{qualityScores} function can be used to quickly retrieve and display the quality score data extracted from a read file.

\begin{scriptsize}
<<>>=
x <- qualityScores(filename=reads,offset=64,nreads=1000)
x[1:10,1:10]
@
\end{scriptsize}

\section{GC content}

The \R{atgcContent} function returns fractions of A, T, G and C bases at each base location of reads or in the entire dataset.

%\begin{scriptsize}
%<<>>=
%## Fraction of A,T,G and C bases in all the reads
%x <- atgcContent(filename=reads)
%x
%## Fraction of A,T,G and C bases at each read location
%xb <- atgcContent(filename=reads,basewise=TRUE)
%xb[,1:5]
%@
%\end{scriptsize}

\section{Mapping percentage}

Function \R{propmapped} returns the proportion of mapped reads included in a SAM/BAM file.
For paired end reads, it can return the proportion of mapped fragments (ie. read pairs).

\begin{scriptsize}
<<>>=
propmapped("alignResults.BAM")
@
\end{scriptsize}

\section{Citation}
Yang Liao, Gordon K Smyth and Wei Shi (2013). The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Research, 41(10):e108.\\

\noindent{Yang Liao, Gordon K Smyth and Wei Shi (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7):923-30}

\section{Authors}
Wei Shi and Yang Liao \\
Bioinformatics Division \\
The Walter and Eliza Hall Institute of Medical Research \\
1G Royal Parade, Parkville, Victoria 3052 \\
Australia \\

\section{Contact}
Please post to Bioconductor mailing list (\url{http://bioconductor.org/}) if you find any bugs and have any inquires.
Or, you may contact Wei Shi (shi at wehi dot edu dot au) directly.

\end{document}