---
title: "Explore"
author: "Wolfgang Huber"
output: 
  BiocStyle::html_document:
    css: wh.css
---

<link rel="stylesheet" type="text/css" href="//fonts.googleapis.com/css?family=Istok+Web" />

```{r global_options, include = FALSE}
knitr::opts_chunk$set(cache = TRUE, autodep = TRUE)
```

# Load packages and the data we read in `load.Rmd`

```{r libs, message = FALSE, warning = FALSE}
library("magrittr")
library("dplyr")
library("reshape2")
library("ggplot2")
library("SummarizedExperiment")
library("pheatmap")
library("ggsci")
library("RColorBrewer")
library("umap")
library("DESeq2")
```

```{r load}
load("ses.RData")
ses <- ses1500
ses
```

`cells` is a `tibble` that contains annotation for each cell. It is derived from the `colData` of `ses` and will later also store the classification of the cells.

```{r markerhist} 
cells <- colData(ses)[, c("perturbation", "barcode")] %>% 
  as("data.frame") %>% as_tibble %>%
  mutate(id = paste(perturbation, barcode, sep = ":"),
         perturbation = as.factor(perturbation))
```
 
# Explore the tissue-specific markers

`markermat` is a matrix of dimensions (number of cells) x (number of tissues), and for each cell it contains the number of reads for *any* marker gene for the given tissue. `markerdf` is a long form version of the same data, suitable for `ggplot2`.

```{r markermat} 
markermat <- matrix(NA_real_, nrow = nrow(cells), ncol = ncol(efroni),
          dimnames = list(cell = cells$id, tissue = colnames(efroni)))

for (i in seq_len(ncol(efroni))) {
  genes = rowData(efroni)$Locus[assay(efroni)[, i]]
  genes = intersect(genes, rownames(ses))
  markermat[, i] = colSums(assay(ses)[genes,, drop = FALSE])
}
markerdf <- reshape2::melt(markermat)
```

Let's look at the histograms of these numbers. Generally, they are quite small. We will struggle to make reliable assignments.

```{r ggmarkerhist, fig.wide = TRUE, fig.width = 8, fig.height = 5}
ggplot(markerdf, aes(x = value)) + 
  geom_histogram(bins = 50) + facet_wrap( ~ tissue, ncol = 4)
ggplot(dplyr::filter(markerdf, value > 0), aes(x = log2(value))) +  
  geom_histogram(bins = 50) + facet_wrap( ~ tissue, ncol = 4)
```


```{r markermapheatmap, fig.width = 6, fig.height = 10, fig.wide = TRUE, fig.cap = "Heatmap of the matrix `markermat`."}
library("pheatmap")
pheatmap(markermat,
    cluster_rows = TRUE, cluster_cols = TRUE, show_rownames = FALSE,
    breaks = sinh(seq(0, 1, length.out = 101) * asinh(max(markermat))), 
    color = colorRampPalette(brewer.pal(9, "Blues"))(100))
```

# Normalize, transform

`DESeq2::estimateSizeFactorsForMatrix` is not applicable since there are no genes which have >0 counts in all cells. We use the `sqrt` rather than a more fancy log/asinh-like transformation since the counts are anyway so small. We'll cluster based on the top `topgenes` genes.

```{r estimateSizeFactors}
if (!file.exists("dds.RData")) {

  dds <- DESeqDataSet(ses, design = ~1)
  # debug(DESeq2:::estimateSizeFactorsIterate)
  # debug(DESeq2:::estimateDispersions.DESeqDataSet)
  # debg(DESeq2:::linearModelMu)
  # Rprof(  "Rprof.log", interval = 0.1)
  dds %<>% estimateSizeFactors(type = "iterate")
  save(dds, file = "dds.RData")
} else {
  load("dds.RData")
}
```

```{r sf, fig.width = 6, fig.height = 6, fig.wide = TRUE, fig.cap = "Size factors."}
sumNorm <- function(x) { stopifnot(is.matrix(x)); cs <- colSums(x); x / (cs / mean(cs)) }
ggplot(tibble(sum = sumNorm(counts(dds)), DESeq2 = sizeFactors(dds)), 
       aes(x = sum, y = DESeq2)) + geom_point()
```

```{r trsf}
sesT <- varianceStabilizingTransformation(dds)
```

# Subset to top genes

```{r topgenes}
topgenes <- 1000
ktop <- order(rowMeans(sesT), decreasing = TRUE)[seq_len(topgenes)]

sesT <- ses[itop,]
```

# UMAP

```{r umap}
set.seed(0xdada)
um <- umap(t(assay(sesT)), n_neighbors = 10)
```
Seems to look better this way
```{r flip}
um$layout <- (-um$layout)
```
```{r plotumap1, fig.wide = TRUE, fig.width = 8, fig.height = 8, fig.cap = sprintf("UMAP layout of %d cells, colored by treatment.", ncol(ses))}
ggplot(as_tibble(um$layout) %>% cbind(cells),
       aes(x = V1, y = V2, col = perturbation)) + geom_point(size = 1.5) +
       coord_fixed() + xlab("UMAP 1") + ylab("UMAP 2") + 
       scale_colour_manual(values = 
         setNames(brewer.pal(8, "Accent")[c(8, 5)], c("control", "mTOR"))) 
```

# Classify cells by marker expression

Questions we would like to ask: do the cells (i.e. their data) form meaningful clusters? How do the  corresponding cells from the mTOR-treated root differ (or diverge) from the untreated?

```{r class}
min_count_in_any_marker <- 8
min_frac <- 0.5
cells$class <- ordered("unassigned", 
      levels = c("unassigned", colnames(markermat)))
rs <- rowSums(markermat)
for (wr in which(rs >= min_count_in_any_marker)) {
  wc <- which.max(markermat[wr, ])
  if (markermat[wr, wc] >= min_frac * rs[wr])
    cells$class[wr] <- colnames(markermat)[wc]
}
table(cells$class)
```


The `xtfrm2` function will be useful for ordering the cells so that the grey points are plotted first, and the colored ones later (and thus, on top.)

```{r umap2a}
xtfrm2 <- function(x) ifelse(is.na(x), 1, xtfrm(x) + 1)
```

Remove unused classes and provide color map.

```{r colores}
cells$class %<>% ordered
colores = setNames(
  c("#d0d0d0", ggsci::pal_ucscgb()(nlevels(cells$class)-1)), 
  levels(cells$class))
```


```{r plotumap2, fig.wide = TRUE, fig.width = 12, fig.height = 8}
as_tibble(um$layout) %>% cbind(cells) %>% arrange(xtfrm2(class)) %>%
ggplot(aes(x = V1, y = V2, shape = perturbation, col = class)) + 
  geom_point(size = 1.5) +
  scale_colour_manual(values = colores) +
  scale_shape_manual(values=c(control=16, mTOR=17)) +
  coord_fixed() + xlab("UMAP 1") + ylab("UMAP 2") 
```


# Session info

```{r sesssion}
devtools::session_info()
```