Exploring Multiomic Cell Village Dataset

Authors: Evie Pless and Max Horlbeck
Published: Aug 4, 2026.

This annotated notebook shows an overview of how we explored a multiomic dataset (10x paired scRNA-seq and scATAC-seq data) from a cell village comprised of human-induced pluripotent stem cells (hiPSCs) derived from people with autism and unaffected controls. In particular, we demo how to use PRINT/seq2PRINT (Hu, Horlbeck, Zhang et al. 2025) to examine chromatin architecture and DNA-binding proteins. For more information about cell villages and the donor assignment pipeline, please see the accompanying blog post.

Thank you to Dr. Ralda Nehme and her team at the Broad Institute for providing the cell village and to the Cell Discovery Network for funding the sequencing.

import os
os.environ['SCPRINTER_DATA'] = '/lab-share/Gene-Horlbeck-e2/Public/.cache/scprinter/'
os.environ['POOCH_CACHE_DIR'] = '/lab-share/Gene-Horlbeck-e2/Public/pooch'
os.environ['XDG_CACHE_HOME'] = '/lab-share/Gene-Horlbeck-e2/Public/.cache'
import scanpy as sc
import scprinter as scp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

base_path = '/lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome'
pilot_path = f'{base_path}/SCBB-2493_NPC/pilot'

/lab-share/Gene-Horlbeck-e2/Public/scprinter_env/lib/python3.11/site-packages/sorted_nearest/__init__.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources

Explore RNA-seq dataset

First we load the filtered RNA-seq dataset and plot UMAPs that show lack of clustering by autism status or donor. However, UMAPs show strong gradients for genes related to neuronal differentiation state (e.g. MAP2) and proliferation (e.g. MKI67)

#Load filtered, processed RNA dataset 
adata = sc.read(f'{base_path}/rna_processed_full.h5ad')

# Plot UMAP by clustering
sc.pl.umap(adata, color='leiden', legend_loc='on data',
           title='Leiden clustering (res=0.5)')

# Plot UMAP by condition
sc.pl.umap(adata, color='condition', palette=['#1f77b4', '#ff7f0e'],
           title='ASD vs Control')

# Plot UMAP by donor
sc.pl.umap(adata, color='donor', legend_loc='right margin')

neural_progenitor_genes = ['MKI67','MAP2']

# Plot gene expression on UMAP
sc.pl.umap(
    adata, 
    color=neural_progenitor_genes,
    ncols=5,
    frameon=False,
    cmap='Reds'
)

Differentially accessible regions between high-MAP2 and low-MAP2

We defined an early commitment neuron-like group of cells (high MAP2, corresponding to Leiden clusters 1, 6, and 8), and then found differentially accessible regions between this group and the rest of the cells.

# Define 2 groups from Leiden clusters 
group1_clusters = {'1', '6', '8'} #High MAP2, neuron-like group
all_clusters = set(adata.obs['leiden'].unique())
group2_clusters = all_clusters - group1_clusters
print("Group 1 (mature neuron-like):", sorted(group1_clusters))
print("Group 2 (other):", sorted(group2_clusters))

adata.obs['atac_group'] = np.where(
    adata.obs['leiden'].isin(group1_clusters), 'group1',
    np.where(adata.obs['leiden'].isin(group2_clusters), 'group2', np.nan)
)
print(adata.obs['atac_group'].value_counts())

barcode_to_group = adata.obs['atac_group'].to_dict()

Group 1 (mature neuron-like): ['1', '6', '8']
Group 2 (other): ['0', '2', '3', '4', '5', '7']
atac_group
group2    2048
group1     447
Name: count, dtype: int64

# Load ATAC printer object
printer = scp.load_printer(f'{pilot_path}/autism_village_scprinter.h5ad', scp.genome.hg38)

# Build per-cell peak matrix
adata_peaks = scp.pp.make_peak_matrix(
    printer,
    regions=f'{pilot_path}/regions.bed',
    region_width=300,
    cell_grouping=None,   # per-cell, not pseudobulked
    group_names=None,
    sparse=True
)

# Attach group labels, keep only matched + labeled cells
adata_peaks.obs['atac_group'] = adata_peaks.obs_names.map(barcode_to_group)
adata_da = adata_peaks[adata_peaks.obs['atac_group'].isin(['group1', 'group2'])].copy()
print(adata_da.obs['atac_group'].value_counts())

# Filter low-count peaks, normalize 
sc.pp.filter_genes(adata_da, min_cells=10)   # "genes" = peaks here
sc.pp.normalize_total(adata_da)
sc.pp.log1p(adata_da)

# Differential accessibility 
sc.tl.rank_genes_groups(
    adata_da,
    groupby='atac_group',
    groups=['group1'],
    reference='group2',
    method='wilcoxon'
)

results = sc.get.rank_genes_groups_df(adata_da, group='group1')
results = results.sort_values('pvals_adj')
print(results.head(20))

loading da_group1_vs_group2_footprints /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/autism_village_scprinter_supp/da_group1_vs_group2_footprints.h5ad
loading control_vs_asd /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/autism_village_scprinter_supp/control_vs_asd.h5ad
Loading insertion profiles

Making peak matrix:   0%|          | 0/202497 [00:00<?, ?it/s]

atac_group
group2    1952
group1     429
Name: count, dtype: int64

... storing 'atac_group' as categorical

                           names    scores  logfoldchanges         pvals  \
202494     chr19:2540823-2541123 -6.427169       -2.082102  1.300024e-10   
202493   chr22:36329323-36329623 -5.883612       -1.841186  4.014093e-09   
202492    chr1:39646405-39646705 -5.816326       -1.278404  6.015503e-09   
202491     chr17:1714587-1714887 -5.549744       -1.553163  2.860878e-08   
202490   chr17:57899008-57899308 -5.455157       -0.982024  4.892962e-08   
202488  chr1:153566510-153566810 -5.315932       -1.327388  1.061125e-07   
202489   chr15:66701022-66701322 -5.349013       -0.973778  8.843530e-08   
202487     chr16:4299140-4299440 -5.312946       -1.191927  1.078669e-07   
202486   chr17:38451877-38452177 -5.279284       -1.409291  1.296895e-07   
202485      chr7:1537155-1537455 -5.191561       -1.692418  2.085381e-07   
202484   chr17:78324537-78324837 -5.155262       -1.568701  2.532767e-07   
202482   chr14:68793192-68793492 -5.122259       -1.150349  3.018966e-07   
202483   chr20:31968435-31968735 -5.131450       -1.581394  2.875182e-07   
0         chr8:79611880-79612180  5.055051        2.365459  4.302751e-07   
202480   chr15:39580765-39581065 -5.056486       -1.557807  4.270517e-07   
202481   chr12:53371260-53371560 -5.065212       -1.080108  4.079461e-07   
202479   chr17:43530730-43531030 -5.041555       -1.068174  4.617630e-07   
202477      chr4:1011349-1011649 -4.946076       -0.931118  7.572435e-07   
202476  chr2:127064882-127065182 -4.938242       -0.932245  7.882984e-07   
202478  chr9:137305033-137305333 -4.955733       -1.307658  7.205824e-07   

        pvals_adj  
202494   0.000026  
202493   0.000406  
202492   0.000406  
202491   0.001448  
202490   0.001982  
202488   0.002730  
202489   0.002730  
202487   0.002730  
202486   0.002918  
202485   0.004223  
202484   0.004662  
202482   0.004703  
202483   0.004703  
0        0.005446  
202480   0.005446  
202481   0.005446  
202479   0.005500  
202477   0.007981  
202476   0.007981  
202478   0.007981  

results.to_csv(f'{pilot_path}/da_peaks_group1_vs_group2_wilcoxon.csv', index=False)

Footprint the most differentially accessible peaks

Although we are able to observe differences in enhancer organization among conditions with equal total accessibility -- a key strength of PRINT/seq2PRINT -- we chose to focus on differntially accessible peaks in order to easily identify loci likely to exhibit differential footprints.

Interestingly, most of differentially accessible regions are less accessible in the high-MAP2 groups. We used PRINT to visualize the footprints from the top 7 most differentially accessible regions as well as the region that was most accessible in the high-MAP2 group compared to other cells (chr8:79611880-7961218).

The X-axis on the footprint plots shows genomic coordinate, and the Y-axis shows size of the footprint (0-200bp), where transcription factors are generally around 40bp, and nucleosomes are around 150bp.

Multiscale footprints (calculated with the scPrinter get_footprint_score function) are generated for 1000bp regions centered on each of the differentially accessible regions (300bp in length), hence the recentering required before plotting.

# Calculate multiscale footprints

# Config
N_TOP_REGIONS = 7
DA_RESULTS_CSV = f'{pilot_path}/da_peaks_group1_vs_group2_wilcoxon.csv'
SAVE_KEY = 'da_group1_vs_group2_footprints'

# ── 1. Load DA results and pick top 7 regions by adjusted p-value ───
da_results = pd.read_csv(DA_RESULTS_CSV)
da_results = da_results.sort_values('pvals_adj')
top_regions = da_results.head(N_TOP_REGIONS)

# 'names' column from rank_genes_groups_df holds the peak/var name,
# expected format "chr:start-end" (matching make_peak_matrix var_names)
# add the accessible region from high-MAP2 group
region_list = top_regions['names'].tolist()
region_list.append('chr8:79611880-79612180')

# ── 2. Build group1 / group2 barcode groupings ───────────────────────
# Rebuild the same RNA-derived group labels used for the DA test
adata_de = sc.read(f'{base_path}/rna_processed_full.h5ad')
group1_clusters = {'1', '6', '8'}
all_clusters = set(adata_de.obs['leiden'].unique())
group2_clusters = all_clusters - group1_clusters

adata_de.obs['atac_group'] = np.where(
    adata_de.obs['leiden'].isin(group1_clusters), 'group1',
    np.where(adata_de.obs['leiden'].isin(group2_clusters), 'group2', np.nan)
)
barcode_to_group = adata_de.obs['atac_group'].dropna().to_dict()

# Restrict to barcodes actually present in the printer object
printer_barcodes = set(printer.obs_names)
barcode_groups = pd.DataFrame(
    [(bc, grp) for bc, grp in barcode_to_group.items() if bc in printer_barcodes],
    columns=['barcode', 'group']
)

cell_grouping, group_names = scp.utils.df2cell_grouping(printer, barcode_groups)

# ── 3. Compute multiscale footprints for group1 vs group2 ────────────
printer.load_disp_model()
scp.tl.get_footprint_score(
    printer=printer,
    cell_grouping=cell_grouping,
    group_names=group_names,
    regions=region_list,
    region_width=1000,
    footprintRadius=None,
    flankRadius=None,
    modes=np.arange(2, 101),
    n_jobs=2,
    save_key=SAVE_KEY,
    backed=True,
    overwrite=True
)

estimated file size: 0.01 GB
Creating da_group1_vs_group2_footprints in printer.footprintsadata
obs=groups, var=regions
8000 100

Submitting jobs:   0%|          | 0/792 [00:00<?, ?it/s]

collecting multi-scale footprints:   0%|          | 0/792 [00:00<?, ?it/s]

import matplotlib.pyplot as plt
save_key = 'da_group1_vs_group2_footprints'

def recenter_region(region_str, width=1000):
    """Expand a peak string to `width` bp, keeping the same midpoint."""
    chrom, coords = region_str.split(':')
    start, end = map(int, coords.split('-'))
    center = (start + end) // 2
    new_start = center - width // 2
    new_end = new_start + width
    return f'{chrom}:{new_start}-{new_end}'

# ── Load top DA regions (original narrow peak coordinates) ──────────────────
da_results = pd.read_csv(f'{pilot_path}/da_peaks_group1_vs_group2_wilcoxon.csv')
da_results = da_results.sort_values('pvals_adj')
top_regions = da_results.head(7)

# Original peak strings (for labeling) and recentered 1000bp strings (for lookup)
orig_region_list = top_regions['names'].tolist()
orig_region_list.append('chr8:79611880-79612180')
region_list = [recenter_region(r, width=1000) for r in orig_region_list]

# ── Rebuild group1/group2 names as used during scoring ──────────────────────
group_names = ['group1', 'group2']  # confirm this matches what get_footprint_score used

# ── Plot ──────────────────────────────────────────────────────────────────
fig, axes = plt.subplots(len(region_list), 2, figsize=(12, 4 * len(region_list)))

stats_lookup = top_regions.set_index('names')[['pvals_adj', 'logfoldchanges']].to_dict('index')

for i, (orig_region, region) in enumerate(zip(orig_region_list, region_list)):
    for j, grp in enumerate(group_names):
        ax = axes[i, j]
        scp.pl.plot_footprints(
            printer,
            save_key=save_key,
            group_names=[grp],
            region=region,
            edge_mode='zeros',
            scales=None,
            ax=ax,
            cmap='Blues',
            vmin=0.5,
            vmax=2.0
        )        
plt.tight_layout()
plt.show()

It's satisfying to see that footprints in these differentially accessible regions are visibly different between the high-MAP2 group and the rest of the cells. We were particularly interested in the bottom figure (chr8:79611530-79612530), which corresponds to the top hit in terms of regions that are more accessible in the high-MAP2 group. This region falls within an enhancer in the pan-neuronal gene STMN2 – a gene that is significantly upregulated in the MAP2-high group. We can see group 2 has a strong nucleosome footprint on the left side of the region which is not present in group 1, potentially facilitating transcription factor binding in the locus.

# Clean up STMN2-enhancer footprint plot
regions_of_interest =  {'chr8:79611530-79612530'}

def recenter_region(region_str, width=1000):
    chrom, coords = region_str.split(':')
    start, end = map(int, coords.split('-'))
    center = (start + end) // 2
    new_start = center - width // 2
    new_end = new_start + width
    return chrom, new_start, new_end

# ─── Load footprint AnnData once ────────────────────────────────────────────
adata = printer.footprintsadata[save_key]
scales_all = np.array(adata.uns['scales'])
index = scales_all * 2  # bp, y-axis

edge_trim = 100  # set >0 if you want to zero out matrix edges like the example

# ─── Plot: rows = regions, cols = groups ────────────────────────────────────
fig, axes = plt.subplots(len(regions_of_interest), len(group_names),
                          figsize=(6 * len(group_names), 3.5 * len(regions_of_interest)))

for i, orig_region in enumerate(regions_of_interest):
    chrom, start, end = recenter_region(orig_region, width=1000)
    region_identifier = f'{chrom}:{start}-{end}'
    x_positions = np.arange(start, end)

    for j, group_name in enumerate(group_names):
        ax = axes[i, j] if len(regions_of_interest) > 1 else axes[j]

        try:
            select_group = adata.obs_ix(np.array([group_name]).astype('str'))
        except Exception:
            select_group = adata.obs.loc[[group_name]]['id']

        fp = adata.obsm[region_identifier][select_group][0].copy()
        if edge_trim > 0:
            fp[:, :edge_trim] = 0
            fp[:, -edge_trim:] = 0

        X, Y = np.meshgrid(x_positions, index)
        ax.pcolormesh(X, Y, fp, cmap='Blues', vmin=0.5, vmax=2.0,
                      shading='auto', rasterized=True)
        ax.set_ylabel('Window\n(bp)', fontsize=9)
        ax.set_title(f'chr8:79611630-79612430\n{group_name}', fontsize=10, fontweight='bold', loc='left')
        ax.set_ylim(index.min(), index.max())
        ax.set_xlim(79611530+100, 79612530-100)
        ax.set_xlabel('Position (bp)', fontsize=9)
        for spine in ax.spines.values():
            spine.set_visible(False)

plt.tight_layout()
plt.show()

Train seq2print models

To continue investigating the STMN2-enhancer region, we could try to line-up the above footprint plot with JASPAR transcription factor motifs, but it would still be hard to figure out exactly which transcription factors were binding in the region. Instead of relying on motif similarity alone, we could learn which sequences are actually driving the footprint by training and querying seq2PRINT models. Then we could use de novo motif calling to find out which transcription factor motifs are most similar to those important regions.

# Generate the config files for training seq2PRINT
import json
import os

group1_cells = adata_de.obs_names[adata_de.obs['atac_group'] == 'group1'].tolist()
group2_cells = adata_de.obs_names[adata_de.obs['atac_group'] == 'group2'].tolist()

group1_in_printer = [bc for bc in group1_cells if bc in printer.obs_names]
group2_in_printer = [bc for bc in group2_cells if bc in printer.obs_names]

os.makedirs(os.path.join(pilot_path, 'configs'), exist_ok=True)
os.makedirs(os.path.join(pilot_path, 'temp'), exist_ok=True)
os.makedirs(os.path.join(pilot_path, 'model'), exist_ok=True)

fold = 0
model_configs = []

# Peaks were called previously using scprinter call_peaks function with 'seq2print' default settings
for group_name, cell_list in [('highMAP2', group1_in_printer), ('lowMAP2', group2_in_printer)]:
    model_config = scp.tl.seq_model_config(
        printer,
        region_path=f'{pilot_path}/seq2print_cleaned_narrowPeak.bed',
        cell_grouping=cell_list,
        group_names=group_name,
        genome=printer.genome,
        fold=fold,
        overwrite_bigwig=False,
        model_name=group_name,
        additional_config={
            "notes": "MAP2_groups",
            "tags": ["ASD_multiome", group_name, f"fold{fold}"]
        },
        path_swap=(pilot_path, ''),
        config_save_path=f'{pilot_path}/configs/{group_name}_fold{fold}.JSON'
    )
    model_configs.append(model_config)
    print(f"✓ Config saved: {group_name}_fold{fold}.JSON")

Creating bigwig for highMAP2

  0%|          | 0/24 [00:00<?, ?it/s]

✓ Config saved: highMAP2_fold0.JSON
Creating bigwig for lowMAP2

  0%|          | 0/24 [00:00<?, ?it/s]

✓ Config saved: lowMAP2_fold0.JSON

#Generate launch command to train seq2PRINT

for group_name in ['highMAP2', 'lowMAP2']:
    scp.tl.launch_seq2print(
        model_config_path=f'{pilot_path}/configs/{group_name}_fold{fold}.JSON',
        temp_dir=f'{pilot_path}/temp',
        model_dir=f'{pilot_path}/model',
        data_dir=pilot_path,
        gpus=0,
        verbose=False,
        launch=False  # prints commands only — copy and run manually in GPU environment or submit via slurm
    )

CUDA_VISIBLE_DEVICES=0 seq2print_train --config /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/configs/highMAP2_fold0.JSON --temp_dir /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/temp --model_dir /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/model --data_dir /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot --project None
CUDA_VISIBLE_DEVICES=0 seq2print_train --config /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/configs/lowMAP2_fold0.JSON --temp_dir /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/temp --model_dir /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/model --data_dir /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot --project None

Generate sequence attribution scores

After locating the trained models, we can generate sequence attribution scores

MAP2_models = ['highMAP2_fold0-4sbo7apm.pt',
                  'lowMAP2_fold0-mfqnhw92.pt']
model_path = [os.path.join(pilot_path, "model", m) for m in MAP2_models]

# This will generate the sequence attribution scores for the footprint; you can input multiple models, and it will iterate over them
scp.tl.seq_attr_seq2print(
    genome=printer.genome,
    region_path=f'{pilot_path}/regions.bed',
    model_type='seq2print',
    model_path=model_path,
    gpus=[0,1],
    preset='footprint',
    overwrite=False,
    verbose=True,
    launch=False)

Using preset, the following parameters would be overwritten
using wrapper: just_sum
using nth_output: 0-30
using decay: 0.85
Please copy the following command in your terminal and run it to start the job
seq2print_attr --pt /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/model/highMAP2_fold0-4sbo7apm.pt --peaks /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/regions.bed --method shap_hypo --wrapper just_sum --nth_output 0-30 --gpus 0 1 --genome hg38 --decay 0.85 --save_key deepshap --model_norm footprint
Using preset, the following parameters would be overwritten
using wrapper: just_sum
using nth_output: 0-30
using decay: 0.85
Please copy the following command in your terminal and run it to start the job
seq2print_attr --pt /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/model/lowMAP2_fold0-mfqnhw92.pt --peaks /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/regions.bed --method shap_hypo --wrapper just_sum --nth_output 0-30 --gpus 0 1 --genome hg38 --decay 0.85 --save_key deepshap --model_norm footprint

# This will generate the sequence attribution scores for the count head
scp.tl.seq_attr_seq2print(
    genome=printer.genome,
    region_path=f'{pilot_path}/regions.bed',
    model_type='seq2print',
    model_path=model_path,
    gpus=[0,1],
    preset='count',
    overwrite=False,
    verbose=True,
    launch=False)

Using preset, the following parameters would be overwritten
using wrapper: count
using nth_output: 0
using decay: 0.85
Please copy the following command in your terminal and run it to start the job
seq2print_attr --pt /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/model/highMAP2_fold0-4sbo7apm.pt --peaks /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/regions.bed --method shap_hypo --wrapper count --nth_output 0 --gpus 0 1 --genome hg38 --decay 0.85 --save_key deepshap --model_norm count
Using preset, the following parameters would be overwritten
using wrapper: count
using nth_output: 0
using decay: 0.85
Please copy the following command in your terminal and run it to start the job
seq2print_attr --pt /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/model/lowMAP2_fold0-mfqnhw92.pt --peaks /lab-share/Gene-Horlbeck-e2/Public/autism_village_multiome/SCBB-2493_NPC/pilot/regions.bed --method shap_hypo --wrapper count --nth_output 0 --gpus 0 1 --genome hg38 --decay 0.85 --save_key deepshap --model_norm count

printer.close()