---
title: "Abundance Normalization"
sidebar_position: 3
format: pixelgen-docusaurus-md
---

# Abundance Normalization

As with any high-throughput techniques, normalization and denoising of Proximity Network abundance data can be a useful tool to make
samples and datasets more comparable to each other as well as help with the biological interpretation of your data.
Therefore in this tutorial we will introduce and discuss [*DSB*](https://www.nature.com/articles/s41467-022-29356-8),
our method of choice for the normalization and denoising of mixed cell population data as well as a more general
approach in the Centered-Log-Ratio normalization([*CLR*](https://pubmed.ncbi.nlm.nih.gov/28759029/)).

We will begin with a short introduction to *DSB* and *CLR* and the algorithms behind them. Then, you will learn how to apply these
techniques to a Proximity Network Assay dataset. We will finish with some words of caution when applying these -or any- type of normalization.

After completing this tutorial, you should be able to:

- Understand the motivation for and the algorithm behind *DSB* and *CLR* normalizations
- Apply these normalization methods to a Proximity Network Assay dataset
- Evaluate the effect of normalization on marker abundance
- Understand the potential pitfalls of normalization

## Why denoise and normalize? {#why-denoise-and-normalize}

By their very nature antibody-based assays can be affected by variable amounts of nonspecific binding. In leukocytes,
this nonspecific binding is partly a result of nonspecific binding of antibodies to Fc-receptors present on the surface
of many white blood cell types. These background levels of antibody binding (or "noise") can fluctuate between samples
and cells as a result of sample type and quality, preprocessing conditions, antibody isotype and other issues.
Ultimately, this can lead to difficulties distinguishing *negative* and *positive* subpopulations of cells expressing a
certain marker, complicating cell type annotation and differential abundance analysis.

Therefore, during the sample preprocessing of Proximity Network Assay samples, several blocking steps are included to reduce the possibility
of nonspecific binding. In addition, computational methods such as the one presented here can, to some extent, be
utilized to reduce the bias of nonspecific binding.

## The *DSB* algorithm {#the-dsb-algorithm}

The *DSB* method was originally developed to denoise and normalize antibody-based assays on droplet-based platforms such
as [CITE-seq](https://www.nature.com/articles/nmeth.4380). It is based on the detection of background marker levels in
so-called "empty" droplets (droplets encapsulating biological material but not a corresponding cell). Proximity Network Assay does not need
cell compartmentalization and therefore does not produce empty droplets as such. However, it has been shown that the
background abundance levels of a marker can be inferred from its abundance distribution across all cells in a mixed cell
population sample (see [DSB paper](https://www.nature.com/articles/s41467-022-29356-8)). By calculating the mean
abundance in this background population per marker and subtracting these values from the marker's abundance, each
marker's background population becomes effectively zero centered. In a subsequent step, the isotype control markers and
a cell-specific background mean are combined in a "noise component" that is regressed out of the zero centered data
corrected for background nonspecific levels.

To summarize these different steps:

1. Transform the count data using $log1p$ transformation.
2. Fit a finite mixed gaussian model with two components (i.e. two normal distributions) to each marker across all
   cells. And as a background correction step, subtract the mean of the first component on the logged data. Note: if no
   two component fit was possible, background correction is skipped.
3. Next, fit a finite mixed gaussian model with two components to each cell across all markers. The mean of the first
   component is the cell-wise background mean. This value is combined with the corrected counts of the isotype controls
   to form a "noise matrix". The first principal component of this noise matrix is regressed out.
4. The resulting residuals are returned as normalized and denoised values.

## The *CLR* algorithm {#clr}

When dealing with pure cell populations such as cell lines or FACS-sorted populations, *DSB* can sometimes result in
unexpected values as it can not reliably detect positive and negative populations for certain markers in these uniform
samples. Therefore, we also include the *CLR* method as an alternative normalization approach. This method was
originally developed for CITE-seq data ([Stoeckius *et al.*, 2017](https://www.nature.com/articles/nmeth.4380)) but
quickly gained traction in other antibody-based assays. It relies on the assumption that the geometric mean of marker
counts per cell is a good reference for normalizing the counts. Here we implemented it as the $log1p$ of the marker
counts minus the geometric mean of the $log1p$ counts per cell.

## Setup {#setup}

To run this tutorial, several python libraries and functions need to be imported.

```{python}
#| label: Abundance Normalization - Python - Setup
import numpy as np
import pandas as pd
import scanpy as sc
from pathlib import Path
from pixelator import read_pna as read
from pixelator.mpx.plot import density_scatter_plot
from pixelator.common.statistics import clr_transformation, dsb_normalize
import seaborn as sns
import matplotlib.pyplot as plt
sns.set_style("whitegrid")
```

```{python}
#| label: Define DATA_DIR - hidden
#| include: false

import os
from pathlib import Path

DATA_DIR = Path(os.environ["QUARTO_DATADIR"])
BASEURL = "https://pixelgen-technologies-datasets.s3.eu-north-1.amazonaws.com/pna-datasets/v1"
```


```{python}
#| label: Download data - hidden
#| include: false

!curl -L -O -C - --create-dirs --output-dir {DATA_DIR} "{baseurl}/PNA062_unstim_PBMCs_1000cells_S02_S2.layout.pxl"
!curl -L -O -C - --create-dirs --output-dir {DATA_DIR} "{baseurl}/PNA062_PHA_PBMCs_1000cells_S04_S4.layout.pxl"
```

```{python}
#| label: Define DATA_DIR
#| eval: false

DATA_DIR = Path("<path to the directory to save datasets to>")
```

## Load Data {#load-data}

To illustrate the normalization procedure we continue with the quality checked and filtered cells we identified in the previous tutorial. We recall that the filtered cells come from combination of a resting PBMC sample and a PHA stimulated PBMC sample.

```{python}
#| label: Quality Control - Python - Load data

pg_data_combined_pxl_object =  read([
    DATA_DIR/ "PNA062_unstim_PBMCs_1000cells_S02_S2.layout.pxl",
    DATA_DIR/ "PNA062_PHA_PBMCs_1000cells_S04_S4.layout.pxl",
])
filtered_cells = list(pd.read_csv(DATA_DIR / "qc_filtered_cells.csv")["component"])
pg_data_combined_pxl_object = pg_data_combined_pxl_object.filter(components=filtered_cells)
# In this tutorial we will mostly work with the AnnData object
# so we will start by selecting that from the pixel data object
pg_data_combined = pg_data_combined_pxl_object.adata()
```

## Normalizations {#dsb_normalize}

*DSB* and *CLR* procedures have been wrapped into single functions in **pixelator**, respectively called `dsb_normalize` and `clr_transformation`. Since *DSB* uses
the antibody isotype controls to perform the normalization, we need to specify these to `dsb_normalize`.
To facilitate
later comparisons, we also construct a layer with $log1p$ transformed data.

```{python}
#| label: Abundance Normalization - Python - Normalization
# Add dsb values as a layer to the object
pg_data_combined.layers["dsb"] = dsb_normalize(
    pg_data_combined.to_df(), 
    isotype_controls = ["mIgG1", "mIgG2a", "mIgG2b"]
    )
# Add clr values as a layer to the object
pg_data_combined.layers["clr"] = clr_transformation(
    pg_data_combined.to_df(), axis=1, non_negative=False
)
# Create log1p layer for later comparisons
pg_data_combined.layers['log1p'] = np.log1p(pg_data_combined.to_df())
```

## Evaluate normalization {#evaluate-normalization}

After having run the normalization step it is a good idea to visualize the overall distribution of some markers of
interest with different normalizations. This will allow you to assess in more detail how each procedure has modified
the abundance levels of markers. Let's have a look at CD19 as an example.

```{python}
#| label: Abundance Normalization - Python - Evaluation
def marker_ridgeplot(adata, layer, marker, title):
    marker_data = adata.to_df(layer)["CD19"]
    conditions = adata.obs['condition']
    plot_data_log = pd.DataFrame({'CD19' : marker_data, 'Conditions' : conditions})
    # Initialize the FacetGrid object
    sns.set_theme(style="white", rc={"axes.facecolor": (0, 0, 0, 0), 'axes.linewidth':2})
    pal = sns.color_palette("Set2", 12)
    g = sns.FacetGrid(plot_data_log, row="Conditions", hue="Conditions", aspect=9, height=1.2, palette=pal)
    # Draw the densities in a few steps
    g.map(sns.kdeplot, "CD19", bw_adjust=.5, clip_on=False, fill=True, alpha=1, linewidth=1.5)
    # Define and use a simple function to label the plot in axes coordinates
    def label(x, color, label):
        ax = plt.gca()
        ax.text(0, .2, label, fontweight="bold", color=color,
                ha="left", va="center", transform=ax.transAxes)
    g.map(label, "CD19")
    # Set the subplots to overlap
    g.figure.subplots_adjust(hspace=-.5)
    # Clean up plot
    g.set_titles("")
    g.set(yticks=[], ylabel="")
    g.despine(bottom=True, left=True)
    g.fig.suptitle(title)
    plt.tight_layout()

pg_data_combined.obs["condition"] = "Resting"
pg_data_combined.obs.loc[
    pg_data_combined.obs["sample"].str.contains("PHA"),
    "condition"
] = "PHA stimulated"
# log1p ridgeplot data
g = marker_ridgeplot(pg_data_combined, layer="log1p", marker="CD19", title='Log1p Transformation')

# DSB ridgeplot
marker_ridgeplot(pg_data_combined, layer="dsb", marker="CD19", title='DSB Transformation')

# CLR ridgeplot
marker_ridgeplot(pg_data_combined, layer="clr", marker="CD19", title='CLR Transformation')
```

CD19 is a cell type marker that is specific for B cells. The population of B cells in these samples should be relatively
small and the majority of cells should therefore not express CD19 at a meaningful level. However, as you can see in the
upper ridge plot prior to normalization, the majority of cells do express low levels of CD19 as a result of nonspecific
binding. In addition, there is clearly some variability between the levels of this nonspecific background noise. As you
can see in the middle ridge plot, the *DSB* normalization centers this negative background population around 0, while
retaining a small CD19+ population of B cells. This normalization and denoising procedure will thus make it easier to
annotate cell clusters and gate for cell population as we will show in more detail in the next section. Similarly in the lower ridge plot, with CLR transformation, the background population has mostly moved over to the negative region.

## Zero-centering background population simplifies downstream analysis {#zero-centering-background-population-simplifies-downstream-analysis}

One of the benefits of this type of normalization is that it facilitates detecting and annotating specific cell
populations in heterogeneous cell type samples. As an example we can create a UMAP of the two samples and visualize
some well-known cell-type specific markers.

```{python, fig-show = 'hold', fig-width = '49%'}
#| label: Abundance Normalization - Python - UMAP
#| layout-ncol: 2

sc.pp.highly_variable_genes(
    pg_data_combined, flavor="seurat_v3", n_top_genes=50
)
pg_data_combined.obsm["pca"] = sc.tl.pca(
    pg_data_combined.layers["dsb"], random_state=42
)
sc.pp.neighbors(
    pg_data_combined,
    n_neighbors=30,
    n_pcs=10,
    use_rep="pca",
    metric="cosine",
)
sc.tl.umap(
    pg_data_combined,
    min_dist=0.3,
)
markers_of_interest = ["CD4", "CD8", "CD19", "CD16"]
```

### DSB normalized abundance across clusters
```{python, fig-show = 'hold', fig-width = '49%'}
#| label: Abundance Normalization - Python - UMAP colored with DSB abundance
#| layout-ncol: 2
sc.pl.umap(pg_data_combined, color=markers_of_interest, ncols=4, layer = "dsb")
```

### Log1p abundance across clusters
```{python, fig-show = 'hold', fig-width = '49%'}
#| label: Abundance Normalization - Python - UMAP colored with Log1p abundance
#| layout-ncol: 2
sc.pl.umap(pg_data_combined, color=markers_of_interest, ncols=4, layer = "log1p")
```

Here, we assess the distribution of abundance of a few cell specific markers: CD4 is positive for T helper cells and
monocytes, CD8 for cytotoxic T cells, CD19 for B cells, and CD16 mainly for monocytes. We can see that the protein abundance is
restricted to one or few clusters for the *DSB* transformed counts in the top row, while the log transformed counts in the
bottom row show an "ambient" level across most clusters. As you can appreciate from these plots the *DSB* normalized
values present a much clearer distinction between positive and negative populations than the $log1p$ transformed values.
Besides these benefits for visualization and annotation, *DSB* and *CLR*-transformed values allow for applying more uniform gating
strategies across markers.

```{python, fig-show = 'hold', fig-width = '49%'}
#| label: Abundance Normalization - Python - Gating
#| layout-ncol: 2
fig, _ = density_scatter_plot(pg_data_combined, "CD3e", "CD4", layer="log1p")
fig.axes[0].set_title("Log1p transformation")

fig, _ = density_scatter_plot(pg_data_combined, "CD3e", "CD4", layer="dsb")
fig.axes[0].set_title("DSB transformation")

fig, _ = density_scatter_plot(pg_data_combined, "CD3e", "CD4", layer="clr")
fig.axes[0].set_title("CLR transformation")
```

## Remarks {#remarks}

As any normalization method, *DSB* has some potential pitfalls that can be avoided by a thorough understanding of your
experimental design and analysis plan. Here we list a few recommendations and good practices:

- The *DSB* normalization method works optimally for markers that have a positive and a negative population in the
  experimental setup; i.e. in datasets consisting of different cell populations (such as PBMCs). This maximizes the
  chances that the GMM fitting will reliably detect the negative background population. Consequently, markers such as
  HLA-ABC that have a relatively consistent high expression across all cell types might display some unexpected
  behaviour (i.e. these can show a low *DSB* value despite being highly expressed in all cells). If your dataset mainly
  consists of very similar cells, such as for example pure cell lines, it might be prudent to compare the *DSB* output
  to log- or *CLR*-normalization.
- Do not apply normalization blindly and always be aware of the effect of a transformation when using normalized marker
  counts for downstream analysis. A good practice is to visualize the distribution of markers you deem important in your
  analysis before and after normalization with, for example, violin plots. Does this correspond to what you expect and
  if not, ask yourself why could this be the case?

::: {.callout-warning}

## Differential Abundance testing on normalized values {#differential-abundance-testing-on-normalized-values}

Make sure the normalized values produced by *DSB* and/or *CLR* are fitting for the test you would like to perform! Some
downstream differential abundance tests demand strictly positive values as input. For example, the `rank_genes_group`
function from Scanpy requires that the data provided to them is strictly non-negative and do not provide Log Fold
Changes for markers that are on average negative in one or more groups. If you would like to use these function as part
of your downstream analysis, we suggest to shift the abundance of markers with negative values to strictly positive.
This would leave the relative differences between groups unaffected and allow to proceed with the testing procedure.
This procedure has been described in the following [Scanpy issue](https://github.com/scverse/scanpy/issues/653).

```{python}
#| eval: false
#| label: Abundance Normalization - Python - Shifted Normalized Values

# Example of code to add a layer of shifted normalized values to the adata object
# Calculate minimum value per marker
mins = pg_data_combined.to_df("dsb").min()
# If minimum value is positive, set to zero
mins[mins > 0] = 0
# Add layer to use in Differential Abundance test
pg_data_combined.layers["shifted_dsb"] = pg_data_combined.to_df("dsb").sub(mins)
```

:::
