Demo Notebook: Using scBoolSeq with scanpy (Binarisation)¶

In this notebook we showcase:

• How to use scanpy, to pre-process scRNA-seq data in order to be used with scBoolSeq. scanpy is a popular package for scRNA-seq data analysis in Python (part of the scverse). This tutorial builds on scanpy's preprocessing and clustering tutorial.
• The integration of scBoolSeq with scanpy.
• How obtain per-cell, coarse-grained gene activities using scBoolSeq.
• How to aggregate these coarse-grained cells to obtain Boolean signatures for cell clusters automatically detected with the Louvain algorithm (using scanpy).

This notebook can be seen as a complement to the presented Case-Study on Early-born Retinal Neurons.

The data can be obtained from its corresponding GEO accession GSE122466.

The reference of the original study which published the data:

Giudice, Q. Lo, Leleu, M., Manno, G. La, & Fabre, P. J. (2019).
Single-cell transcriptional logic of cell-fate specification and axon guidance in early-born retinal neurons. 
Development (Cambridge), 146(17). https://doi.org/10.1242/dev.178103

Read-in the raw dataset¶

Standard Single-Cell RNA-seq Quality Control and Preprocessig Pipeline using scanpy¶

scanpy provides numerous functions for scRNA-seq data analysis. Here we perform a standard pipeline which consists of:

0. Preprocessing:
   1. Analysing the percentage of mitochondrial genes.
   2. Doublet detection.
   3. Normalisation and log-transformation.
   4. Selecting Highly Variable Genes.
1. Performing a Principal Component Analysis (PCA) on top of previously selected Highly Variable Genes.
2. Computing a nearest neighbour distance matrix plus a neighbourhood graph of observations.
3. Performing a UMAP embedding.
4. Unsupervised community detection using the Louvain algorithm.

And finally:

6. The binarisation of scRNA-seq data with scBoolSeq

Basic QC: Counts and Percentage of Mitochondrial Genes¶

Doublet detection¶

Normalisation¶

Filter out low-quality cells¶

We are happy with these results because the previous normalisation and quality control were effective, we can judge this because:

1. The two batches (Lanes) have similar forms and are superimposed on one another.
2. Common technical biases do not fully account for the observed topologies on the projected space:
   • There is no cluster determined by the doublet_score or by the percentage of mitochondrial counts.
   • Some clusters have a higher sequencing depth but this does not represent a major axis of variation.
3. The projection coincides with that of the original article (see section Inferring Boolean signatures

Binarisation of Highly Variable Genes with scBoolSeq¶

As explained in the main text, scBoolSeq should not be applied on the whole dataset but rather on the subset of "Highly Variable Genes" which are deemed to be explicative of the phenomenon of interest.

These can be accessed using the "highly variable" annotation of the var attribute of an AnnData object as follows:

Create a scBoolSeq instance and fit it.¶

By overriding the parameter zeroes_are=0 of scBoolSeq, cells with zero counts for genes classified as "ZeroInflated" will have their coarse-grained activity assigned to be inactive rather than unknown. This is a modelling choice that can be changed depending on the application.

As shown by the ipython magic command %% time, estimating all parameters needed by scBoolSeq to coarse-grain experimental data and generate synthetic counts from Boolean dynamics is relatively fast.

Treating a reference dataset comprised of 4717 cells and 1980 genes took a little over a minute Wall time: 1min 15s on modern laptop with an Intel i7-10875H (16) @ 5.100GHz processor.

It is important to note that scBoolSeq objects are pickable. Even though fitting the instance (computing the classification parameters) is relatively fast, there is no reason to compute them more than once:

• Save a scBoolSeq instance:

with open("scbool_GSE122466.pkl", "wb") as _out_f:
    pickle.dump(scbool, _out_f)

• Reload a scBoolSeq instance for use in another session:

with open("scbool_GSE122466.pkl", "rb") as _in_f:
    scbool = pickle.load(scbool, _in_f)

• Furthermore, scBoolSeq instances require little disk space:

$ du -sh scbool_GSE122466.pkl
# 808K scbool_GSE122466.pkl

Some genes selected as being Highly Variable by scanpy have been rejected for binarisation by our classification scheme.

Coarse-grain/Binarise expression data¶

As seen here, coarse-graining expression data is almost instant.

One might notice that numerous genes have NaN everywhere. This is expected as 387 genes were classified as being Discarded.

Inferring the Boolean signature of Louvain clusters¶

One may want to pool groups of cells together to perform analyses on some sort of consensus Boolean signature, rather than a single cells. For this purpose, there is a module scboolseq.meta containing classes and functions facilitating the agreggation of Boolean cells. By defining these meta-observations of Boolean values, one can reason in terms of the common/consensus activation state of cell clusters.

scboolseq.meta.CellAggregator has two main arguments:

• binary_df, the first argument with coarse-grained cells (the output of scboolseq.scBoolSeq().fit_transform(log_expression)).
• groups, the second argument, which should be either a pandas GroupBy object or a dictionnary whose keys are the group names and values are pandas.Index.

Combining these two, the aggregator has numerous properties representing different aggregation strategies:

Here, we demonstrate what we call binary strategies. These use the following statistics to yield a consensus signature:

We can asssess if these strategies permit identifying the groups (in this case cell clusters determined via the Louvain algorithm: scanpy.tl.louvain(adata)

The mode strategy corresponds to taking the by-group most frequent coarse-grained value, for each gene. It seems like this strategy allows yields a high average number of positive markers (genes being classified as active within the group and in no other group).

To better assess the quality of these markers, we can add them to the adata.obs DataFrame. We can subsequently visualise the number of markers with the corresponding clusters on the UMAP projection.

We observe that some clusters have few Boolean markers. Extremes representing the starting point of the differentiation process (Early Retinogenesis) and the final phenotypes have a higher number of markers, while transient clusters have less extreme expression values and will therefore have less genes whose coarse-grained activity enables their identification.

These can be compared with the figure published in the original article:

Giudice, Q. Lo, Leleu, M., Manno, G. La, & Fabre, P. J. (2019).
Single-cell transcriptional logic of cell-fate specification and axon guidance in early-born retinal neurons. 
Development (Cambridge), 146(17). https://doi.org/10.1242/dev.178103