---
title: "trisect.Rd"
output: html_document
---

### Set up

```{r, message=FALSE}
library('gbm')
library('ROCR')
library("kohonen")
library('gplots')
library('RColorBrewer')
library("biomaRt")
library("ChIPpeakAnno")
library("org.Hs.eg.db")
```


### Loading source files and data

In this vignette, we provide an example dataset for ZNF143 transcription factor in the form of a `list`, which stores data `list` of four cell types: `Gm12878`, `H1hesc`, `Helas3`, and `K562`. For each cell type, the `list` includes four objects: `trainseq`, `testseq`, `training`, and `testing`. The `training` and `testing` matrix each stores sequences in the rows and features in the columns, with an additional column of sequence labels. 


Load EMT functions.

```{r}
source('R/emt.r') 
source('R/clustering.r')
source('R/coherence.r')
```

Load ZNF143 data.

```{r}
tf <- 'ZNF143'
data("ZNF143")
cells <- names(datalist)
print(cells)
```

The `trainseq` columns are : 

1. `sequence`: 100 base pair sequence.
2. `seqnames`: Name of the chromosome
3. `start`: The starting position of the binding site in the chromosome
4. `end`: The ending position of the binding site in the chromosome
5. `score`: Signal value
6. `value`: The overall measurement of the enrichment in the region; for non-binding sites it is the value of associated binding site.
7. `peakid`: Id.
8. `label`: Indicator of binding or non-binding.


```{r}
head(datalist$Gm12878$trainseq)
```

The `train` columns are 

1. `label`: Indicator of binding or non-binding.
2. `M00001`: TRANSFAC id (i.e. feature 1)
3. `M00002`: TRANSFAC id (i.e. feature 2)
4. `M00003`: TRANSFAC id (i.e. feature 3)
5. `M00004`: TRANSFAC id (i.e. feature 4)


```{r}
datalist$Gm12878$train[1:5, 1:5]
```


### Building model for each cell

Note that this steps takes considerable amount of time unless multiple cores are used.

```{r, eval = FALSE}
models <- lapply(datalist, function(x) build_emt(x))
names(models) <- cells
```

For the purpose of this vignette, we pre-computed the model and stored the results in an R object which can be easily loaded as follows.

```{r}
data("models")
sapply(cells, function(x) models[[x]]$auc) 
```

The number of sub-models for each of the four cell-type:

```{r}
sapply(cells, function(x) models[[x]]$model$n.tree)
```

Let's create an object that stores the sub-model results only.

```{r}
submodels <- lapply(models, function(x) x$model)
```


### Clustering the submodels across cells

Summarize the cluster membership matrix and plot the constituent cell-specific sub-model counts. In the figure, each row denotes a cluster number and each column denotes a cell line. Each value of the matrix gives the number of sub-models coming from a cell line and belonging to a cluster.


```{r, message=FALSE, fig.width = 8, fig.height = 8}
ld <- get_cluster_membership(submodels, clen = 16)
fit <- ld$fit
cluster.membership <- ld$cluster.membership
plot_cluster_membership(cluster.membership, tf)
```


Make a new ensemble object for each cluster.

```{r}
newEnsembles <- make_cluster_ensembles(submodels, fit, clen = 16)
head(sapply(newEnsembles, class))
```


### Functional assesment of the target genes

Get the target genes (i.e. determined by the nearest gene from a binding site which belong to the cluster) of each cell type from each cluster. 

```{r, message=FALSE}
targets <- get_targets(datalist, newEnsembles, clen = 16)
```


How many clusters have expression coherence?

```{r}
data("expression")
coherence <- setup_exp_coherence(targets, exprsn, exptheK = 1, verbose = F)
df <- exp_coherence(coherence)
sum(df$odds > 1 & df$p.value < 0.05, na.rm = T)
head(df) #clusters with single-celll and/or insufficient data have NA in the corresponding row
```


How many clusters have pathway coherence?

```{r}
data("keggmat")
coherence2 <- setup_pathway_coherence(targets$entrezmat, pathwaymat, verbose=F)
df <- pathway_coherence(coherence2)
sum(df$odds > 1 & df$p.value < 0.05, na.rm = T) 
head(df) #clusters with single-celll and/or insufficient data have NA in the corresponding row

```



### Session information


```{r}
sessionInfo()
```
