Single-cell ATAC-sequencing analysis

In this use case, we reproduce part of the analysis from the study “HOX13-dependent chromatin accessibility underlies the transition towards the digit development program” (Desanlis et al., 2020, Nature Communications) to demonstrate how Drylab can replicate published computational workflows.

Objective:
The transition from proximal to distal limb identity and the formation of digits require precise chromatin remodeling. 
HOX13 proteins (HOXA13 and HOXD13) are known to be essential for digit development, but their exact role at the chromatin level remained incompletely understood. 

Methods Summary

• Data source: GSE145657 — two .snap files (GSM4323464, KO; GSM4323465, WT), extracted from a 2.5 GB tarball (Desanlis et al., Nat Commun (2020)) to provide a comprehensive, reproducible view of HOX13-dependent chromatin accessibility changes in the developing mouse forelimb bud.

• QC: Three-threshold evaluation (lenient / standard / strict) on log10(total fragments), TSS-region signal (SE), and FRIP (PL/PP). Standard threshold selected.

• Dimensionality reduction: Feature selection (top 50,000 variable bins), spectral embedding (50 components, SnapATAC2 tl.spectral), UMAP at k = 15, 30, 50.

• Clustering: Leiden algorithm grid search over k in {15, 30, 50} and resolution in {0.2, 0.4, 0.6, 0.8, 1.0, 1.5} (18 combinations). Stability assessed via silhouette score and 5-fold 80% subsampling ARI.

• Statistical testing: Global chi-squared test; per-cluster two-sided Fisher's exact test with Benjamini-Hochberg FDR correction; log2 odds ratio with 1,000-iteration bootstrap 95% CI.

Key Findings

1. QC Filtering

The standard threshold was selected, retaining 5,581 WT and 3,856 KO cells with median FRIP approximately 0.92, indicating high library quality.

2. Dimensionality Reduction

• Spectral embedding yielded 46 non-trivial components; elbow detection placed the dominant biological signal at component 5.

• UMAP topology was stable across k = 15, 30, 50, confirming robustness of the manifold structure.

Clustering Stability Grid Search

The optimal parameters (k = 50, resolution = 0.6) achieved the highest silhouette score and subsampling ARI simultaneously, yielding 12 clusters. [,]

3. WT vs KO Distribution

Global chi-squared test: statistic = 2,631.6, df = 11, p approximately 0 — highly significant genotype-cluster association.

Per-cluster Fisher's exact test (BH-FDR corrected):

Interpretation

The scATAC-seq landscape resolves 12 stable chromatin accessibility states in the developing limb bud , consistent with the multiple mesenchymal clusters reported in the original publication. The moderate overall silhouette score (0.172) is expected for scATAC-seq data with continuous developmental transitions, yet the high subsampling ARI (0.802) confirms that these 12 clusters represent reproducible chromatin states rather than noise-driven partitions.

Cluster C4 is the most HOX13-dependent population. C4 contains 859 KO cells and 0 WT cells, representing a completely KO-exclusive chromatin state (log2 OR = -11.64, FDR = 0) [,]. This corresponds to the KO-specific cluster (cluster 3 in the original numbering) described by Desanlis et al., which displayed a proximal-like but distinct chromatin signature that emerges only in the absence of HOX13 Desanlis et al., Nat Commun 2020.

Cluster C1 represents the HOX13-dependent digit program. C1 contains 1,645 WT cells and only 1 KO cell (log2 OR = +10.07, FDR = 0), identifying it as the chromatin state that requires HOX13 for establishment.

This aligns with the original study's distal-specific clusters (clusters 2 and 6) that were enriched for the HOXA13 motif and associated with digit morphogenesis gene ontology terms, and that were absent in Hox13-/- samples.

The extreme bimodal genotype segregation — with KO cells completely excluded from C1 and shunted into C4 — demonstrates that HOX13 loss does not uniformly perturb all cell states but instead prevents cells from accessing the digit-specific chromatin program and diverts them into an aberrant KO-specific configuration . The original paper explicitly noted that "chromatin accessibility did not simply switch to a proximal limb bud profile," consistent with the observation here that C4 is a distinct state rather than a simple reversion.

These data support a model in which HOX13 acts as a master regulator — and likely a pioneer transcription factor — of the digit-specific chromatin accessibility program, with its absence preventing the opening of distal limb-specific loci. The moderately affected clusters (C0: 55% KO, FDR = 3.9e-39; C5: 49.5% KO, FDR = 1.5e-6; C7: 78.5% WT, FDR = 1.5e-18) reflect intermediate progenitor states where HOX13 activity modulates but does not completely determine chromatin accessibility. The KO-exclusive C4 population may represent an alternative, partially proximal-like chromatin configuration that is normally suppressed by HOX13 activity during the critical window of digit morphogenesis.

The cluster numbering in this reanalysis (C1, C4) does not directly correspond to the original study's labels (clusters 2/6 for WT-distal, cluster 3 for KO-specific), as expected given differences in preprocessing parameters, bin size, and clustering resolution. Nonetheless, the biological phenomena — a WT-exclusive digit chromatin state and a KO-exclusive aberrant state — are fully reproduced.

Recommended Next Steps

• Differential accessibility analysis: Apply findDAR or diff_test between C1 (digit program) and C4 (KO-aberrant) to identify the specific genomic loci whose accessibility is HOX13-dependent.

• Motif enrichment: Perform motif enrichment analysis on differentially accessible regions to confirm HOX13 motif prevalence at C1-specific peaks and identify compensatory TF motifs enriched in C4.

• Integration with scRNA-seq: Link chromatin accessibility changes to gene expression (e.g., Bmp2, Shh pathway genes) via gene activity scores or multi-omic integration.

• Pseudotime trajectory: Model the developmental trajectory from shared progenitor states toward the digit program (C1) versus the KO-aberrant state (C4) to identify the chromatin bifurcation point.