---
title: "CNE identification and visualisation"
author: "Ge Tan"
date: "`r doc_date()`"
package: "`r pkg_ver('CNEr')`"
abstract: >
  Identification of Conserved Nocoding Elements (CNEs) from _Axt_ alignment
  and visualisation of CNE density in horizon plot.
  
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 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 for novel genes that are not annotated yet.

# Pipeline of the package
This section will briefly demonstrate the pipeline of CNE identification 
and visualisation.
More detailed usage of each step will be described in following sections 
with a concreate example of CNE identification and visualisation 
for Human (hg19) and Zebrafish (danRer7).

## CNE identification

1. axtNets: The axtNet files can be downloaded from UCSC 
  or generated by itself
2. scan alignments: The regions with minimal __I__ identities 
  over __C__ columns.
3. remove elements: Elements that overlap with annotated exons/repeats.
4. merge elements to get CNEs: Elements that overlap on both genomes.

## CNE visualisation
1. display parameters: Chromosome, start, end, smooth window size.
2. horizon plot: Visualise CNE densities.

# Input
`r Biocpkg("CNEr")` starts from 
[axt](http://genome.ucsc.edu/goldenPath/help/axt.html)
net files of two genomes pairwise alignments and _bed_ files for filtering.
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.

## axt file
By the time I wrote this package, there was no suitable class to store the axt files in Bioconductor.
Hence we created a __S4__ class `Axt` to hold the content from _axt_ files.
Basically, it utilises `GRanges` from `r Biocpkg("GenomicRanges")` package
and `DNAStringSet` from `r Biocpkg("Biostrings")` package.

To read axt file into R, 
`r Biocpkg("CNEr")` provides `readAxt` function
for highly efficient reading, 
which heavily depends on Kent's utilities source code [@Kent:2002bw].
The axt file can be gzippped or in plain text file.
The alignments between two genomes can also be in one 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)
axtFilesHg19DanRer7 <- file.path(system.file("extdata", package="CNEr"), 
                                "hg19.danRer7.net.axt")
axtFilesDanRer7Hg19 <- file.path(system.file("extdata", package="CNEr"), 
                                "danRer7.hg19.net.axt")
```

```{r readAxt, eval=FALSE, echo=TRUE}
axtHg19DanRer7 <- readAxt(axtFilesHg19DanRer7)
axtDanRer7Hg19 <- readAxt(axtFilesDanRer7Hg19)
```

```{r loadAxt, eval=TRUE, echo=FALSE}
data(axtHg19DanRer7)
data(axtDanRer7Hg19)
```

```{r showAxt, eval=TRUE, echo=TRUE}
axtHg19DanRer7
axtDanRer7Hg19
```

## filtering information
The gene annotation information, including exons coordinates, 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 
hg19        |  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 
danRer7     |  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 rmak table from UCSC
library(rtracklayer)
mySession <- browserSession("UCSC")
genome(mySession) <- "hg19"
hg19.rmsk <- getTable(ucscTableQuery(mySession, track="RepeatMasker",
               table="rmsk"))
hg19.rmskGRanges <- GRanges(seqnames=hg19.rmsk$genoName,
                            ## The UCSC coordinate is 0-based.
                            ranges=IRanges(start=hg19.rmsk$genoStart+1,
                                           end=hg19.rmsk$genoEnd),
                            strand=hg19.rmsk$strand)
## To fetch ensembl gene exons from BioMart
library(biomaRt)
ensembl <- useMart("ensembl", host="feb2014.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, "+", "-")
                       )
## Use the existing Bioconductor annotation package
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
exonsRanges <- exons(TxDb.Hsapiens.UCSC.hg19.knownGene)
```

Of course, more customised gene annotation can be included in the bed file.
To import bed file into `r Biocpkg("GRanges")` in `R`, 
`r Biocpkg("rtracklayer")` provides a general function 
`import.bed` to do that.
However, if your bed file only contains 3 columns: chromosome, start and end,
that are the information needed by `r Biocpkg("CNEr")`,
then it is better to use the `readBed` function that is much faster with C implementation.

```{r Bed, eval=TRUE, echo=TRUE}
bedHg19Fn <- file.path(system.file("extdata", package="CNEr"), 
                      "filter_regions.hg19.bed")
bedHg19 <- readBed(bedHg19Fn)
bedHg19
bedDanRer7Fn <- file.path(system.file("extdata", package="CNEr"), 
                         "filter_regions.danRer7.bed")
bedDanRer7 <- readBed(bedDanRer7Fn)
bedDanRer7
```

## Chromosome size of query assembly
The chromosome size of the query assembly is necessary 
when the filters for query assembly is provided.

To facilitate the fetch of chromosome sizes,
we prepared a function `fetchChromSizes` to automate the procedure.
Currently the assemblies from UCSC are supported,
and more assemblies from other sources will be implemented in the future.
A object of `Seqinfo` is returned 
which suits the input to downstream analysis.

```{r chromSizes, eval=FALSE, echo=TRUE}
qSizesHg19 <- fetchChromSizes("hg19")
qSizesDanRer7 <- fetchChromSizes("danRer7")
```

```{r chromSizesData, eval=TRUE, echo=FALSE}
data(qSizesHg19)
data(qSizesDanRer7)
```

```{r showchromSizesData, eval=TRUE, echo=TRUE}
qSizesHg19
qSizesDanRer7
```

`Seqinfo` can also be generated from local two bit file 
with `seqinfo` from 
`r Biocpkg("rtracklayer")`.

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

## net alignments scan
Detecting CNEs highly relies on the whole-genome pairwise net alignments.
To correct the bias of a chosen genome 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}
## axt, GRanges objects as input
CNEHg19DanRer7 <- ceScan(axts=axtHg19DanRer7, tFilter=bedHg19,
                         qFilter=bedDanRer7, qSizes=qSizesDanRer7,
                         thresholds=c("45_50", "48_50", "49_50"))
CNEDanRer7Hg19 <- ceScan(axts=axtDanRer7Hg19, tFilter=bedDanRer7,
                         qFilter=bedHg19, qSizes=qSizesHg19,
                         thresholds=c("45_50", "48_50", "49_50"))

## axt and bed files as input
CNEHg19DanRer7 <- ceScan(axts=axtFilesHg19DanRer7, tFilter=bedHg19Fn,
                         qFilter=bedDanRer7Fn, qSizes=qSizesDanRer7,
                         thresholds=c("45_50", "48_50", "49_50"))
CNEDanRer7Hg19 <- ceScan(axts=axtFilesDanRer7Hg19, tFilter=bedDanRer7Fn,
                         qFilter=bedHg19Fn, qSizes=qSizesHg19,
                         thresholds=c("45_50", "48_50", "49_50"))
```

At this stage, a list of `data.frame` is returned 
from `ceScan`,
which contains the preliminary CNEs.
Here is some exemple output:

```{r CNEScanHead, eval=TRUE, echo=TRUE}
lapply(CNEHg19DanRer7, head)
```

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

## Merge CNEs
Because we do two rounds of CNE detections with each genome as reference,
some conserved elements overlap on both genomes and are supposed to be removed.
But elements that only overlap on one of the genomes are kept, 
so that duplicated elements remain distinct.

```{r CNEMerge, eval=TRUE, echo=TRUE}
cneMergedDanRer7Hg19 <- mapply(cneMerge, CNEDanRer7Hg19, CNEHg19DanRer7, 
                               SIMPLIFY=FALSE)
lapply(cneMergedDanRer7Hg19, head)
```

## 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 exceed certain threshold, for instance, 8,
that CNE will be discarded.

This step can be very time-consuming when the number of CNEs are large.
Other alignment method can also be considerd, such as Bowtie2, BWA.
The two bit file for each assembly is required.

```{r CNERealignment, eval=FALSE, echo=TRUE}
assemblyHg19Twobit <- "/Users/gtan/CSC/CNEr/2bit/hg19.2bit"
assemblyDanRer7Twobit <- "/Users/gtan/CSC/CNEr/2bit/danRer7.2bit"
cneBlatedDanRer7Hg19 <- list()
for(i in 1:length(cneMergedDanRer7Hg19)){
  cneBlatedDanRer7Hg19[[names(cneMergedDanRer7Hg19)[i]]] <-
    blatCNE(cneMergedDanRer7Hg19[[i]], 
            as.integer(sub(".+_.+_\\d+_", "", names(cneMergedDanRer7Hg19)[i])), 
            cutoffs1=8L, cutoffs2=8L, 
            assembly1Twobit=assemblyDanRer7Twobit, 
            assembly2Twobit=assemblyHg19Twobit,
            blatBinary="blat")
}
```

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 to store these information.

## CNE storage and query
Since the CNEs `data.frame` is just a table and 
can be imported into a SQL table naturally.
To speed up the query from the SQL database,
the bin indexing system is acquired.
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(cneBlatedDanRer7Hg19)
for(i in 1:length(cneBlatedDanRer7Hg19)){
  tableName <- paste("danRer7_hg19", names(cneBlatedDanRer7Hg19)[i],
                      sep="_")
	saveCNEToSQLite(cneBlatedDanRer7Hg19[[i]], dbName, tableName, overwrite=TRUE)
}

## on CNE class
data(finalCNE)
saveCNEToSQLite(finalCNE, dbName=dbName, overwrite=TRUE)
```

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

```{r queryCNE, eval=TRUE, echo=TRUE}
chr <- "chr11"
start <- 31000000L
end <-  33000000L
minLength <- 50L
tableName <- "danRer7_hg19_45_50"
fetchedCNERanges <- readCNERangesFromSQLite(dbName, tableName, chr, 
                                           start, end, whichAssembly="L",
                                           minLength=minLength)
```


# CNEs visualisation
To visualise the CNEs together with 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 of hg19 vs danRer7 in this vignette, 
we choose hg19 as the reference and 
show the range of developmental gene __PAX6__.

```{r queryUCSC, eval=FALSE, echo=TRUE}
library(Gviz)
genome <- "hg19"
chr <- "chr11"
start <- 31000000L
end <-  33000000L
axisTrack <- GenomeAxisTrack()
ideoTrack <- IdeogramTrack(genome=genome, chromosome=chr)
cpgIslands <- UcscTrack(genome=genome, chromosome=chr,
                       track="cpgIslandExt", from=start, to=end,
                       trackType="AnnotationTrack", start="chromStart",
                       end="chromEnd", id="name", shape="box",
                       fill="#006400", name="CpG Islands")
refGenes <- UcscTrack(genome="hg19", 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)
biomTrack <- BiomartGeneRegionTrack(genome="hg19", chromosome=chr, 
                                   start=start , end=end, name="Ensembl")
```

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

```{r plotAnnotation, eval=TRUE, echo=TRUE, fig.height=10, fig.width=7.5}
library(Gviz)
plotTracks(list(axisTrack, ideoTrack, cpgIslands, refGenes), 
           collapseTranscripts=TRUE, shape="arrow",
           showId=TRUE, 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"),
                    "cne.sqlite")
windowSize <- 300L
minLength <- 50L
cneHg19DanRer7_45_50 <- 
  CNEDensity(dbName=dbName, 
             tableName="danRer7_hg19_45_50", 
             assembly1="hg19", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
cneHg19DanRer7_48_50 <-
  CNEDensity(dbName=dbName, 
             tableName="danRer7_hg19_48_50", 
             assembly1="hg19", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
cneHg19DanRer7_49_50 <-
  CNEDensity(dbName=dbName, 
             tableName="danRer7_hg19_49_50", 
             assembly1="hg19", chr=chr, start=start,
             end=end, windowSize=windowSize, 
             minLength=minLength)
```

```{r GvizDataTrack, eval=TRUE, echo=TRUE}
data(cneHg19DanRer7_45_50)
data(cneHg19DanRer7_48_50)
data(cneHg19DanRer7_49_50)
genome <- "hg19"
chr <- "chr11"
start <- 31000000L
end <-  33000000L
#axisTrack = GenomeAxisTrack()
#ideoTrack = IdeogramTrack(genome=genome, chromosome=chr)
strand <- "+"
dataMatrix <- cneHg19DanRer7_45_50
dTrack1 <- DataTrack(start=dataMatrix[ ,1], end=dataMatrix[ ,1], 
                     data=dataMatrix[ ,2], chromosome=chr, strand=strand, 
                     genome=genome, type="horiz", horizon.scale=0.1, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="danRer7 45/50")
dataMatrix <- cneHg19DanRer7_48_50
dTrack2 <- DataTrack(start=dataMatrix[ ,1], end=dataMatrix[ ,1], 
                     data=dataMatrix[ ,2], chromosome=chr, strand=strand, 
                     genome=genome, type="horiz", horizon.scale=0.1, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="danRer7 48/50")
dataMatrix <- cneHg19DanRer7_49_50
dTrack3 <- DataTrack(start=dataMatrix[ ,1], end=dataMatrix[ ,1], 
                     data=dataMatrix[ ,2], chromosome=chr, strand=strand, 
                     genome=genome, type="horiz", horizon.scale=0.1, 
                     fill.horizon=c("#B41414", "#E03231", "#F7A99C", 
                                    "yellow", "orange", "red"), 
                     name="danRer7 49/50")
```

```{r plotCNE, eval=TRUE, echo=TRUE, fig.height=15, fig.width=7.5}
plotTracks(list(axisTrack, ideoTrack, cpgIslands, 
                refGenes, 
                dTrack1, dTrack2, dTrack3),
           collapseTranscripts=TRUE, shape="arrow",
           showId=TRUE, transcriptAnnotation="symbol")
```

From this horizon plot compared with Zebrafish with Human as reference genome, 
the developmental gene PAX6 was surrounded by the density peaks of CNEs.

# Conclusion
With this package, we are able to identify CNEs efficiently 
and handle 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
