# Differential expression and abundance between conditions

```{r multiSplComp_setup, include=FALSE, echo=FALSE}
# First, set some variables:
require(knitr)

opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
opts_chunk$set(echo = TRUE)
opts_chunk$set(eval = TRUE) 
options(stringsAsFactors = FALSE)
opts_chunk$set(fig.width=7, fig.height=7)
set.seed(123) # for reproducibility
```

```{r}
splSetToGet <- "PBMMC,ETV6-RUNX1"
splSetVec <- unlist(strsplit(splSetToGet, ","))
splSetToGet2 <- gsub(",", "_", splSetToGet)
nbPcToComp <- 50
figSize <- 7
```

```{r, message=FALSE, warning=FALSE}
library(scater)
library(scran)
library(batchelor)
library(edgeR)
library(tidyverse)
library(patchwork)
library(DT)
fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
```

Source: [Multi-sample comparisons](https://osca.bioconductor.org/multi-sample-comparisons.html) of the OSCA book.

## Exercise 1 - differential expression

Identify label-specific DE genes that are significant in 'c10' yet not DE in any other label.

Plot the top-ranked gene for inspection.

## Setting up the data

Load the SCE object (with 1200 cells per sample):

```{r}
# Read object in:
merged <- readRDS("../Robjects/caron_sce_nz_postDeconv_1p2kcps_dsi_PBMMC_ETV6-RUNX1_merged.Rds")
# also get raw counts that were written to a separate file
# (to help file sharing)
merged_counts <- readRDS("../Robjects/caron_sce_nz_postDeconv_1p2kcps_dsi_PBMMC_ETV6-RUNX1_merged_counts.Rds")
# put raw counts back:
counts(merged) <- merged_counts
# tidy:
rm(merged_counts)
```

A brief inspection of the results shows clusters contain varying contributions from samples:

```{r}
colLabels(merged) <- merged$clusters.mnn
tab <- table(colLabels(merged), merged$SampleName)
tab
```

On the t-SNE plots below, cells are coloured by type or sample ('batch of origin'). Cluster numbers are superimposed based on the median coordinate of cells assigned to that cluster. 

```{r, fig.height=0.5*figSize, fig.width=1.2*figSize}
p1 <- plotTSNE(merged, colour_by="SampleGroup", text_by="label", point_size=0.3)
p2 <- plotTSNE(merged, colour_by="SampleName", point_size=0.3) +
  facet_wrap(~colData(merged)$SampleName)
p1
p2
```

## Creating pseudo-bulk samples

Sum counts together for all cells with the same combination of label and sample,
with `aggregateAcrossCells`.

```{r}
# Using 'label' and 'sample' as our two factors; each column of the output
# corresponds to one unique combination of these two factors.
columnsToUse <- c("batch", "SampleName", "SampleGroup", "clusters.mnn")
colData(merged) <- colData(merged) %>% data.frame() %>% select(all_of(columnsToUse)) %>% DataFrame
summed <- aggregateAcrossCells(merged, 
    				id = DataFrame(
    					label = merged$clusters.mnn,
    					sample = merged$SampleName
					)
)
colData(summed) %>% head(3)
```

## Performing the DE analysis

Filter out all sample-label combinations with insufficient cells, 20 cells or less.

```{r}
summed.filt <- summed[,summed$ncells >= 20]
```

Construct a common design matrix that will be used in the analysis for each label.

```{r}
# Pulling out a sample-level 'targets' data.frame:
targets <- colData(merged)[!duplicated(merged$SampleName),] %>%
  data.frame() %>%
  select(-clusters.mnn)

# Constructing the design matrix:
design <- model.matrix(~factor(SampleGroup), data=targets)
rownames(design) <- targets$SampleName
```

Apply the `pseudoBulkDGE` function to obtain a list of DE genes for each label.

```{r}
summed.filt$SampleGroup <- factor(summed.filt$SampleGroup)

de.results <- pseudoBulkDGE(summed.filt, 
    label = summed.filt$label,
    design = ~SampleGroup,
    coef = "SampleGroupPBMMC",
    condition = summed.filt$SampleGroup
)
```

Examine the numbers of DEGs at a FDR of 5% for each label using the `decideTestsPerLabel` function.

```{r}
is.de <- decideTestsPerLabel(de.results, threshold=0.05)
summarizeTestsPerLabel(is.de)
```

For each gene, we compute the percentage of cell types in which that gene is upregulated or downregulated. (Here, we consider a gene to be non-DE if it is not retained after filtering.).

<!-- TODO: add gene names -->

```{r}
# Upregulated across most cell types.
up.de <- is.de > 0 & !is.na(is.de)
head(sort(rowMeans(up.de), decreasing=TRUE), 10)
```

```{r}
# Downregulated across cell types.
down.de <- is.de < 0 & !is.na(is.de)
head(sort(rowMeans(down.de), decreasing=TRUE), 10)
```

Identify label-specific DE genes that are significant in 'c10' yet not DE in any other label (FDR threshold of 50%).

Plot the top-ranked gene for inspection.

```{r}
remotely.de <- decideTestsPerLabel(de.results, threshold=0.5)
not.de <- remotely.de==0 | is.na(remotely.de)
```

```{r}
# get c10's 'unique.degs':






```

```{r}
# Choosing the top-ranked gene for inspection:





# plot expression of top gene
# use plotExpression()

```

