Real data example

Vignette on Larsson 2019 data can be found here, which has allelic single-cell RNA-seq with 4 cell states.

Simulated data example I

The airpart package takes input data of counts from each of two alleles across genes (rows) and cells (columns) from a single-cell RNA-seq experiment.

For demonstration in the package vignette, we will simulate some data using makeSimulatedData function provided within the airpart package. We will examine the allelic counts and then perform QC steps before analyzing the data for allelic imbalance across groups of cells.

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)
## 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
## 
## ct1 ct2 ct3 ct4 
##  20  20  20  20
##       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.

## [1] "a1" "a2"
##  [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

Create allelic ratio matrix

In the preprocess step, we add a pseudo-count for gene clustering and visualization (not used for inference later on allelic imbalance though, which uses original allelic counts). From the heatmap, we can clearly identify the three gene clusters (across rows), and we also see cell type differences (across columns). Within each cell type, there are some cells with noisier estimates (lower total count) than others. Again, the allelic ratio tells us how much more of the a1 allele is expressed, with 1 indicating all of the expression coming from the a1 allele and 0 indicating all of the expression coming from the a2 allele.

Quality control steps

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.

## 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.

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.

Gene clustering

airpart provides a function to cluster genes by their allelic imbalance profile across cells (not using cell grouping information, e.g. sce$x). We then recommend providing genes within a cluster to the partition function. Clustering genes increases power for detecting cell type partitions, and improves speed as it reduces the number of times the partition must be estimated.

We provide two methods for gene clustering.

  1. Gaussian Mixture modeling

Gaussian mixture modeling is the default method for gene clustering. The scatter plot is shown based on top 2 PCs of the smoothed allelic ratio data. The argument plot=FALSE can be used to avoid showing the plot.

## model-based optimal number of clusters: 3

## [1] 25 25 25
  1. Hierarchical clustering