% \VignetteIndexEntry{ChIPseeker: an R package for ChIP peak Annotation, Comparison, and Visualization}
% \VignettePackage{ChIPseeker}
% \VignetteEngine{knitr::knitr}

% To compile this document, run the commands within R
% R CMD Sweave --engine=knitr::knitr --pdf ChIPseeker.Rnw

\documentclass[12pt]{article}

<<knitr, echo=FALSE, result="hide">>=
knitr::opts_chunk$set(tidy         = FALSE,
                      out.truncate = 80,
                      out.lines    = 6,
                      dev          = 'pdf',
                      include      = TRUE,
                      fig.width    = 6,
                      fig.height   = 6,
                      resolution   = 100,
                      fig.show     = 'asis',
                      warning      = FALSE,
                      message      = FALSE)
@

<<style-knitr, eval=TRUE, echo=FALSE, results="asis">>=
BiocStyle::latex()
@ 

<<loadPKG, echo=FALSE>>=
library(GenomicFeatures)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
library(clusterProfiler)
library(ChIPseeker)
@ 


\author{Guangchuang Yu \\[1em]
  \small{School of Public Health} \\
  \small{The University of Hong Kong} \\
  \small{\email{guangchuangyu@gmail.com}}
}

\title{ChIPseeker: an R package for ChIP peak Annotation, Comparision and Visualization}


\begin{document}
\maketitle

<<options,results='hide',echo=FALSE>>=
options(digits=3, width=80, prompt=" ", continue=" ")
@

\begin{abstract}
ChIPseeker is an R package for annotating ChIP-seq data analysis. It supports annotating ChIP peaks and provides functions to visualize ChIP peaks coverage over chromosomes and profiles of peaks binding to TSS regions. Comparison of ChIP peak profiles and annotation are also supported. Moreover, it supports evaluating significant overlap among ChIP-seq datasets. Currently, ChIPseeker contains 17,000 bed file information from GEO database. These datasets can be downloaded and compare with user's own data to explore significant overlap datasets for inferring co-regulation or transcription factor complex for further investigation.

  \begin{center}
    
    \vspace{1em}
    \textbf{\Biocpkg{ChIPseeker} version:} \Sexpr{packageVersion("ChIPseeker")}
    \vspace{1em}
    
    \begin{tabular}{ | l | }
      \hline
      If you use \Biocpkg{ChIPseeker} in published research, please cite: \\
      \\
      G Yu, LG Wang, QY He. \textbf{ChIPseeker: an R/Bioconductor package for} \\
      \textbf{ChIP peak annotation, comparison and visualization}. \\
      \emph{Bioinformatics} 2015. 31(14): 2382-2383.\\
      \url{http://dx.doi.org/10.1093/bioinformatics/btv145} \\
      \hline
    \end{tabular}
    
  \end{center}
\end{abstract}


\newpage
\tableofcontents
\newpage

\section{Introduction}
Chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) has become standard technologies for genome wide identification of DNA-binding protein target sites. After read mappings and peak callings, the peak should be annotated to answer the biological questions. Annotation also create the possibility of integrate expression profile data to predict gene expression regulation. \Biocpkg{ChIPseeker} was developed for annotating nearest genes and genomic features to peaks. 
\\
\\
ChIP peak data set comparison is also very important. We can use it as an index to estimate how well biological replications are. Even more important is applying to infer cooperative regulation. If two ChIP seq data, obtained by two different binding proteins, overlap significantly, these two proteins may form a complex or have interaction in regulation chromosome remodelling or gene expression. \Biocpkg{ChIPseeker} support statistical testing of significant overlap among ChIP seq data sets, and incorporate open access database GEO for users to compare their own dataset to those deposited in database. Protein interaction hypothesis can be generated by mining data deposited in database. Converting genome coordinations from one genome version to another is also supported, making this comparison available for different genome version and different species.
\\
\\
Several visualization functions are implemented to visualize the coverage of the ChIP seq data, peak annotation, average profile and heatmap of peaks binding to TSS region. 
\\
\\
Functional enrichment analysis of the peaks can be performed by my Bioconductor packages \Biocpkg{DOSE} \cite{yu_dose_2015}, \Biocpkg{ReactomePA}, \Biocpkg{clusterProfiler} \cite{yu_clusterprofiler_2012} . 

<<loadPkgs>>=
## loading packages
require(ChIPseeker)
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
require(clusterProfiler)
@ 
\section{ChIP profiling}
The datasets CBX6 and CBX7 in this vignettes were downloaded from GEO (GSE40740) \cite{pemberton_genome-wide_2014} while ARmo\_0M, ARmo\_1nM and ARmo\_100nM were downloaded from GEO (GSE48308) \cite{urbanucci_overexpression_2012} . \Biocpkg{ChIPseeker} provides \Rfunction{readPeakFile} to load the peak and store in \Robject{GRanges} object. 

<<samplefiles>>=
files <- getSampleFiles()
print(files)
peak <- readPeakFile(files[[4]])
peak
@ 
\subsection{ChIP peaks coverage plot}
After peak calling, we would like to know the peak locations over the whole genome, \Rfunction{covplot} function calculates the coverage of peak regions over chromosomes and generate a figure to visualize.
<<chromosomeCoverage, fig.height=8, fig.width=7, out.width="0.8\\textwidth">>=
covplot(peak, weightCol="V5")
@ 


<<covplot, fig.height=3, fig.width=7, out.width="0.6\\textwidth">>=
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4.5e7, 5e7))
@ 

\subsection{Profile of ChIP peaks binding to TSS regions}
First of all, for calculate the profile of ChIP peaks binding to TSS regions, we should prepare the TSS regions, which are defined as the flanking sequence of the TSS sites. Then align the peaks that are mapping to these regions, and generate the tagMatrix.

<<promoter>>= 
## promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
## tagMatrix <- getTagMatrix(peak, windows=promoter)
##
## to speed up the compilation of this vignettes, we use a precalculated tagMatrix
data("tagMatrixList")
tagMatrix <- tagMatrixList[[4]]
@

In the above code, you should notice that tagMatrix is not restricted to TSS regions. The regions can be other types that defined by the user.

\subsubsection{Heatmap of ChIP binding to TSS regions}

<<heatmap, eval=TRUE, fig.cap="Heatmap of ChIP peaks binding to TSS regions", fig.align="center", fig.height=12, fig.width=4, out.width="0.2\\textwidth", fig.pos="h">>=
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
@ 

\Biocpkg{ChIPseeker} provide a one step function to generate this figure from bed file. The following function will generate the same figure as above.
<<heatmap2, eval=FALSE>>=
peakHeatmap(files[[4]], TxDb=txdb, upstream=3000, downstream=3000, color="red")
@ 

\subsubsection{Average Profile of ChIP peaks binding to TSS region}
<<avgprof, eval=TRUE, fig.cap="Average Profile of ChIP peaks binding to TSS region", fig.align="center", fig.height=4, fig.width=7, out.width="0.6\\textwidth", fig.pos="h">>=
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
@ 


The function \Rfunction{plotAvgProf2} provide a one step from bed file to average profile plot. The following command will generate the same figure as shown above.
<<avgprof2, eval=FALSE>>=
plotAvgProf2(files[[4]], TxDb=txdb, upstream=3000, downstream=3000, 
             xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
@

Confidence interval estimated by bootstrap method is also supported for characterizing ChIP binding profiles.

<<avgprofCI, fig.cap="Average Profile of ChIP peaks binding to TSS region", fig.align="center", fig.height=4, fig.width=7, out.width="0.6\\textwidth", results='hide', fig.pos="h">>=
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 500)
@ 

\section{Peak Annotation}
<<peak annotation>>=
peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Hs.eg.db")
@ 

Peak Annotation is performed by \Rfunction{annotatePeak}. User can define TSS (transcription start site) region, by default TSS is defined from -3kb to +3kb. The output of \Rfunction{annotatePeak} is \Robject{csAnno} instance. \Biocpkg{ChIPseeker} provides \Rfunction{as.GRanges} to convert \Robject{csAnno} to \Robject{GRanges} instance, and \Rfunction{as.data.frame} to convert \Robject{csAnno} to \Robject{data.frame} which can be exported to file by \Rfunction{write.table}.
\\
\\
\Robject{TxDb} object contained transcript-related features of a particular genome. Bioconductor provides several package that containing \Robject{TxDb} object of model organisms with multiple commonly used genome version, for instance \Biocannopkg{TxDb.Hsapiens.UCSC.hg38.knownGene}, \Biocannopkg{TxDb.Hsapiens.UCSC.hg19.knownGene} for human genome hg38 and hg19, \Biocannopkg{TxDb.Mmusculus.UCSC.mm10.knownGene} and

\Biocannopkg{TxDb.Mmusculus.UCSC.mm9.knownGene} for mouse genome mm10 and mm9, etc. User can also prepare their own \Robject{TxDb} object by retrieving information from UCSC Genome Bioinformatics and BioMart data resources by R function

\Rfunction{makeTranscriptDbFromBiomart} and \Rfunction{makeTranscriptDbFromUCSC}. \Robject{TxDb} object should be passed for peak annotation. 
\\
\\
All the peak information contained in peakfile will be retained in the output of \Rfunction{annotatePeak}. The position and strand information of nearest genes are reported. The distance from peak to the TSS of its nearest gene is also reported. The genomic region of the peak is reported in annotation column. Since some annotation may overlap, \Biocpkg{ChIPseeker} adopted the following priority in genomic annotation.
\begin{itemize}
  \item Promoter
  \item 5' UTR
  \item 3' UTR
  \item Exon
  \item Intron
  \item Downstream
  \item Intergenic 
\end{itemize}

Downstream is defined as the downstream of gene end.
\\
\Biocpkg{ChIPseeker} also provides parameter \Robject{genomicAnnotationPriority} for user to prioritize this hierachy.
\\
\\
\Rfunction{annotatePeak} report detail information when the annotation is Exon or Intron, for instance "Exon (uc002sbe.3/9736, exon 69 of 80)", means that the peak is overlap with an Exon of transcript uc002sbe.3, and the corresponding Entrez gene ID is 9736 (Transcripts that belong to the same gene ID may differ in splice events), and this overlaped exon is the 69th exon of the 80 exons that this transcript uc002sbe.3 prossess. 
\\
\\
Parameter annoDb is optional, if provided, extra columns including SYMBOL, GENENAME, ENSEMBL/ENTREZID will be added. The geneId column in annotation output will be consistent with the geneID in TxDb. If it is ENTREZID, ENSEMBL will be added if annoDb is provided, while if it is ENSEMBL ID, ENTREZID will be added.

\subsection{Visualize Genomic Annotation}
To annotate the location of a given peak in terms of genomic features, \Rfunction{annotatePeak} assigns peaks to genomic annotation in "annotation" column of the output, which includes whether a peak is in the TSS, Exon, 5' UTR, 3' UTR, Intronic or Intergenic. Many researchers are very interesting in these annotations. TSS region can be defined by user and \Rfunction{annotatePeak} output in details of which exon/intron of which genes as illustrated in previous section.
\\
\\
Pie and Bar plot are supported to visualize the genomic annotation.
<<pieplot, fig.cap="Genomic Annotation by pieplot", fig.align="center", fig.height=6, fig.width=8, out.width="0.6\\textwidth", fig.pos="h">>= 
plotAnnoPie(peakAnno)
@ 

<<barplot_anno, fig.cap="Genomic Annotation by barplot", fig.align="center", fig.height=4, fig.width=10, out.width="0.7\\textwidth", fig.pos="h">>= 
plotAnnoBar(peakAnno)
@ 

Since some annotation overlap, user may interested to view the full annotation with their overlap, which can be partially resolved by \Rfunction{vennpie} function.

<<vennpie, fig.cap="Genomic Annotation by vennpie", fig.align="center", fig.height=8, fig.width=11, out.width="0.8\\textwidth", fig.pos="h">>= 
vennpie(peakAnno)
@ 


\subsection{Visualize distribution of TF-binding loci relative to TSS}
The distance from the peak (binding site) to the TSS of the nearest gene is calculated by \Rfunction{annotatePeak} and reported in the output. We provide \Rfunction{plotDistToTSS} to calculate the percentage of binding sites upstream and downstream from the TSS of the nearest genes, and visualize the distribution. 
<<barplot_tss, fig.cap="Distribution of Binding Sites", fig.align="center", fig.height=2, fig.width=6, out.width="0.8\\textwidth", fig.pos="h">>= 
plotDistToTSS(peakAnno, 
              title="Distribution of transcription factor-binding loci\nrelative to TSS")
@

\section{Functional enrichment analysis}
Once we have obtained the annotated nearest genes, we can perform functional enrichment analysis to identify predominant biological themes among these genes by incorporating biological knowledge provided by biological ontologies. For instance, Gene Ontology (GO) \cite{ashburner_gene_2000} annotates genes to biological processes, molecular functions, and cellular components in a directed acyclic graph structure, Kyoto Encyclopedia of Genes and Genomes (KEGG) \cite{kanehisa_kegg_2004} annotates genes to pathways, Disease Ontology (DO) \cite{schriml_disease_2011} annotates genes with human disease association, and Reactome \cite{croft_reactome_2013} annotates gene to pathways and reactions. 
\\
\\
Enrichment analysis is a widely used approach to identify biological themes. I have developed several Bioconductor packages for investigating whether the number of selected genes associated with a particular biological term is larger than expected, including \Biocpkg{DOSE} \cite{yu_dose_2015} for Disease Ontology, \Biocpkg{ReactomePA} for reactome pathway, \Biocpkg{clusterProfiler} \cite{yu_clusterprofiler:_2012} for Gene Ontology and KEGG enrichment analysis. 

<<enrichment>>= 
require(clusterProfiler)
bp <- enrichGO(as.data.frame(peakAnno)$geneId, ont="BP", readable=TRUE)
head(summary(bp), n=3)
@ 

More information can be found in the vignettes of Bioconductor packages \Biocpkg{DOSE} \cite{yu_dose_2015}, \Biocpkg{ReactomePA}, \Biocpkg{clusterProfiler} \cite{yu_clusterprofiler_2012}, which also provide several methods to visualize enrichment results. The \Biocpkg{clusterProfiler} \cite{yu_clusterprofiler_2012} is designed for comparing and visualizing functional profiles among gene clusters, and can directly applied to compare biological themes at GO, DO, KEGG, Reactome perspective.
\\
\\
\section{ChIP peak data set comparison}
\subsection{Profile of several ChIP peak data binding to TSS region}
Function \Rfunction{plotAvgProf} and \Rfunction{tagHeatmap} can accept a list of \Robject{tagMatrix} and visualize profile or heatmap among several ChIP experiments, while \Rfunction{plotAvgProf2} and \Rfunction{peakHeatmap} can accept a list of bed files and perform the same task in one step.

\subsubsection{Average profiles}
<<avgprof3, eval=TRUE, fig.cap="Average Profiles of ChIP peaks among different experiments", fig.align="center", fig.height=4, fig.width=6, out.width="0.7\\textwidth", fig.pos='H', fig.show='asis'>>=
## promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
## tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
## 
## to speed up the compilation of this vigenette, we load a precaculated tagMatrixList
data("tagMatrixList")
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
@ 

<<avgprof4, eval=TRUE, fig.cap="Average Profiles of ChIP peaks among different experiments", fig.align="center", fig.height=7, fig.width=6, out.width="0.7\\textwidth", fig.pos='H', results='hide', message=FALSE, fig.show='asis'>>=
## resample = 500 by default, here use 100 to speed up the compilation of this vignette.
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=100, facet="row")
@ 


\subsubsection{Peak heatmaps}
<<heatmap3, eval=TRUE, fig.cap="Heatmap of ChIP peaks among different experiments", fig.align="center", fig.height=8, fig.width=7, out.width="0.8\\textwidth", fig.pos="h", fig.show='asis'>>=
tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)
@ 


\subsection{ChIP peak annotation comparision}
The \Rfunction{plotAnnoBar} and \Rfunction{plotDistToTSS} can also accept input of a named list of annotated peaks (output of \Rfunction{annotatePeak}).


<<List_of_peakAnno>>=
peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb, 
                       tssRegion=c(-3000, 3000), verbose=FALSE)
@

We can use \Rfunction{plotAnnoBar} to comparing their genomic annotation.
<<barplot_anno_list, fig.cap="Genomic Annotation among different ChIPseq data", fig.align="center", fig.height=4, fig.width=6, out.width="0.6\\textwidth", fig.pos="h", fig.show='asis'>>= 
plotAnnoBar(peakAnnoList)
@

R function \Rfunction{plotDistToTSS} can use to comparing distance to TSS profiles among ChIPseq data.
<<barplot_tss_list, fig.cap="Distribution of Binding Sites among different ChIPseq data", fig.align="center", fig.height=5, fig.width=8, out.width="0.6\\textwidth", fig.pos="h">>= 
plotDistToTSS(peakAnnoList)
@

\subsection{Functional profiles comparison}
As shown in section 4, the annotated genes can analyzed by \Biocpkg{clusterProfiler} \cite{yu_clusterprofiler:_2012}, \Biocpkg{DOSE} and \Biocpkg{ReactomePA} for Gene Ontology, KEGG, Disease Ontology and Reactome Pathway enrichment analysis. 
\\
\\
The \Biocpkg{clusterProfiler} \cite{yu_clusterprofiler:_2012} package provide \Rfunction{compareCluster} function for comparing biological themes among gene clusters, and can be easily adopted to compare different ChIP peak experiments.

<<compareCluster, eval=FALSE>>=
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
names(genes) = sub("_", "\n", names(genes))
compGO <- compareCluster(geneCluster   = genes, 
                         fun           = "enrichGO", 
                         ont           = "MF", 
                         organism      = "human", 
                         pvalueCutoff  = 0.05, 
                         pAdjustMethod = "BH")
plot(compGO, showCategory = 20, title = "Molecular Function Enrichment")
@

\begin{figure}[h]
  \begin{center}
    \includegraphics[width=.8\linewidth]{figures/mf.png}
    \caption{Compare Biological themes among different experiments}
    \label{Fig:MF}
  \end{center}
\end{figure}


\subsection{Overlap of peaks and annotated genes}
User may want to compare the overlap peaks of replicate experiments or from different experiments. \Biocpkg{ChIPseeker} provides \Rfunction{peak2GRanges} that can read peak file and stored in GRanges object. Several files can be read simultaneously using lapply, and then passed to \Rfunction{vennplot} to calculate their overlap and draw venn plot.
\\
\\
\Rfunction{vennplot} accept a list of object, can be a list of GRanges or a list of vector. Here, I will demonstrate using \Rfunction{vennplot} to visualize the overlap of the nearest genes stored in peakAnnoList.

<<vennplot, fig.cap="Overlap of annotated genes", fig.align="center", fig.height=7, fig.width=7, out.width="0.6\\textwidth", fig.pos="h">>= 
genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)
@ 

\section{Statistical testing of ChIP seq overlap} 
Overlap is very important, if two ChIP experiment by two different proteins overlap in a large fraction of their peaks, they may cooperative in regulation. Calculating the overlap is only touch the surface. \Biocpkg{ChIPseeker} implemented statistical methods to measure the significance of the overlap.

\subsection{Shuffle genome coordination}

<<shuffle>>=
p <- GRanges(seqnames=c("chr1", "chr3"), 
             ranges=IRanges(start=c(1, 100), end=c(50, 130)))
shuffle(p, TxDb=txdb)
@ 

We implement the \Rfunction{shuffle} function to randomly permute the genomic locations of ChIP peaks defined in a genome which stored in \Robject{TxDb} object.

\subsection{Peak overlap enrichment analysis}
With the ease of this \Rfunction{shuffle} method, we can generate thousands of random ChIP data and calculate the background null distribution of the overlap among ChIP data sets.

<<peakEnrichment>>=
enrichPeakOverlap(queryPeak     = files[[5]], 
                  targetPeak    = unlist(files[1:4]), 
                  TxDb          = txdb, 
                  pAdjustMethod = "BH", 
                  nShuffle      = 50, 
                  chainFile     = NULL,
                  verbose       = FALSE)
@ 

Parameter \Robject{queryPeak} is the query ChIP data, while \Robject{targetPeak} is bed file name or a vector of bed file names from comparison; \Robject{nShuffle} is the number to shuffle the peaks in \Robject{targetPeak}. To speed up the compilation of this vignettes, we only set \Robject{nShuffle} to 50 as an example for only demonstration. User should set the number to 1000 or above for more robust result. Parameter \Robject{chainFile} are chain file name for mapping the \Robject{targetPeak} to the genome version consistent with \Robject{queryPeak} when their genome version are different. This creat the possibility of comparison among different genome version and cross species.

In the output, \Robject{qSample} is the name of \Robject{queryPeak} and \Robject{qLen} is the the number of peaks in \Robject{queryPeak}. \Robject{N\_OL} is the number of overlap between \Robject{queryPeak} and \Robject{targetPeak}. 


\section{Data Mining with ChIP seq data deposited in GEO}
There are many ChIP seq data sets that have been published and deposited in GEO database. We can compare our own dataset to those deposited in GEO to search for significant overlap data. Significant overlap of ChIP seq data by different binding proteins may be used to infer cooperative regulation and thus can be used to generate hypotheses.
\\
\\
We collect about 15,000 bed files deposited in GEO, user can use \Rfunction{getGEOspecies} to get a summary based on speices.

\subsection{GEO data collection}
<<getGEOspecies>>=
getGEOspecies()
@ 


The summary can also based on genome version as illustrated below:
<<getGEOgenomeVersion>>=
getGEOgenomeVersion()
@ 

User can access the detail information by \Rfunction{getGEOInfo}, for each genome version.
<<getGEOInfo>>=
hg19 <- getGEOInfo(genome="hg19", simplify=TRUE)
head(hg19)
@ 

If \Robject{simplify} is set to \Robject{FALSE}, extra information including \Robject{source\_name}, \Robject{extract\_protocol}, \Robject{description}, \Robject{data\_processing}, and \Robject{submission\_date} will be incorporated.

\subsection{Download GEO ChIP data sets}
\Biocpkg{ChIPseeker} provide function \Rfunction{downloadGEObedFiles} to download all the bed files of a particular genome.

<<downloadGEO, eval=FALSE>>=
downloadGEObedFiles(genome="hg19", destDir="hg19")
@ 

Or a vector of GSM accession number by \Rfunction{downloadGSMbedFiles}.
<<downloadGSM, eval=FALSE>>=
gsm <- hg19$gsm[sample(nrow(hg19), 10)]
downloadGSMbedFiles(gsm, destDir="hg19")
@ 


\subsection{Overlap significant testing}
After download the bed files from GEO, we can pass them to \Rfunction{enrichPeakOverlap} for testing the significant of overlap. Parameter \Robject{targetPeak} can be the folder, e.g. hg19, that containing bed files. \Rfunction{enrichPeakOverlap} will parse the folder and compare all the bed files. It is possible to test the overlap with bed files that are mapping to different genome or different genome versions, \Rfunction{enrichPeakOverlap} provide a parameter \Robject{chainFile} that can pass a chain file and liftOver the \Robject{targetPeak} to the genome version consistent with \Robject{queryPeak}. Signifcant overlap can be use to generate hypothesis of cooperative regulation.By mining the data deposited in GEO, we can identify some putative complex or interacted regulators in gene expression regulation or chromsome remodelling for further validation.  

\section{External documents}

\begin{itemize}
  \item \href{http://ygc.name/2014/01/14/bug-of-r-package-chippeakanno/}{Bug of R package ChIPpeakAnno}
  \item \href{http://ygc.name/2014/04/13/chipseeker-for-chip-peak-annotation/}{ChIPseeker for ChIP peak annotation}
  \item \href{http://ygc.name/2014/04/30/visualization-methods-in-chipseeker/}{visualization methods in ChIPseeker}
  \item \href{http://ygc.name/2014/10/01/multiple-annotation-in-chipseeker/}{multiple annotation in ChIPseeker}
\end{itemize}
    
\section{Session Information}

Here is the output of \Rcode{sessionInfo()} on the system on which this document was compiled:

<<sessInfo, results='asis', echo=FALSE>>=
toLatex(sessionInfo())
@


\bibliography{ChIPseeker}
\end{document}