---
title: "BSEQ-sc: Deconvolution of Bulk Sequencing Experiments using Single Cell Data"
shorttitle: "BSEQ-sc: Deconvolution Using scRNA-seq"
author: 
  - {name: Renaud Gaujoux, affiliation: "Systems Immunology Lab, Technion, Haifa, Israel"}
  - {name: Maayan Baron, affiliation: "Institute for Computational Medicine, NYU School of Medicine, New York, USA"}
output: 
  BiocStyle::html_document:
    self_contained: false
    fig_caption: true
  pkgdown:
    css: style.css
date: ""
package: bseqsc
bibliography: ~/Documents/articles/library.bib
abstract: >
  This vignette provides a walk through tutorial on how to use Bseq-sc to 
  estimate cell type proportions from bulk sequencing data based 
  on prior single cell data, and integrate them into a cell type-specific
  differential analysis.
  It serves as a companion document to the analysis performed in the original publication (Baron et al. 2016 -- Cell Systems).
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{BSEQ-sc: Deconvolution of Bulk Sequencing Experiments using Single Cell Data}
  %\VignetteEncoding{UTF-8} 
  \usepackage[utf8]{inputenc}
---

```{r setup, include = FALSE}
# knitr options
library(knitr)
knit_hooks$set(backspace = pkgmaker::hook_backspace())
opts_chunk$set(dpi=300, dev="png")#, dev.args=list(pointsize=15))

IS_PKGDOWN <- nzchar(Sys.getenv('_R_PKGDOWN_RENDER'))
# packages
library(BiocStyle)
# pkgmaker::load_project('pancreas/bseq-sc')
library(bseqsc)
library(ggplot2)
library(reshape2)
library(plyr)

# load color theme
COLOR_THEME <- bseqsc:::.ggtheme()

BSEQSC <- 'BSEQ-sc'
BSEQSCpkg <- 'bseqsc'
```

```{r local_data, include = FALSE}
if( local_data <- dir.exists(data_dir <- file.path(pkgmaker::find_devpackage('pancreas'), 'docs/data')) ){
  eset <- readRDS(file.path(data_dir, 'GSE50244.rds'))
  eislet <- readRDS(file.path(data_dir, 'islet-eset.rds'))
}
```

```{r cite, eval = !IS_PKGDOWN, results = 'asis', echo = FALSE}
cat("
## How to cite `r BSEQSC` {-}
When using this work in any way, please cite the following work:

> *A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure*<br />
> <small>M. Baron, A. Veres, S.L. Wolock, A. L. Faust, R. Gaujoux, A. Vetere, J. Hyoje Ryu, B. K. Wagner, S. Shen-Orr, A. M. Klein, D. A. Melton, I. Yanai<br /></small>
> Cell Systems. 2016 Oct 26 [10.1016/j.cels.2016.08.011](https://www.ncbi.nlm.nih.gov/pubmed/27667365)
") 
```

# Requirements
## Software

In order to run the entire deconvolution pipeline, users need to install the `r Githubpkg('shenorrlab/bseq-sc', 'bseqsc')` R package from GitHub.
See the [installtion instructions](pages/installation.html) for details.

## Data
`r BSEQSC` uses two types of input data:

  * bulk data obtained from RNA sequencing of mixed tissue samples: these are the data we want to deconvolve; 
  * data from single cell RNA sequencing: these serve to compute cell type-specific reference profiles that are used 
  to estimate cell type proportions in the bulk data. 

*****
# `r BSEQSC` analysis pipeline

Single cell RNA-seq information is used to deconvolve bulk RNA-seq samples by estimating the cell-type proportion of key cell-types. Statistical deconvolution is then used to leverage the variation in cell type frequencies between individuals to estimate average cell-type specific expression in 2 groups of individuals and to compute cell type-specific differential expression.

We used this approach to investigate gene expression differences between diabetics and healthy individuals, using single cell RNA-seq data obtained from an independent set of pancreatic islets (Figure 1).
However, it is applicable to other type of tissue, for which both bulk and single cell gene expression data are available.

```{r pipeline, echo = FALSE, fig.cap = "**Figure 1:** The BSEQ-sc pipeline applied to pancreatic islets.", results = 'asis'}
cat("![](images/pipeline-800x600.png)")
```

*****
# Sample analysis

We reproduce here the analysis detailed in [@Baron2016].

  * Total gene expression data were obtained from GEO dataset [GSE50244](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50244), which contains data on gobal transcriptomic analysis of bulk human pancreatic islets to study glucose metabolism in healthy and hyper-glycemic conditions [@Fadista2014].
  * Single cell RNA-seq data were generated by us from pancreatic islets from an independent set of 3 healthy indiviuals.

## Data preparation

Bioconductor base package provides the `ExpressionSet` class, which is a convenient data structure to hold expression data along with sample/feature annotation.
Here we use two `ExpressionSet` objects to handle the bulk and single cell data respectively. 

### Bulk data

The dataset's [GEO entry (GSE50244)](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50244) contains raw RNA-seq and sample annotation data. 
For the purpose of this vignette, we will use data that we pre-processed and
made available on the [data
download page](pages/data.html).

After pre-processing, the data were filtered to remove unmapped tags, or tags that mapped to pseudogenes, or known ribosomal/mithocondrial genes (list available in file `data/genes_mitrib.txt.gz` or via function `getMITRIB()` provided in the package).
Moreover, counts from transcripts that mapped to a same gene were averaged.
Sample phenotypic annotations were also processed to define glycemic groups based on *hba1c* measurements, as per the original publication [@Fadista2014].

```{r, eval = !local_data}
# download GEO dataset from Github
eset <- readRDS('https://shenorrlab.github.io/bseqsc/data/GSE50244.rds')
```
```{r, eset}
# for this analysis we only look at samples with hba1c data
eset <- droplevels(eset[, !is.na(eset$hba1c)])
eset
```

### Single cell data
We generated single cell RNA-seq data of pancreatic islets from 3 healthy individuals, which provided gene expression profiles for 17434 genes in 7729 cells.
The raw counts for these data are available on the [data download
page](pages/data.html), in the form of an
`ExpressionSet`: 
```{r islet_eset, eval = !local_data}
eislet <- readRDS('https://shenorrlab.github.io/bseqsc/data/islet-eset.rds')
eislet
```

After identification of the cell types by clustering, we extract marker genes for each cell type, i.e. sets of genes that are mostly highly expressed in each cell type, which can be done in multiple ways, see for example details in [@Baron2016].

The `r BSEQSCpkg` package provides the list of selected marker genes for the pancreatic islets, stored in data object `pancreasMarkers`:

```{r p_markers}
data(pancreasMarkers)
str(pancreasMarkers)
```

*****
## Estimation of cell type proportions
Once cell type-specific markers have been identified, we build a reference basis matrix from their single cell gene expression, which we use to infer the composition of bulk tissue samples from their total gene expression.

### Building the reference basis matrix
The purpose of this matrix is to provide reference expression profiles that together discriminate the cell types of interest.
We compute it by averaging, within each cell type, the gene expression of each marker gene across all the cells in the cell type.

Notably, the computation is performed on scaled single cell count data, to better reflect cell type-specific variations in average transcript counts.
Indeed, looking at the average counts within each cell type shows how cell types have varying levels of RNA transcripts:
```{r cell_totals, out.width = '800px'}
# average counts computed within each cell type in each sample
plotCellTotals(eislet, 'cellType', 'sampleID')
```

Because we want the basis matrix of reference profiles to reflect these biological differences, we use these average counts to re-scale the CPM data before computing the average gene expression.
```{r basis_matrix, cache = TRUE}
B <- bseqsc_basis(eislet, pancreasMarkers, clusters = 'cellType', samples = 'sampleID', ct.scale = TRUE)
```

The resulting basis matrix should show a clear cell type-specific expression pattern (Figure 2):
```{r basis_heatmap, out.width = '800px', fig.cap='**Figure 2:** Basis matrix of gene expression reference profiles from pancreatic cell type derived from single cell RNA-seq data.'}
plotBasis(B, pancreasMarkers, Colv = NA, Rowv = NA, layout = '_', col = 'Blues')
```

### Estimating proportions
Cell type proportions are estimated using CIBERSORT, a method that was developed to investigate immune infiltrating landscape in tumour tissue [@Newman2015; @Gentles2015] (See section Requirements for setup instructions).

```{r bseqsc_proportions, cache = TRUE, collapse = TRUE}
fit <- bseqsc_proportions(eset, B, verbose = TRUE)
```

To facilitate plotting and model fitting in downstream analysis, we add the estimated proportions as extra phenotypic variables to the `ExpressionSet`:
```{r emix}
pData(eset) <- cbind(pData(eset), t(coef(fit)))
```

Estimated proportions can then be compared between the two HBA1C groups, within each cell type (Figure 3).
```{r prop_plot, out.width = '800px', fig.width = 16, fig.height = 3, echo = FALSE, fig.cap = '**Figure 3:** Estimated cell type proportions in each Normoglycemic (blue) and hyperglycemic (red) individuals.'}
local({
  df <- melt(pData(eset), measure.vars = rownames(coef(fit)))
  ggplot(df, aes(x = hba1c_class2, y = value, color = hba1c_class2)) + geom_boxplot(width = .3) + facet_wrap( ~ variable, nrow = 1, scales = 'free') + scale_x_discrete(breaks = NULL) + scale_color_brewer(type = 'qual', palette = 'Set1', direction = -1) + xlab('') +   ylab('Frequency') + guides(color = guide_legend(''))
})
```

*****
## Gene expression differences adjusted for proportions
Now that we have estimates for the main cell types, we can use `r Biocpkg('EdgeR')` to derive a set of genes that are likely to be differentially regulated at the cell type level.

To do so, we fit 2 models of gene expression differences between HBA1C groups (IGT, Diabetic, Normal):
  * a base model only includes the HBA1C group variable, as well as gender and age covariates;
  * an extended model that additionally includes the proportions of alpha, beta, ductal and acinar cells.

Genes that are likely to be differentially regulated at the cell type level are those whose p-values improve (i.e. decrease) when adjusting for cell type proportions (red dots on Figure 4).
Here we used a p-value threshold of 0.01.

To fit these models we use the utility function `fitEdgeR` included in the package:
```{r fitEdgeR}
fitEdgeR
```

```{r edgeR_models, cache = TRUE}
# base
fit_edger <- fitEdgeR(eset, ~ gender + ageN + hba1c_class, coef = c('hba1c_classIGT', 'hba1c_classDiabetic'))
# extended
fit_edger_ext <- fitEdgeR(eset, ~ gender + ageN + alpha + beta + ductal + acinar + hba1c_class
                          , coef = c('hba1c_classIGT', 'hba1c_classDiabetic'))
```

We then visualize 
```{r edgeR_agg, out.width = '800px', fig.cap = '**Figure 4:** P-values comparison between the 2 EdgeR modesl: Base model (y-axis) and extended model (x-axis).'}
# gather P-values from both models
edger_pvals <- cbind(fit_edger$table[, 'Symbol', drop = FALSE], Base = fit_edger$table$PValue
               , Adjusted = fit_edger_ext$table$PValue[match(fit_edger$table$Symbol, fit_edger_ext$table$Symbol)]
               , stringsAsFactors = FALSE)

edger_pvals <- mutate(edger_pvals, Regulated = Adjusted <= 0.01 & Adjusted <= Base)

# plot Base vs Adjusted
pvalueScatter(Base ~ Adjusted, edger_pvals, pval.th = 0.01, label.th = 3.5)
```

*****
## Cell type-specific differential expression

Now, we run a cell type-specific differential analysis using `csSAM` [@Shen-Orr2010], on the subset of selected genes, to which we add a set of 6 ER stress-related genes, that were identified in the single cell data analysis (see details in [@Baron2016]).

`csSAM` fits a linear model that estimates interaction terms between cell proportions and a group variable.
This imposes a constraint on the sample size within each group.
For this reason we only inlude into the model the 4 cell types (alpha, beta, ductal and acinar cells), either due to their variation (Figure 3) or biological relevance.

```{r bseqsc_csdiff, cache = TRUE, backspace = TRUE}
# ER genes
genes_ER <- c('HSPA5', 'MAFA', 'HERPUD1', 'DDIT3', 'UCN3', 'NEUROD1')
# Fit on ER stress genes and genes regulated beyond cell type proportion differences
genes <- union(genes_ER, subset(edger_pvals, Regulated)$Symbol)
# NB: we fix seed for reproducibility
csfit <- bseqsc_csdiff(eset[genes, ] ~ gender + ageN + hba1c_class2 | alpha + beta + ductal + acinar,  verbose = 2, nperms = 5000, .rng = 12345)
```

We observe a strong signal for up-regulation on beta cells (at FDR 0.05), and milder signals of up and down-regulation on acinar and beta cells respectively (at looser FDR 0.25).
```{r fdrPlot, out.width = '800px', warning = FALSE}
p <- plot(csfit, alt = c('great', 'less'), by = 'alt', xlab = 'Number of significant genes', ylab = 'FDR cutoff')
p$data$Cell.type <- gsub(":Hyper", "", p$data$Cell.type)
p$data$Alternative <- c('greater' = 'Up-regulated', less = 'Down-regulated')[as.character(p$data$Alternative)]
p + theme_get()
```

We filter these results for genes with the strongest signal in alpha, beta and acinar cells:
```{r csfilter, out.width = '800px', fig.height = 12}
cst <- csTopTable(csfit, threshold = c(alpha = 0.05, beta = 0.25, acinar = 0.25), alt = 'min', merge = TRUE, n = Inf, nodups = TRUE)
plotEffectSize(cst)
```

# Session Info
```{r sessionInfo, echo = FALSE}
sessionInfo()
```

# References
