---
title: "CNE identification and visualisation"
author: "Ge Tan"
date: "`r doc_date()`"
package: "`r pkg_ver('CNEr')`"
abstract: >
  Comparative genomics has revealed noncoding DNA regions with extremely high conservation across large evolutionary distances. These regions have been termed conserved noncoding elements (CNEs). Despite of various resources of CNEs, there is no publicly available tool to identify these elements from scratch. We describe CNEr, an R package for large-scale identification and advanced visualisation of sets of CNEs from _Axt_ alignments. Given whole genome pairwise alignments of two species as input, our pipeline facilitates alignment screening, elements merging, calculation of CNE density and visualisation of CNEs in the form of horizon plots. Furthermore, it provides efficient scalable data structures for representing paired ranges on the genome, especially the support for pair alignments in addition to functions potentially useful for exploratory analysis of the identified elements. 

vignette: >
  %\VignetteIndexEntry{CNE identification and visualisation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output: BiocStyle::html_document
bibliography: CNEr.bib
---
```{r code, echo = FALSE}
date = "`r doc_date()`"
pkg = "`r pkg_ver('CNEr')`"
```

```{r global_options, echo=FALSE}
short=TRUE #if short==TRUE, do not echo code chunks
debug=FALSE
knitr::opts_chunk$set(echo=!short, warning=debug, message=debug, error=FALSE,
               cache.path = "cache/", 
               fig.path = "figures/")
```

# Introduction
Conserved noncoding elements (CNEs) are a pervasive class of elements 
clustering around genes with roles in development and differentiation 
in Metazoa  [@Woolfe:2004ur]. 
While many have been shown to act as long-range developmental enhancers 
[@Sandelin:2004bd], 
the source of their extreme conservation remains unexplained. 
To study the evolutionary dynamics of these elements 
and their relationship to the genes around which they cluster, 
it is essential to be able to produce genome-wide sets of these elements 
for a large number of species comparisons, 
each with multiple size and conservation thresholds.

This `r Biocpkg("CNEr")` package aims to detect CNEs and 
visualise them along the genome.
For performance reasons, the implementation of CNEs detection 
and corresponding I/O functions are primarily written as C extensions to R. 
We have used `r Biocpkg("CNEr")` to produce sets of CNEs by scanning pairwise whole-genome net alignments 
with multiple reference species, 
each with two different window sizes and a range of minimum identity thresholds. 
Then, to pinpoint the boundaries of CNE arrays, 
we compute the CNE densities as the percentages of length 
covered by CNEs within a user specified window size. 
Finally, we describe a novel visualisation method using horizon plot tracks 
that shows a superior dynamic range to the standard density plots, 
simultaneously revealing CNE clusters characterized 
by vastly different levels of sequence conservation. 
Such CNE density plots generated using precise locations of CNEs 
can be used to identify genes involved in developmental regulation, 
even novel genes that are not annotated yet.

# Workflow of the package
This section briefly demonstrates the pipeline of CNE identification 
and visualisation.
More detailed usage of each step is described in following sections with a concise example of CNE identification and visualisation for the "barhl2" and "sox14" (chr6:24,000,000..27,000,000) loci in Zebrafish (danRer10) genome against Human (hg38).

## CNE identification

1. Preparation of axtNets: the axtNet files can be downloaded from UCSC 
  or generated by this package.
2. Identification of conserved noncoding regions: scan the axt alignments for the regions, with minimal __I__ identities over __C__ columns, which do not overlap with annotated exons/repeats.
3. Merging overlapped elements : scanned elements from two pairs of alignments (one from each genome as reference) that overlap on both genomes are merged.
4. Removal of unannotated repetitive sequences: realign the elements back to respective genome with _BLAT_ and discard the elements which exceed a certain number of hits on either genome.

## CNE visualisation
1. Display parameters: chromosome, start, end, smooth window size.
2. Horizon plot: visualise CNE densities with the package `r Biocpkg("Gviz")`.

# Input
The minimal input for `r Biocpkg("CNEr")` includes the whole genome pairwise alignment of two assemblies, [axt](http://genome.ucsc.edu/goldenPath/help/axt.html)
net files.
UCSC already provides a set of precomputed axt files on http://hgdownload.soe.ucsc.edu/downloads.html for most of popular genomes.
In case the _axt_ net files are not available from UCSC, you can always generate the axt net files by following another vignette "Pairwise whole genome alignment" in this `r Biocpkg("CNEr")` package.

Another essential information is the annotation of exons and repeats, which could be retrieved from various resources.

## axt alignment file
During the development of this package, there was no suitable class to store the axt alignments in Bioconductor.
Hence, we created two new __S4__ classes, `Axt` and `GRangePairs`, to easily manipulate the axt alignment files.

`GRangePairs` is designed to hold a pair of `GRanges` objects, which have the same length.
It builds on the `Pairs` class of Bioconductor and inherits many useful methods from it.
This `Axt` class inherits from `GRangePairs` with no extra slots to hold the content from _axt_ files, but many specific methods are created for `Axt`.
The ranges of the target and the query organism are stored in two `GRanges` objects with alignment sequences as metadata columns.
The Blastz scores and the widths of the alignments are stored in metadata columns of `GRangePairs`.
For more information on the usage of these two classes, please refer to the documentation.

To read axt file into R, 
`r Biocpkg("CNEr")` provides `readAxt` function for highly efficient reading. 
This function is built on a backend C code of Kent's utilities [@Kent:2002bw].
The axt alignment files can be either gzippped or in plain text file.
The alignments between two genomes can also be in one big file or in several files, such as "chr1.hg19.mm10.net.axt.gz", "chr2.hg19.mm10.net.axt.gz", etc.

```{r Axt, eval=TRUE, echo=TRUE}
library(CNEr)
## These axt files are specially prepared for the region
## (chr6:24,000,000..27,000,000)
axtFilesHg38DanRer10 <- file.path(system.file("extdata", package="CNEr"), 
                                  "hg38.danRer10.net.axt")
axtFilesDanRer10Hg38 <- file.path(system.file("extdata", package="CNEr"), 
                                  "danRer10.hg38.net.axt")
```

```{r readAxt, eval=TRUE, echo=TRUE}
axtHg38DanRer10 <- readAxt(axtFilesHg38DanRer10)
axtDanRer10Hg38 <- readAxt(axtFilesDanRer10Hg38)
```

```{r showAxt, eval=TRUE, echo=TRUE}
## Axt class is shown in UCSC axt format
axtHg38DanRer10
axtDanRer10Hg38
```

```{r matchDistribution, eval=TRUE, echo=TRUE}
## Distribution of matched alignments; Given an Axt alignment, plot a heatmap with percentage of each matched alignment
matchDistribution(axtHg38DanRer10)
matchDistribution(axtDanRer10Hg38)
```

```{r syntenyDotplot, eval=TRUE, echo=TRUE}
## Example of chr4 on hg19 and galGal3
## The synteny of human and zebrafish is not quite obvious on the dotplot.
library(BSgenome.Hsapiens.UCSC.hg19)
library(BSgenome.Ggallus.UCSC.galGal3)
fn <- file.path(system.file("extdata", package="CNEr"),
                "chr4.hg19.galGal3.net.axt.gz")
axt <- readAxt(fn, 
               tAssemblyFn=file.path(system.file("extdata",
                                     package="BSgenome.Hsapiens.UCSC.hg19"),
                                     "single_sequences.2bit"),
               qAssemblyFn=file.path(system.file("extdata",
                                     package="BSgenome.Ggallus.UCSC.galGal3"),
                                     "single_sequences.2bit"))
library(GenomeInfoDb)
syntenicDotplot(axt, firstChrs=c("chr4"), secondChrs="chr4", type="dot")
```

There are methods defined for handling `Axt` objects, including subsetting, output to axt files.
More details can be found in the man page.

## Filtering information
The gene annotation information, including exons and repeats, is used to filter out the undesired regions.
Here we summarise a table of filtering information we used:

Assembly    |  Name     |  Exon    |   Repeat 
----------- | ----------- |  ----------------------------------------- | ---------
hg38        |  Human    |   RefSeq Genes, Ensembl Genes, UCSC Known Genes | RepeatMasker 
mm10        |  Mouse    |   RefSeq Genes, Ensembl Genes, UCSC Known Genes | RepeatMasker  
xenTro3     |  Frog     |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
tetNig2     |  Tetraodon |  Ensembl Genes   |   | 
canFam3     |  Dog      |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
galGal4     |  Chicken  |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
danRer10     |  Zebrafish |  RefSeq Genes, Ensembl Genes  | RepeatMasker 
fr3         |  Fugu     |   RefSeq Genes                 | RepeatMasker 
anoCar2     |  Lizard   |   Ensembl Genes                | RepeatMasker 
equCab2     |  Horse    |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
oryLat2     |  Medaka   |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
monDom5     |  Opossum  |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
gasAcu1     |  Stickleback | RefSeq Genes, Ensembl Genes | RepeatMasker 
rn5         |  Rat      |   RefSeq Genes, Ensembl Genes  | RepeatMasker 
dm3     |  D. melanogaster | RefSeq Genes, Ensembl Genes | RepeatMasker 
droAna2     |  D. ananassae  |    | RepeatMasker 
dp3         |  D. pseudoobscura  |     | RepeatMasker 
ce4         |  C. elegans        |     RefSeq Genes  |    RepeatMasker 
cb3         |  C. briggsae       |    |             RepeatMasker 
caeRem2     |  C. remanei        |     |  RepeatMasker  
caePb1      |  C. brenneri       |    |  RepeatMasker  

For the sake of simplicity, all the information listed above can be fetched easily
with Bioconductor package `r Biocpkg("rtracklayer")`, `r Biocpkg("biomaRt")` and
precompiled Bioconductor annotation packages.
A few examples are given here:

```{r UCSC, eval=FALSE, echo=TRUE}
## To fetch rmsk table from UCSC
library(rtracklayer)
mySession <- browserSession("UCSC")
genome(mySession) <- "hg38"
hg38.rmsk <- getTable(ucscTableQuery(mySession, track="RepeatMasker",
                                     table="rmsk"))
hg38.rmskGRanges <- GRanges(seqnames=hg38.rmsk$genoName,
                            ## The UCSC coordinate is 0-based.
                            ranges=IRanges(start=hg38.rmsk$genoStart+1,
                                           end=hg38.rmsk$genoEnd),
                            strand=hg38.rmsk$strand)
## To fetch ensembl gene exons from BioMart
library(biomaRt)
ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL",
                   host="dec2015.archive.ensembl.org")
ensembl <-  useDataset("hsapiens_gene_ensembl",mart=ensembl)
attributes <- listAttributes(ensembl)
exons <- getBM(attributes=c("chromosome_name", "exon_chrom_start",
                            "exon_chrom_end", "strand"), 
               mart=ensembl)
exonsRanges <- GRanges(seqnames=exons$chromosome_name,
                       ranges=IRanges(start=exons$exon_chrom_start,
                                      end=exons$exon_chrom_end),
                       strand=ifelse(exons$strand==1L, "+", "-")
                       )
seqlevelsStyle(exonsRanges) <- "UCSC"
## Use the existing Bioconductor annotation package for hg38
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
exonsRanges <- exons(TxDb.Hsapiens.UCSC.hg38.knownGene)
```

The regions to filter out can also be provided in a bed file.
To import the bed file into `r Biocpkg("GRanges")` in `R`, 
`r Biocpkg("rtracklayer")` provides a general function 
`import.bed` to do that.
Since only the chromosome names, start and end coordinates are used in `r Biocpkg("CNEr")`,
we provide a more efficient `readBed` function.

```{r Bed, eval=TRUE, echo=TRUE}
## Existing bed file for chr6:24,000,000..27,000,000 of Zebrafish danRer10
bedDanRer10Fn <- file.path(system.file("extdata", package="CNEr"), 
                           "filter_regions.danRer10.bed")
danRer10Filter <- readBed(bedDanRer10Fn)
danRer10Filter

## Existing bed file for alignment region in Human hg38 against
## chr6:24,000,000..27,000,000 of danRer10
bedHg38Fn <- file.path(system.file("extdata", package="CNEr"), 
                       "filter_regions.hg38.bed")
hg38Filter <- readBed(bedHg38Fn)
hg38Filter
```

## Creating a `CNE` class
We designed a `CNE` class to store all metadata of running the pipeline for identifying a set of CNEs between two species, including the intermediate and final results.
`CNE` can be created by providing the paths of the twoBit files of two assemblies, and
the paths of axt files, with each assembly as reference.

```{r CNE, eval=TRUE, echo=TRUE}
## Here we have the twoBit files from Bioconductor package
## BSgenome.Drerio.UCSC.danRer10 and BSgenome.Hsapiens.UCSC.hg38
cneDanRer10Hg38 <- CNE(
  assembly1Fn=file.path(system.file("extdata",
                                    package="BSgenome.Drerio.UCSC.danRer10"),
                        "single_sequences.2bit"),
  assembly2Fn=file.path(system.file("extdata",
                                    package="BSgenome.Hsapiens.UCSC.hg38"),
                        "single_sequences.2bit"),
  axt12Fn=axtFilesDanRer10Hg38, axt21Fn=axtFilesHg38DanRer10,
  cutoffs1=8L, cutoffs2=4L)
cneDanRer10Hg38
```

__Note__: the order of assemblies when creating CNE object is important.
Here we have danRer10 as assembly1 and hg38 as assembly2.
Then the `axt12Fn` contains the axt alignment with assembly1 danRer10 as reference and `axt21Fn` contains the alignment with assembly2 hg38 as reference.
The `cutoffs1` and `cutoffs2` are the maximal number of hits during the realignment in later steps.
Because zebrafish has undergone additional whole genome duplication compared to human, the cutoffs of zebrafish also doubles the cutoffs of human.

# CNE identification
In this section, we will go through the details of CNE identification.

## Scan axt alignments
Detecting CNEs highly relies on the whole-genome pairwise net alignments.
To correct the bias of a chosen genome (which bias?) and 
capture the duplicated CNEs during genome evolution,
we scan two sets of nets for each pairwise comparison,
one as reference from each of the genomes.

We identify CNEs by scanning the alignments for regions with 
at least __I__ identities over __C__ alignment columns.
Because different genes and loci may favor various similarity scores,
we usually scan at two diffrent window sizes 30 and 50 with 
several similarity criterias (__I/C__) range from 70% to 100%.

```{r CNEScan, eval=TRUE, echo=TRUE}
identities <- c(45L, 48L, 49L)
windows <- c(50L, 50L, 50L)
## Here danRer10Filter is tFilter since danRer10 is assembly1
cneListDanRer10Hg38 <- ceScan(x=cneDanRer10Hg38, tFilter=danRer10Filter,
                              qFilter=hg38Filter,
                              window=windows, identity=identities)
```

At this stage, a list of `CNE` is returned from `ceScan`, which contains the preliminary two sets of CNEs from a pair of axt alignments.
We can examine the intermediate CNEs by

```{r CNEScanHead, eval=TRUE, echo=TRUE}
## CNEs from the alignments with danRer10 as reference
CNE12(cneListDanRer10Hg38[["45_50"]])
## CNEs from the alignments with hg38 as reference
CNE21(cneListDanRer10Hg38[["45_50"]])
```

In the result table, even though the strand for query element can be negative,
the coordinate for that query element is already on the positive strand.

It is essential to scan two sets of pairwise net alignments with each assembly as reference, in order not to miss any duplicated elements in either lineage.
This is particularly important for the comparison between teleost fishes and other vertebrates, because one (for instance, the case of zebrafish) or two (the case of common carp) extra whole genome duplications occured.

## Merge CNEs
As we perform two rounds of CNE detection with each genome as reference,
some conserved elements overlap on both genomes and should be removed.
Elements, however, that overlap only on one of the genomes are kept, 
so that duplicated elements remain distinct.

```{r CNEMerge, eval=TRUE, echo=TRUE}
cneMergedListDanRer10Hg38 <- lapply(cneListDanRer10Hg38, cneMerge)
```

## Realignment of CNEs
Some CNEs might be unannotated repeats.
To remove them, currently we use __blat__ [@Kent:2002jd] to realign
each sequence of CNEs against the respective genomes.
When the number of matches exceeds a certain threshold, for instance 8,
that CNE will be discarded.

This step can be very time-consuming when the number of CNEs is large.
Other alignment methods could also be considered, for example Bowtie2 or BWA (provided that they are installed on the user's machine)

```{r CNERealignment, eval=FALSE, echo=TRUE}
cneFinalListDanRer10Hg38 <- lapply(cneMergedListDanRer10Hg38, blatCNE)
```

<!-- Now at this stage, these elements are the final CNEs. -->
<!-- We also prepare a one step function `ceScanOneStep` -->
<!-- and it returns a `CNE` object directly -->
<!-- which wrapps all the necessary information. -->
<!-- This one-step function is highly recommended to avoid the tedious steps above. -->

<!-- ```{r ceScanOneStep, eval=FALSE, echo=TRUE} -->
<!-- assemblyHg19Twobit <- "/Users/gtan/CSC/CNEr/2bit/hg19.2bit" -->
<!-- assemblyDanRer7Twobit <- "/Users/gtan/CSC/CNEr/2bit/danRer7.2bit" -->
<!-- finalCNE <- ceScanOneStep(axt1=axtHg19DanRer7, filter1=bedHg19, -->
<!--                           sizes1=qSizesHg19, assembly1="hg19",  -->
<!--                           twoBit1=assemblyHg19Twobit, -->
<!--                           axt2=axtDanRer7Hg19, filter2=bedDanRer7, -->
<!--                           sizes2=qSizesDanRer7, assembly2="danRer7", -->
<!--                           twoBit2=assemblyDanRer7Twobit, -->
<!--                           thresholds=c("45_50", "48_50", "49_50"), -->
<!--                           blatBinary="blat", blatCutoff1=8L, blatCutoff2=8L) -->
<!-- ``` -->

## CNE storage and query
As the computation of CNEs from the whole pipeline and 
the preparation of annotation package can be very time-consuming,
for a smoother visualisation experience, 
we decided to use a local SQLite database in order to store CNEs.

Since the CNEs `data.frame` is a table, it can be imported into a SQL table.
To speed up the query from the SQL database,
the bin indexing system is adopted.
For more information, please refer to the paper [@Kent:2002bw]
and [genomewiki](http://genomewiki.ucsc.edu/index.php/Bin_indexing_system).

```{r saveCNE, eval=TRUE, echo=TRUE}
## on individual tables
dbName <- tempfile()
data(cneFinalListDanRer10Hg38)
tableNames <- paste("danRer10", "hg38", names(cneFinalListDanRer10Hg38),
                    sep="_")
for(i in 1:length(cneFinalListDanRer10Hg38)){
	saveCNEToSQLite(cneFinalListDanRer10Hg38[[i]], dbName, tableNames[i],
	                overwrite=TRUE)
}
```

When querying results from the local SQLite database based on the chr, 
coordinates and other criterias, 
a `GRanges` object is returned.

```{r queryCNE, eval=TRUE, echo=TRUE}
chr <- "chr6"
start <- 24000000L
end <-  27000000L
minLength <- 50L
tableName <- "danRer10_hg38_45_50"
fetchedCNERanges <- readCNERangesFromSQLite(dbName, tableName, chr, 
                                            start, end, whichAssembly="first",
                                            minLength=minLength)
fetchedCNERanges
```

## CNE length distribution
As the lengths of CNEs [@Salerno:2006sc] and the distances between consecutive elements [@Polychronopoulos:2014pl] exhibit power-law distributions, we implemented a function that might be useful in showing interesting patterns in the distribution of the identified elements.
```{r CNEWidthDistribution, eval=TRUE, echo=TRUE}
dbName <- file.path(system.file("extdata", package="CNEr"),
                    "danRer10CNE.sqlite")
tAssemblyFn <- file.path(system.file("extdata",
                         package="BSgenome.Drerio.UCSC.danRer10"),
                         "single_sequences.2bit")
qAssemblyFn <- file.path(system.file("extdata",
                         package="BSgenome.Hsapiens.UCSC.hg38"),
                         "single_sequences.2bit")
cneGRangePairs <- readCNERangesFromSQLite(dbName=dbName, 
                                          tableName="danRer10_hg38_45_50",
                                          tAssemblyFn=tAssemblyFn,
                                          qAssemblyFn=qAssemblyFn)
plotCNEWidth(cneGRangePairs)
```

## Genomic distribution of CNEs along the chromosome
CNEs tend to form clusters. A quick check of the genomic distribution of CNEs
is available.

```{r plotCNEDistribution, eval=TRUE, echo=TRUE}
plotCNEDistribution(first(cneGRangePairs)) 
```

## Output of bed and bedGraph files
For visualisation in other Genome Browser, we provide functions to generate the CNE in bed files and CNE density in bedGraph files.
 For example, to get the first 1000 coordinates of CNEs:

```{r outputBedBW, eval=FALSE, echo=TRUE}
makeCNEDensity(cneGRangePairs[1:1000])
```

# CNEs visualisation
To visualise CNEs alongside other gene annotations, 
we choose to use the Bioconductor package `r Biocpkg("Gviz")` in this vignette.
`r Biocpkg("Gviz")`, based on the `r CRANpkg("grid")` graphics scheme, 
is a very powerful package for plotting data and annotation information
along genomic coordinates.
The functionality of integrating publicly available genome annotation data, 
such as UCSC or Ensembl,
significantly reduced the burden of preparing annotations for common assemblies.
Since the Bioconductor release 2.13 of `r Biocpkg("Gviz")`, 
it provides the data track in horizon plot,
which exactly meets our needs for visualisation of CNEs density plots.
For more detailed usage, please check the vignette or manual of `r Biocpkg("Gviz")`.

Another option for visualisation is the package `r Biocpkg("ggbio")`, 
which is based on `r CRANpkg("ggplot2")`.
The advantage of `r Biocpkg("ggbio")` is the simplicity of 
adding any customised `r CRANpkg("ggplot2")` style track
into the plot without tuning the coordinate systems.
The densities generated in the following section can be easily plot in the 
horizon plot.
A short straightforward tutorial regarding horizon plot 
in `ggplot2` format
is available from 
http://timelyportfolio.blogspot.co.uk/2012/08/horizon-on-ggplot2.html.

## Gene annotation visualisation
For the example case of hg38 vs danRer10 in this vignette, 
we choose danRer10 as the reference and 
show the range of developmental gene __barhl2__ and __sox14__. 

```{r queryUCSC, eval=FALSE, echo=TRUE, cache=TRUE}
library(Gviz)
library(biomaRt)
genome <- "danRer10"
axisTrack <- GenomeAxisTrack()
cpgIslands <- UcscTrack(genome=genome, chromosome=chr,
                        track="cpgIslandExt", from=start, to=end,
                        trackType="AnnotationTrack", start="chromStart",
                        end="chromEnd", id="name", shape="box",
                        showId=FALSE,
                        fill="#006400", name="CpG",
                        background.title="brown")
refGenes <- UcscTrack(genome=genome, chromosome=chr,
                      track="refGene", from=start, to=end,
                      trackType="GeneRegionTrack", rstarts="exonStarts",
                      rends="exonEnds", gene="name2", symbol="name2",
                      transcript="name", strand="strand", fill="#8282d2",
                      name="refSeq Genes", collapseTranscripts=TRUE,
                      showId=TRUE, background.title="brown")
ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL",
                   host="dec2015.archive.ensembl.org")
ensembl <-  useDataset("drerio_gene_ensembl",mart=ensembl)
biomTrack <- BiomartGeneRegionTrack(genome=genome, chromosome=chr, 
                                    biomart=ensembl,
                                    start=start , end=end, name="Ensembl Genes")
```

```{r loadAnnotation, eval=TRUE, echo=FALSE}
data(axisTrack)
data(cpgIslands)
data(refGenes)
```

```{r plotAnnotation, eval=TRUE, echo=TRUE, fig.height=5, fig.width=7}
library(Gviz)
plotTracks(list(axisTrack, cpgIslands, refGenes), 
           collapseTranscripts=TRUE, shape="arrow",
           transcriptAnnotation="symbol")
```

It is also possible to plot the annotation from an ordinary `R` object,
such as `data.frame`, `GRanges`, `IRanges`
or even from a local file.
Usually the __gff__ file containing the gene annotation can be processed by
`r Biocpkg("Gviz")` directly.
For more details, please look into the vignette of `r Biocpkg("Gviz")`.

## CNEs horizon plot

```{r CNEDensity, eval=TRUE, echo=TRUE}
dbName <- file.path(system.file("extdata", package="CNEr"),
                    "danRer10CNE.sqlite")
genome <- "danRer10"
windowSize <- 200L
minLength <- 50L
cneDanRer10Hg38_21_30 <- 
  CNEDensity(dbName=dbName, 
             tableName="danRer10_hg38_21_30",
             whichAssembly="first", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
cneDanRer10Hg38_45_50 <-
  CNEDensity(dbName=dbName, 
             tableName="danRer10_hg38_45_50", 
             whichAssembly="first", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
cneDanRer10Hg38_49_50 <-
  CNEDensity(dbName=dbName, 
             tableName="danRer10_hg38_49_50", 
             whichAssembly="first", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
cneDanRer10AstMex102_48_50 <-
  CNEDensity(dbName=dbName, 
             tableName="AstMex102_danRer10_48_50",
             whichAssembly="second", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
cneDanRer10CteIde1_75_75 <-
  CNEDensity(dbName=dbName, 
             tableName="cteIde1_danRer10_75_75", 
             whichAssembly="second", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
```

```{r GvizDataTrack, eval=TRUE, echo=TRUE}
dTrack1 <- DataTrack(range=cneDanRer10Hg38_21_30,
                     genome=genome, type="horiz", 
                     horizon.scale=max(cneDanRer10Hg38_21_30$score)/3, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="human 21/30", background.title="brown")
dTrack2 <- DataTrack(range=cneDanRer10Hg38_45_50,
                     genome=genome, type="horiz", 
                     horizon.scale=max(cneDanRer10Hg38_45_50$score)/2, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="human 45/50", background.title="brown")
dTrack3 <- DataTrack(range=cneDanRer10Hg38_49_50,
                     genome=genome, type="horiz", 
                     horizon.scale=max(cneDanRer10Hg38_21_30$score)/3, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="human 49/50", background.title="brown")
dTrack4 <- DataTrack(range=cneDanRer10AstMex102_48_50,
                     genome=genome, type="horiz", 
                     horizon.scale=max(cneDanRer10Hg38_21_30$score)/3, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="blind cave fish 48/50", background.title="brown")
dTrack5 <- DataTrack(range=cneDanRer10CteIde1_75_75,
                     genome=genome, type="horiz", 
                     horizon.scale=max(cneDanRer10CteIde1_75_75$score)/3, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="grass carp 75/75", background.title="brown")
```

```{r plotCNE, eval=TRUE, echo=TRUE, fig.height=10, fig.width=8}
ht <- HighlightTrack(trackList=list(refGenes, dTrack5, dTrack4, 
                                    dTrack1, dTrack2, dTrack3), 
                     start=c(24200000, 25200000, 26200000), 
                     end=c(25100000, 26150000, 27000000),
                     chromosome =chr)
plotTracks(list(axisTrack, cpgIslands, ht),
           collapseTranscripts=TRUE, shape="arrow",
           transcriptAnnotation="symbol",
           from=24000000, to=27000000)
```

From this horizon plot and when comparing Zebrafish with Human as a reference genome, 
we notice that the genes barhl2, lmo4b and sox14 are surrounded by the density peaks of CNEs.

# Conclusions
`r Biocpkg("CNEr")` efficiently identifies CNEs
and handles the corresponding objects conveniently in R.  
Horizon plot shows a superior dynamic range to the standard density plots, 
simultaneously revealing CNE clusters characterized 
by vastly different levels of sequence conservation. 
Such CNE density plots generated using precise locations of CNEs 
can be used to identify genes involved in developmental regulation, 
even for novel genes that are not yet annotated.

The following is the session info that generated this vignette:
```{r sessionInfo, eval=TRUE, echo=TRUE}
  sessionInfo()
```


# References