Preserve sparse cells for downstream analysis

[synatkeamsk]

Dear StereoPy Team,

Thank you for your recent assistance with my Stereo-seq analysis. I was able to successfully implement the batch effect correction pipeline on my dataset. However, I encountered an issue where T cells and B cells—being relatively sparse—were filtered out during the merging step.

Since these cell populations are biologically important for my study, I would appreciate any guidance on how to retain them for downstream analysis. I have attached the CD3D gene visualization from StereoMap for reference. My goal is to preserve these low-abundance cell populations throughout the analysis workflow.

For reference, I followed this batch effect pipeline:
https://stereopy.readthedocs.io/en/v1.5.1/Tutorials%28Multi-sample%29/Batch_Effect.html#

Thank you very much for your support.

[wanruiwen-genomics-cn]

Dear @synatkeamsk,

Thank you for reaching out and sharing the visualization!

Short answer: T cells and B cells are likely being removed by the preprocessing filtering steps (e.g., filter_cells, filter_genes, highly_variable_genes) rather than the merge itself. You can retain these rare populations by relaxing filtering thresholds and ensuring their marker genes are preserved.

---

Assessment

• Type: usage/question
• Severity: low
• Affected module: stereo/preprocess/filter.py, stereo/tools/highly_variable_genes.py

What Is Happening

In the batch effect correction pipeline, there are several steps that can reduce the representation of sparse cell populations:

1. filter_cells(min_genes=...) — Cells with very few detected genes (common for some immune cell subsets) will be removed if min_genes is set too high.
2. filter_genes(min_cells=...) — Genes expressed in very few cells (such as T cell or B cell markers like CD3D, CD19, CD79A) may be dropped if min_cells is too stringent.
3. highly_variable_genes(n_top_genes=2000) — Only the top 2000 most variable genes are selected by default. Marker genes for rare populations may not be among the most variable genes and could be excluded.
4. Gene intersection during merge — MSData.integrate() uses var_type='intersect' by default, keeping only genes shared across all samples. If a marker gene is absent from one sample, it will be dropped from all.

These are all working as designed — the filtering is intentionally aggressive to improve statistical power for the majority of cells, but it can inadvertently remove sparse populations.

Is This Expected?

Yes. The default filtering parameters are optimized for common cell types. Rare populations like T cells and B cells in spatial transcriptomics data (where immune cells can be very sparse) require adjusted thresholds.

Recommended Workflow

Step 1: Relax cell filtering thresholds

When calling filter_cells, use lower thresholds to retain cells with fewer detected genes:

# Use a lower min_genes threshold to keep sparse immune cells
ms_data.tl.filter_cells(min_genes=50)  # or even lower, depending on your data

Step 2: Relax gene filtering thresholds

Use a lower min_cells to retain genes expressed in few cells:

# Keep genes expressed in at least 3 cells (instead of higher defaults)
ms_data.tl.filter_genes(min_cells=3)

Step 3: Use gene union instead of intersection during merge

Set var_type to 'union' so that genes unique to some samples are retained (filled with zeros in other samples):

ms_data._var_type = 'union'
ms_data.integrate()

Step 4: Ensure marker genes survive HVG selection

After running highly_variable_genes, verify that key immune marker genes (CD3D, CD3E, CD19, MS4A1/CD20, etc.) are in the selected gene set. If not, consider increasing n_top_genes:

ms_data.tl.highly_variable_genes(n_top_genes=4000)

Alternatively, you can skip filter_by_hvgs and use all genes for downstream clustering, though this will be slower.

Alternative Approaches

If adjusting thresholds alone is not sufficient:

• Option A: Subset before filtering — Extract T/B cell regions using spatial coordinates (filter_coordinates) or known marker gene expression before running the full pipeline, then process them separately.
• Option B: Skip HVG filtering for clustering — Run clustering on all genes instead of only HVGs. This retains all signal but increases computation time.
• Option C: Manual gene list — Use filter_genes(gene_list=[...]) with a curated list of immune marker genes to ensure they are not dropped.

Useful References

Resource	Purpose
Batch Effect Tutorial	The pipeline you followed
Multi-sample Clustering Tutorial	Alternative multi-sample workflow

Notes

• When using var_type='union', missing genes in some samples will be filled with zeros. This is generally safe for downstream analysis but can affect normalization — re-normalize after merging.
• Consider inspecting how many T/B cells remain after each filtering step to identify which step causes the most loss.
• The CD3D visualization from StereoMap confirms the cells exist in your data — the issue is purely about filtering thresholds.

Please let us know if you have further questions!

Best regards,
Stereopy Maintainer

[wanruiwen-genomics-cn]

Dear @synatkeamsk,

Thank you for reporting this issue and for following the batch effect pipeline!

Short answer: The merging step in Stereopy does not inherently remove cells — sparse T cells and B cells are most likely lost during the explicit filter_cells QC step or become indistinguishable due to gene intersection reducing marker resolution. Adjusting filtering thresholds and gene merge strategy should help retain them.

---

Assessment

• Type: usage/question
• Severity: medium
• Affected module: stereo/preprocess/filter.py, stereo/utils/data_helper.py

What Is Happening

When running the batch effect pipeline, two mechanisms can cause low-abundance cell populations to appear "filtered out":

1. Explicit cell filtering (filter_cells): The filter_cells step removes cells below user-specified thresholds (min_counts, min_genes, pct_counts_mt, etc.). If thresholds are tuned for the majority cell types, sparse populations like T and B cells — which may have lower total UMI counts or fewer detected genes — can fall below the cutoff.

2. Gene intersection during merge: By default, MSData uses var_type='intersect', which keeps only genes present in all samples. If key marker genes (e.g., CD3D, CD3E for T cells; CD19, MS4A1 for B cells) are missing or lowly expressed in one sample, they may be dropped, making these populations harder to detect during clustering.

Neither integrate() nor batches_integrate() (Harmony) removes cells — they concatenate all cells and adjust PCA embeddings, respectively.

Is This Expected?

Yes — this is expected behavior when using strict QC thresholds or gene intersection. It is not a bug but a parameter tuning consideration.

Recommended Workflow

Step 1: Review and relax your filter_cells thresholds

Before merging, check the distribution of total counts and gene counts in your samples. For sparse immune cell populations, consider using more permissive thresholds:

# Example: more permissive filtering
data.tl.cal_qc()
data.tl.filter_cells(min_counts=50, min_genes=20, pct_counts_mt=20)

Compare the number of cells before and after filtering with data.shape to verify sparse populations are retained.

Step 2: Consider using var_type='union' instead of 'intersect'

When constructing the MSData object, you can use gene union to keep all genes across batches (missing genes are filled with 0):

ms_data = st.io.read_h5ms(h5ms_path=None, _data_list=[data1, data2], _var_type='union')
# Or set on an existing MSData:
ms_data._var_type = 'union'
ms_data.integrate()

This preserves batch-specific marker genes that are critical for identifying sparse populations.

Step 3: Verify cell retention after each step

Add checkpoints to monitor cell counts:

print(f"Before merge: {[d.shape[0] for d in [data1, data2]]}")
ms_data.integrate()
print(f"After merge: {ms_data.merged_data.shape[0]}")
ms_data.tl.filter_cells(min_gene=100, mode='integrate')
print(f"After filter: {ms_data.merged_data.shape[0]}")

Step 4 (optional): Use batch-aware HVG selection

When computing highly variable genes, pass the groups='batch' parameter to select HVGs per batch, which can help retain markers of sparse populations that are enriched in specific samples:

ms_data.tl.highly_variable_genes(groups='batch', mode='integrate')

Alternative Approaches

If relaxing thresholds is not sufficient:

• Option A: Pre-label sparse populations — Run clustering per sample before merging (mode='isolated'), annotate T/B cells, then merge. After merging, use the saved annotations rather than re-clustering.
• Option B: Isolated-mode filtering — Use mode='isolated' for filter_cells so each sample gets individually tuned thresholds, preventing one high-complexity sample from dominating the cutoff.

Notes

• The Harmony step (batches_integrate) only modifies the PCA embedding and never removes cells.
• After merge, per-sample annotations in cells.obs (except batch) are cleared — save important per-sample labels separately before integrating.
• If you see CD3D expression in StereoMap but cells disappear after the pipeline, the most likely cause is the QC filtering step, not the merge.

Please let us know if you have further questions!

Best regards,
Stereopy Maintainer