Ambient RNA removal

Droplets capture cell-free mRNA from the cell suspension together with their cell. This ambient RNA, the soup, comes from cells that lysed before or during loading, so every droplet holds some of it, including droplets without a cell. A cell then shows counts of genes it does not express, such as haemoglobin from red blood cells or immunoglobulins from plasma and B cells. These counts add apparent expression of other cell types’ markers and can blur differential expression.

CyteArc follows SoupX as Young and Behjati describe it Young & Behjati, 2020 in three steps:

1. qc.ambient_profile estimates the composition of the soup from droplets that hold no cell.

2. qc.ambient_contamination measures the fraction of each cell’s counts that come from the soup, using genes that some cells cannot express.

3. qc.correct_ambient removes each cell’s share of soup counts and writes the result as a new assay of integer counts.

The first two steps save immutable results with provenance. The third writes a new RNA assay that every analysis can use by name. The method needs the unfiltered droplets of each channel, not only the barcodes that Cell Ranger called as cells.

Download a raw and filtered matrix pair¶

We use the 10x Genomics dataset “1k PBMCs from a Healthy Donor (v2 chemistry)”, published under CC BY 4.0. The raw matrix (18 MB) holds every barcode with at least one count, and the filtered matrix (3 MB) holds the barcodes that Cell Ranger called as cells. We use the filtered matrix only for its list of cell barcodes.

The example uses a temporary directory. Use a project directory instead if you want to keep the downloaded files and Zarr store after closing the notebook.

Import every droplet¶

The soup profile needs droplets without a cell, so import the raw matrix rather than the filtered one. CrH5Reader(is_filtered=False) keeps the barcodes whose total count exceeds filtering_cutoff; a cutoff of 1 keeps every droplet with at least two counts, the smallest droplets that the profile uses. Most barcodes of a raw matrix hold only a few counts. The default feature floor of the first open would remove most of them, so it removes none and warns instead; min_features_per_cell=-1 keeps them all active without the warning.

The import writes the counts of tens of thousands of droplets and takes a few minutes.

The barcode rank plot orders the droplets by their total count. Cells form the plateau at the top. Below the steep drop lie droplets with soup alone, and the shaded band marks the 2 to 10 counts that the paper uses for the soup profile.

Estimate the soup profile¶

qc.ambient_profile sums the counts of the chosen droplets in one bounded pass and stores each gene’s share of them. Choose the droplets with total_counts=(low, high), an inclusive range of total counts, or with empty_droplets=, a cell selection of empty barcodes. The paper uses droplets with 2 to 10 counts and leaves out droplets with one count, whose barcodes may be sequencing errors; profiles from droplets with up to about 100 counts agreed well with it.

Abundant transcripts of every cell type, such as MALAT1, B2M, and ribosomal protein genes, dominate the soup. Transcripts of abundant cell types follow, so the soup resembles the pooled expression of the cells that lysed.

Cluster the cells¶

The contamination estimate leaves out every cell of a cluster that expresses a marker gene set, so it needs cluster labels for the cells. An RNA run on the called cells provides them; its cells= names the column that marks them. Filtering is switched off so that every called cell gets a label, and the run’s marker search names the clusters.

Estimate the contamination fraction¶

The contamination fraction is measured with genes that some cells cannot express. In a cell without endogenous expression of a gene set, every count of the set comes from the soup, so the ratio of its observed counts to the counts that a droplet of pure soup with the same total would hold is the contamination fraction. Haemoglobin genes, outside red blood cells, and immunoglobulin genes, outside B and plasma cells, are the usual choices for blood.

For each gene set, a cell expresses the set when its count is significantly larger than pure soup could give it: a Poisson test at the largest plausible contamination, max_contamination (1 by default), with a Benjamini-Hochberg false discovery rate of fdr (0.05). The estimate leaves out those cells and every other cell of their clusters, then pools the remaining cells and gene sets.

About one percent of each cell’s counts come from the soup in this dataset, and the two gene sets agree. The evidence shows which cells measured it. The B-cell cluster 4 expresses the immunoglobulin genes, so it is left out of that set. Cluster 7 has five cells that express the haemoglobin genes and five that express the immunoglobulin genes, so it is left out of both sets. That is why the per-cluster estimate below fails for it.

Leaving out whole clusters matters. Without clusters=, only the cells that pass the test are left out, and cells whose expression is real but too weak to pass it inflate the estimate:

level="cluster" estimates each cluster from its own cells instead. Each cluster then needs a gene set that it does not express, and a cluster that expresses every set raises an error that names it. The paper finds a global estimate sufficient in most cases.

Supply a known contamination fraction¶

Pass fraction= to assign a known fraction to every selected cell, as SoupX’s setContaminationFraction does. The example uses 0.01 (1%); choose a value supported by your experiment.

Pass supplied_fraction to qc.correct_ambient to apply that value. The correction below uses estimate, the fraction measured from the haemoglobin and immunoglobulin gene sets.

Remove the ambient counts¶

qc.correct_ambient removes each cell’s contamination fraction of its counts. The removal is the most likely one under a multinomial model with the soup’s gene fractions that never removes more of a gene than the cell holds: genes are emptied in increasing order of their count divided by their soup fraction, and the rest is taken from the other genes in proportion to their soup fraction. Genes absent from the soup keep every count.

Most genes are due far less than one count in a cell, so the removal must be rounded to whole counts. As in the paper, each gene’s removal is rounded up at random with a probability equal to its fractional part: a gene due 0.2 counts in each cell loses one count in about a fifth of the cells, so summed over the cells every gene loses its share of the soup, however small. The draws follow systematic sampling with one draw per cell, so each cell loses its share rounded down or up, and the corrected counts never exceed the raw counts. random_seed, 0 by default, sets the draws, so a repeated correction gives the same counts.

The new assay has a row for every barcode but measures only the corrected cells, which its membership column RNA_soupx_I records; the other rows hold no counts.

Compare markers before and after¶

Ambient counts appear as low expression of one cell type’s markers in other cell types. The monocyte markers S100A9 and LYZ are among the most abundant genes of the soup, while the immunoglobulin IGKC holds a smaller share of it.

The table shows the percentage of each cluster’s cells that hold any count of these markers, in the raw and the corrected counts.

In the lymphocyte clusters, fewer cells keep S100A9 and LYZ counts after the correction. Cells that express these genes at high levels, such as classical monocytes, keep them, while a cell with only a few counts can lose them, because the correction cannot tell a cell’s own few counts from the soup’s. IGKC changes little: only three cells lose it, one in cluster 1 and two in cluster 8. The correction takes counts in proportion to the soup, and with about one percent contamination a cell’s expected soup count of IGKC is far below one: a cell loses its IGKC count with that small probability, so a single IGKC count in a T cell usually stays. Summed over the cells, the correction removes the soup’s expected share of each gene, not every count of a marker in the wrong cell type.

Use the corrected counts¶

The corrected assay is an RNA assay with countsT, so every analysis takes it by name, such as ds.pipeline.rna(assay="RNA_soupx", cells=run), which analyzes the same cells as the run above, from_assay="RNA_soupx" in plots and marker searches, or ds.features.universe(from_assay="RNA_soupx"). Select the corrected cells: the assay did not measure the other barcodes, and analyses refuse them.

Corrected counts from other tools¶

CellBender removes ambient RNA with a generative model and writes its corrected counts in the Cell Ranger HDF5 layout. Import that file with cytearc.CrToZarr(cytearc.CrH5Reader("cellbender_output_filtered.h5"), zarr_loc=...) as you would a Cell Ranger matrix. Its counts are then ordinary counts of a store, without the provenance of the steps above.

Common mistakes and limitations¶

• A soup profile from droplets that hold cells is not a soup profile. Keep the count range below the cells, as the barcode rank plot shows, and import the raw matrix, not the filtered one.

• A gene set that every cluster expresses measures nothing: every cell is left out. Choose genes that are specific to a few cell types and abundant in the soup.

• Without clusters, weakly expressing cells raise the estimate. Pass a clustering of the cells.

• The correction removes counts in proportion to the soup’s composition; it cannot tell which molecule of a gene came from the soup. Genes that are abundant in the soup lose counts in every cell, including cells that express them.

• The rounding is random: another random_seed takes a gene’s rare soup counts from other cells. A conclusion that rests on a few counts of a gene should hold for another seed.

• Estimate each channel separately: each channel has its own soup. One profile and estimate correct one set of cells into one assay.

• The paper’s automatic estimate, its hierarchical model of cell-level fractions, and its correction of counts aggregated within clusters are not implemented.