Simulation set-up

The simulated example dataset has 3 gene clusters with differential allelic imbalance (DAI):

• the first cluster has pairs of cell types with same allelic ratio with 0.2 and 0.8 (larger DAI)
• the second cluster has balanced allelic ratio
• the third cluster has pairs of cell types with same allelic ratio with 0.7 and 0.9 (smaller DAI)

Below we specify a number of simulation settings as arguments to the simulation function:

• the “noisy” cell count is 2
• the normal cell count is 10
• 4 cell types
• 20 cells within each cell type
• 25 genes within each gene cluster
• overdispersion parameter theta in rbetabinom is 20 (higher is less dispersion)

library(airpart)
suppressPackageStartupMessages(library(SingleCellExperiment))
p.vec <- rep(c(0.2, 0.8, 0.5, 0.5, 0.7, 0.9), each = 2)
set.seed(2021)
sce <- makeSimulatedData(
  mu1 = 2, mu2 = 10, nct = 4, n = 20,
  ngenecl = 25, theta = 20, ncl = 3,
  p.vec = p.vec
)

unique(rowData(sce)) # the true underlying allelic ratios

## DataFrame with 3 rows and 4 columns
##              ct1       ct2       ct3       ct4
##        <numeric> <numeric> <numeric> <numeric>
## gene1        0.2       0.2       0.8       0.8
## gene26       0.5       0.5       0.5       0.5
## gene51       0.7       0.7       0.9       0.9

table(sce$x) # counts of each cell type

## 
## ct1 ct2 ct3 ct4 
##  20  20  20  20

assays(sce)[["a1"]][1:5, 1:5] # allelic counts for the effect allele

##       cell1 cell2 cell3 cell4 cell5
## gene1     0     2     0     1     0
## gene2     0     1     0     0     2
## gene3     0     0     1     0     0
## gene4     0     2     0     0     0
## gene5     0     1     0     1     0

Required input data

In summary, airpart expects a SingleCellExperiment object with:

• discrete cell types recorded as a variable x in the colData(sce)
• effect and non-effect allelic counts as assays a1 and a2

The allelic ratio is calculated as a1 / (a1 + a2).

Note: We assume that the cell types have been either provided by the experiment, or identified based on total count. We assume the allelic ratio was not used in determining the cell groupings in x.

assayNames(sce)

## [1] "a1" "a2"

sce$x

##  [1] ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1 ct1
## [20] ct1 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2 ct2
## [39] ct2 ct2 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3 ct3
## [58] ct3 ct3 ct3 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4 ct4
## [77] ct4 ct4 ct4 ct4
## Levels: ct1 ct2 ct3 ct4

QC on cells

We recommend both QC on cells and on genes. We begin with cell allelic ratio quality control. For details on these metrics, see ?cellQC.

cellQCmetrics <- cellQC(sce, mad_detected = 4)
cellQCmetrics

## DataFrame with 80 rows and 7 columns
##               x       sum  detected spikePercent filter_sum filter_detected
##        <factor> <numeric> <numeric>    <numeric>  <logical>       <logical>
## cell1       ct1   2.19590        61            0       TRUE            TRUE
## cell2       ct1   2.23045        64            0       TRUE            TRUE
## cell3       ct1   2.09342        57            0       TRUE            TRUE
## cell4       ct1   2.18184        59            0       TRUE            TRUE
## cell5       ct1   2.19590        63            0       TRUE            TRUE
## ...         ...       ...       ...          ...        ...             ...
## cell76      ct4   2.88309        75            0       TRUE            TRUE
## cell77      ct4   2.83187        73            0       TRUE            TRUE
## cell78      ct4   2.86153        75            0       TRUE            TRUE
## cell79      ct4   2.89763        74            0       TRUE            TRUE
## cell80      ct4   2.83123        73            0       TRUE            TRUE
##        filter_spike
##           <logical>
## cell1          TRUE
## cell2          TRUE
## cell3          TRUE
## cell4          TRUE
## cell5          TRUE
## ...             ...
## cell76         TRUE
## cell77         TRUE
## cell78         TRUE
## cell79         TRUE
## cell80         TRUE

Now define cell filtering automatically or users can manually filter out based on sum,detected and spikePercent.

keep_cell <- (
  cellQCmetrics$filter_sum | # sufficient features (genes)
    cellQCmetrics$filter_detected | # sufficient molecules counted
    # sufficient features expressed compared to spike genes,
    # high quality cells
    cellQCmetrics$filter_spike
)
sce <- sce[, keep_cell]

QC on genes

We also recommend QC on genes for allelic ratio analysis. Note that we require genes to be expressed in at least 25% of cells within each cell type and the genes to have high allelic imbalance variation across cell types. The following code chunk is recommended (not evaluated here though). If users want to estimate homogeneous cell type allelic imbalance, they can set sd = 0 and examine the below summary step to find interesting gene clusters with weighted mean deviating from 0.5.

featureQCmetric <- featureQC(sce)
keep_feature <- (featureQCmetric$filter_celltype &
  featureQCmetric$filter_sd &
  featureQCmetric$filter_spike)
sce <- sce[keep_feature, ]

Simple summary table of allelic ratio

We first quickly look at the weighted mean of allelic ratio for each gene cluster. From this step we will identify the interesting gene clusters. The mean is calculated, weighting the information from each gene x cell element of the matrices by the total count.

summary <- summaryAllelicRatio(sce)
summary

## $`gene cluster 1 with 25 genes`
##     x weighted.mean      mean        var
## 1 ct1     0.2084022 0.2056936 0.06355708
## 2 ct2     0.1981464 0.1879503 0.05634768
## 3 ct3     0.8133540 0.8090369 0.05572735
## 4 ct4     0.8054237 0.8136062 0.05143953
## 
## $`gene cluster 2 with 25 genes`
##     x weighted.mean      mean        var
## 1 ct1     0.5119248 0.5096426 0.09555925
## 2 ct2     0.4949174 0.5000794 0.09566643
## 3 ct3     0.4961165 0.5135395 0.09642713
## 4 ct4     0.5048356 0.4883957 0.08759873
## 
## $`gene cluster 3 with 25 genes`
##     x weighted.mean      mean        var
## 1 ct1     0.7085658 0.7216087 0.07258102
## 2 ct2     0.6942398 0.6803795 0.08580046
## 3 ct3     0.9117748 0.8952586 0.03821103
## 4 ct4     0.8956159 0.9017706 0.02942587

The following step is a complement of the QC on genes step. We recommend users only run airpart when the largest ordered allelic ratio difference > 0.05 for speed concerns. We find that the allelic ratio of most of the gene clusters in such cases (small absolute allelic ratio differences) won’t provide enough evidence to detect differential allelic imbalance.

sapply(1:length(summary), function(i) {
  inst <- summary[[i]]
  inst_order <- inst[order(inst$weighted.mean), ]
  max(diff(inst_order$weighted.mean)) > 0.05
})

## [1]  TRUE FALSE  TRUE

Experiment-wide beta-binomial over-dispersion

We recommend examining the experiment-wide beta-binomial over-dispersion, which helps to inform whether to use a binomial likelihood or a nonparametric approach to partitioning the cell types by allelic imbalance.

We focus on the first gene cluster (if a gene cluster is not provided, estDisp will choose the largest cluster).

The blue trend line gives the typical values of over-dispersion across all the genes in the cluster, and across all the cell types (accounting for differences across the cell types in the expected ratio).

estDisp(sce, genecluster = 1)

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Modeling using fused lasso with binomial likelihood

airpart offers a method for partitioning cell types using the generalized fused lasso with binomial likelihood, as implemented in the smurf package. Cell types are merged based on their similarity of allelic ratios, accounting for excess variance on the ratio from low counts. The penalization is determined using deviance on held-out data, with a 1 SE cross-validation rule for favoring smaller models (more fused cell types).

The fusion step can also taken into account both cell-level and gene-level baseline effects, through the use of a formula (see ?fusedLasso for example).

sce_sub <- fusedLasso(sce,
  model = "binomial",
  genecluster = 1, ncores = 1,
  niter = 2
)

The partition groups and the penalty λ from the fused lasso are stored in the metadata:

knitr::kable(metadata(sce_sub)$partition, row.names = FALSE)

metadata(sce_sub)$lambda

##     part1     part2 
## 0.2053138 0.2053138

Above, ncores is the number of CPU used for parallelization. As a guide, one can specify niter=5 when the cts weighted allelic ratio difference is smaller than 0.1, in order to provide additional estimator robustness.

sce_sub <- consensusPart(sce_sub)
knitr::kable(metadata(sce_sub)$partition, row.names = FALSE)

Modeling using pairwise Mann-Whitney-Wilcoxon extension

An alternative to the fused lasso with binomial likelihood is an extension we have implemented wherein all pairs cell types are compared with Mann-Whitney-Wilcoxon rank sum tests. In practice, we find that when the allelic counts deviates strongly from a binomial (e.g. large over-dispersion, small values of theta), the wilcoxExt function can offer improved performance, in terms of recovery of the true partition of cell types by allelic imbalance. The partition is decided based on a loss function motivated by the Bayesian Information Criteria.

thrs <- 10^seq(from = -2, to = -0.4, by = 0.2)
sce_sub_w <- wilcoxExt(sce, genecluster = 1, threshold = thrs)
knitr::kable(metadata(sce_sub_w)$partition, row.names = FALSE)

metadata(sce_sub_w)$threshold

## [1] 0.01

Calculating allelic ratio estimates via beta-binomial model

After airpart determines a partition of cell types either by the fused lasso with binomial likelihood or the nonparametric approach described above, it uses those fused lasso estimates or weighted means as the center of a Cauchy prior for posterior estimation of allelic ratios per cell type and per gene. Posterior mean and credible intervals are provided. The posterior inference makes use of a beta-binomial likelihood, and a moderated estimate of the over-dispersion. The prior from the partition and the moderated estimate of over-dispersion are provided to the apeglm function from the Bioconductor package of the same name.

Note that the estimates and credible intervals are not equal for cell types in the same partition and for genes, because in this step we re-estimate the conditional cell type means per cell type (not per partition) and account for each gene’s moderated estimate of over-dispersion.

sce_sub <- allelicRatio(sce_sub, DAItest = TRUE)
makeForest(sce_sub, showtext = TRUE)

## svalue shown in columns per cell type

Allelic ratio estimates (ar) as well as svalue and credible interval (lower and upper) are stored in rowData. Can use extractResult function to derive them.

genepoi <- paste0("gene", seq_len(5))
ar <- extractResult(sce_sub)
knitr::kable(ar[genepoi,])

makeStep(sce_sub[genepoi,])

s <- extractResult(sce_sub, "svalue")
apply(s[genepoi,],2, function(s){s<0.005})

##        ct1  ct2  ct3  ct4
## gene1 TRUE TRUE TRUE TRUE
## gene2 TRUE TRUE TRUE TRUE
## gene3 TRUE TRUE TRUE TRUE
## gene4 TRUE TRUE TRUE TRUE
## gene5 TRUE TRUE TRUE TRUE

adj.p <- mcols(sce_sub)$adj.p.value
adj.p < 0.05

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE