feat: add converter of reference from seurat to h5 and parquet#114
feat: add converter of reference from seurat to h5 and parquet#114mariiabilous wants to merge 2 commits into
Conversation
WalkthroughThis PR adds a Seurat reference ChangesReference RDS to H5/Parquet Conversion Pipeline
Estimated code review effort: 3 (Moderate) | ~25 minutes Sequence Diagram(s)sequenceDiagram
participant Snakefile
participant runConvertSeuratReferenceToH5
participant convertScript as convert_seurat_2_h5.R
participant CellTypeAnnotation as cell_type_annotation.smk
Snakefile->>Snakefile: _collect_reference_rds_paths()
Snakefile->>Snakefile: build REFERENCE_PATH_TO_RDS, REFERENCE_STEM
Snakefile->>runConvertSeuratReferenceToH5: resolve reference_stem input
runConvertSeuratReferenceToH5->>convertScript: invoke with assay/layer/genome params
convertScript->>convertScript: extract counts, write HDF5
convertScript->>convertScript: build metadata, write Parquet
convertScript-->>runConvertSeuratReferenceToH5: outputs (_counts.h5, _metadata.parquet)
CellTypeAnnotation->>CellTypeAnnotation: get_h5path4reference_based_annotation
CellTypeAnnotation->>CellTypeAnnotation: get_parquetpath4reference_based_annotation
CellTypeAnnotation->>Snakefile: lookup via REFERENCE_PATH_TO_RDS
Possibly related issues
Poem
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
workflow/Snakefile (1)
1406-1422: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick winDefault
rule allnow excludes several pipeline outputs.get_input2all()still defaultsreports,coexpression,joint_analysis,count_correction, andpost_count_correction_cell_type_annotationtoTrue, butrule alloverrides them toFalse, so a plainsnakemakerun no longer builds the same target set. Leave those arguments unset, or restore the previous defaults, unless this narrowing is intentional.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@workflow/Snakefile` around lines 1406 - 1422, The default target set in rule all is narrower than get_input2all()’s defaults, so a plain snakemake run no longer builds the full expected pipeline outputs. Update rule all to use get_input2all() with the same default behavior for reports, coexpression, joint_analysis, count_correction, and post_count_correction_cell_type_annotation, or remove those overrides if narrowing is not intended. Use the rule all block and the get_input2all() call as the places to align the target selection.
🧹 Nitpick comments (5)
workflow/rules/_reference_preprocessing/convert_seurat_2_h5.smk (1)
5-14: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winStale header comment doesn't match the actual
script:path.The comment references
"../../../scripts/utils/seurat_to_h5_parquet.R"(3 levels up,utils/dir), but the actual directive at Line 35 is"../../scripts/_reference_preprocessing/convert_seurat_2_h5.R"(2 levels up,_reference_preprocessing/dir). Left as-is, this will mislead anyone relocating the rule file later.📝 Suggested fix
-# NOTE: `script:` path is relative to *this* rule file's location — as -# written it assumes the same nesting depth as seurat.R -# ("../../../scripts/utils/seurat_to_h5_parquet.R"); adjust the leading -# "../" count if you place this .smk file elsewhere. +# NOTE: `script:` path is relative to *this* rule file's location +# ("../../scripts/_reference_preprocessing/convert_seurat_2_h5.R"); +# adjust the leading "../" count if you place this .smk file elsewhere.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@workflow/rules/_reference_preprocessing/convert_seurat_2_h5.smk` around lines 5 - 14, The header comment in convert_seurat_2_h5.smk is stale and no longer matches the rule’s actual script target. Update the NOTE block near the top to reference the real script path used by the rule’s script directive, and make sure the nesting-depth guidance matches the current location of the `convert_seurat_2_h5` rule so future moves of this file are not guided by the old `seurat_to_h5_parquet.R` path.workflow/scripts/_reference_preprocessing/convert_seurat_2_h5.R (2)
21-26: 🩺 Stability & Availability | 🔵 Trivial | ⚡ Quick winWrap conversion in
tryCatch/on.exitto guarantee sink cleanup and clean error surfacing.If
readRDS(),stopifnot(), orstop()insideseurat_to_h5_parquet()fail, execution jumps straight past Lines 185-187, leavingsink()diverted andlog_conunclosed for the remainder of the R session (harmless on process exit, but fragile if this script is ever sourced/re-run interactively, and means the log file may not be flushed consistently on failure).♻️ Suggested refactor
log_file <- snakemake@log[[1]] log_con <- file(log_file, open = "wt") sink(log_con, type = "output") sink(log_con, type = "message") +on.exit({ + sink(type = "message") + sink(type = "output") + close(log_con) +}, add = TRUE)Also applies to: 185-187
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@workflow/scripts/_reference_preprocessing/convert_seurat_2_h5.R` around lines 21 - 26, Wrap the script’s logging and conversion flow in a guarded cleanup so sinks are always restored even on failure. In the `log setup` block and around `seurat_to_h5_parquet()`, use `tryCatch` with an `on.exit` cleanup to close `log_con` and undo both `sink()` calls, and make sure errors from `readRDS()`, `stopifnot()`, or `stop()` are rethrown or logged clearly. Use the `log_file`, `log_con`, and `seurat_to_h5_parquet()` symbols to locate the setup and cleanup points.
113-129: 🗄️ Data Integrity & Integration | 🔵 Trivial | ⚡ Quick winNo collision guard when merging derived columns into metadata.
meta[shared, colnames(emb)] <- emb[shared, , drop = FALSE](Line 123) will silently overwrite an existing metadata column if a reduction's generated name (e.g.pca_1) already exists inmeta.data. Similarly,df[["cell_id"]] <- rownames(df)(Line 137) silently overwrites any pre-existingcell_idmetadata column. Both are edge cases, but failures would be silent data corruption in the exported Parquet rather than a visible error.🛡️ Suggested guard
colnames(emb) <- paste0(red_name, "_", seq_len(n_keep)) + clobbered <- intersect(colnames(emb), colnames(meta)) + if (length(clobbered)) { + warning(sprintf("Reduction '%s' overwrites existing metadata columns: %s", + red_name, paste(clobbered, collapse = ", "))) + } shared <- intersect(rownames(meta), rownames(emb))df <- as.data.frame(meta) + if ("cell_id" %in% colnames(df)) { + warning("Existing 'cell_id' metadata column will be overwritten with cell barcodes.") + } df[["cell_id"]] <- rownames(df)Also applies to: 136-138
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@workflow/scripts/_reference_preprocessing/convert_seurat_2_h5.R` around lines 113 - 129, The metadata merge in convert_seurat_2_h5.R can silently overwrite existing columns when adding reduction-derived fields or the cell_id field. Update the preprocessing flow around meta, the extra_reductions loop, and the cell_id assignment to detect name collisions before writing into df/meta.data, and either skip, rename, or stop with a clear error if any generated column already exists. Use the existing symbols Embeddings, extra_reductions, and cell_id handling to place the guard where the derived columns are added.workflow/rules/cell_type_annotation.smk (1)
120-138: 🎯 Functional Correctness | 🔵 Trivial | ⚡ Quick winWire the reference-based outputs through the new helpers The
_cell_type_annotation/_reference_based/*.smkrules still useget_path2reference4reference_based_annotation;get_h5path4reference_based_annotationandget_parquetpath4reference_based_annotationare currently unused, so the shared stem logic and.rdscheck never apply.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@workflow/rules/cell_type_annotation.smk` around lines 120 - 138, Wire the reference-based rules to the new path helpers so the shared stem logic is actually used. Update the `_cell_type_annotation/_reference_based/*.smk` rule inputs/outputs to call get_h5path4reference_based_annotation and get_parquetpath4reference_based_annotation instead of get_path2reference4reference_based_annotation, and keep _reference_path4reference_based_annotation as the single source of truth for deriving the reference key and validating the .rds-backed stem.workflow/rules/reference_preprocessing.smk (1)
6-12: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick winRemove the unused
get_path2reference_by_reference_pathhelper
workflow/rules/_reference_preprocessing/convert_seurat_2_h5.smkusesreference_stemdirectly, and there are no repo-wide call sites for this helper. Dropping it would remove dead code.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@workflow/rules/reference_preprocessing.smk` around lines 6 - 12, Remove the dead `get_path2reference_by_reference_path` helper from `reference_preprocessing.smk`; it is not referenced anywhere and `convert_seurat_2_h5.smk` already uses `reference_stem` directly. Delete the unused function and keep the surrounding Snakefile logic intact, making sure no remaining rules or imports depend on `REFERENCE_PATH_TO_RDS` through this helper.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@workflow/Snakefile`:
- Around line 126-159: The debug print in the reference path collection block
should be removed or replaced with a properly available logger, since `sys` is
not defined where `_collect_reference_rds_paths()` and `REFERENCE_STEM` are
built. Also harden `_collect_reference_rds_paths()` so it skips or normalizes
null `references` entries before calling `.items()`, using the existing config
access around `get_dict_value` and the `reference_files_by_condition` loop to
avoid crashes on explicit nulls.
---
Outside diff comments:
In `@workflow/Snakefile`:
- Around line 1406-1422: The default target set in rule all is narrower than
get_input2all()’s defaults, so a plain snakemake run no longer builds the full
expected pipeline outputs. Update rule all to use get_input2all() with the same
default behavior for reports, coexpression, joint_analysis, count_correction,
and post_count_correction_cell_type_annotation, or remove those overrides if
narrowing is not intended. Use the rule all block and the get_input2all() call
as the places to align the target selection.
---
Nitpick comments:
In `@workflow/rules/_reference_preprocessing/convert_seurat_2_h5.smk`:
- Around line 5-14: The header comment in convert_seurat_2_h5.smk is stale and
no longer matches the rule’s actual script target. Update the NOTE block near
the top to reference the real script path used by the rule’s script directive,
and make sure the nesting-depth guidance matches the current location of the
`convert_seurat_2_h5` rule so future moves of this file are not guided by the
old `seurat_to_h5_parquet.R` path.
In `@workflow/rules/cell_type_annotation.smk`:
- Around line 120-138: Wire the reference-based rules to the new path helpers so
the shared stem logic is actually used. Update the
`_cell_type_annotation/_reference_based/*.smk` rule inputs/outputs to call
get_h5path4reference_based_annotation and
get_parquetpath4reference_based_annotation instead of
get_path2reference4reference_based_annotation, and keep
_reference_path4reference_based_annotation as the single source of truth for
deriving the reference key and validating the .rds-backed stem.
In `@workflow/rules/reference_preprocessing.smk`:
- Around line 6-12: Remove the dead `get_path2reference_by_reference_path`
helper from `reference_preprocessing.smk`; it is not referenced anywhere and
`convert_seurat_2_h5.smk` already uses `reference_stem` directly. Delete the
unused function and keep the surrounding Snakefile logic intact, making sure no
remaining rules or imports depend on `REFERENCE_PATH_TO_RDS` through this
helper.
In `@workflow/scripts/_reference_preprocessing/convert_seurat_2_h5.R`:
- Around line 21-26: Wrap the script’s logging and conversion flow in a guarded
cleanup so sinks are always restored even on failure. In the `log setup` block
and around `seurat_to_h5_parquet()`, use `tryCatch` with an `on.exit` cleanup to
close `log_con` and undo both `sink()` calls, and make sure errors from
`readRDS()`, `stopifnot()`, or `stop()` are rethrown or logged clearly. Use the
`log_file`, `log_con`, and `seurat_to_h5_parquet()` symbols to locate the setup
and cleanup points.
- Around line 113-129: The metadata merge in convert_seurat_2_h5.R can silently
overwrite existing columns when adding reduction-derived fields or the cell_id
field. Update the preprocessing flow around meta, the extra_reductions loop, and
the cell_id assignment to detect name collisions before writing into
df/meta.data, and either skip, rename, or stop with a clear error if any
generated column already exists. Use the existing symbols Embeddings,
extra_reductions, and cell_id handling to place the guard where the derived
columns are added.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Organization UI
Review profile: CHILL
Plan: Pro
Run ID: 7623f6ee-bb80-4d46-9df4-1b959427e56d
📒 Files selected for processing (5)
workflow/Snakefileworkflow/rules/_reference_preprocessing/convert_seurat_2_h5.smkworkflow/rules/cell_type_annotation.smkworkflow/rules/reference_preprocessing.smkworkflow/scripts/_reference_preprocessing/convert_seurat_2_h5.R
| def _collect_reference_rds_paths() -> list[str]: | ||
| """ | ||
| All unique .rds reference paths configured across conditions/reference | ||
| names under experiments.<_cell_type_annotation>.<reference_files>. | ||
| """ | ||
| mapping = {} | ||
| reference_files_by_condition = get_dict_value( | ||
| config, | ||
| "experiments", | ||
| cc.EXPERIMENTS_CELL_TYPE_ANNOTATION_NAME, | ||
| cc.EXPERIMENTS_CELL_TYPE_ANNOTATION_REFERENCE_FILES_NAME, | ||
| replace_none={}, | ||
| ) | ||
| for condition_id, references in reference_files_by_condition.items(): | ||
| for reference_name, path in references.items(): | ||
| if path: | ||
| mapping[f"{condition_id}/{reference_name}"] = path | ||
| return mapping | ||
|
|
||
| REFERENCE_PATH_TO_RDS = _collect_reference_rds_paths() | ||
|
|
||
| # reference_stem is the .rds path with the extension stripped — the | ||
| # conversion rule writes <reference_stem>_counts.h5 and | ||
| # <reference_stem>_metadata.parquet right next to the source .rds. | ||
| REFERENCE_STEM = uniquify_elements_in_list( | ||
| [ | ||
| path[: -len(".rds")] | ||
| for path in REFERENCE_PATH_TO_RDS.values() | ||
| if path.endswith(".rds") | ||
| ] | ||
| ) | ||
|
|
||
| print(f"[DEBUG] REFERENCE_STEM: {REFERENCE_STEM}", file=sys.stderr) | ||
|
|
There was a problem hiding this comment.
🩺 Stability & Availability | 🟠 Major | ⚡ Quick win
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
rg -n '^import sys|^import re' workflow/scripts/utils/config_utils.py workflow/scripts/utils/raw_data_utils.py workflow/scripts/utils/run_time_utils.py 2>/dev/null
rg -n 'import sys' workflow/SnakefileRepository: bdsc-tds/xenium_analysis_pipeline
Length of output: 328
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "== Snakefile imports and sys usages =="
rg -n '^(import|from) ' workflow/Snakefile
rg -n '\bsys\b' workflow/Snakefile
echo "== Snakefile includes =="
rg -n '^include:' workflow/Snakefile
echo "== Surrounding lines for the reviewed block =="
nl -ba workflow/Snakefile | sed -n '110,170p'
echo "== Search for the referenced config helper behavior =="
rg -n 'def get_dict_value|replace_none|None' workflow/scripts/utils/config_utils.py
nl -ba workflow/scripts/utils/config_utils.py | sed -n '1,220p'Repository: bdsc-tds/xenium_analysis_pipeline
Length of output: 1183
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "== All sys imports in repo =="
rg -n '^import sys$|^from sys import|\\bsys\\b' workflow -g '!**/.snakemake/**' || true
echo "== get_dict_value definition =="
rg -n 'def get_dict_value|replace_none|replace_' workflow/scripts/utils/config_utils.py
sed -n '1,220p' workflow/scripts/utils/config_utils.py
echo "== included utility files (top of file) =="
for f in workflow/scripts/utils/config_utils.py workflow/scripts/utils/raw_data_utils.py workflow/scripts/utils/run_time_utils.py; do
echo "--- $f ---"
sed -n '1,40p' "$f"
done
echo "== Snakefile block around reference stem =="
sed -n '120,165p' workflow/SnakefileRepository: bdsc-tds/xenium_analysis_pipeline
Length of output: 13631
Import sys or remove this debug print, and guard null reference entries
print(..., file=sys.stderr) will raise NameError here because workflow/Snakefile and the included utility scripts don’t define sys. Remove or gate the debug line before merge. Also handle references before calling .items(), since an explicit null config entry would currently crash the loop.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@workflow/Snakefile` around lines 126 - 159, The debug print in the reference
path collection block should be removed or replaced with a properly available
logger, since `sys` is not defined where `_collect_reference_rds_paths()` and
`REFERENCE_STEM` are built. Also harden `_collect_reference_rds_paths()` so it
skips or normalizes null `references` entries before calling `.items()`, using
the existing config access around `get_dict_value` and the
`reference_files_by_condition` loop to avoid crashes on explicit nulls.
mariiabilous
changed the title
1 participant
Add Seurat → HDF5/Parquet reference conversion
Summary
Adds a rule that converts scRNA-seq reference
.rdsobjects into HDF5 counts + Parquet metadata, and wires those files as inputs to the reference-based annotation rule. Lets Python tooling (e.g. py-RCTD) read the reference without going through R/rhdf5.Changes
Snakefile:REFERENCE_STEMwildcard +REFERENCE_PATH_TO_RDSlookup, built fromexperiments.<_cell_type_annotation>.<reference_files>.rules/_reference_preprocessing/convert_seurat_2_h5.smk(new):runConvertSeuratReferenceToH5—{reference_stem}.rds→_counts.h5+_metadata.parquet.scripts/_reference_preprocessing/convert_seurat_2_h5.R(new): does the actual conversion viawrite10xCounts()+arrow::write_parquet().rules/reference_preprocessing.smk: helpers to resolve an annotation rule's wildcards to the corresponding.h5/.parquetpaths.Notes
mem_mbmultiplier after an OOM onexternal_CRC.rds.rules/seurat_to_h5.smk,scripts/utils/seurat_to_h5.R) are unused, safe to delete.Testing
rules/_cell_type_annotation/_reference_based/rctd.smk./run.sh -nconfirms DAG pulls in conversion automatically when outputs are missing.external_CRC.rds; succeeded after memory fix.