result_AD.rmd

---
title: "Counterfactual Confounder Adjustment for Single-cell Differential Expression Analysis"
author: "Yongjin Park"
---

```{r include=FALSE}
library(tidyverse)
library(data.table)
library(patchwork)
library(R.utils)
source("R/Util-rmd.R")

# Suppress summarise info
options(dplyr.summarise.inform = FALSE)

## document options
fig.dir = 'Fig/AD/'
dir.create(fig.dir, recursive = TRUE, showWarnings = FALSE)
knitr::opts_chunk$set(fig.path = fig.dir, fig.width = 6, fig.height = 4)
knitr::opts_chunk$set(dev.args = list(bg = "transparent"), echo = FALSE)
options(stringsAsFactors = FALSE)

## helper functions
.fread <- function(...) fread(..., header=FALSE)
.fwrite <- function(...) fwrite(..., sep="\t", col.names=FALSE, row.names=FALSE)
.read.mat <- function(...) .fread(...) %>% as.matrix
```

```{r}
.colors <- fread("data/brain_colors.txt") %>%
    mutate(hex = as.character(hex))
```

```{r include=FALSE}
fit.umap <- function(V, .file, rank = 10, refit = FALSE) {
    .mkdir(dirname(.file))

    if(!file.exists(.file) || refit) {
        .mat <- as.matrix(V)
        .mat <- .mat[, 1:min(ncol(.mat), rank), drop = FALSE]
        .umap <- uwot::umap(.mat,
                            fast_sgd = TRUE,
                            n_threads = 10,
                            spread = 5,
                            verbose = TRUE) %>%
            as.data.table
        .fwrite(.umap, .file)
    }
    .ret <- .fread(.file)
    colnames(.ret) <- "UMAP " %&% 1:ncol(.ret)
    return(.ret)
}

.meta <-
    fread("data/brain_2018-05-03/filtered_column_metadata.txt.gz") %>%
    mutate(temp = TAG) %>%
    separate("temp", c("barcode","batch"), sep = "[.]")

.read.annot <- function(...) {
    .fread(...,
           col.names=c("TAG", "ct", "pr", "llik")) %>%
        left_join(.meta, by = "TAG") %>% 
        na.omit %>% 
        as.data.table
}

read.annot.consensus <- function(.annot.1, .annot.2) {
    .annot.dt <- .read.annot(.annot.1)
    .annot.2.dt <- .read.annot(.annot.2)

    .valid <- rbind(.annot.dt[, .(TAG, ct, pr)], .annot.2.dt[, .(TAG, ct, pr)]) %>%
        (function(x) { ret <- x[, .(n = .N, min.pr = min(pr)), by = .(TAG, ct)]; ret[n > 1] }) %>%
        select(TAG, ct, min.pr)

    return(left_join(.valid, .annot.dt, by = c("TAG", "ct")))
}
```

## Cell type annotation

```{r include=FALSE}
plot.five.umap <- function(.dt) {

    nn <- nrow(.dt)

    p1 <-
        .gg.plot(.dt[sample(nn)], aes(x = `UMAP 1`, y = `UMAP 2`, colour = factor(msex, c(0, 1), c("F","M")))) +
        geom_point(pch = 19, stroke = 0, size = .5) +
        scale_colour_manual("Sex", values = c("#ef8a62","#67a9cf"))

    p2 <-
        .gg.plot(.dt[sample(nn)], aes(x = `UMAP 1`, y = `UMAP 2`, colour = factor(age.death, 0:2, c("(, 80]","(80, 90]","(90, )")))) +
        geom_point(pch = 19, stroke = 0, size = .5) +
        scale_colour_brewer("Age", direction=1, palette="Greys")

    p3 <- 
        .gg.plot(.dt[sample(nn)], aes(x = `UMAP 1`, y = `UMAP 2`, colour = factor(cog, c(0,1), c("Decline", "Incline")))) +
        geom_point(pch = 19, stroke = 0, size = .5) +
        scale_colour_manual("Cognition", values=c("#2c7fb8","#edf8b1"))

    p4 <- 
        .gg.plot(.dt[sample(nn)], aes(x = `UMAP 1`, y = `UMAP 2`, colour = factor(apoe.e4, c(0,1), c("WT", "E4")))) +
        geom_point(pch = 19, stroke = 0, size = .5) +
        scale_colour_manual("APOe", values = c("#e7e1ef","#dd1c77"))

    p5 <- 
        .gg.plot(.dt[sample(nn)][pathoAD >= 0], aes(x = `UMAP 1`, y = `UMAP 2`, colour = factor(pathoAD, c(0,1), c("HC", "AD")))) +
        geom_point(pch = 19, stroke = 0, size = .5) +
        scale_colour_manual("pathoAD", values = c("#e0ecf4","#8856a7"))

    plt <- p1/p2/p3/p4/p5
}
```

### Broad cell types

```{r fig.width=4, fig.height=8}
.annot.dt <- read.annot.consensus("result/celltype.annot.gz", "result/bbknn_celltype.annot.gz")
.pheno.dt <- fread("result/phenotyped.txt.gz")

tot.umap.dt <- .fread("result/bbknn/Total.factors.gz") %>%
    fit.umap("result/umap/Total.umap.gz") %>%
    cbind(.fread("result/bbknn/Total.cols.gz", col.names="TAG")) %>%
    left_join(.annot.dt, by = "TAG") %>%
    left_join(.pheno.dt, by = c("TAG", "projid")) %>% 
    na.omit

.hex <- tot.umap.dt %>%
    select(ct) %>%
    unique %>%
    left_join(.colors, by = "ct") %>%
    arrange(ct)

.dt <- tot.umap.dt %>%
    mutate(ct = factor(ct, .hex$ct, .hex$ct.name))

p1 <-
    .gg.plot(.dt, aes(x = `UMAP 1`, y = `UMAP 2`, colour = ct)) +
    scale_colour_manual(values = .hex$hex) +
    geom_point(pch = 19, stroke = 0, size = .5)

plt <- (p1 / plot.five.umap(.dt)) + plot_layout(heights=c(1,5))
print(plt)
.file <- fig.dir %&% "/Fig_umap_broad.pdf"
.gg.save(.file, plot=plt, width=4, height=8)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

```{r Fig_Marker_Weighted, fig.width=6, fig.height=2}
suppressMessages(library(cba))
suppressMessages(library(proxy))

.labels <- .fread("result/bbknn_celltype.label_names.gz") %>% unlist
.genes <- .fread("result/bbknn_celltype.marker_names.gz") %>% unlist

.dt <-
    .fread("result/bbknn_celltype.marker_profile.gz", col.names=.labels) %>%
    mutate(geneSymbol = .genes) %>%
    melt(id.vars = "geneSymbol", value.name = "mu", variable.name = "celltype") %>%
    filter(celltype %in% unique(.annot.dt$ct)) %>% 
    mutate(col = geneSymbol, row = celltype, weight = mu) %>% 
    order.pair(ret.tab=TRUE)

.marker.order <- sort(unique(.dt$geneSymbol))

plt <- 
    .gg.plot(.dt, aes(y=row, x=col, fill = sqrt(pmax(mu, 0)))) +
    theme(legend.position = "top") +
    geom_tile() +
    scale_fill_distiller("", palette="RdPu", direction=1) +
    xlab(num.int(length(unique(.dt$geneSymbol))) %&% " genes present in snRNA-seq") +
    ylab("cell types") + 
    theme(axis.ticks.x=element_blank()) +
    theme(axis.text.x=element_blank())

print(plt)
.file <- fig.dir %&% "/Fig_Marker_Weighted.pdf"
.gg.save(.file, plot=plt, width=6, height=2)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

#### Basic statistics for cell type annotation

```{r include = FALSE}
.take.annot.freq <- function(.annot.dt) {
    .tot.i.freq.dt <-
        .annot.dt[, .(n = .N), by = .(projid)]

    .tot.ct.freq.dt <-
        .annot.dt[, .(n = .N), by = .(ct)]

    .ct.order <-
        .tot.ct.freq.dt %>%
        arrange(n) %>%
        select(ct) %>%
        left_join(.colors, by = "ct")

    .ct.freq.dt <- .annot.dt[, .(n = .N), by = .(ct, projid)]

    .max.ct <- .ct.freq.dt[which.max(.ct.freq.dt$n), .(ct)] %>% unlist

    .projid <-
        .ct.freq.dt %>%
        filter(ct == .max.ct) %>%
        arrange(desc(n)) %>%
        select(projid) %>%
        unlist

    .ct.freq.dt <- .ct.freq.dt %>%
        mutate(projid = factor(projid, .projid)) %>%
        mutate(ct = factor(ct, .ct.order$ct, .ct.order$ct.name)) %>%
        as.data.table

    .ct.freq.dt[, pr := n/sum(n), by = .(projid)]

    .avg.ct.dt <-
        .ct.freq.dt[,
                    .(.mean = mean(n), .sd = sd(n), .max = max(n), .min = min(n),
                      .mean.f = mean(pr), .sd.f = sd(pr), .max.f = max(pr), .min.f = min(pr)),
                    by = .(ct)] %>%
        mutate(str = num.round(.mean.f*100)  %&% "% (N=" %&% num.int(round(.mean)) %&% " $\\pm$ " %&% round(.sd) %&% ") cells for " %&% ct)

    list(tot = .tot.i.freq.dt, avg.ct = .avg.ct.dt, ct.freq.dt = .ct.freq.dt, nind = length(.projid), ct.order = .ct.order)
}
```

```{r Fig_CT_Freq, fig.width=5, fig.height=4}
freq <- .take.annot.freq(.annot.dt)

p2 <- 
    .gg.plot(freq$ct.freq.dt, aes(x = projid, y = pr, fill = ct)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks.x = element_blank()) +
    xlab(freq$nind %&% " individuals") +
    ylab("Cell type fraction") +
    scale_fill_manual(values = freq$ct.order$hex)

p1 <- 
    .gg.plot(freq$ct.freq.dt, aes(x = projid, y = n, fill = ct)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks.x = element_blank()) +
    xlab(freq$nind %&% " individuals") +
    ylab("#Cells") +
    scale_fill_manual(values = freq$ct.order$hex)

plt <- p1 / p2

print(plt)
.file <- fig.dir %&% "/Fig_CT_Freq.pdf"
.gg.save(.file, plot=plt, width = 5, height = 4)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

* `r num.int(round(mean(freq$tot$n)))` cells per individual: `r paste(freq$avg.ct$str, collapse=", ")` on average.

#### Correlation with the annotations of Mathys _et al._

```{r Fig_cor_with_mathys, fig.width=4.5, fig.height=3}
.cor.freq.dt <- .annot.dt[, .(n = .N), by = .(ct, broad.cell.type)] %>% rename(mathys = broad.cell.type)
.cor.freq.dt[, p := n/sum(n), by = .(ct)]
.cor.freq.dt[, q := n/sum(n), by = .(mathys)]
.cor.freq.dt[, ct := factor(ct, rev(c("Ex", "In", "Oligo", "OPC", "Microglia", "Astro", "Per", "Endo")))]
.cor.freq.dt[, mathys := factor(mathys, c("Ex", "In", "Oli", "Opc", "Mic", "Ast", "Per", "End"))]

plt <-
    .gg.plot(.cor.freq.dt, aes(y = ct, x = mathys, fill = p)) +
    scale_x_discrete(position = "top") +
    scale_fill_distiller("P(Mathys|vMF)", palette = "Blues", direction=1) +
    geom_tile() +
    geom_text(aes(label = num.int(n)), size = 2) +
    xlab("Mathys et al.") + ylab("This work (von Mises Fisher)")

print(plt)
.file <- fig.dir %&% "/Fig_cor_with_mathys.pdf"
.gg.save(.file, plot=plt, width = 4.5, height = 3)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```


### The cortical layers of excitatory neurons

```{r fig.width=3.2, fig.height=8}
.Ex.annot.dt <- .read.annot("result/subtype/bbknn_Ex.annot.gz")

Ex.umap.dt <- .fread("result/bbknn/Ex.factors.gz") %>%
    fit.umap("result/umap/Ex.umap.gz") %>%
    cbind(.fread("result/bbknn/Ex.cols.gz", col.names="TAG")) %>%
    left_join(.Ex.annot.dt, by = "TAG") %>%
    left_join(.pheno.dt, by = c("TAG", "projid")) %>% 
    na.omit

nn <- nrow(Ex.umap.dt)

p1 <-
    .gg.plot(Ex.umap.dt[sample(nn)], aes(x = `UMAP 1`, y = `UMAP 2`, fill = ct)) +
    geom_point(pch = 21, stroke = .05, size = 1) +
    scale_fill_brewer(palette="PiYG")

plt <- (p1 / plot.five.umap(Ex.umap.dt)) + plot_layout(heights=c(1,5))
print(plt)
.file <- fig.dir %&% "/Fig_umap_Ex.pdf"
.gg.save(.file, plot=plt, width = 3.2, height = 8)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```


```{r Fig_Ex_CT_Freq}
freq <- .take.annot.freq(.Ex.annot.dt)

p2 <- 
    .gg.plot(freq$ct.freq.dt, aes(x = projid, y = pr, fill = ct)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks.x = element_blank()) +
    xlab(freq$nind %&% " individuals") +
    ylab("Cell type fraction") +
    scale_fill_manual(values = freq$ct.order$hex)

p1 <- 
    .gg.plot(freq$ct.freq.dt, aes(x = projid, y = n, fill = ct)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks.x = element_blank()) +
    xlab(freq$nind %&% " individuals") +
    ylab("#Cells") +
    scale_fill_manual(values = freq$ct.order$hex)

plt <- p1 / p2

print(plt)
.file <- fig.dir %&% "/Fig_Ex_CT_Freq.pdf"
.gg.save(.file, plot=plt, width = 5, height = 4)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

* `r num.int(round(mean(freq$tot$n)))` cells per individual: `r paste(freq$avg.ct$str, collapse=", ")` on average.

### Inhibitory neuron subtypes

```{r fig.width=3, fig.height=8}
.In.annot.dt <- .read.annot("result/subtype/bbknn_In.annot.gz")

In.umap.dt <- .fread("result/bbknn/In.factors.gz") %>%
    fit.umap("result/umap/In.umap.gz") %>%
    cbind(.fread("result/bbknn/In.cols.gz", col.names="TAG")) %>%
    left_join(.In.annot.dt, by = "TAG") %>%
    left_join(.pheno.dt, by = c("TAG", "projid")) %>% 
    na.omit

nn <- nrow(In.umap.dt)

p1 <- 
    .gg.plot(In.umap.dt[sample(nn)], aes(x = `UMAP 1`, y = `UMAP 2`, fill=ct)) +
    geom_point(pch = 21, stroke = .05, size = 1) +
    scale_fill_brewer(palette="YlGnBu")

plt <- (p1 / plot.five.umap(In.umap.dt)) + plot_layout(heights=c(1,5))
print(plt)
.file <- fig.dir %&% "/Fig_umap_In.pdf"
.gg.save(.file, plot=plt, width = 3, height = 8)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```


```{r Fig_In_CT_Freq}
freq <- .take.annot.freq(.In.annot.dt)

p2 <- 
    .gg.plot(freq$ct.freq.dt, aes(x = projid, y = pr, fill = ct)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks.x = element_blank()) +
    xlab(freq$nind %&% " individuals") +
    ylab("Cell type fraction") +
    scale_fill_manual(values = freq$ct.order$hex)

p1 <- 
    .gg.plot(freq$ct.freq.dt, aes(x = projid, y = n, fill = ct)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks.x = element_blank()) +
    xlab(freq$nind %&% " individuals") +
    ylab("#Cells") +
    scale_fill_manual(values = freq$ct.order$hex)

plt <- p1 / p2

print(plt)
.file <- fig.dir %&% "/Fig_In_CT_Freq.pdf"
.gg.save(.file, plot=plt, width = 5, height = 4)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

* `r num.int(round(mean(freq$tot$n)))` cells per individual: `r paste(freq$avg.ct$str, collapse=", ")` on average.

## Differentially Expressed Genes

```{r}
.stat.files <- list.files("result/glob_stat/",
                          pattern=".stat.gz",
                          full.names=TRUE)
```

### Association results

```{r results="asis"}
.names <- names(fread("result/phenotyped.txt.gz",header=TRUE))[-(1:2)]
.idx <- sapply(.stat.files, basename) %>% str_remove_all(".stat.gz")
.names <- .names[as.integer(.idx)]
cat("* [Tab for the associations with " %&% .names %&% "](" %&% .stat.files %&% ")\n\n")
```

```{r read_full_stat}
stat.dt <-
    lapply(.stat.files, fread, header=TRUE) %>%
    do.call(what=rbind) %>%
    mutate(celltype = str_remove(celltype, "subtype-")) %>%
    filter(pheno != "age.death") %>% 
    as.data.table

fit.ashr <- function(.dt) {
    .beta <- .dt$ADE.beta
    .se <- .dt$ADE.se
    .ash <- ashr::ash(.beta, .se, "normal")
    list(.ash$result$lfsr, .ash$result$svalue,
         .ash$result$lfdr, .ash$result$qvalue)
}

stat.dt[tot.1 >= 100 | tot.0 >= 100,
        c("lfsr", "svalue", "lfdr", "qvalue") := fit.ashr(.SD),
        by = .(pheno, celltype)]

stat.dt[,
        fwer.cocoa := p.adjust(pv.cocoa, "holm"),
        by = .(pheno, celltype)]

stat.dt[,
        fdr.cocoa := p.adjust(pv.cocoa, "fdr"),
        by = .(pheno, celltype)]
```

```{r read_combined_stat}
stat.combined.dt <- fread("result/glob_stat_combined/1.stat.gz")
stat.combined.dt[tot.1 >= 100 | tot.0 >= 100,
                 c("lfsr", "svalue", "lfdr", "qvalue") := fit.ashr(.SD),
                 by = .(pheno, celltype)]

stat.combined.dt[,
                 fwer.cocoa := p.adjust(pv.cocoa, "holm"),
                 by = .(pheno, celltype)]

stat.combined.dt[,
                 fdr.cocoa := p.adjust(pv.cocoa, "fdr"),
                 by = .(pheno, celltype)]
```

#### The number of genes associated with PathoAD by non-parametric ranksum test

```{r}
n.fdr01.genes <- stat.dt[pheno == "pathoAD" & fdr.cocoa < .01, .(gene)] %>%
    unique %>% nrow

n.fwer01.genes <- stat.dt[pheno == "pathoAD" & fwer.cocoa < .01, .(gene)] %>%
    unique %>% nrow

n.tot.genes <- stat.dt[, .(gene)] %>% unique %>% nrow

n.fdr01.celltype.genes <- stat.dt[pheno == "pathoAD" & fdr.cocoa < .01,
                            .(n=length(unique(gene))), by=.(celltype)] %>%
    mutate(str = num.int(n) %&% " found in " %&% celltype) %>%
    select(str) %>%
    unlist %>%
    paste(collapse =", ")

n.fwer01.celltype.genes <- stat.dt[pheno == "pathoAD" & fwer.cocoa < .01,
                             .(n=length(unique(gene))), by=.(celltype)] %>%
    mutate(str = num.int(n) %&% " found in " %&% celltype) %>%
    select(str) %>%
    unlist %>% 
    paste(collapse =", ")
```

* DEGs with FDR controlled at 1%: `r num.int(n.fdr01.genes)` genes (`r num.round(100 * n.fdr01.genes / n.tot.genes)`%) consist of `r paste(n.fdr01.celltype.genes, collapse=", ")`

* DEGs with FWER controlled at 1%: `r num.int(n.fwer01.genes)` genes (`r num.round(100 * n.fwer01.genes / n.tot.genes)`%) consist of `r paste(n.fwer01.celltype.genes, collapse=", ")`

* `r num.int(n.tot.genes)` genes in total.


#### Trending the number of discoveries with different p-value cutoffs

```{r}
pv.cutoff <- 10^(-seq(2,6,by=.25))

.fun <- function(.cutoff, .dt) {
    ret <-
        .dt[pv.cocoa < .cutoff, .(n = .N), by = .(celltype, pheno)] %>%
        mutate(cutoff = .cutoff, method = "cocoa")
    ret2 <-
        .dt[pv.tot < .cutoff, .(n = .N), by = .(celltype, pheno)] %>%
        mutate(cutoff = .cutoff, method = "tot")
    rbind(ret, ret2)
}

.fun2 <- function(.cutoff, .dt) {
    ret <-
        .dt[pv.cocoa < .cutoff, .(n = length(unique(gene))), by = .(pheno)] %>%
        mutate(cutoff = .cutoff, method = "cocoa")
    ret2 <-
        .dt[pv.tot < .cutoff, .(n = length(unique(gene))), by = .(pheno)] %>%
        mutate(cutoff = .cutoff, method = "tot")
    rbind(ret, ret2)
}

.dt1 <- lapply(pv.cutoff, .fun, .dt = stat.dt) %>% do.call(what = rbind)
.dt1.tot <- lapply(pv.cutoff, .fun2, .dt = stat.dt) %>% do.call(what = rbind) %>%
    mutate(celltype = "stratified")
.dt2 <- lapply(pv.cutoff, .fun, .dt = stat.combined.dt) %>% do.call(what = rbind)
```

```{r Fig_pvalue_cutoff_celltype, fig.width=3, fig.height=3}
.ct.order <-
    .dt1[, head(.SD, 1), by = .(celltype)] %>%
    arrange(n) %>%
    select(celltype) %>%
    unlist

.temp <- .colors[match(.ct.order, .colors$ct), ]

.dt1 <- .dt1 %>%
    mutate(celltype = factor(celltype, .ct.order, .temp$ct.name))

.fill <- .temp[, "hex"] %>%
    unlist %>%
    as.character

plt <- 
    .gg.plot(.dt1[pheno == "pathoAD"], aes(x=log10(cutoff), y=n, fill=celltype)) +
    ggtitle("pathoAD") +
    facet_grid(method~., scales="free", space="free") +
    theme(legend.position=c(0,1), legend.justification=c(0,1)) +
    geom_bar(stat="identity") +
    ylab("number of genes discovered in each cell type") +
    xlab("p-value cutoff") +
    scale_fill_manual(values = .fill) +
    scale_colour_manual(values = .fill)

print(plt)
.file <- fig.dir %&% "/Fig_pvalue_cutoff_celltype.pdf"
.gg.save(.file, plot=plt, width=3, height=3)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

```{r Fig_pvalue_cutoff_discovery, fig.width=3, fig.height=3}
.df <- rbind(.dt2, .dt1.tot) %>%
    mutate(celltype = factor(celltype, c("stratified", "combined")))

plt <-
    .gg.plot(.df[pheno == "pathoAD"], aes(x = log10(cutoff), y = n, fill = celltype)) +
    ggtitle("pathoAD") +
    facet_grid(method~., scales="free", space="free") +
    scale_x_continuous(labels = function(x) num.sci(10^x)) +
    theme(legend.position=c(0,1), legend.justification=c(0,1)) +
    geom_line() +
    geom_point(pch=21) +
    scale_fill_manual(values=c("gray20", "white")) +
    ylab("number of unqiue genes discovered") +
    xlab("p-value cutoff")

print(plt)
.file <- fig.dir %&% "/Fig_pvalue_cutoff_discovery.pdf"
.gg.save(.file, plot=plt, width=3, height=3)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

### Check: Histogram of p-values

```{r Fig_histogram, fig.width=5, fig.height=7}
.hist.dt <- stat.dt[!is.na(qvalue),
                    .(n = .N),
                    by = .(pv = round(20*pv.cocoa)/20, celltype, pheno)]

plt <- .gg.plot(.hist.dt, aes(x = pv, y = n)) +
    facet_grid(celltype ~ pheno, scale="free_y") +
    geom_bar(stat = "identity") +
    ylab("Frequency") + xlab("p-value (Wilcoxon)")

print(plt)
.file <- fig.dir %&% "/Fig_histogram.pdf"
.gg.save(.file, plot=plt, width=5, height=7)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

### Check: Will CoCoA manipulate variance asymmetrically?

```{r Fig_variance, fig.width=5, fig.height=6}
sig.dt <- stat.dt %>% 
    select(celltype, pheno, starts_with("sig"))

p1 <-
    .gg.plot(sig.dt, aes(x=celltype, y=sig.cocoa.1 - sig.cocoa.0)) +
    ggtitle("CoCoA (delta)") +
    ylab("SD in case - SD in control") +
    facet_grid(pheno~., scales="free") +
    geom_violin(draw_quantiles=.5, size=.2, fill="gray") +
    geom_hline(yintercept=0, lty=2, col=2, size=.2)

p2 <- 
    .gg.plot(sig.dt, aes(x=celltype, y=sig.tot.1 - sig.tot.0)) +
    ggtitle("log(1+ pseudo-bulk)") +
    ylab("SD in case - SD in control") +
    facet_grid(pheno~., scales="free") +
    geom_violin(draw_quantiles=.5, size=.2, fill="gray") +
    geom_hline(yintercept=0, lty=2, col=2, size = .2)

plt <- p1/p2
print(plt)
.file <- fig.dir %&% "/Fig_variance.pdf"
.gg.save(.file, plot=plt, width=5, height=6)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```


### Check: Can confounding factors associate with biological variables?


```{r Fig_correlation, fig.width=7, fig.height=4}

.aes <- aes(x = ADE.beta, y = ADE.beta.cf)

p1 <-
    .gg.plot(stat.dt, .aes) +
    xlab("Average disease effect") +
    ylab("Average\nconfounding effect") +
    scale_y_continuous(limits = c(-.8, .8)) +
    scale_x_continuous(limits = c(-.8, .8)) +
    facet_wrap(~ pheno, scales="free", nrow=1) +    
    geom_hex(aes(fill = log10(..count..), colour = log10(..count..))) +
    scale_fill_distiller(palette="RdPu", direction = 1) +
    scale_colour_distiller(palette="RdPu", direction = 1) +
    geom_hline(yintercept=0, size=.5, lty=2) +
    geom_vline(xintercept=0, size=.5, lty=2)

.aes <- aes(x = ADD.beta, y = ADC.beta)

p2 <-
    .gg.plot(stat.dt, .aes) +
    xlab("Average disease effect on the case") +
    ylab("Average disease effect\n on the control") +
    scale_y_continuous(limits = c(-.8, .8)) +
    scale_x_continuous(limits = c(-.8, .8)) +
    facet_wrap(~ pheno, scales="free", nrow=1) +    
    geom_hex(aes(fill = log10(..count..), colour = log10(..count..))) +
    scale_fill_distiller(palette="RdPu", direction = 1) +
    scale_colour_distiller(palette="RdPu", direction = 1) +
    geom_hline(yintercept=0, size=.5, lty=2) +
    geom_vline(xintercept=0, size=.5, lty=2)

plt <- p1/p2
print(plt)
.file <- fig.dir %&% "/Fig_correlation.pdf"
.gg.save(.file, plot=plt, width=7, height=4)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

### The number of genes identified by significant average disease effect (ADE) (local FSR 5%) and non-parametric rank test (Wilcoxon FDR 5%): `r num.int(nrow(stat.dt[lfsr < .05 & fdr.cocoa < .05, .(gene)] %>% unique))`

```{r Fig_count_5, fig.width=3, fig.height=3}
.count.df <-
    stat.dt[lfsr < .05 & fdr.cocoa < .05,
            .(n = .N), by = .(celltype, pheno)] %>%
    group_by(celltype) %>%
    mutate(freq = n / sum(n)) %>% 
    ungroup() %>% 
    mutate(row = pheno, col = celltype, weight = log10(freq)) %>%
    order.pair(ret.tab=TRUE)

plt <- 
    .gg.plot(.count.df, aes(x=row, y=col, fill=weight, label=num.int(n))) +
    geom_tile(size=.1, colour="gray20") +
    geom_text(size=2, colour="white") +
    xlab("phenotype") + ylab("celltype") +
    scale_fill_distiller("Pr(assoc|celltype)",
                         direction=1, palette = "Blues",
                         labels=function(x) num.round(10^x))

print(plt)
.file <- fig.dir %&% "/Fig_count_5.pdf"
.gg.save(.file, plot=plt, width=3, height=3)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```

### The number of genes identified by significant average disease effect (ADE) (local FSR 1%) and non-parametric rank test (Wilcoxon FDR 1%): `r num.int(nrow(stat.dt[lfsr < .01 & fdr.cocoa < .01, .(gene)] %>% unique))`

```{r Fig_count, fig.width=3, fig.height=3}
.count.df <-
    stat.dt[lfsr < .01 & fdr.cocoa < .01,
            .(n = .N), by = .(celltype, pheno)] %>%
    group_by(celltype) %>%
    mutate(freq = n / sum(n)) %>% 
    ungroup() %>% 
    mutate(row = pheno, col = celltype, weight = log10(freq)) %>%
    order.pair(ret.tab=TRUE)

plt <- 
    .gg.plot(.count.df, aes(x=row, y=col, fill=weight, label=num.int(n))) +
    geom_tile(size=.1, colour="gray20") +
    geom_text(size=2, colour="white") +
    xlab("phenotype") + ylab("celltype") +
    scale_fill_distiller("Pr(assoc|celltype)",
                         direction=1, palette = "Blues",
                         labels=function(x) num.round(10^x))

print(plt)
.file <- fig.dir %&% "/Fig_count.pdf"
.gg.save(.file, plot=plt, width=3, height=3)
```

```{r results="asis"}
cat("[PDF](" %&% .file %&% ")\n\n")
```



### Show top genes for each cell type and phenotype


```{r read_top_data}
.data.files <- list.files("result/glob_stat/",
                          pattern=".stat_data.gz",
                          full.names=TRUE)

data.dt <-
    lapply(.data.files, fread, header=TRUE) %>%
    do.call(what=rbind) %>%
    mutate(celltype = str_remove(celltype, "subtype-")) %>%
    filter(pheno != "age.death") %>% 
    as.data.table
```

```{r}
show.gene.plot <- function(g, .pheno, .lab, .ct = NULL) {

    .stat <- stat.dt[gene == g & pheno == .pheno]

    .dt <- data.dt[gene == g & pheno == .pheno]

    .dir <- 0

    if(!is.null(.ct)) {
        .stat <- .stat[celltype == .ct]
        .dt <- .dt[celltype == .ct]
        .dir <- sign(.stat$ADE.beta[1])
    }

    .dt <- .dt %>% 
        filter(tot >= 5) %>% 
        mutate(y = ln.mu, y.se = ln.mu.sd) %>% 
        mutate(cc = ln.cf, cc.se = ln.cf.sd) %>% 
        mutate(x = factor(label, c(0, 1), .lab))

    if(.dir >= 0) {
        .dt <- .dt %>%
            as_tibble %>% 
            group_by(celltype) %>% 
            arrange(label, ln.mu) %>% 
            mutate(idx = 1:n()) %>% 
            ungroup

    } else {
        .dt <- .dt %>% 
            as_tibble %>% 
            group_by(celltype) %>% 
            arrange(label, desc(ln.mu)) %>% 
            mutate(idx = 1:n()) %>% 
            ungroup
    }

    .dt <- .dt %>% 
        as.data.table

    .scale <- scale_size_continuous(range=c(.1, 2),
                                    labels = function(x) round(x*x))

    .fill <- scale_fill_manual("", values=c("#e0ecf4","#8856a7"))
    .colour <- scale_colour_manual("", values=c("gray20", "#8856a7"))

    .show.violin <- function(.aes) {
        .gg.plot(.dt, .aes) +
            facet_wrap(~celltype, scales="free", nrow = 1) +
            xlab("label") + .scale +
            geom_violin(size=0, scale="width", fill="gray80", trim=FALSE) +
            geom_jitter(aes(fill=x), width=.1, pch=21, stroke=.2) +
            theme(legend.position = "none") + .fill + 
            ggpubr::stat_compare_means(aes(label = ..p.format..),
                                       method = "wilcox",
                                       size = 2,
                                       colour = "magenta")
    }

    .show.linebar <- function(.aes) {
        .gg.plot(.dt, .aes) +
            facet_wrap(~celltype, scales="free", nrow = 1) +
            geom_linerange(size=.1) +
            geom_point(pch=21, stroke=.1) +
            .scale + .fill + .colour + xlab("individuals")
    }

    .aes <- aes(y = tot, x = x, size = sqrt(tot))
    p1 <- .show.violin(.aes) + ylab("pseudo-bulk") + ggtitle(g)


    .aes <- aes(y = y, x = x, size = sqrt(tot))
    p2 <- .show.violin(.aes) + ylab("disease effect") + ggtitle(g)

    .aes <- aes(y = y, x = idx, size = sqrt(tot),
                fill = x, colour = x, ymin = y - y.se, ymax = y + y.se)

    p2.1 <- .show.linebar(.aes) + ylab("disease effect")

    .aes <- aes(y = cc, x = x, size = sqrt(tot))
    p3 <- .show.violin(.aes) + ylab("confounders") + ggtitle(g)

    .aes <- aes(y = cc, x = idx, size = sqrt(tot), colour = x,
                fill = x, ymin = cc - cc.se, ymax = cc + cc.se)
    p3.1 <- .show.linebar(.aes) + ylab("confounders") +
        theme(legend.position = "none")

    p1/p2/p2.1/p3/p3.1
}
```


```{r}
.pheno <- "pathoAD"
.lab <- c("HC","AD")

.genes <- data.dt[pheno == .pheno, .(gene)] %>% unlist

.temp <-
    stat.dt[pheno == .pheno &
            fdr.cocoa < .01 &
            lfsr < .01,
            .(gene, celltype)] %>%
    unique

n.g.ct <- .temp %>% nrow
n.g <- .temp[, .(gene)] %>% unique %>% nrow
n.g.ct.stat <- .temp[, .(n = .N), by = .(celltype)]

.gene.ct <-
    stat.dt[pheno == .pheno &
            (gene %in% .genes) &
            fdr.cocoa < .01 &
            lfsr < .01,
            .(gene, celltype)] %>%
    unique %>%
    mutate(ct = celltype) %>% 
    left_join(.colors) %>% 
    select(-ct) %>% 
    as.data.frame

.save.example <- function(gi) {
    g <- .gene.ct[gi, 1]
    ct <- .gene.ct[gi, 2]
    plt <- show.gene.plot(g, .pheno, .lab, ct)

    .file <- fig.dir %&% "/examples/" %&%
        .pheno %&% "/Fig_" %&%
        .pheno %&% "_" %&% g %&% "_" %&%
        ct %&% ".pdf"
    
    .gg.save(.file, plot = plt, width=2, height=6)
    return(.file)
}

.files <- sapply(1:nrow(.gene.ct), .save.example)
```

### Controlling Wilcoxon FDR 1% and local FSR 1%, we identified `r num.int(n.g.ct)` gene and cell type pairs, `r num.int(n.g)` AD genes.

* Average: `r num.round(mean(n.g.ct.stat$n))` $\pm$ `r num.round(sd(n.g.ct.stat$n))`

#### A genomic view combining all the cell types, highlighting `r nrow(.gene.ct)` gene-cell-type pairs and `r length(unique(.gene.ct$gene))` genes (at most top 20 genes per cell type).

```{r}
.onto <- read.ontology()
gene.info <- .onto$coding
```

```{r Fig_genomic_combined, fig.width=8, fig.height=4}
.dt <- stat.dt[!is.na(ADE.beta)]

.dt <- .dt[order(stat.dt$pv.cocoa), head(.SD, 1), by = .(pheno, gene)] %>%
    filter(pheno == .pheno) %>%
    mutate(hgnc_symbol = gene) %>% 
    left_join(gene.info) %>% 
    select(-hgnc_symbol) %>% 
    mutate(chr = factor(chromosome_name, c(1:22,"X","Y","MT"))) %>% 
    filter(!is.na(chr)) %>% 
    as.data.table

.dt.highlight <- .gene.ct %>% left_join(.dt) %>% na.omit %>% as.data.table
.dt.lab <- .dt.highlight %>% filter(celltype == "Microglia")

.ct.order <- sort(unique(.dt.highlight$celltype))

.temp <- .colors[match(.ct.order, .colors$ct), ]

.dt.highlight <- .dt.highlight %>%
    mutate(celltype = factor(celltype, .ct.order, .temp$ct.name)) %>%
    as.data.table

.fill <- .temp[, "hex"] %>%
    unlist %>%
    as.character

plt <- 
    .gg.plot(anti_join(.dt, .gene.ct), aes(x=transcription_start_site, y=ADE.beta)) +
    ylab("Average Disese Effect") +
    theme(axis.ticks.x = element_blank()) +
    theme(axis.text.x = element_blank()) +
    theme(legend.position = "bottom") +
    facet_grid(.~chr, space="free", scales="free") +    
    geom_point(colour = "gray", stroke=0, size=.5) +
    geom_point(aes(fill = celltype), data = .dt.highlight, pch=21, stroke = .1, size = 1.5) +
    geom_text_repel(aes(label = gene), data = .dt.lab, size=2) +
    scale_fill_manual(values = .fill)

print(plt)

.file <- fig.dir %&% "/Fig_genomic_combined.pdf"
.gg.save(.file, plot=plt, width=8, height=4)
```

```{r echo=FALSE, results="asis"}
cat("\n\n[Combined genomic view](" %&% .file %&% ")\n\n")
```


```{r results="asis"}
as_tibble(.gene.ct) %>%
    mutate(file = .files) %>% 
    left_join(stat.dt[pheno == .pheno], by = c("gene","celltype")) %>%
    mutate(celltype = ct.name) %>% 
    mutate(ADE = num.round(ADE.beta) %&% " \\pm " %&% num.sci(ADE.se)) %>% 
    mutate(`p-value` = num.sci(pv.cocoa)) %>% 
    select(pheno, gene, celltype, ADE, lfdr, lfsr, `p-value`, file) %>%
    mutate(lfdr = num.sci(lfdr),
           lfsr = num.sci(lfsr),
           file) %>%
    mutate(file = "[PDF](" %&% file %&% ")") %>% 
    arrange(gene, celltype) %>% 
    knitr::kable(format = "html", escape = FALSE) %>%
    kableExtra::collapse_rows(columns = c(1, 2), valign = "top") %>%
    kableExtra::kable_styling(c("striped", "condensed"), full_width=FALSE, font_size = 10)
```