## Introduction

We investigate
a dataset of 102 bulk [RNA-seq](https://en.wikipedia.org/wiki/RNA-Seq) samples
extracted from the prefrontal cortex
of 102 deceased individuals
([details](https://github.com/numpde/sta426_deconvo/tree/main/data/20201128-FGCZ)).
Roughly half of the samples are
from patients diagnosed with 
[ALS](https://en.wikipedia.org/wiki/Amyotrophic_lateral_sclerosis),
the remaining ones are referred to as the control group.
We attempt to infer 
cell type proportions
hidden in each bulk RNA-seq sample.
This procedure,
known as cell type deconvolution
([overview](https://bit.ly/sta426_deconvo_sketch)),
should be integral to
inference from the sample to disease condition.
We focus here on reference-based deconvolution:
in essence,
if $X$ is a bulk RNA-seq count vector,
and $a, b, \ldots$
are single-cell RNA-seq count vectors of
representatives of relevant cells types,
one aims to find non-negative coefficients
$\alpha, \beta, \ldots$
such that
$X \approx \alpha a + \beta b + \ldots$.
Unfortunately,
this reconstruction
is confounded by 
batch effects
such as difference in sequencing platforms,
in particular
different statistical properties of 
the technical noise 
of bulk versus single-cell RNA-seq.

As for cell type reference datasets we resort to

- the single-cell RNA-seq dataset
[Darmanis et al](https://github.com/numpde/sta426_deconvo/tree/main/data/20201206-RefDatasets/2015-Darmanis)
(~300 labeled cells used),
and to

- the single-nucleus RNA-seq dataset
[Allen Brain M1](https://github.com/numpde/sta426_deconvo/tree/main/data/20201206-RefDatasets/2019-AllenBrain-M1)
(>70k labeled cells).

The two datasets are compared 
[here](https://github.com/numpde/sta426_deconvo/tree/main/code/sandbox_ra/20201220-Darmanis_vs_Allen).

We apply the following packages/algorithms:

- **Bisque** (R)
prior to estimating 
cell type proportions transforms 
the bulk gene expression gene-wise
to match the distribution from
the reference dataset.
It expects the reference single cells
to come from several individuals.
Links:
[paper](https://www.nature.com/articles/s41467-020-15816-6),
[vignette](https://cran.r-project.org/web/packages/BisqueRNA/vignettes/bisque.html),
[manual](https://cran.r-project.org/web/packages/BisqueRNA/BisqueRNA.pdf),
[github](https://github.com/cozygene/bisque),
[ReferenceBasedDecomposition](https://rdrr.io/cran/BisqueRNA/src/R/reference_based.R)
(uses [limSolve](https://cran.r-project.org/web/packages/limSolve/limSolve.pdf)).

- **MuSiC** (R)
"up-weighs genes with low cross-subject variance" and vice versa,
and
recursively estimates the proportion of clusters of similar cell types
using genes of low within-cluster variance.
Links:
[paper](https://www.nature.com/articles/s41467-018-08023-x),
[tutorial](https://xuranw.github.io/MuSiC/articles/MuSiC.html),
[github](https://github.com/xuranw/MuSiC),
[music_prop](https://rdrr.io/github/xuranw/MuSiC/src/R/utils.R).

- **NNLS**
is a barebones Python implementation 
of cell type proportion reconstruction
using non-negative least squares.
Links: 
[code](https://github.com/numpde/sta426_deconvo/blob/main/code/sandbox_ra/20201206-NNLS_Darmanis).

- **RNA-Sieve** (Python/Mathematica)
incorporates 
the dependence
of the observation noise on cell type composition
of the bulk sample 
into the maximum likelihood estimator,
and can also provide
confidence intervals.
Links:
[preprint](https://doi.org/10.1101/2020.10.01.322867),
[github](https://github.com/songlab-cal/rna-sieve),
[code](https://github.com/numpde/sta426_deconvo/blob/main/code/sandbox_ra/20201226-RNASieve/a_deconvolution.py).

- **DWLS** (R)
identifies differentially expressed genes as markers
and averages those 
for each cell type.
The absolute square residual 
in the usual NLLS
is replaced by a relative one,
requiring an iterative procedure to find a minimizer.
Links:
[paper](https://www.nature.com/articles/s41467-019-10802-z),
[github](https://github.com/dtsoucas/DWLS),
[bitbucket](https://bitbucket.org/yuanlab/dwls/src/master/),
[quickstart](https://github.com/dtsoucas/DWLS/blob/master/Manual.docx),
[manual](https://bitbucket.org/yuanlab/dwls/src/master/manual.pdf).

- We were unable to run 
[SCDC](https://meichendong.github.io/SCDC/articles/SCDC.html).


This script:

- [source code](https://github.com/numpde/sta426_deconvo/tree/main/code/sandbox_ra/20201226-AllDeco),

- [output folder](https://github.com/numpde/sta426_deconvo/tree/main/code/sandbox_ra/20201226-AllDeco/output).


## Preliminaries {.tabset .tabset-fade .tabset-pills}

### Setup    

```{r setup, include=TRUE}
suppressPackageStartupMessages({
  # install.packages("kableExtra")
  requireNamespace("kableExtra")

  # install.packages("dplyr")
  library(dplyr)
  # install.packages("ggplot2")
  library(ggplot2)

  # install.packages("pheatmap")
  requireNamespace("pheatmap")
  # install.packages("stringr")
  requireNamespace("stringr")
  # install.packages("reshape2")
  requireNamespace("reshape2")
  # install.packages("pbapply")
  requireNamespace("pbapply")
  # install.packages("pathlibr")
  requireNamespace("pathlibr")

  # avoid importing `exprs` that leads to clashes
  requireNamespace("rlang")

  # https://vroom.r-lib.org/articles/vroom.html
  # install.packages("vroom")
  requireNamespace("vroom")

  # Provides `ExpressionSet` data structure
  # install.packages("BiocManager")
  # BiocManager::install("Biobase")
  requireNamespace("Biobase")

  #
  # Deconvolution packages follow
  #

  # install.packages("BisqueRNA")
  requireNamespace("BisqueRNA")

  # remotes::install_github("renozao/xbioc@1354168")
  requireNamespace("xbioc") # For MuSiC and such

  # https://xuranw.github.io/MuSiC/articles/MuSiC.html
  # remotes::install_github("xuranw/MuSiC@01e51ba")
  requireNamespace("MuSiC")

  # We need to expose some functions from those packages for DWLS:
  # BiocManager::install("ROCR")
  # BiocManager::install("MAST")
  requireNamespace("ROCR")
  requireNamespace("MAST")

  # If installed afresh, this takes ages and then some:
  # remotes::install_bitbucket("yuanlab/dwls@f13dcf9")
  requireNamespace("DWLS")
})
```


### Session info

```{r, echo=FALSE}
sessionInfo()
```


### Paths

```{r}
BASEPATH <- pathlibr::Path$new(".")
stopifnot("deconvolution.Rmd" %in% names(BASEPATH$dir))
```

```{r}
out_file <- (function(.) BASEPATH$join("output")$join(.)$show)
dir.create(out_file(""), showWarnings = FALSE)
```

```{r}
path_to <- (function(.) Sys.glob(BASEPATH$join("../../..")$join(.)$show))
```

### Random-seed

```{r}
set.seed(43)
```

### Row names, etc.

```{r}
# Use first column as index
by_col1 <- (function(.) tibble::column_to_rownames(., colnames(.)[1]))
# Use index as new column `name`
ind2col <- (function(., name) tibble::rownames_to_column(., var = name))
```

```{r}
group_by_and_sum <-
  function(., column_name) {
    (.) %>%
      dplyr::group_by(!!rlang::sym(column_name)) %>%
      dplyr::summarise_all(sum) %>%
      tibble::column_to_rownames(column_name)
  }
```

```{r}
fix_names <- function(.) {
  (.) %>%
    # https://stackoverflow.com/a/55184433
    rlang::set_names(stringr::str_replace(names(.), "cell type", "celltype"))
}
```

```{r}
col_norm2 <- (function(.) t(t(.) / sqrt(colSums((.)**2))))
```

### File I/O

The function `utils::read.delim` is too slow
to read wide tables (Allen Brain M1 has >70k cells).

```{r}
Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 32)
from_csv <- function(.) {
  vroom::vroom(., del = "\t") %>%
    by_col1() %>%
    suppressMessages()
}
```

A call to `to_csv` should be preceded by
`ind2col("name for index")`
to write the row names.

```{r}
to_csv <- function(., f) {
  vroom::vroom_write(., out_file(f), delim = "\t")
}
```

### Plotting

```{r}
hush <- base::invisible
```

```{r}
ggplot2::theme_set(theme_light(base_size = 15))
```

```{r}
kable <- function(.) {
  kableExtra::kbl(., align = "c") %>%
    kableExtra::kable_paper("hover", full_width = F)
}
```

```{r}
# Wide plots (inches?)
WIDTH1 <- 20
HEIGHT1 <- 3
```

```{r}
hist_theme <-
  theme(
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.title.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_blank(),
  )
```

```{r}
save_barplot <-
  function(., filename) {
    as.data.frame(.) %>%
      ind2col("celltype") %>%
      reshape2::melt(
        id = "celltype",
        var = "sample",
        value.name = "y"
      ) %>%
      {
        ggplot(., aes(x = sample, y = y, fill = celltype)) +
          geom_bar(stat = "identity") +
          scale_fill_brewer(palette = "Paired") +
          ylim(0, 1) +
          hist_theme +
          theme(
            legend.title = element_blank(),

            axis.text.y = element_blank(),

            axis.title.x = element_blank(),
            axis.text.x = element_text(
              a = -90,
              vjust = 0.5,
              hjust = 0,
              size = 10,
            )
          ) +
          ggsave(
            filename = out_file(filename),
            width = WIDTH1,
            height = HEIGHT1,
            device = "png"
          )
      } %>%
      suppressWarnings()
    return(.)
  }
```

```{r}
save_heatmap <-
  function(., filename) {
    as.data.frame(.) %>%
      pheatmap::pheatmap(
        filename = out_file(filename),
        cluster_rows = FALSE,
        cluster_cols = FALSE,
        fontsize_row = 14,
        fontsize_col = 10,
        width = WIDTH1,
        height = HEIGHT1
      )
    return(.)
  }
```





## Data sources {.tabset .tabset-fade .tabset-pills}


### Bulk expression: FGCZ


```{r, cache=TRUE}
fgcz_data <-
  path_to("data/20201128-FGCZ/*count.zip") %>%
  unz(., unzip(., list = TRUE)$Name) %>%
  from_csv() %>%
  # Collapse ENSG IDs by gene_name:
  group_by_and_sum("gene_name")

fgcz_meta <-
  from_csv(path_to("data/20201128-FGCZ/*infos.tsv"))
```

This will be used later.

```{r}
assert_fgcz_meta_order <-
  function(.) {
    stopifnot(all(rownames(.) == rownames(fgcz_meta)))
    return(.)
  }
```

```{r}
t(fgcz_data) %>%
  assert_fgcz_meta_order() %>%
  hush()
```

### Ref scRNA: Darmanis

```{r, cache=TRUE}
darm_data <- from_csv(path_to("data/*/2015-Darmanis/b*/data.csv.gz"))
darm_meta <- from_csv(path_to("data/*/2015-Darmanis/b*/meta.csv.gz")) %>%
  fix_names()
```

```{r}
exclude_celltypes <- c("fetal_quiescent", "fetal_replicating", "hybrid")
```

Plot cell type counts.

```{r}
darm_meta$celltype %>%
  data.frame(x = .) %>%
  ggplot(aes(x = x, fill = if_else(x %in% exclude_celltypes, "drop", "keep"))) +
  geom_bar() +
  ggtitle("Cell types in the 'Darmanis' reference dataset") +
  scale_y_log10() +
  hist_theme +
  theme(axis.title.x = element_blank()) +
  theme(axis.text.x = element_text(a = 45, hjust = 1)) +
  theme(legend.title = element_blank())
```

Remove unnecessary cell types
from the scRNA reference dataset.

```{r, cache=TRUE}
darm_meta <-
  darm_meta %>%
  dplyr::filter(!(celltype %in% exclude_celltypes))

darm_data <-
  darm_data %>%
  dplyr::select_if(names(.) %in% rownames(darm_meta))

stopifnot(285 == nrow(darm_meta))
stopifnot(285 == ncol(darm_data))
```


### Ref snRNA: Allen Brain M1

```{r, cache=TRUE}
abm1_meta <- from_csv(path_to("data/*/2019-AllenBrain-M1/b*/meta.csv*"))
abm1_data <- from_csv(path_to("data/*/2019-AllenBrain-M1/b*/data.csv*"))
```

Some sanity checks.

```{r}
# Number of single cells
stopifnot(ncol(abm1_data) > 70000)
stopifnot(nrow(abm1_meta) > 70000)
# Number of genes and some examples
stopifnot(nrow(abm1_data) == 141)
stopifnot(all(c("PIK3CD", "WNT4", "LDLRAP1") %in% rownames(abm1_data)))
```

Drop samples with too little expression,
otherwise Bisque doesn't go through.

```{r, cache=TRUE}
abm1_data <-
  abm1_data %>%
  # Keep only genes common with the FGCZ dataset
  {
    (.)[rownames(.) %in% rownames(fgcz_data), ]
  } %>%
  # Drop zero columns (dplyr variant is too slow here)
  {
    (.)[, colSums(.) != 0]
  }

abm1_meta <-
  abm1_meta %>%
  dplyr::filter(rownames(.) %in% colnames(abm1_data))
```

Plot cell type counts.

```{r}
data.frame(
  x = abm1_meta$celltype,
  donor = abm1_meta$donor
) %>%
  ggplot(aes(x = x, fill = donor)) +
  geom_bar(position = "dodge") +
  ggtitle("Cell types in the 'Allen Brain M1' reference dataset") +
  scale_fill_brewer(palette = "Set1") +
  scale_y_log10() +
  hist_theme +
  theme(axis.title.x = element_blank()) +
  theme(axis.title.y = element_blank())
```


### Marker genes

These marker genes 
[were extracted](https://github.com/numpde/sta426_deconvo/tree/main/data/20201206-RefDatasets/2019-AllenBrain-M1)
from the Allen Brain M1 dataset.
Technically, these should be the same as 
in `abm1_data` before filtering.
Just in case, we subset them the maximal common subset.


```{r, cache=TRUE}
abm1_markergenes <-
  path_to("data/*/*AllenBrain-M1/b*/marker_genes.csv") %>%
  from_csv() %>%
  rownames() %>%
  {
    (.)[(.) %in% rownames(fgcz_data)]
  } %>%
  {
    (.)[(.) %in% rownames(abm1_data)]
  }

stopifnot(117 == length(abm1_markergenes))
```


## Exploratory {.tabset .tabset-fade .tabset-pills}

### Initial

Some cross-sections of 
the bulk dataset
can be found
[here](https://github.com/numpde/sta426_deconvo/tree/main/code/sandbox_ra/20201130-FirstLook).

### More

Of particular interest for the following 
are 
[RIN](https://en.wikipedia.org/wiki/RNA_integrity_number)
and post-mortem delay.
They exhibit a dependence on the source
but appear to be uncorrelated.

```{r}
fgcz_meta %>%
  ggplot(aes(x = pmDelay, y = RIN, shape = Condition, color = Source, size = LibConc_100_800bp)) +
  scale_x_log10() +
  labs(
    x = "Post-mortem delay (hours)",
    y = "RIN",
    color = "Sample source",
    shape = "Condition",
    size = "Lib. conc. (a.u.)"
  ) +
  geom_point(alpha = 0.8) +
  scale_color_brewer(palette = "Set1")
```

When we order the samples
by post-mortem delay or RIN,
we can see short
runs of samples
where the condition ALS/Control does not change.
Thus if we introduce a quality cut-off
(e.g. low RIN),
we should take care not to bias
the number of ALS/Control samples.
For example,
dropping 14 samples of lowest RIN could be a good choice.

```{r}
fgcz_meta %>%
  {
    lapply(
      c("RIN", "pmDelay"),
      function(field) {
        (.) %>%
          mutate(value = !!rlang::sym(field)) %>%
          arrange(value) %>%
          mutate(rank = (1:nrow(.))) %>%
          mutate(cdf_als = cumsum(Condition == "ALS") - (rank / 2)) %>%
          mutate(ranked_by = field) %>%
          mutate(median_value = median(value))
      }
    )
  } %>%
  bind_rows() %>%
  ggplot(aes(x = rank, y = cdf_als, shape = Condition, color = Source)) +
  labs(x = "Sample number", y = "(Cumulative number of ALS cases) - Expected") +
  facet_wrap(~ranked_by, nrow = 1) +
  geom_point(size = 2, alpha = 0.7) +
  scale_color_brewer(palette = "Set1") +
  scale_x_continuous(breaks = seq(0, 100, by = 10))
```


## Deconvolution


### From Darmanis {.tabset .tabset-fade .tabset-pills}

In this section we use
**Darmanis et al**
as the reference dataset.

Please select one of the methods:

#### Bisque {#fgcz_darm_bisque}

First, repackage the data
as ExpressionSet
following the 
[Bisque vignette](https://cran.r-project.org/web/packages/BisqueRNA/vignettes/bisque.html).
This will also be used in other methods.

```{r}
fgcz_eset <- Biobase::ExpressionSet(
  # Subset to marker genes for residual norm computation
  assayData = as.matrix(fgcz_data[abm1_markergenes, ])
)
```

```{r}
darm_eset <- Biobase::ExpressionSet(
  # Expression data
  assayData = as.matrix(darm_data),
  # Metadata
  phenoData = Biobase::AnnotatedDataFrame(
    data = data.frame(
      row.names = rownames(darm_meta),
      cellType = darm_meta$celltype,
      SubjectName = darm_meta$experiment_sample_name,
      check.names = FALSE,
      check.rows = FALSE,
      stringsAsFactors = FALSE
    ),
    varMetadata = data.frame(
      row.names = c("cellType", "SubjectName"),
      labelDescription = c("cellType", "SubjectName")
    )
  )
)
```

Deconvolution.

```{r, cache=TRUE}
bisque_report <-
  BisqueRNA::ReferenceBasedDecomposition(
    # BULK DATA
    bulk.eset = fgcz_eset,
    # REFERENCE
    sc.eset = darm_eset,
    #
    use.overlap = FALSE,
    markers = abm1_markergenes,
    verbose = FALSE
  )
```

```{r}
# Fields of the Bisque result
bisque_report %>%
  summary() %>%
  kable()
```


Bisque returns proportions 
that sum to one:

```{r}
stopifnot(max(abs(1 - colSums(bisque_report$bulk.props))) <= 1e-10)
```

```{r}
fgcz_darm_bisque <- bisque_report$bulk.props
```

Cluster bulk samples by composition.

```{r}
samples_order <-
  fgcz_darm_bisque %>%
  {
    (.)[, hclust(dist(t(.)))$order]
  } %>%
  colnames()
```

Save to disk
and
visualize inferred cell type composition
by bulk sample.

```{r}
save_heat_bars_data <-
  function(., prefix = deparse(substitute(.))) {
    sorted.df <- (.)[order(tolower(rownames(.))), samples_order]

    filenames <-
      (function(ext) paste(prefix, ext, sep = "_")) %>% {
        list(
          data = (.)("data.csv"),
          heat = (.)("heat.png"),
          bars = (.)("bars.png"),
          path = out_file("")
        )
      }

    # Attach the clinical Condition to sample IDs
    sorted.df.tagged <-
      sorted.df %>%
      t() %>%
      as.data.frame() %>%
      ind2col("Name") %>%
      mutate(suffix = if_else(fgcz_meta[Name, "Condition"] == "ALS", " (ALS)", "")) %>%
      mutate(Name = paste0(Name, suffix)) %>%
      select(-suffix) %>%
      by_col1() %>%
      t()

    # Now plot as heatmap and as barplot
    sorted.df.tagged %>%
      save_heatmap(filenames$heat) %>%
      save_barplot(filenames$bars)

    # Save table to disk
    sorted.df %>%
      as.data.frame() %>%
      ind2col("celltype") %>%
      to_csv(filenames$data)

    return(filenames)
  }
```

```{r}
save_heat_bars_data(fgcz_darm_bisque) %>% hush()
```

![](`r out_file("fgcz_darm_bisque_heat.png")`)
![](`r out_file("fgcz_darm_bisque_bars.png")`)


The cell type proportions [returned by Bisque](#fgcz_darm_bisque)
vary rather more than expected.
The samples come from four
different hospitals 
(the field `fgcz_meta$Source`),
which could introduce the strongest batch effect.
We can also compare that 
to the [RIN](https://en.wikipedia.org/wiki/RNA_integrity_number)
(RNA integrity number).
Age could be another important factor.

```{r}
PC1_and_friends <-
  function(.) {
    plots <-
      prcomp(t(.))$x %>%
      as.data.frame() %>%
      assert_fgcz_meta_order() %>%
      cbind(fgcz_meta) %>%
      cbind(y = (2 * rank((.)$PC1) / length((.)$PC1) - 1)) %>%
      {
        # Make several plots silently
        lapply(
          c("PC2", "RIN", "Age"),
          function(field) {
            cbind((.), x = (.)[[field]]) %>% {
              ggplot(., aes(x = x, y = y, shape = Condition, color = Source)) +
                labs(
                  x = field,
                  y = "Deviation from the average composition (~PC1)",
                  color = "Sample source",
                  shape = "Condition"
                ) +
                geom_point(size = 5, alpha = 0.2 + abs((.)$y)) +
                scale_color_brewer(palette = "Set1")
            }
          }
        )
      }

    t <- theme(
      legend.position = "none",
      axis.title.y = element_blank(),
      axis.ticks.y = element_blank(),
      axis.text.y = element_blank(),
    )

    filename <- out_file(paste(
      deparse(substitute(.)),
      paste(lapply(plots, function(plot) plot$labels$x), collapse = "_"),
      "vs__PC1.png",
      sep = "__"
    ))

    plots <- gridExtra::arrangeGrob(
      grobs = c(
        # Remove redundant stuff
        lapply(plots, function(.) (. + t)),
        # but keep one legend
        list(cowplot::get_legend(plots[[1]]))
      ),
      # All grobs in one row
      nrow = 1,
      # Y-axis label
      left = grid::grid.text(
        plots[[1]]$labels$y,
        rot = 90,
        draw = FALSE,
        gp = grid::gpar(fontsize = ggplot2::theme_get()$text$size)
      )
    )

    ggsave(
      filename = filename,
      plot = plots,
      width = WIDTH1,
      height = WIDTH1 / (length(plots) - 1)
    )

    return(filename)
  }
```

Principal components of cell type proportions, etc.

```{r}
PC1_and_friends(fgcz_darm_bisque) %>% hush()
```

![](`r out_file("fgcz_darm_bisque__PC2_RIN_Age__vs__PC1.png")`)


#### MuSiC

```{r, cache=TRUE}
# Prevent errors from MuSiC
exprs <- xbioc::exprs
pVar <- xbioc::pVar

music_report.darm <-
  MuSiC::music_prop(
    bulk.eset = fgcz_eset,
    sc.eset = darm_eset,
    clusters = darm_meta$celltype,
    samples = names(darm_data)
  ) %>%
  suppressMessages()
```

```{r}
music_report.darm %>%
  summary() %>%
  kable()
```

Estimate cell type proportions (celltype x sample).

```{r}
fgcz_darm_music <-
  t(music_report.darm$Est.prop.weighted) %>% {
    (.)[rownames(fgcz_darm_bisque), ]
  }
```

```{r}
save_heat_bars_data(fgcz_darm_music) %>% hush()
PC1_and_friends(fgcz_darm_music) %>% hush()
```

![](`r out_file("fgcz_darm_music_heat.png")`)
![](`r out_file("fgcz_darm_music_bars.png")`)

![](`r out_file("fgcz_darm_music__PC2_RIN_Age__vs__PC1.png")`)


We would like to know
whether there is 
the cell type composition differs 
in ALS vs Control samples.
First,
the following figure shows the distribution
of the first principal component
(of cell type proportions)
separated by sample source.

```{r}
fgcz_darm_music.pca <-
  fgcz_darm_music %>%
  {
    as.data.frame(prcomp(t(.))$x)
  } %>%
  assert_fgcz_meta_order() %>%
  cbind(fgcz_meta)
```

```{r, eval=TRUE}
fgcz_darm_music.pca %>% {
  ggplot(., aes(x = PC1, fill = Source, linetype = Source)) +
    labs(fill = "Sample source") +
    guides(linetype = FALSE) +
    geom_density(
      position = "identity", alpha = 0.4
    ) +
    geom_histogram(
      aes(y = stat(density) / 5),
      position = "dodge", bins = 20, color = 1, linetype = 1
    ) +
    hist_theme +
    theme(axis.text.y = element_blank()) +
    theme(
      panel.grid.major.y = element_blank(),
      panel.grid.minor.y = element_blank(),
    ) +
    scale_fill_brewer(palette = "Set1")
}
```

The separation by **Source**
seen in the above figure 
suggests looking at 
the **Condition** (ALS vs Control)
for each **Source** individually.

```{r}
fgcz_darm_music.pca %>%
  {
    ggplot(., aes(x = PC1, fill = Condition)) +
      geom_density(
        position = "identity", alpha = 0.4
      ) +
      geom_histogram(
        aes(y = stat(density) / 5),
        position = "dodge", bins = 20, color = 1, linetype = 1
      ) +
      hist_theme +
      theme(
        axis.text.x = element_blank(),
        axis.text.y = element_blank()
      ) +
      scale_fill_brewer(palette = "Set1") +
      facet_grid(cols = vars(Source))
  } %>%
  ggsave(
    filename = out_file("fgcz_darm_music_PC1.png"),
    width = WIDTH1,
    height = WIDTH1 / 4
  ) %>%
  hush()
```

![](`r out_file("fgcz_darm_music_PC1.png")`)


Among those,
samples from **Lo**
show a conspicuous further separation 
by **Condition**,
supported by the following 
two-sided two-sample t-test.

```{r}
fgcz_darm_music.pca %>%
  filter(Source != "PD") %>%
  group_by(Source) %>%
  group_map(~ {
    t.test(
      x = filter((.x), Condition == "ALS")$PC1,
      y = filter((.x), Condition != "ALS")$PC1,
      alternative = "two"
    ) %>% {
      list(
        `Sample source` = dplyr::pull((.y)),
        `p-value` = round((.)$p.value, 3),
        `t-statistic` = round((.)$statistic, 2),
        `#dof` = round((.)$parameter, 1)
      )
    }
  }) %>%
  data.table::rbindlist() %>%
  kable()
```


Ordering 
the samples
by the first principal component 
reveals its meaning
as a measure of 
the content of astrocytes vs neurons.
The samples from **Lo** suggest
that high proportion of neurons
is associated with ALS.

```{r}
fgcz_darm_music %>%
  t() %>%
  as.data.frame() %>%
  assert_fgcz_meta_order() %>%
  cbind(Source = fgcz_meta$Source) %>%
  cbind(PC1 = fgcz_darm_music.pca$PC1) %>%
  tibble::rownames_to_column("ID") %>%
  mutate(ID = paste(ID, if_else(fgcz_meta$Condition == "ALS", "(ALS)", ""))) %>%
  arrange(PC1) %>%
  select(-PC1) %>%
  group_by(Source) %>%
  group_map(~ {
    (.x) %>%
      tibble::column_to_rownames("ID") %>%
      t() %>%
      save_barplot(paste("fgcz_darm_music_", pull(.y), "_bars.png", sep = ""))
  }) %>%
  hush()
```

Lo:
![](`r out_file("fgcz_darm_music_Lo_bars.png")`)

NB:
![](`r out_file("fgcz_darm_music_NB_bars.png")`)

Ox:
![](`r out_file("fgcz_darm_music_Ox_bars.png")`)

PD:
![](`r out_file("fgcz_darm_music_PD_bars.png")`)


Let's have a look again at the quality metrics,
this time annotated by 
the estimated fraction of astrocytes.

```{r}
fgcz_darm_music %>%
  t() %>%
  as.data.frame() %>%
  assert_fgcz_meta_order() %>%
  cbind(fgcz_meta) %>%
  ggplot(aes(x = pmDelay, y = RIN, shape = Condition, color = Source, alpha = 100 * astrocytes)) +
  scale_x_log10() +
  labs(
    x = "Post-mortem delay (hours)",
    y = "RIN",
    color = "Sample source",
    shape = "Condition",
    alpha = "Astrocytes, %"
  ) +
  geom_point(size = 5) +
  scale_color_brewer(palette = "Set1")
```

A trend
towards a higher content of (intact) astrocytes
at lower RIN
is now visible:

```{r}
fgcz_darm_music %>%
  t() %>%
  as.data.frame() %>%
  assert_fgcz_meta_order() %>%
  cbind(fgcz_meta) %>%
  ggplot(aes(x = RIN, y = 100 * astrocytes, shape = Condition, color = Source)) +
  scale_x_log10() +
  labs(
    y = "Estimated astrocyte content, %",
    color = "Sample source",
    shape = "Condition"
  ) +
  geom_point(size = 5, alpha = 0.8) +
  scale_color_brewer(palette = "Set1")
```

Notably,
among **Lo** samples,
samples with lower RIN are predominantly 
from the control group.
Let's have a second look at
the p-values, now
as a function of RIN cut-off
(we drop $n$ samples with lowest RIN from each Source).

```{r}
fgcz_darm_music.pca %>%
  filter(Source != "PD") %>%
  group_by(Source) %>%
  group_map(~ {
    lapply(
      0:10,
      function(ndrop) {
        (.x) %>%
          arrange(RIN) %>%
          {
            (.)[(1 + ndrop):nrow(.), ]
          } %>%
          {
            t.test(
              x = filter((.), Condition == "ALS")$PC1,
              y = filter((.), Condition != "ALS")$PC1,
              alternative = "two"
            ) %>% {
              list(
                `Sample source` = dplyr::pull((.y)),
                `p-value` = round((.)$p.value, 3),
                `t-statistic` = round((.)$statistic, 2),
                `#dof` = round((.)$parameter, 1),
                `Drop` = ndrop
              )
            }
          }
      }
    ) %>%
      data.table::rbindlist()
  }) %>%
  data.table::rbindlist() %>%
  ggplot(aes(x = Drop, y = `p-value`, color = `Sample source`)) +
  labs(x = "Number of lowest-RIN samples to drop (from each Source)") +
  geom_point(size = 3) +
  geom_line(linetype = "dashed") +
  scale_x_continuous(breaks = 0:10) +
  scale_y_log10(
    breaks = c(0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1),
    labels = c("1%", "2%", "5%", "10%", "20%", "50%", "100%")
  )
```



#### Baseline NNLS

For the code and additional analysis see 
[here](https://github.com/numpde/sta426_deconvo/blob/main/code/sandbox_ra/20201206-NNLS_Darmanis).

```{r}
fgcz_darm_nnls <-
  path_to("code/*/*-NNLS_Darmanis/a_*/celltypes.csv") %>%
  # The first column is not unique
  # Keep header as is, e.g. "H-0C083"
  utils::read.delim(., check.names = F) %>%
  # `cell type` ~> `celltype`
  fix_names() %>%
  # Collapse cell types:
  group_by_and_sum("celltype")
```

```{r}
save_heat_bars_data(fgcz_darm_nnls) %>% hush()
PC1_and_friends(fgcz_darm_nnls) %>% hush()
```

![](`r out_file("fgcz_darm_nnls_heat.png")`)
![](`r out_file("fgcz_darm_nnls_bars.png")`)

![](`r out_file("fgcz_darm_nnls__PC2_RIN_Age__vs__PC1.png")`)



These results are similar enough to MuSiC's
to warrant a closer look.
For each sample we compute the cosine similarity
between the two methods
and plot a histogram of those
separated by sample source.

```{r}
similarity_plot <- function(X1, X2) {
  X1 <- X1[order(rownames(X1)), order(colnames(X1))]
  X2 <- X2[order(rownames(X2)), order(colnames(X2))]

  # Check that samples are arranged consistently
  stopifnot(all(colnames(X1) == colnames(X2)))
  stopifnot(all(rownames(X1) == rownames(X2)))

  colSums(col_norm2(X1) * col_norm2(X2)) %>%
    data.frame(
      x = .,
      s = fgcz_meta$Source
    ) %>%
    {
      ggplot(., aes(x = x, fill = s, linetype = s)) +
        labs(fill = "Sample source") +
        guides(linetype = FALSE) +
        geom_density(
          position = "identity", alpha = 0.4
        ) +
        geom_histogram(
          aes(y = stat(density) / 5),
          position = "dodge", bins = 30, color = 1, linetype = 1,
        ) +
        hist_theme +
        theme(axis.text.y = element_blank()) +
        theme(
          panel.grid.major.y = element_blank(),
          panel.grid.minor.y = element_blank(),
        ) +
        scale_fill_brewer(palette = "Set1")
    }
}
```

```{r}
similarity_plot(fgcz_darm_music, fgcz_darm_nnls) +
  labs(x = "Cosine similarity of a sample (MuSiC vs NNLS)")
```


#### RNA-Sieve

For the code see 
[here](https://github.com/numpde/sta426_deconvo/blob/main/code/sandbox_ra/20201226-RNASieve/a_deconvolution.py).


```{r}
fgcz_darm_sieve <-
  path_to("code/*/*-RNASieve/a_*/celltypes.csv") %>%
  from_csv() %>%
  fix_names()
```

```{r}
save_heat_bars_data(fgcz_darm_sieve) %>% hush()
PC1_and_friends(fgcz_darm_sieve) %>% hush()
```

![](`r out_file("fgcz_darm_sieve_heat.png")`)
![](`r out_file("fgcz_darm_sieve_bars.png")`)

![](`r out_file("fgcz_darm_sieve__PC2_RIN_Age__vs__PC1.png")`)


```{r}
similarity_plot(fgcz_darm_music, fgcz_darm_sieve) +
  labs(x = "Cosine similarity of a sample (MuSiC vs RNA-Sieve)")
```



#### DWLS

Preprocesses the reference dataset.

```{r, cache=TRUE, results='hide', message=FALSE}
# This takes some 20min

prediction <- ROCR::prediction
performance <- ROCR::performance
zlm <- MAST::zlm

dwls_darm_signature <-
  DWLS::buildSignatureMatrixMAST(
    scdata = darm_data,
    id = setNames(as.list(darm_meta[["celltype"]]), rownames(darm_meta)),
    path = out_file("dwls_darm_signature/")
  )
```

Apply DWLS deconvolution
to each bulk sample.

```{r, cache=TRUE}
fgcz_darm_dwls <-
  fgcz_data %>%
  {
    genes <- intersect(rownames(dwls_darm_signature), rownames(.))

    pbapply::pblapply(
      (.)[genes, ],
      function(sample) {
        # Silence silly printout
        purrr::quietly(DWLS::solveDampenedWLS)(
          as.matrix(dwls_darm_signature[genes, ]), as.matrix(sample)
        )$result
      },
      # Parallel processes -- requires (cl x RAM)
      cl = 1
    )
  } %>%
  as.data.frame(check.names = F)
```

```{r}
fgcz_darm_dwls <-
  fgcz_darm_dwls %>% {
    fgcz_darm_dwls[order(tolower(rownames(.))), ]
  }
```

```{r}
fgcz_darm_dwls %>%
  save_heat_bars_data("fgcz_darm_dwls") %>%
  hush()
```

![](`r out_file("fgcz_darm_dwls_heat.png")`)
![](`r out_file("fgcz_darm_dwls_bars.png")`)


```{r}
similarity_plot(fgcz_darm_music, fgcz_darm_dwls) +
  labs(x = "Cosine similarity of a sample (MuSiC vs DWLS)")
```



### From Allen Brain M1 {.tabset .tabset-fade .tabset-pills}

In this section we use
**Allen Brain M1**
as the reference dataset.

#### Bisque

Repackage:

```{r, eval=FALSE, echo=FALSE}
# # HACK TO SPEED UP
# a <- unique(abm1_meta$donor)[1]
# b <- unique(abm1_meta$donor)[2]
# a <- rownames(abm1_meta[abm1_meta$donor == a, ])[1:10]
# b <- rownames(abm1_meta[abm1_meta$donor == b, ])[1:10]
# s <- c(a, b)
# abm1_meta <- abm1_meta[s, ]
# abm1_data <- abm1_data[, s]
```

```{r}
abm1_eset <- Biobase::ExpressionSet(
  # Expression data
  assayData = as.matrix(abm1_data),
  # Metadata
  phenoData = Biobase::AnnotatedDataFrame(
    data = data.frame(
      row.names = rownames(abm1_meta),
      cellType = abm1_meta$celltype,
      SubjectName = abm1_meta$donor,
      check.names = FALSE,
      check.rows = FALSE
    ),
    varMetadata = data.frame(
      row.names = c("cellType", "SubjectName"),
      labelDescription = c("cellType", "SubjectName")
    )
  )
)
```

Deconvolution:

```{r, cache=TRUE}
bisque_report <-
  BisqueRNA::ReferenceBasedDecomposition(
    # BULK DATA
    bulk.eset = fgcz_eset,
    # REFERENCE
    sc.eset = abm1_eset,
    #
    use.overlap = FALSE,
    markers = abm1_markergenes,
    verbose = FALSE
  )
```

```{r}
bisque_report$bulk.props %>%
  save_heat_bars_data("fgcz_abm1_bisque") %>%
  hush()
```

![](`r out_file("fgcz_abm1_bisque_heat.png")`)
![](`r out_file("fgcz_abm1_bisque_bars.png")`)


#### MuSiC

```{r, cache=TRUE}
# Read from disk if already available
# otherwise it takes ages
if (!file.exists(out_file("fgcz_abm1_music_data.csv"))) {
  # MuSiC with ABM1
  music_report.abm1 <-
    MuSiC::music_prop(
      bulk.eset = fgcz_eset,
      sc.eset = abm1_eset,
      clusters = abm1_meta$celltype,
      samples = names(abm1_data)
    ) %>%
    suppressMessages()

  fgcz_abm1_music <-
    t(music_report.abm1$Est.prop.weighted)
} else {
  fgcz_abm1_music <-
    from_csv(out_file("fgcz_abm1_music_data.csv"))
}

fgcz_abm1_music %>%
  save_heat_bars_data("fgcz_abm1_music") %>%
  hush()
```

![](`r out_file("fgcz_abm1_music_heat.png")`)
![](`r out_file("fgcz_abm1_music_bars.png")`)



## Discussion

While many deconvolution tools have been developed, 
it appears a priori
entirely unclear 
which one/s will give 
the most reliable results. 
We have restricted ourselves to a few tools
that estimate cell type proportions
in the bulk sample
from a reference single-cell dataset.
There are several sources of uncertainty 
that should be taken into account,
including
technical and biological variability of the bulk sample,
technical noise of the reference dataset,
difference in sequencing platforms,
and composition of the reference dataset.
Yet,
perhaps the weakest point of the present mini-study
is insufficient *coverage*:
the reference dataset does not necessarily
cover the aberrant cell states that may
appear in ALS,
and not all tools acknowledge
the possibility.
We have focused on the results from MuSiC 
(with **Darmanis** as reference dataset)
as they seemed more reasonable
and more consistent with two other methods.
However, 
MuSiC inferred a higher proportion of neurons
in ALS patients vs Control
for a subset of samples,
which is unexpected
for a neurodegenerative disease
(unless the brain attempts to compensate for lack in functionality). 


## Acknowledgements

This is a project under the aegis 
of the course
[STA426](https://www.math.uzh.ch/index.php?id=1402&key1=0&key2=3968&key3=352&semId=41)/[CBB](https://cbb.ethz.ch/).
Many thanks to:

- [K Dobler](https://github.com/KalvinDobler)
for contributions to 
the [project sketch](https://bit.ly/sta426_deconvo_sketch),

- [PL Germain](https://schrattlab.ethz.ch/the-group/people/person-detail.MjUwODk0.TGlzdC8yOTEzLDExMTE5MzUxNzY=.html) (and mystery colleague)
for useful pointers and comments,

- [B Jew](https://brandonjew.com/)
for clarifications on Bisque,

- [S Lewandowski](https://staff.ki.se/people/seblew)
for an in-depth discussion and suggestions,

- [H Rehrauer](https://fgcz.ch/the-center/people/rehrauer.html) 
for the dataset and helpful comments.


## Cite as

R Andreev, 
Cell type deconvolution 
from bulk RNA-seq of the human prefrontal cortex,
`r format(Sys.time(), '%Y-%m-%d')`,
[http://bit.ly/deco_ra](http://bit.ly/deco_ra).


***


```{r include=FALSE}
knitr::knit_exit()
```



#### SCDC -- couldn't get it to work

```{r}
# https://github.com/meichendong/SCDC
# install.packages("devtools")
# devtools::install_github("meichendong/SCDC")
requireNamespace("SCDC")

# SCDC fails with:
#
# Error in y[y < q15] <- q15[y < q15] :
#   NAs are not allowed in subscripted assignments
# In addition: There were 50 or more warnings:
# ...
# 50: In FUN(newX[, i], ...) : no non-missing arguments to max; returning -Inf

# https://meichendong.github.io/SCDC/articles/SCDC.html
SCDC::SCDC_ENSEMBLE(
  bulk.eset = fgcz_eset,
  sc.eset.list = list(darmanis = darm_eset),
  sample = "SubjectName",
  ct.varname = "cellType",
  ct.sub = c("astrocytes", "neurons")
)
```



[//]: # XX 

```{r, eval=FALSE, echo=FALSE}
library(SingleCellExperiment)
library(scater)
factor_name <- "Source"
fgcz_data %>%
  edgeR::cpm(.) %>%
  {
    (.)[abm1_markergenes, ]
  } %>%
  {
    list(logcounts = log(1 + .))
  } %>%
  SingleCellExperiment(assays = .) %>%
  scater::runMDS() %>%
  SingleCellExperiment::reducedDim() %>%
  as.data.frame() %>%
  setNames(c("x", "y")) %>%
  ggplot(aes(x = x, y = y, color = fgcz_meta[[factor_name]])) +
  labs(color = factor_name) +
  geom_point(size = 5)
```

```{r, eval=FALSE, echo=FALSE}
fgcz_abm1 %>%
  t() %>%
  as.data.frame() %>%
  dplyr::mutate(condition = fgcz_meta$Condition) %>%
  {
    MASS::lda(formula = condition ~ ., data = (.))
  }
```