
Feature linkages are defined as pairs of genomic features, such as peaks and genes, that have significant correlation in signals across cells. Feature linkages provide a basis for inferring enhancer-gene targeting relationships and constructing transcriptional networks. The features with strong linkages are considered to be “co-expressed” and enrich for a shared regulatory mechanism. For example, the accessibility of an enhancer and the expression of its target gene can display a very synchronized differential pattern across a heterogeneous population of cells. A highly accessible enhancer increases transcription factor (TF) binding, which in turn activates (or represses) gene expression. Conversely, when the enhancer is inaccessible, no TF can bind to the enhancer and impact transcription, which leads to reduced target gene expression (Figure 1).

<Image src={"https://cdn.10xgenomics.com/image/upload/v1694206492/software-support/Cell%20Ranger%20ARC/feature_linkage_cartoon.png"} />

Figure 1. A schematic example of a feature linkage. The arc height denotes the magnitude of linkage correlation. The accessibility of the blue peak is much more correlated with the gene expression level than the green peak, hence the stronger feature linkage.

The computation of feature linkage consists of 2 steps (Figure 2)

1. Normalize feature-barcode matrix and perform matrix smoothing.
2. Compute feature linkage correlations and significance.

<Image src={"https://cdn.10xgenomics.com/image/upload/v1694206914/software-support/Cell%20Ranger%20ARC/algo_schema.png"} />

Figure 2. A schematic representation of the feature linkage algorithm Matrix normalization and smoothing (A) and computation of feature linkage matrices (B) is shown in the schematic.

<LinkableSection id={"matrix-normalization"} title={"Matrix normalization and smoothing"}>

The feature-barcode matrix is normalized using a method based on the depth-adaptive negative binomial normalization (DANB) that was originally proposed by Andrews et al. (Bioinformatics 2019). The ATAC peak-barcode matrix and GEX gene-barcode matrix are normalized independently to accommodate the different signal and sparsity profile of the two data types. The DANB method which was originally proposed for single cell gene expression data has been modified at the dispersion estimation step to make it compatible with both gene expression and ATAC data. Instead of the method of moments estimation of the dispersion in the original DANB method, we computed the maximum likelihood estimator for the dispersion, which has a better fitting performance for ATAC peak features.

Due to the sparsity of single cell data, especially the cut site count in ATAC peaks, signals of a peak and a gene will likely not be detected at the same time in one cell when both are expected to have high expression levels. As a result, a direct computation of correlation or other measure of dependence on the raw counts between two features is unable to distinguish highly co-expressed features from noise (Figure 3).

To clarify with an example, GZMB is a known NK and CD8 T cell marker in PBMC and its expression is highly enriched in NK and CD8 cytotoxic T cells. As the GZMB promoter displays significant NK and CD8 T cell-specific accessibility, the GZMB promoter and the gene are expected to display strong feature linkage (Figure 3A-B). However, because of the sparsity in single cell data, more than 80% of the cells (3094 out of 3799) have zero count in both GZMB gene expression and GZMB promoter. Moreover, more than 40% (355 out of 877) of annotated CD8/NK cells have zeros in both features (Figure 3A). The correlation of raw single-cell counts between GZMB gene expression and promoter accessibility is 0.285, which is lower than what we would expect based on known biology.

To overcome this barrier, we perform matrix smoothing. The value of a feature in a given cell is enhanced by “borrowing” values of the same feature from “neighboring” cells. Here, neighboring cells describe a population of cells whose GEX or ATAC profiles have high similarity, which is determined by applying K-nearest neighbor (KNN) algorithm on the gene expression PCA reduced dimensions. The information “borrowing” is achieved by a weighted sum of signals from 30 closest neighboring cells. The process of the information “borrowing” is called smoothing. After smoothing, putative co-expressed features show a much stronger correlation pattern compared to randomly selected pairs of features. In the GZMB example, the correlation of the smoothed values of gene expression and promoter accessibility is 0.87, with a clear visual pattern of correlation. The promoter of KHNYN, a nearby gene, does not have NK or CD8 T cell-specific differential accessibility. Instead, it has comparable accessibility in all major PBMC cell types. Hence, the linkage between KGNYN promoter and GZMB gene is weak, as evident from the low correlation value (Figure 3D).

The smoothing weight is defined as the Gaussian kernel transformation of the euclidean distance between two cells in PCA space. The Gaussian kernel transformation guarantees that the smoothing weight is high only when two cells have a highly similar gene expression profile and quickly decays to zero when the similarity between cells decreases. This is necessary to prevent over-smoothing in KNN-based methods (Andrews and Hemberg, F1000 Research 2019).

<Image src={"https://cdn.10xgenomics.com/image/upload/v1694207056/software-support/Cell%20Ranger%20ARC/matrix_smoothing.png"} />

Figure 3. Matrix smoothing improves interpretability of the linkage correlation. Examples of strong feature linkage (A, B) and weak feature linkage (C, D) from a PBMC dataset with 3,799 cells. Raw counts are colored based on the number of barcodes that share the same raw counts (A, C). Note that more than 3,000 cells have raw counts of (0, 0). Smoothed values are colored as the density in the 2-D scatter plot (C, D).

</LinkableSection>

<LinkableSection id={"compute-feature-linkages"} title={"Compute feature linkages"}>

Feature linkage correlation is determined by computing the Pearson correlation coefficient between two features after smoothing. Linkage correlation is a direct and interpretable measure of the linkage strength, and ranges from -1 to 1 (Figure 4). The sign of the correlation indicates a positive or negative association.

Feature linkage significance provides a measure of statistical uncertainty for feature linkage inference and offers another way to contrast strong linkages relative to weak linkages. A modified version of the Hotspot algorithm (DeTomaso et al. 2020) is applied to compute “local correlation” for a feature linkage and the value is transformed to a Gaussian random variable to enable hypothesis testing.

<Image src={"https://cdn.10xgenomics.com/image/upload/v1694207138/software-support/Cell%20Ranger%20ARC/linkage_dist.png"} />

Figure 4. A typical distribution of feature linkage correlation and significance.

</LinkableSection>

<LinkableSection id={"sparsity-generation"} title={"Sparsity generation of the feature linkage matrix"}>

Because the number of computable linkages increases quadratically with the number of features, it is expected that the majority of computable linkages are not biologically significant, and sparsity in the inference of feature linkage is expected. We applied a thresholding method in which linkages with significance < 5 are removed from the linkage matrix (Figure 4). The threshold is determined based on an analysis of continuously down-sampling reads and comparing the decay of linkage significance and correlation. A significance threshold of 5 strikes a good balance between retaining linkages with lower correlation scores and stability against downsampling.

</LinkableSection>

<LinkableSection id={"discussion"} title={"Discussion"}>

Feature linkages describe associations between genomic features. The biological, causal relationship between interactive enhancers or between enhancers and their targeting gene expressions are dependencies generating such associations. Several alternative molecular assays exist that examine these dependencies. For example, H3K27Ac ChIP-seq measures active enhancers, Hi-C measures chromosome conformation, and CRISPR-based enhancer screen tests whether an enhancer is required for target gene expression. Such alternative assays, individually or in combination, characterize enhancers and generate feature linkage candidates with different throughput and different levels of fidelity to the biological, causal interactions. However, when restricting the data source to transcriptome and chromatin accessibility profiling, the feature linkages we infer are candidates that are enriched for causal relationships but do not directly imply causality.

</LinkableSection>

<LinkableSection id={"references"} title={"References"}>

Andrews, T. S., & Hemberg, M. False signals induced by single-cell imputation. *F1000Research*, 7, 1–36 (2019).

Andrews, T. S., & Hemberg, M. M3Drop: dropout-based feature selection for scRNASeq. *Bioinformatics*, 35(16), 2865–2867 (2019).

DeTomaso, D., & Yosef, N. Identifying Informative Gene Modules Across Modalities of Single Cell Genomics. *BioRxiv*, 2020.02.06.937805 (2020).

</LinkableSection>


Algorithms for computation of feature linkages

Cell Ranger ARC Algorithm Overview


This step associates a subset of barcodes observed in the library to the nuclei loaded from the sample. Identification of these cell barcodes enables data analysis at single cell resolution. For each barcode, we have an associated count of transposition events in peaks (from the chromatin accessibility library), as well as gene expression UMIs (from the gene expression library). The cell-calling algorithm uses this paired information to identify the cell population from the non-cell background.

Prior to performing cell calling, barcodes must pass three filtering criteria.

**First filter - ATAC "low targeting" barcodes:** Identifies "low-targeting" barcodes in which the fraction of fragments overlapping called peaks is lower than the fraction of genome in peaks (for this calculation, peaks are padded by 2,000 bp on both sides to account for fragment length). These barcodes typically have their cut sites randomly distributed over the genome, are not targeted to be enriched near functional regions, and do not exhibit the typical ATAC-seq nucleosome phasing signal.

**Second filter - gel bead doublets:** Identifies putative gel bead doublets where a partition contains one cell and two barcoded gel beads. These cells then manifest as two barcodes of the same cell type in the dataset. The rate of such doublets is low and the presence of these few extra barcodes doesn't significantly affect secondary analysis such as clustering or differential analysis, although it can potentially inflate abundance measurements of very rare cell types. The chromatin accessibility library is used to identify a minor-major pair of barcodes (B1, B2) that are part of a putative gel bead doublet by observing if the pair of barcodes shares more genomically adjoining "linked" fragments (fragments sharing a transposition event) with each other (B1-B2) as opposed to themselves (B1-B1 or B2-B2). The minor barcode is identified as the one with fewer fragments.

**Third filter - greater than 1 count observed in each library:** Identifies barcodes that have been observed in both chromatin accessibility and gene expression libraries. The minimum threshold is set at >1 count in each library.

The first two filters rely solely on data from the chromatin accessibility library, and are often referred to as ATAC exclusions. A barcode that fails any of the above three filtering criteria is masked from the total set of barcodes prior to cell calling.

We then perform cell calling on the remainder barcodes. The algorithm has four key steps:

1. **De-duplication**: The cell-calling algorithm places each barcode into a 2-D space defined by their individual ATAC and gene expression counts. Barcodes with identical coordinates i.e. ATAC and gene expression counts are collapsed into a single measurement to generate a more uniform density of points across the 2-D space. This de-duplication step de-emphasizes over-represented low-count barcodes and allows suppression of noise without using thresholds or making assumptions about the count distribution profiles.
2. **Ordmag-derived initial grouping**: The de-duplicated count data is partitioned using thresholds derived from “ordmag” into cell and non-cell populations. Ordmag is a published algorithm that finds a threshold that is 10-fold less than the maximum value after removing outliers (Zheng, et al., 2017). A threshold is defined independently for each dimension and barcodes above both ATAC and gene expression thresholds are labeled as cells, with the remainder labeled as non-cells.
3. **K-means boundary refinement**: Since “ordmag” is insufficient to capture the spread of paired ATAC and gene expression data, we refine the boundaries of these initial set of cells using K-means. The K-means algorithm is initialized using centroids calculated from the ordmag-defined cell and non-cell groups. The K-means algorithm is implemented using this single initialization and set to K=2.
4. **Map classification to de-duplicated barcodes**: The K-means assignment of the de-duplicated count data into cells and non-cells are used to classify the full set of non-excluded barcodes. Barcodes that were masked during the de-duplication step are assigned to the cell or non-cell group based on the K-means classification of their counterpart with identical ATAC and gene expression counts.

Cell calling is limited to produce < 20k cells as the assay is currently designed to support 500-10k cells. To override the cell-caller, two additional parameters need to be provided to `cellranger-arc count`: `--min-atac-count=N` and `--min-gex-count=N`, the former specifies the minimum number of transposition events in peaks for a cell barcode and the latter specifies the minimum number of UMI counts for a cell barcode. Any barcode that satisfies *both* thresholds will be classified as a cell, unless it is excluded from the cell calling algorithm because it fails the above specified filtering criteria.

In this example, the cell caller identifies 3009 cells as shown in the [web summary](/support/software/cell-ranger-arc/latest/analysis/web-summary#joint-view-cfe384). The toggle spike lines feature is useful in determining appropriate thresholds to define a cell. With override parameters set to values `--min-atac-count=1690` and `--min-gex-count=570`, barcodes in the shaded yellow area are selected as cells, while barcodes in the shaded blue region are treated as background.

<Image src={"https://cdn.10xgenomics.com/image/upload/v1694205925/software-support/Cell%20Ranger%20ARC/force-cells.png"} />

#### References

Zheng, G.X., Terry, J.M., Belgrader, P., Ryvkin, P., Bent, Z.W., Wilson, R., Ziraldo, S.B., Wheeler, T.D., McDermott, G.P., Zhu, J. & Gregory, M.T. Massively parallel digital transcriptional profiling of single cells. *Nature communications* 8, 1-12 (2017).


Calling Cell Barcodes

Cell Ranger ARC is an advanced analytical suite designed for the Chromium Single Cell Multiome ATAC + Gene Expression sequencing. It provides in-depth analysis of gene expression and chromatin accessibility at a single cell level, uniquely linking these aspects for enhanced genomic understanding.

Overview

Getting Started

Downloads & Installation

Analysis

Tutorials

Algorithms

Advanced

Troubleshooting

Release Notes

Datasets

Loupe Browser

Cell Ranger ARC

A tutorial on using Cell Ranger ARC mkref to create a reference for Rattus norvegicus

Custom Reference for Rat

Learn more about input files required to run Cell Ranger ARC

Inputs

Generating FASTQs (mkfastq)

Specifying FASTQs

FASTQs with Illumina Software

Custom Reference (mkref)

Learn more about how to run Cell Ranger ARC

Running Pipelines

Analyzing GEX + ATAC (count)

Aggregating Multiple Samples (aggr)

Custom Analysis (reanalyze)

Learn more about the output files generated by Cell Ranger

Outputs

Web Summary

ATAC + Gene Expression Analysis

Feature-Barcode Matrices

Per Barcode Metrics

Per Barcode QC Metrics

ATAC Barcoded BAM

ATAC Fragments

ATAC Peaks

ATAC Peak Annotations

ATAC Transposition Counts

GEX Barcoded BAM

GEX Molecule Info (H5)

HDF5 Feature Barcode Matrix Format

What is Cell Ranger ARC?

We provide resources from protocols, user guides, and free-to-use analysis tools

Algorithms overview

Advanced pipeline content for Cell Ranger ARC.

View the latest updates to Cell Ranger ARC.

Browse Q&As from the 10x Genomics support team.

Access free public Chromium data from 10x Genomics.

Loupe Browser is a desktop application that allows you to interactively visualize data.

Algorithms Overview