%
%
%
%\VignetteIndexEntry{An R package for visualization and analysis of Cufflinks high-throughput sequencing data}
%\VignetteKeywords{cummeRbund,visualization,sequencing,cufflinks,cuffdiff}
%\VignettePackage{cummeRbund}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[10pt]{article}
\usepackage{amsmath}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}
\hypersetup{
    colorlinks,
    citecolor=black,
    filecolor=black,
    linkcolor=red,
    urlcolor=black
}
\usepackage{theorem}
\usepackage{float}
\usepackage{ifthen}
\usepackage[OT1]{fontenc}

\newcommand{\R}{{\textsf{R}}}
\newcommand{\code}[1]{{\texttt{#1}}}
\newcommand{\term}[1]{{\emph{#1}}}
\newcommand{\Rpackage}[1]{\textsf{#1}}
\newcommand{\Rfunction}[1]{\texttt{#1}}
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rclass}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\textit{#1}}}
\newcommand{\Rfunarg}[1]{{\textit{#1}}}

\bibliographystyle{plainnat}
\title{cummeRbund: Visualization and Exploration of Cufflinks High-throughput Sequencing Data}

\author{Loyal A. Goff, Cole Trapnell}
\date{1 April, 2011}

\begin{document}
<<init, echo=FALSE>>=
options(width=65)
@ 

\maketitle

\tableofcontents

\section{Requirements}
\begin{itemize}
	\item Cufflinks $\ge$ v1.1.0 (Note: as of the release of this package, the current version of cufflinks is 1.3.3.  While this package will work with >v1.1, we recommend updating your cufflinks prior to using cummeRbund)
	\item R $\ge$ v2.7.0
	\item Packages:
	\begin{itemize}
		\item \Rpackage{RSQLite}
		\item \Rpackage{ggplot2}
		\item \Rpackage{reshape}
		\item \Rpackage{plyr}
	\end{itemize}
\end{itemize}
	

\section{Introduction}
	\Rpackage{cummeRbund} is a visualization package for Cufflinks high-throughput sequencing data. The base class, \Rclass{cuffSet} 
	is a 'pointer' to cufflinks data that are stored out-of-memory in a sqlite database.
	
\section{CummeRbund Classes}

\subsection{CuffSet Class}
	A pointer class to control access to the sqlite tables holding the Cufflinks data. The primary slot is DB which contains the RSQLite connection object. This can be accessed using the \Rmethod{DB()} accessor.
	The additional slots (genes, isoforms, TSS, and CDS) are each instances of the \Rclass{CuffData} class and are pointers to sets of tables for each data subtype. They can be accessed with similar accessor wrappers.
	This is the default class created by \Rmethod{readCufflinks}.  By default, \Rclass{CuffData} accessor methods applied to a \Rclass{CuffSet} class will operate on the 'genes' slot.

\subsection{CuffData Class}
	The \Rclass{CuffData} class is also a pointer class to the SQL backend, but each instance is specific for a data subtype (genes, isoforms, TSS, CDS). Again, there is an DB slot (accessible using \Rmethod{DB()}) that contains the RSQLite connection object.
	There are several accessor, setter, and plotting methods that allow for global analysis of all features within a \Rmethod{CuffData} class.Subsetting is currently being re-written, however, it is primarily done through the 'gene\_id' field.
	Available slots for the CuffData class are: 
	\begin{itemize}
		\item DB: RSQLite connection object
		\item tables: A \Rclass{list} of tables in the SQLite DB that contain the cufflinks data.
		\item filters: A \Rclass{list} of filters for subsetting (not implemented yet).
		\item type: A \Rclass{character} field describing the data (ie. 'genes','isoforms','TSS','CDS','other')
		\item idField: The name of the identifying index field for this object (eg. 'gene\_id' for type='gene', or 'isoform\_id' for type='isoform')
	\end{itemize}
	Making the best use of either the CuffSet or CuffData classes will enable you to keep the entire dataset out of memory and significantly improve performance for large cufflinks datasets.

\subsection{CuffFeatureSet Class}
	The \Rclass{CuffFeatureSet} class is a data-storage container that holds all available data for a pre-determined list of features. Slots for FPKM data, differential regulation data, and feature-level annotation are all available. Unlike the previous classes, this class contains no connection information to the SQL database, but
	rather contains several slots with \Rclass{data.frame} objects storing multiple-features worth of information.  There are available accessors, and plotting methods that are designed to present multiple-features worth of information (eg. heatmaps, scatterplots, etc)
	Available slots for a \Rclass{CuffFeatureSet} object include:
	\begin{itemize}
		\item annotation: Holds all feature-level annotation information for all features in object.
		\item fpkm: A data frame of FPKM data across all samples, for all features in object.
		\item diff: A data frame of differential expression/regulation data for all features in object.
	\end{itemize}
	
	A specialized sub-class of \Rclass{CuffFeatureSet} is the \Rclass{CuffGeneSet} class. This subclass adds additional slots to contain all isoforms, TSS, and CDS information for a given set of gene\_ids.  The \Rclass{CuffGeneSet} class is designed to aggregate all relevant
	information for a set of genes into one object for easy analysis and/or manipulation.
	The \Rclass{CuffGeneSet} object adds the following slots:
	\begin{itemize}
		\item ids: A 'character' list of all gene\_ids used in object.
		\item isoforms: A \Rclass{CuffFeatureSet} object for all isoforms of genes in object.
		\item TSS: A \Rclass{CuffFeatureSet} object for all TSS of genes in object.
		\item CDS: A \Rclass{CuffFeatureSet} object for all CDS of genes in object.
	\end{itemize}

\subsection{CuffFeature Class}
	The \Rclass{CuffFeature} class is designed for single-feature-level data analysis and plotting.  The methods available for this object are designed to analyze or visualize information about a specific feature.
	This is a 'data' object, as opposed to a 'pointer' object to the database backend. There is a validity requirement that a \Rclass{CuffFeature} object only point to data from a single feature.
	Available slots for a \Rclass{CuffFeature} object include:
	\begin{itemize}
		\item annotation: Holds feature-level annotation information for a given feature.
		\item fpkm: A data frame of FPKM data across all samples for a given feature.
		\item diff: A data frame of differential expression/regulation data for a given feature.
	\end{itemize}
	
	A specialized sub-class of \Rclass{CuffFeature} is the \Rclass{CuffGene} class. This subclass adds additional slots to contain all isoform, TSS, and CDS information for a given gene.
	The \Rclass{CuffGene} object adds the following slots:
	\begin{itemize}
		\item id: The common 'gene\_id' for all data in object
		\item isoforms: A \Rclass{CuffFeature} object for all isoforms of a given gene.
		\item TSS: A \Rclass{CuffFeature} object for all TSS of a given gene.
		\item CDS: A \Rclass{CuffFeature} object for all CDS of a given gene.
	\end{itemize}
	
Note: Future versions of cummeRbund may try to collapse the redundant functionality of the CuffFeature and CuffFeatureSet classes.

\section{Reading cuffdiff output}
\Rpackage{cummeRbund} was designed to process the multi-file output format for a 'cuffdiff' differential expression analysis.  In this type of analysis, a user will use a reference .gtf file (either known annotation or a .gtf file created from a cufflinks assembly or merge of assemblies) and quantitate the expression values and differential regulation of the annotation(s) in the .gtf file across two or more SAM/BAM files.
By design, cuffdiff produces a number of output files that contain test results for changes in expression at the level of transcripts, primary transcripts, and genes. It also tracks changes in the relative abundance of transcripts sharing a common transcription start site, and in the relative abundances of the primary transcripts of each gene. Tracking the former allows one to see changes in splicing, and the latter lets one see changes in relative promoter use within a gene. \\

Note: Cuffdiff requires that transcripts in the input GTF be annotated with certain attributes in order to look for changes in primary transcript expression, splicing, coding output, and promoter use. \\

These attributes are:
\begin{itemize}
	\item tss\_id: The ID of this transcript's inferred start site. Determines which primary transcript this processed transcript is believed to come from. Cuffcompare appends this attribute to every transcript reported in the .combined.gtf file.
	\item p\_id	The ID of the coding sequence this transcript contains. This attribute is attached by Cuffcompare to the .combined.gtf records only when it is run with a reference annotation that include CDS records. Further, differential CDS analysis is only performed when all isoforms of a gene have p\_id attributes, because neither Cufflinks nor Cuffcompare attempt to assign an open reading frame to transcripts.
\end{itemize}

cuffdiff calculates the FPKM of each transcript, primary transcript, and gene in each sample. Primary transcript and gene FPKMs are computed by summing the FPKMs of transcripts in each primary transcript group or gene group. The results are output in FPKM tracking files, the structure of which can be found in the cufflinks manual.\\

There are four FPKM tracking files:
\begin{itemize}
	\item \emph{isoforms.fpkm\_tracking}	Transcript FPKMs
	\item \emph{genes.fpkm\_tracking}	Gene FPKMs. Tracks the summed FPKM of transcripts sharing each gene\_id
	\item \emph{cds.fpkm\_tracking}	Coding sequence FPKMs. Tracks the summed FPKM of transcripts sharing each p\_id, independent of tss\_id
	\item \emph{tss\_groups.fpkm\_tracking}	Primary transcript FPKMs. Tracks the summed FPKM of transcripts sharing each tss\_id
\end{itemize}

cuffdiff also performs differential expression tests between supplied conditions. This tab delimited file lists the results of differential expression testing between samples for spliced transcripts, primary transcripts, genes, and coding sequences. For detailed file structure see cufflinks manual. \\

Four .diff files are created:
\begin{itemize}
	\item \emph{isoform\_exp.diff}	Transcript differential FPKM.
	\item \emph{gene\_exp.diff}	Gene differential FPKM. Tests difference sin the summed FPKM of transcripts sharing each gene\_id
	\item \emph{tss\_group\_exp.diff}	Primary transcript differential FPKM. Tests differences in the summed FPKM of transcripts sharing each tss\_id
	\item \emph{cds\_exp.diff}	Coding sequence differential FPKM. Tests differences in the summed FPKM of transcripts sharing each p\_id independent of tss\_id
\end{itemize}

In addition, cuffdiff also performs differential splicing, CDS usage, and promoter usage tests for each gene across conditions:

\begin{itemize}
	\item \emph{splicing.diff}	Differential splicing tests.
	\item \emph{CDS.diff}	Differential coding output.
	\item \emph{promoters.diff}	Differential promoter use.
\end{itemize}

All of these output files are related to each other through their various tracking\_ids, but parsing through individual files to query for important result information requires both a good deal of patience and a strong grasp of command-line text manipulation. Enter cummeRbund, an R solution to aggregate, organize, and help visualize this multi-layered dataset. \\
One of the principle benefits of using cummeRbund is that data are stored in a SQLite database.  This allows for out-of-memory analysis of data, quick retrieval, and only a one-time cost to setup the tables. By default, cummeRbund assumes that all output files from cuffdiff are in the current working directory.
To read these files, populate the 'cuffData.db' database backend, and return the \Rclass{CuffSet} pointer object, you can do the following. 

<<loadLib>>=
library(cummeRbund)
@
%%fileDir<-("../../extdata/")
<<read>>=
cuff <- readCufflinks(system.file("extdata", package="cummeRbund"))
cuff
@
Again, by default $dir$ is assumed to be the current working directory and \code{cuff<-readCufflinks()} should work if all appropriate files are in the current working directory. Should you need to rebuild the SQLite backend for any reason, you can add the option \Rfunarg{rebuild=T} to \Rmethod{readCufflinks}.  Once the database is created, \Rmethod{readCufflinks} will default to using the SQL backend and should not need to rebuild this database.
Each R session should begin with a call to \Rmethod{readCufflinks} so as to initialize the database connection and create an object with the appropriate RSQLite connection information. 

\subsection{Adding additional feature annotation}
Gene- or feature-level annotation can be permanently added to the database tables for future querying. If you have a data.frame where the first column contains the 'tracking\_id' (eg. 'gene\_id' for genes, 'isoform\_id' for isoforms, etc). You can easily add feature level annotation using the \Rfunction{addFeatures()} function:

<<add_features>>=
#annot<-read.table("gene_annotation.tab",sep="\t",header=T,na.string="-")
#addFeatures(cuff,annot,level="genes")
@
By default, features added to a \Rclass{CuffSet} object are assumed to be gene-level annotations, but the level can selected using the argument \Rfunarg{level}. Features added to a \Rclass{CuffData} object are assumed to be of the same type as the 'type' value for that given object (e.g. gene-level features for 'genes', isoform-level features for isoforms, etc.)


\section{Global statistics}
Several plotting methods are available that allow for quality-control or global analysis of cufflinks data. For example, to assess the
distributions of FPKM scores across samples, you can use the \Rmethod{csDensity} plot (Figure 1).
<<global_plots_1,include=FALSE>>=
dens<-csDensity(genes(cuff))
dens
@


<<label=global_plots_dens,fig=TRUE,echo=FALSE>>=
<<global_plots_1>>
print(dens)
@

Boxplots can be visualized using the \Rmethod{csBoxplot} method (Figure 2).
<<global_plots_2,include=FALSE>>=
b<-csBoxplot(genes(cuff))
b
@


<<label=global_plots_box,fig=TRUE,echo=FALSE>>=
<<global_plots_2>>
print(b)
@

Pairwise comparisons can be made by using \Rmethod{csScatter}. You must specify the sample names to use for the $x$ and $y$ axes:
<<global_plots_3,include=FALSE>>=
s<-csScatter(genes(cuff),"hESC","Fibroblasts",smooth=T)
s
@

<<label=global_plots_scatter,fig=TRUE,echo=FALSE>>=
<<global_plots_3>>
print(s)
@

Volcano plots are also available for the \Rclass{CuffData} objects. Again, you must specify the comparisons by sample name.
<<global_plots_4,include=FALSE>>=
v<-csVolcano(genes(cuff),"hESC","Fibroblasts")
v
@

<<label=global_plots_volcano,fig=TRUE,echo=FALSE>>=
<<global_plots_4>>
print(v)
@

\clearpage

\section{Accessing Data}
Feature-level information can be accessed directly from a \Rclass{CuffData} object using the \Rmethod{fpkm}, \Rmethod{diffData}, or \Rmethod{features} methods:

<<data_access_1>>=
gene.features<-features(genes(cuff))
head(gene.features)
gene.fpkm<-fpkm(genes(cuff))
head(gene.fpkm)
isoform.fpkm<-fpkm(isoforms(cuff))
head(isoform.fpkm)
gene.diff<-diffData(genes(cuff))
head(gene.diff)
@

Vectors of sample names and feature names are available by using the \Rmethod{samples} and \Rmethod{featureNames} methods:

<<data_access_2>>=
sample.names<-samples(genes(cuff))
head(sample.names)
gene.featurenames<-featureNames(genes(cuff))
head(gene.featurenames)
@

To facilitate Bioconductor-like operations, an 'FPKM-matrix' can be returned easily using the \Rmethod{fpkmMatrix} method:
<<data_access_3>>=
gene.matrix<-fpkmMatrix(genes(cuff))
head(gene.matrix)
@

\subsection{Writing your own SQL accessors}
Since the cufflinks is a SQLite database backend, if you are familiar with SQL and/or RSQLite query construction, you can simply design your own SQL queries to access the data that you are after. 

\begin{figure}[h]
\centering

\includegraphics[angle=90, height=\textheight]{cuffData_schema.pdf}

\end{figure}

\clearpage 

\section{Creating Gene Sets}
Gene Sets (stored in a \Rclass{CuffGeneSet} object) can be created using the \Rmethod{getGenes} method on a CuffSet object.
You must first create a vector of 'gene\_id' or 'gene\_short\_name' values to identify the genes you wish to select:

<<create_geneset_1>>=
data(sampleData)
myGeneIds<-sampleIDs
myGeneIds
myGenes<-getGenes(cuff,myGeneIds)
myGenes
@
The same \Rmethod{fpkm}, \Rmethod{fpkmMatrix}, \Rmethod{features}, \Rmethod{diffData}, \Rmethod{samples}, and \Rmethod{featureNames} are available for instances of the \Rmethod{CuffGeneSet} class.


\subsection{Geneset level plots}
There are several plotting functions available for gene-set-level visualization:

<<geneset_plots_1,include=FALSE>>=
h<-csHeatmap(myGenes,cluster='both')
h
@

<<label=geneset_plots_heatmap,fig=TRUE,echo=FALSE>>=
<<geneset_plots_1>>
print(h)
@

<<geneset_plots_2,include=FALSE>>=
s<-csScatter(myGenes,"Fibroblasts","hESC",smooth=T)
s
@

<<label=geneset_plots_scatter,fig=TRUE,echo=FALSE>>=
<<geneset_plots_2>>
print(s)
@

<<geneset_plots_3,include=FALSE>>=
v<-csVolcano(myGenes,cluster='both')
v
@

<<label=geneset_plots_volcano,fig=TRUE,echo=FALSE>>=
<<geneset_plots_3>>
print(v)
@

Similar plots can be made for all sub-level features of a \Rclass{CuffGeneSet} class by specifying which slot you would like to plot (eg. \Rfunarg{isoforms(myGenes)},\Rfunarg{TSS(myGenes)},\Rfunarg{CDS(myGenes)}).

<<geneset_plots_4,include=FALSE>>=
ih<-csHeatmap(isoforms(myGenes),cluster='both',labRow=F)
ih
@

<<label=geneset_plots_isoform_heatmap,fig=TRUE,echo=FALSE>>=
<<geneset_plots_4>>
print(ih)
@

Rudimentary k-means clustering is implemented as well.

<<geneset_plots_5,include=FALSE>>=
ic<-csCluster(myGenes,k=4)
ic
@

<<label=geneset_plots_cluster,fig=TRUE,echo=FALSE>>=
<<geneset_plots_5>>
print(ic)
@

\clearpage

\section{Individual Genes}
An individual CuffGene object can be created by using the \Rfunction{getGene} function for a given 'gene\_id'.
<<gene_level_1>>=
myGeneId<-"PINK1"
myGene<-getGene(cuff,myGeneId)
myGene
head(fpkm(myGene))
head(fpkm(isoforms(myGene)))
@

\subsection{Gene-level plots}
<<gene_plots_1,include=FALSE>>=
gl<-expressionPlot(myGene)
gl
@

<<label=gene_plots_line,fig=TRUE,echo=FALSE>>=
<<gene_plots_1>>
print(gl)
@

<<gene_plots_2,include=FALSE>>=
gb<-expressionBarplot(myGene)
gb
@

<<label=gene_plots_bar,fig=TRUE,echo=FALSE>>=
<<gene_plots_2>>
print(gb)
@

<<gene_plots_3,include=FALSE>>=
igb<-expressionBarplot(isoforms(myGene))
igb
@

<<label=gene_plots_bar_isoforms,fig=TRUE,echo=FALSE>>=
<<gene_plots_3>>
print(igb)
@

\clearpage

\section{Data Exploration}
The cummeRbund package is more than just a visualization tool as well.  We are working to implement several different means of data exploration from gene and condition clustering, finding features with similar expression profiles, as well as incorporating Gene Ontology analysis.

\subsection{Finding similar genes}
One common question in large-scale gene expression analyses is 'How can I find genes with similar expression profiles to gene $x$?'. We have implemented a method, \Rmethod{findSimilar} to allow you to identify a fixed number of the most similar genes to a given gene of interest.
For example, if you wanted to find the 20 genes most similar to "PINK1", you could do the following:

<<similar_1, include=FALSE>>=
mySimilar<-findSimilar(cuff,"PINK1",n=20)
mySimilar.expression<-expressionPlot(mySimilar,logMode=T,showErrorbars=F)
@

<<label=similar_plots_1,fig=TRUE,echo=FALSE>>=
<<similar_1>>
print(mySimilar.expression)
@

You are also able to provide your own expression profile in lieu of a 'gene\_id'.  The vector provided must match the order and length of \Rmethod{samples()}.

<<similar_2, include=FALSE>>=
myProfile<-c(500,0,400)
mySimilar2<-findSimilar(cuff,myProfile,n=10)
mySimilar2.expression<-expressionPlot(mySimilar2,logMode=T,showErrorbars=F)
@

<<label=similar_plots_2,fig=TRUE,echo=FALSE>>=
<<similar_2>>
print(mySimilar2.expression)
@

\Rmethod{findSimilar()} uses the Jensen-Shannon distance between the probability distributions of each gene across conditions to determine the similarity.  
We have found this to be a more robust way to determine distance between genes using the high dynamic range of FPKM data. Future versions may allow for other dissimilarity measures to be used instead.

\clearpage

\section{Miscellaneous}
\begin{itemize}
	\item All plotting functions return ggplot objects and the resulting objects can be manipulated/faceted/altered using standard ggplot2 methods.
	\item There are occasional DB connectivity issues that arise.  Not entirely sure why yet.  If necessary, just \Rfunction{readCufflinks} again and this should solve connectivity issues with a new
	RSQLite connection object.  If connectivity continues to be a problem, try \Rfunction{cuff<-readCufflinks(rebuild=T)}
	\item I am still working on fully documenting each of the methods.  There are a good number of arguments that exist, but might be hard to find without looking at the source.
\end{itemize}

\clearpage

\section{Known Issues}
\begin{itemize}
	\item You must have at least one p\_id field (see cufflinks manual) in your cuffdiff reference gtf file. Otherwise no results will be populated for the CDS.diff files and nothing will be available for cummeRbund to parse. This is described in more detail in the cuffdiff section of the cufflinks user guide.
	\item Large cuffdiff runs (e.g. $\ge$10 conditions) produce very large results files.  These will take some time to parse and populate the cuffData.db sqlite database.  While this is only a one time cost, the process can take a while. We are working on making the table writes and indexing significantly faster.
	\item Cuffdiff does not 'require' that gene\_ids, isoform\_ids, TSS\_group\_ids, or CDS\_ids be unique in your reference gtf file. In fact, duplicate IDs will be aggregated by cummeRbund in the indexing phase and will produce undesireable effects. Please ensure that all of your IDs are unique prior to running cuffdiff (see cuffmerge for help) to avoid this issue.
\end{itemize}

\clearpage

\section{Session info}
<<label=close_connection,echo=FALSE>>=
end<-sqliteCloseConnection(cuff@DB)
@

<<session>>=
sessionInfo()
@

\end{document}