Make qc consensus floor changeable

Author: cakirb

README.md

@@ -74,6 +74,7 @@ The pipeline requires the following primary inputs, typically configured via com
 
 *   Gather mode (`--gather_mode`): `auto`, `starsolo`, or `cellbender` (default: `cellbender`). In `auto` mode, samples annotated as `cell` use `cellbender`; samples annotated as `nuclei` use `starsolo`.
 *   QC mode (`--qc_mode`): `original`, `multires`, or `combined` (default: `original`).
+*   Consensus floor (`--consensus_floor`): minimum consensus fraction for multires/combined consensus rescue (default: `0.90`).
 *   CellTypist model (`--celltypist_model`): `gut` preset or a path to a custom model.
 *   Integration HVG strategy (`--hvg_strategy`): `cell_nuclei_union` or `global` (default: `cell_nuclei_union`).
 *   Optional batch-aware HVG selection (`--hvg_batch_key`): obs column passed to `scanpy.pp.highly_variable_genes(batch_key=...)`.
@@ -396,8 +397,8 @@ This step requires the output of `gather_matrices` step which is the h5ad object
 
 After optional CellTypist annotation, `run_qc` applies the automatic QC to each sample. Three QC modes are supported and selectable via `qc_mode`:
 * `original` mode uses the main automatic QC logic (summarised [here](https://teichlab.github.io/sctk/notebooks/automatic_qc.html)). It (i) computes the eight QC metrics, (ii) builds a QC embedding (PCA → neighbors → UMAP → clustering), then (iii) loops over mitochondrial upper bounds (20, 50, 80). For each bound it refits Gaussian Mixture Models (GMMs) for four core metrics (n_counts, n_genes, percent_mito, percent_spliced) (default thresholds for these four core metrics differ if the sample is single-nuc `nuclei` or single-cell `cell`) to derive sample‑specific pass ranges, assigns per‑cell pass/fail across all required metrics, and calls clusters good if ≥50% of their cells pass (via `min_frac`). The loop produces cellwise and clusterwise pass flags per mito threshold; final summary columns (`pass_auto_filter`, `good_qc_cluster`) encode the most stringent (lowest mito) threshold each cell/cluster still satisfies (smaller value = stricter pass).
-* `multires` mode keeps the same metric computation and initial QC embedding but fixes the mitochondrial bound at 20 (no mito loop). It runs the GMM cellwise QC once to obtain `cell_passed_qc`, then performs multi‑resolution clustering (grid 0.1–1.0) applying clusterwise QC at each resolution. From these runs it records, for every cell, a `consensus_fraction` (fraction of resolutions where its cluster passes) and selects a consensus cut‑off that maximizes the Jaccard overlap of failing cells versus the cellwise GMM calls (explicit search over unique consensus_fraction values). Cells with consensus_fraction ≥ chosen threshold become `consensus_passed_qc`; a degeneracy guard forces all False if the fractions carry no signal, and an additional `keep_multires` flag rescues high‑agreement cells (≥0.90) or cluster passes. This yields per‑cell QC flags less sensitive to one arbitrary clustering resolution.
-* `combined` mode nests the multi‑resolution procedure inside the mitochondrial loop (20, 50, 80). For each mito bound it refits the GMMs (adjusting percent_mito upper limit), performs cellwise QC (`good_qc_cell_mitoX`), runs multi‑resolution cluster QC (producing `good_qc_cluster_mitoX`, `consensus_fraction_mitoX`), searches a consensus threshold (same Jaccard strategy) to define `consensus_passed_qc_mitoX`, and applies a consensus fraction floor (0.90) plus degeneracy checks. Pass propagation merges stricter passes upward (20 → 50 → 80) and collapses to numeric summaries (`pass_auto_filter`, `consensus_pass_auto_filter`) representing the most stringent mitochondrial window retained. This delivers robust QC calls across both clustering resolution and mitochondrial burden.
+* `multires` mode keeps the same metric computation and initial QC embedding but fixes the mitochondrial bound at 20 (no mito loop). It runs the GMM cellwise QC once to obtain `cell_passed_qc`, then performs multi‑resolution clustering (grid 0.1–1.0) applying clusterwise QC at each resolution. From these runs it records, for every cell, a `consensus_fraction` (fraction of resolutions where its cluster passes) and selects a consensus cut‑off that maximizes the Jaccard overlap of failing cells versus the cellwise GMM calls (explicit search over unique consensus_fraction values). Cells with consensus_fraction ≥ chosen threshold become `consensus_passed_qc`; a degeneracy guard forces all False if the fractions carry no signal, and an additional `keep_multires` flag rescues high‑agreement cells controlled by `--consensus_floor` or cluster passes. This yields per‑cell QC flags less sensitive to one arbitrary clustering resolution.
+* `combined` mode nests the multi‑resolution procedure inside the mitochondrial loop (20, 50, 80). For each mito bound it refits the GMMs (adjusting percent_mito upper limit), performs cellwise QC (`good_qc_cell_mitoX`), runs multi‑resolution cluster QC (producing `good_qc_cluster_mitoX`, `consensus_fraction_mitoX`), searches a consensus threshold (same Jaccard strategy) to define `consensus_passed_qc_mitoX`, and applies the `--consensus_floor` consensus fraction floor plus degeneracy checks. Pass propagation merges stricter passes upward (20 → 50 → 80) and collapses to numeric summaries (`pass_auto_filter`, `consensus_pass_auto_filter`) representing the most stringent mitochondrial window retained. This delivers robust QC calls across both clustering resolution and mitochondrial burden.
 
 A few more details on `run_qc` step:
 - Thresholds for the four core metrics (n_counts, n_genes, percent_mito, percent_spliced) to use in GMM are predefined in this step. Different thresholds are used if a sample is single-nuc `nuclei `or single-cell `cell`.

bin/qc.py

@@ -174,7 +174,7 @@ def run_qc_mito_loop_original(ad, qc_metrics, metrics_custom, threshold):
     for max_mito in reversed(MITO_THRESHOLDS):
         ad.obs.loc[ad.obs[f"good_qc_cluster_mito{max_mito}"], "good_qc_cluster"] = max_mito
 
-def run_qc_multi_res(ad, qc_metrics, metrics_custom, threshold):
+def run_qc_multi_res(ad, qc_metrics, metrics_custom, threshold, consensus_floor):
     """Multi-resolution QC approach with consensus across clusterings."""
 
     sk._pipeline.generate_qc_clusters(ad, metrics=qc_metrics)
@@ -229,14 +229,13 @@ def run_qc_multi_res(ad, qc_metrics, metrics_custom, threshold):
         ad.uns["consensus_threshold_degenerate"] = False
 
     # Optional stricter rescue: require a high raw fraction on top of kp9's boolean.
-    CONS_FLOOR = 0.90   # adjust (0.85–0.95) as you like
     ad.obs["keep_multires"] = (
         ad.obs["cluster_passed_qc"] |
-        (ad.obs["consensus_passed_qc"] & (frac >= CONS_FLOOR))
+        (ad.obs["consensus_passed_qc"] & (frac >= consensus_floor))
     )
 
 
-def run_qc_combined(ad, qc_metrics, metrics_custom, threshold):
+def run_qc_combined(ad, qc_metrics, metrics_custom, threshold, consensus_floor):
     """Combined original + multi-resolution QC across mitochondrial thresholds."""
 
     for max_mito in MITO_THRESHOLDS:
@@ -309,8 +308,6 @@ def run_qc_combined(ad, qc_metrics, metrics_custom, threshold):
 
     # --- Degenerate guard + controlled collapse (combined) ---
     eps = 1e-9
-    CONS_FLOOR = 0.90   # optional: demand high cross-resolution agreement
-
     for max_mito in MITO_THRESHOLDS:
         frac_key = f"consensus_fraction_mito{max_mito}"
         pass_key = f"consensus_passed_qc_mito{max_mito}"
@@ -321,8 +318,7 @@ def run_qc_combined(ad, qc_metrics, metrics_custom, threshold):
             ad.uns.setdefault("consensus_threshold_degenerate", {})[max_mito] = "degenerate_all_zero"
 
         # Optional stricter rescue: require high raw fraction too
-        if CONS_FLOOR is not None:
-            ad.obs[pass_key] &= (ad.obs[frac_key] >= CONS_FLOOR)
+        ad.obs[pass_key] &= (ad.obs[frac_key] >= consensus_floor)
 
     # Optional interpretability: if pass at stricter mito, pass at looser too
     ad.obs["consensus_passed_qc_mito50"] |= ad.obs["consensus_passed_qc_mito20"]
@@ -490,7 +486,14 @@ def generate_qc_plots(ad, qc_metrics, qc_mode, metric_pairs, ctp_models, ctp_nam
         qc_consensus_ufig,
         qc_consensus_venn
     )
-def run_qc(ad, ctp_models=None, qc_mode="original", metrics_custom=None, threshold=0.5):
+def run_qc(
+    ad,
+    ctp_models=None,
+    qc_mode="original",
+    metrics_custom=None,
+    threshold=0.5,
+    consensus_floor=0.90,
+):
     """End-to-end QC run: compute metrics, run QC, and build figures."""
 
     qc_metrics = list(DEFAULT_QC_METRICS)
@@ -513,28 +516,33 @@ def run_qc(ad, ctp_models=None, qc_mode="original", metrics_custom=None, thresho
     logging.info("Calculating QC metrics")
     calculate_qc(ad, run_scrublet=("scrublet_score" in qc_metrics))
 
-    logging.info(f"Starting QC: mode={qc_mode}, cutoff={args.gmm_cutoff}")
+    logging.info(
+        f"Starting QC: mode={qc_mode}, cutoff={args.gmm_cutoff}, "
+        f"consensus_floor={consensus_floor}"
+    )
+    ad.uns["consensus_floor"] = float(consensus_floor)
     if qc_mode == "original":
         run_qc_mito_loop_original(ad, qc_metrics, metrics_custom, threshold)
     elif qc_mode == "multires":
-        run_qc_multi_res(ad, qc_metrics, metrics_custom, threshold)
+        run_qc_multi_res(ad, qc_metrics, metrics_custom, threshold, consensus_floor)
     elif qc_mode == "combined":
-        run_qc_combined(ad, qc_metrics, metrics_custom, threshold)
+        run_qc_combined(ad, qc_metrics, metrics_custom, threshold, consensus_floor)
     else:
         raise ValueError(f"Unknown QC mode: {qc_mode}")
 
     logging.info("Generating QC plots")
     return generate_qc_plots(ad, qc_metrics, qc_mode, metric_pairs, ctp_models, ctp_name)
 
 
-def process_sample(ad, ctp_models, qc_mode, metrics_custom, min_frac):
+def process_sample(ad, ctp_models, qc_mode, metrics_custom, min_frac, consensus_floor):
     """Wrapper around run_qc to keep the original structure."""
     return run_qc(
         ad=ad,
         ctp_models=ctp_models,
         qc_mode=qc_mode,
         metrics_custom=metrics_custom,
         threshold=min_frac,
+        consensus_floor=consensus_floor,
     )
 
 
@@ -571,6 +579,10 @@ def main(args):
     else:
         min_frac = float(0.5)
 
+    consensus_floor = float(args.consensus_floor)
+    if not 0 <= consensus_floor <= 1:
+        raise ValueError("--consensus_floor must be between 0 and 1.")
+
     if args.metrics_csv is not None:
         logging.info(f"Using custom metrics from {args.metrics_csv}")
         metrics_custom = pd.read_csv(Path(args.metrics_csv), index_col=0)
@@ -596,7 +608,9 @@ def main(args):
         qc_cluster_ufig,
         qc_consensus_ufig,
         qc_consensus_venn
-    ) = process_sample(ad, ctp_models, qc_mode, metrics_custom, min_frac)
+    ) = process_sample(
+        ad, ctp_models, qc_mode, metrics_custom, min_frac, consensus_floor
+    )
 
     logging.info("Saving QC plots")
     def _save(fig, filename):
@@ -705,6 +719,7 @@ def _save(fig, filename):
     parser.add_argument("--metrics_csv", type=nullable_string, nargs='?', help="CSV file of metric cutoffs")
     parser.add_argument("--celltypist_model", nargs='?', help="comma-separated <name>:<model.pkl> pairs or a predefined set key (e.g., 'gut')")
     parser.add_argument("--min_frac", default=None, help="min frac of pass_auto_filter for a cluster to be called good [default: 0.5]")
+    parser.add_argument("--consensus_floor", default=0.90, type=float, help="minimum consensus fraction required for multires/combined consensus rescue [default: 0.90]")
     parser.add_argument("--gath_obj", default=None, help="path to the AnnData object from gather_matrices step")
     parser.add_argument("--gmm_cutoff", default=None, help="whether to use closest or furthest point relative to GMM PDF")
     parser.add_argument("--sample_id", default=None, help="sample id")

main.nf

@@ -41,7 +41,7 @@ process run_qc {
   script:
   """
   export BASE_DIR=${projectDir}
-  python ${projectDir}/bin/qc.py --sample_id ${samp} --metrics_csv ${params.metrics_csv} --celltypist_model ${params.celltypist_model} --qc_mode ${params.qc_mode} --gath_obj ${gath_out} --gmm_cutoff ${params.gmm_cutoff}
+  python ${projectDir}/bin/qc.py --sample_id ${samp} --metrics_csv ${params.metrics_csv} --celltypist_model ${params.celltypist_model} --qc_mode ${params.qc_mode} --gath_obj ${gath_out} --gmm_cutoff ${params.gmm_cutoff} --consensus_floor ${params.consensus_floor}
   """
 }
 

nextflow.config

@@ -74,6 +74,7 @@ params {
   gather_mode = "cellbender" // Options: auto, cellbender, starsolo
   gmm_cutoff = "inner"  // Options: inner, outer
   qc_mode = "original"       // Options: original, multires, combined
+  consensus_floor = 0.90     // Minimum consensus fraction for multires/combined rescue
   celltypist_model = ""
   from_scautoqc = true
   n_top_genes = 5000