diff --git a/CHANGELOG.md b/CHANGELOG.md index 4f027ad..2a11c58 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,8 @@ * `workflows/single_cell/process_integrate_annotate`: Perform integration and annotation with multiple methods in parallel rather than sequentially using the `workflows/single_cell/parallel_annotation` and `workflows/single_cell/parallel_integration` workflows (PR #21). +* `workflows/single_cell/parallel_subtyping`: Add a workflow that subtypes each major cell type independently via reference-based label projection, splitting the reference by matching major cell type (`--reference_obs_major_cell_type`) so each type is subtyped against its own reference cells. Query cell types absent from the reference raise an error by default, or are passed through unannotated when `--allow_missing_reference_cell_type` is set. The subtype labels are combined into a single output h5mu (PR #23). + ## MINOR CHANGES * Migration of test resources to the package-specific `s3://openpipelines-bio/openpipeline_composed/resources_test` bucket (PR #24): diff --git a/src/single_cell/parallel_subtyping/config.vsh.yaml b/src/single_cell/parallel_subtyping/config.vsh.yaml new file mode 100644 index 0000000..1d3cd99 --- /dev/null +++ b/src/single_cell/parallel_subtyping/config.vsh.yaml @@ -0,0 +1,527 @@ +name: parallel_subtyping +namespace: single_cell +scope: public +description: | + Cell subtyping of h5mu file annotated with major types. + + The input is expected to already carry a major cell-type annotation in + .obs (--obs_major_cell_type) on the modality to be subtyped (--modality). + Each major cell type is then subtyped independently, so that the annotation + is tailored to the cells within that group rather than the dataset as a + whole. The reference is likewise split by its matching major cell-type + annotation (--reference_obs_major_cell_type) so that each major cell type is + subtyped only against the reference cells of that same type. The resulting + subtype labels are combined back into a single output h5mu, preserving all + original modalities and cells. + + The annotation methods and per-method options are those of the + parallel_annotation workflow and are applied identically to every major + cell type. + + Available annotation methods: + - CellTypist + - ScanVI with scArches + - SingleR + - SCVI integration followed by KNN label transfer + - Harmony integration followed by KNN label transfer + +authors: + - __merge__: /src/authors/dorien_roosen.yaml + roles: [author, maintainer] + +argument_groups: + - name: Inputs + arguments: + - name: "--id" + type: string + description: ID of the sample. + required: true + example: foo + - name: "--input" + type: file + description: | + Input query h5mu file with a major cell-type annotation. Must be + preprocessed in the same way as the reference. + required: true + example: input.h5mu + - name: "--modality" + type: string + description: | + Which modality to subtype. Must match the modality in --reference. + Other modalities present in the input are passed through unaltered + and merged back into the output. + default: "rna" + required: false + - name: "--obs_major_cell_type" + type: string + description: | + The .obs column of --modality holding the major cell-type labels to + split on. The target modality is split by this column and each major + cell type is subtyped independently. + required: true + example: cell_type + - name: "--input_layer" + type: string + description: | + The layer of the input query containing raw counts. Required by + celltypist (which log-normalizes internally), scanvi_scarches and + scvi_knn (which model raw counts directly). If not provided, .X + will be used. + required: false + - name: "--input_layer_lognormalized" + type: string + description: | + The layer of the input query containing log-normalized counts. + Required by harmony_knn and scvi_knn. + default: log_normalized + required: false + - name: "--input_obs_batch_label" + type: string + description: | + .obs column in the input query containing batch labels. Required by + harmony_knn, scanvi_scarches and scvi_knn. + default: sample_id + required: false + - name: "--input_var_gene_names" + type: string + description: | + .var column containing gene names. By default the .var index is used. + required: false + - name: "--input_reference_gene_overlap" + type: integer + default: 100 + min: 1 + description: Minimum number of genes required to overlap between query and reference. + - name: "--sanitize_ensembl_ids" + type: boolean + default: true + description: | + Whether to sanitize ensembl gene ids by removing version numbers. Read by + celltypist, scanvi_scarches, scvi_knn and singler. + - name: "--input_obs_categorical_covariates" + type: string + multiple: true + description: | + .obs column(s) in the input query with categorical covariates for + scanvi_scarches to correct for in addition to the batch label. Treated as + nuisance factors; their order must match --reference_obs_categorical_covariates. + None by default. + required: false + - name: "--input_obs_numerical_covariates" + type: string + multiple: true + description: | + .obs column(s) in the input query with continuous (numerical) covariates for + scanvi_scarches to correct for in addition to the batch label. Treated as + nuisance factors; their order must match --reference_obs_numerical_covariates. + None by default. + required: false + + - name: Methods + arguments: + - name: "--annotation_methods" + type: string + description: Annotation methods to run in parallel on each major cell type. + choices: [celltypist, harmony_knn, scanvi_scarches, scvi_knn, singler] + multiple: true + required: true + example: [celltypist, harmony_knn] + + - name: Reference + arguments: + - name: "--reference" + type: file + description: | + Reference h5mu file with labeled observations. Used by all annotation + methods (celltypist, harmony_knn, scanvi_scarches and scvi_knn). + required: true + example: reference.h5mu + - name: "--reference_layer" + type: string + description: Layer of the reference containing raw counts. If not provided, .X will be used. + required: false + - name: "--reference_layer_lognormalized" + type: string + description: Layer of the reference containing log-normalized counts. + required: false + default: log_normalized + - name: "--reference_obs_target" + type: string + description: .obs column in the reference containing the target (sub)cell-type labels. + required: false + example: cell_type + - name: "--reference_obs_major_cell_type" + type: string + description: | + The .obs column in --reference holding the major cell-type labels, + matching the label set in --obs_major_cell_type. The reference is split + by this column and each query major cell type is subtyped against the + matching subset of the reference. + required: true + example: major_cell_type + - name: "--allow_missing_reference_cell_type" + type: boolean_true + description: | + Controls what happens when a query major cell type has no matching + subset in the reference. When false (default), the workflow errors. + When true, such cell types are passed through unannotated and merged + back into the output alongside the subtyped cell types. + - name: "--reference_obs_batch_label" + type: string + description: .obs column in the reference containing batch labels. + required: false + default: sample_id + - name: "--reference_var_gene_names" + type: string + description: .var column in the reference containing gene names. Defaults to the .var index. + required: false + - name: "--reference_var_input" + type: string + description: | + .var column in the reference flagging highly-variable genes (boolean). + Used by celltypist and scanvi_scarches to subset features if provided. + required: false + default: filter_with_hvg + - name: "--reference_obs_categorical_covariates" + type: string + multiple: true + description: | + .obs column(s) in the reference with categorical covariates for + scanvi_scarches to correct for in addition to the batch label. Treated as + nuisance factors; their order must match --input_obs_categorical_covariates. + None by default. + required: false + - name: "--reference_obs_numerical_covariates" + type: string + multiple: true + description: | + .obs column(s) in the reference with continuous (numerical) covariates for + scanvi_scarches to correct for in addition to the batch label. Treated as + nuisance factors; their order must match --input_obs_numerical_covariates. + None by default. + required: false + - name: "--reference_obs_label_unlabeled_category" + type: string + default: "Unknown" + description: Value in --reference_obs_target that indicates unlabeled observations. + + - name: Clustering options + description: Leiden clustering options shared by harmony_knn, scanvi_scarches, scvi_knn. + arguments: + - name: "--leiden_resolution" + type: double + description: Leiden clustering resolution(s). Higher values lead to more clusters. + default: [1] + multiple: true + + - name: Neighbor classifier options + description: KNN options shared by harmony_knn, scanvi_scarches, scvi_knn. + arguments: + - name: "--knn_weights" + type: string + choices: [uniform, distance] + default: "uniform" + description: Weight function for KNN classifier. + - name: "--knn_n_neighbors" + type: integer + min: 5 + default: 15 + description: Number of neighbors for the KNN classifier. + + - name: scVI / scANVI / scArches training options + description: | + Training arguments shared by scanvi_scarches and scvi_knn. + arguments: + - name: "--scvi_early_stopping" + type: boolean + description: Whether to perform early stopping with respect to the validation set. + required: false + - name: "--scvi_early_stopping_monitor" + type: string + choices: [elbo_validation, reconstruction_loss_validation, kl_local_validation] + default: "elbo_validation" + description: Metric logged during validation set epoch. + - name: "--scvi_early_stopping_patience" + type: integer + min: 1 + default: 45 + description: Number of validation epochs with no improvement after which training stops. + - name: "--scvi_early_stopping_min_delta" + type: double + min: 0 + default: 0.0 + description: Minimum change in the monitored metric to qualify as an improvement. + - name: "--scvi_max_epochs" + type: integer + description: Number of training passes through the dataset. + required: false + - name: "--scvi_reduce_lr_on_plateau" + type: boolean + default: true + description: Reduce learning rate when validation metric plateaus. + - name: "--scvi_lr_factor" + type: double + min: 0 + default: 0.6 + description: Factor by which to reduce learning rate. + - name: "--scvi_lr_patience" + type: double + min: 0 + default: 30 + description: Number of epochs with no improvement after which lr is reduced. + + - name: CellTypist options + arguments: + - name: "--celltypist_majority_voting" + type: boolean + default: false + description: Refine predicted labels by majority voting after over-clustering. + - name: "--celltypist_feature_selection" + type: boolean + default: false + description: Whether to perform feature selection during training. + - name: "--celltypist_C" + type: double + default: 1.0 + description: Inverse of regularization strength in logistic regression. + - name: "--celltypist_max_iter" + type: integer + default: 1000 + description: Maximum iterations before reaching the cost-function minimum. + - name: "--celltypist_use_SGD" + type: boolean + default: false + description: Whether to use stochastic gradient descent. + - name: "--celltypist_min_prop" + type: double + default: 0 + description: | + Minimum proportion of cells in a sub-cluster supporting a label for + majority voting; only relevant if --celltypist_majority_voting is set. + + - name: harmony_knn options + arguments: + - name: "--harmony_knn_n_hvg" + type: integer + default: 2000 + description: Number of highly-variable genes to keep for harmony_knn. + - name: "--harmony_knn_pca_num_components" + type: integer + description: Number of principal components to compute for harmony_knn. + required: false + - name: "--harmony_knn_theta" + type: double + description: Diversity clustering penalty for harmony. + default: [2] + multiple: true + + - name: scvi_knn options + arguments: + - name: "--scvi_knn_n_hvg" + type: integer + default: 2000 + description: Number of highly-variable genes to keep for scvi_knn. + + - name: SingleR options + arguments: + - name: "--singler_input_obs_clusters" + type: string + required: false + description: | + .obs column with cluster identities. If provided, SingleR annotates the + aggregated cluster profiles; otherwise it annotates per observation. + - name: "--singler_de_method" + type: string + choices: [classic, t, wilcox] + default: "classic" + description: Method to detect differentially expressed genes between pairs of labels. + - name: "--singler_de_n_genes" + type: integer + min: 1 + required: false + description: | + Number of differentially expressed genes per label pair to use as markers. + Defaults to a SingleR heuristic based on the number of labels. + - name: "--singler_quantile" + type: double + min: 0 + max: 1 + default: 0.8 + description: Quantile of the correlation distribution used to score each label. + - name: "--singler_fine_tune" + type: boolean + default: true + description: Whether to fine-tune label assignment after the initial classification. + - name: "--singler_fine_tuning_threshold" + type: double + default: 0.05 + description: Maximum difference from the maximum correlation to use in fine-tuning. + - name: "--singler_prune" + type: boolean + default: true + description: | + Whether to prune low-quality labels (replaced with NA in the pruned + predictions output). + + - name: Consensus vote options + description: | + A final consensus_vote step (run inside parallel_annotation on each major + cell type) combines the per-method predictions into a single consensus + label via a (probability-weighted) majority vote. + + As with the per-method slots, the consensus .obs columns default to + consensus_pred_ and + consensus_score_ when not set, and the suffix is + omitted when no --reference_obs_target is set. + arguments: + - name: "--run_consensus" + type: boolean + default: true + description: | + Run a final consensus vote combining the per-method predictions. Only + takes effect when more than one --annotation_methods is selected. + - name: "--consensus_obs_predictions" + type: string + required: false + description: .obs slot for the consensus predicted cell type. Defaults to consensus_pred_. + - name: "--consensus_obs_score" + type: string + required: false + description: | + .obs slot for the consensus score, defined as the fraction of total + weight assigned to the winning cell type. Defaults to + consensus_score_. + - name: "--consensus_use_probabilities" + type: boolean + default: true + description: | + Scale each method's vote by its per-cell prediction probability. When + false, every method votes with equal weight regardless of confidence. + - name: "--consensus_tie_label" + type: string + description: | + Label to assign when two or more cell types receive equal votes. If not + provided, tied cells are assigned a missing value. + required: false + example: "Unknown" + + - name: Output slot naming + description: | + Names of the .obs and .obsm slots each method writes to the subtyped output. + arguments: + - name: "--celltypist_obs_predictions" + type: string + required: false + description: .obs column for CellTypist predicted labels. Defaults to celltypist_pred_. + - name: "--celltypist_obs_probability" + type: string + required: false + description: .obs column for CellTypist probabilities. Defaults to celltypist_probability_. + + - name: "--harmony_knn_obs_predictions" + type: string + required: false + description: .obs column for harmony_knn predicted labels. Defaults to harmony_knn_pred_. + - name: "--harmony_knn_obs_probability" + type: string + required: false + description: .obs column for harmony_knn probabilities. Defaults to harmony_knn_probability_. + - name: "--harmony_knn_obsm_integrated" + type: string + default: "X_integrated_harmony" + + - name: "--scanvi_scarches_obs_predictions" + type: string + required: false + description: .obs column for scanvi_scarches predicted labels. Defaults to scanvi_pred_. + - name: "--scanvi_scarches_obs_probability" + type: string + required: false + description: .obs column for scanvi_scarches probabilities. Defaults to scanvi_probability_. + - name: "--scanvi_scarches_obsm_integrated" + type: string + default: "X_integrated_scanvi" + + - name: "--scvi_knn_obs_predictions" + type: string + required: false + description: .obs column for scvi_knn predicted labels. Defaults to scvi_knn_pred_. + - name: "--scvi_knn_obs_probability" + type: string + required: false + description: .obs column for scvi_knn probabilities. Defaults to scvi_knn_probability_. + - name: "--scvi_knn_obsm_integrated" + type: string + default: "X_integrated_scvi" + + - name: "--singler_obs_predictions" + type: string + required: false + description: .obs column for SingleR predicted labels. Defaults to singler_pred_. + - name: "--singler_obs_probability" + type: string + required: false + description: .obs column for SingleR probabilities. Defaults to singler_probability_. + - name: "--singler_obs_delta_next" + type: string + required: false + description: .obs column for the SingleR delta-next score. Defaults to singler_delta_next_. + - name: "--singler_obs_pruned_predictions" + type: string + required: false + description: .obs column for SingleR pruned predicted labels. Defaults to singler_pruned_labels_. + - name: "--singler_obsm_scores" + type: string + required: false + description: .obsm key for the SingleR per-label score matrix. Defaults to singler_scores_. + + - name: Output + arguments: + - name: "--output" + alternatives: [-o] + type: file + description: | + Output h5mu file — the input query with each selected method's + subtype predictions and embeddings merged into the target modality, + and the other modalities merged back in unaltered. + direction: output + required: true + example: output.h5mu + - name: "--output_compression" + type: string + choices: [gzip, lzf] + required: false + example: gzip + description: Compression format for the output h5mu file. + +dependencies: + - name: single_cell/parallel_annotation + - name: dataflow/split_modalities + repository: openpipeline + - name: dataflow/split_h5mu + repository: openpipeline + - name: dataflow/concatenate_h5mu + repository: openpipeline + - name: dataflow/merge + repository: openpipeline + # Used only by test.nf to inject a major cell-type annotation into the query. + - name: metadata/add_id + repository: openpipeline + +resources: + - type: nextflow_script + path: main.nf + entrypoint: run_wf + +test_resources: + - type: nextflow_script + path: test.nf + entrypoint: test_wf + - type: nextflow_script + path: test.nf + entrypoint: test_wf_2 + - path: /resources_test/pbmc_1k_protein_v3/pbmc_1k_protein_v3_mms.h5mu + - path: /resources_test/annotation_test_data/TS_Blood_filtered.h5mu + +runners: + - type: nextflow diff --git a/src/single_cell/parallel_subtyping/integration_test.sh b/src/single_cell/parallel_subtyping/integration_test.sh new file mode 100644 index 0000000..e475b52 --- /dev/null +++ b/src/single_cell/parallel_subtyping/integration_test.sh @@ -0,0 +1,29 @@ +#!/bin/bash + +set -eo pipefail + +# get the root of the directory +REPO_ROOT=$(git rev-parse --show-toplevel) + +# ensure that the command below is run from the root of the repository +cd "$REPO_ROOT" + +nextflow \ + run . \ + -main-script src/single_cell/parallel_subtyping/test.nf \ + -entry test_wf \ + -resume \ + -profile docker \ + -c src/configs/labels_ci.config \ + -c src/configs/integration_tests.config \ + --publish_dir test + +nextflow \ + run . \ + -main-script src/single_cell/parallel_subtyping/test.nf \ + -entry test_wf_2 \ + -resume \ + -profile docker \ + -c src/configs/labels_ci.config \ + -c src/configs/integration_tests.config \ + --publish_dir test diff --git a/src/single_cell/parallel_subtyping/main.nf b/src/single_cell/parallel_subtyping/main.nf new file mode 100644 index 0000000..2fdbac9 --- /dev/null +++ b/src/single_cell/parallel_subtyping/main.nf @@ -0,0 +1,315 @@ +workflow run_wf { + take: + input_ch + + main: + // === Split modalities === + // Preserve the user-requested target modality and the original sample id, + // then split the input into one h5mu file per modality. + modalities_ch = input_ch + | map { id, state -> + def new_state = state + [ + "target_modality": state.modality, + "original_id": id, + "_meta": ["join_id": id] + ] + [id, new_state] + } + | split_modalities.run( + fromState: [ + "input": "input", + "output_compression": "output_compression" + ], + toState: [ + "split_modalities_output": "output", + "split_modalities_types": "output_types" + ] + ) + // Fan the modality directory out into one event per modality. + | flatMap { id, state -> + def outputDir = state.split_modalities_output + def csv = state.split_modalities_types + .splitCsv(strip: true, sep: ",") + .findAll { !it[0].startsWith("#") } + def header = csv.head() + def rows = csv.tail().collect { row -> [header, row].transpose().collectEntries() } + rows.collect { dat -> + def new_id = "${state.original_id}_${dat.name}" + def new_state = state + [ + "input": outputDir.resolve(dat.filename), + "modality": dat.name + ] + [new_id, new_state] + } + } + | map { id, state -> + def keysToRemove = ["split_modalities_output", "split_modalities_types"] + [id, state.findAll { it.key !in keysToRemove }] + } + + // Only the target modality is subtyped; other modalities pass through + // unaltered and are merged back at the end. + target_ch = modalities_ch + | filter { id, state -> state.modality == state.target_modality } + + passthrough_ch = modalities_ch + | filter { id, state -> state.modality != state.target_modality } + + // === Subtype the target modality === + // Split the target modality by major cell type. + query_split_ch = target_ch + | split_h5mu.run( + fromState: [ + "input": "input", + "modality": "modality", + "obs_feature": "obs_major_cell_type", + "output_compression": "output_compression" + ], + // Cell-type names can collide after sanitizing; keep filenames unique. + args: ["ensure_unique_filenames": true], + toState: [ + "split_h5mu_output": "output", + "split_h5mu_files": "output_files" + ] + ) + + // Split the reference by its matching major cell-type annotation, + // so each major cell type can be subtyped against the reference cells of that same type. + reference_split_ch = target_ch + | split_h5mu.run( + key: "split_reference_h5mu", + fromState: [ + "input": "reference", + "modality": "modality", + "obs_feature": "reference_obs_major_cell_type", + "output_compression": "output_compression" + ], + args: ["ensure_unique_filenames": true], + toState: [ + "split_reference_output": "output", + "split_reference_files": "output_files" + ] + ) + | map { id, state -> [id, state.subMap("split_reference_output", "split_reference_files")] } + + // Recombine the two splits per (sample, modality) id, then fan out into one + // event per major cell type. + // Keep the modality-level id as the group key so the subsets + // can be concatenated back together after annotation. + per_cell_type_ch = query_split_ch.join(reference_split_ch) + | map { id, query_state, reference_state -> [id, query_state + reference_state] } + | flatMap { id, state -> + def outputDir = state.split_h5mu_output + def csv = state.split_h5mu_files + .splitCsv(strip: true, sep: ",") + .findAll { !it[0].startsWith("#") } + def header = csv.head() + def rows = csv.tail().collect { row -> [header, row].transpose().collectEntries() } + + // Map each major cell-type name to its reference subset. + def refByCellType = [:] + def refDir = state.split_reference_output + def refCsv = state.split_reference_files + .splitCsv(strip: true, sep: ",") + .findAll { !it[0].startsWith("#") } + def refHeader = refCsv.head() + refCsv.tail() + .collect { row -> [refHeader, row].transpose().collectEntries() } + .each { refByCellType[it.name] = refDir.resolve(it.filename) } + + rows.collect { dat -> + def new_id = "${id}_${dat.name}" + + // Match the query major cell type to its reference subset. A query + // cell type absent from the reference either errors or, with + // --allow_missing_reference_cell_type, is passed through unannotated. + def reference = state.reference + def annotate = true + if (refByCellType.containsKey(dat.name)) { + reference = refByCellType[dat.name] + } else if (state.allow_missing_reference_cell_type) { + annotate = false + } else { + throw new RuntimeException( + "Major cell type '${dat.name}' has no matching subset in the " + + "reference (--reference_obs_major_cell_type). Set " + + "--allow_missing_reference_cell_type to true to pass such cell " + + "types through without subtyping instead." + ) + } + + def new_state = state + [ + "id": new_id, + "input": outputDir.resolve(dat.filename), + "subtype_group_id": id, + "major_cell_type": dat.name, + "reference": reference, + "annotate": annotate + ] + [new_id, new_state] + } + } + | map { id, state -> + def keysToRemove = ["split_h5mu_output", "split_h5mu_files", "split_reference_output", "split_reference_files"] + [id, state.findAll { it.key !in keysToRemove }] + } + + // Only cell types with a matching subset in the reference are subtyped. + annotated_ch = per_cell_type_ch + | filter { id, state -> state.annotate } + // Derive the predicted-label / probability .obs column names from the + // reference target so subtype predictions are tagged with the label set + // they were transferred from. Explicitly provided names take precedence; + // the suffix is omitted when no reference target is set. + | map { id, state -> + def suffix = state.reference_obs_target ? "_${state.reference_obs_target}" : "" + def slot_defaults = [ + "celltypist_obs_predictions": "celltypist_pred${suffix}", + "celltypist_obs_probability": "celltypist_probability${suffix}", + "harmony_knn_obs_predictions": "harmony_knn_pred${suffix}", + "harmony_knn_obs_probability": "harmony_knn_probability${suffix}", + "scanvi_scarches_obs_predictions": "scanvi_pred${suffix}", + "scanvi_scarches_obs_probability": "scanvi_probability${suffix}", + "scvi_knn_obs_predictions": "scvi_knn_pred${suffix}", + "scvi_knn_obs_probability": "scvi_knn_probability${suffix}", + "singler_obs_predictions": "singler_pred${suffix}", + "singler_obs_probability": "singler_probability${suffix}", + "singler_obs_delta_next": "singler_delta_next${suffix}", + "singler_obs_pruned_predictions": "singler_pruned_labels${suffix}", + "singler_obsm_scores": "singler_scores${suffix}", + "consensus_obs_predictions": "consensus_pred${suffix}", + "consensus_obs_score": "consensus_score${suffix}" + ] + def filled = slot_defaults.collectEntries { k, v -> [(k): state[k] ?: v] } + [id, state + filled] + } + | parallel_annotation.run( + fromState: [ + "id": "id", + "input": "input", + "modality": "modality", + "input_layer": "input_layer", + "input_layer_lognormalized": "input_layer_lognormalized", + "input_obs_batch_label": "input_obs_batch_label", + "input_var_gene_names": "input_var_gene_names", + "input_reference_gene_overlap": "input_reference_gene_overlap", + "sanitize_ensembl_ids": "sanitize_ensembl_ids", + "input_obs_categorical_covariates": "input_obs_categorical_covariates", + "input_obs_numerical_covariates": "input_obs_numerical_covariates", + "annotation_methods": "annotation_methods", + "reference": "reference", + "reference_layer": "reference_layer", + "reference_layer_lognormalized": "reference_layer_lognormalized", + "reference_obs_target": "reference_obs_target", + "reference_obs_batch_label": "reference_obs_batch_label", + "reference_var_gene_names": "reference_var_gene_names", + "reference_var_input": "reference_var_input", + "reference_obs_categorical_covariates": "reference_obs_categorical_covariates", + "reference_obs_numerical_covariates": "reference_obs_numerical_covariates", + "reference_obs_label_unlabeled_category": "reference_obs_label_unlabeled_category", + "leiden_resolution": "leiden_resolution", + "knn_weights": "knn_weights", + "knn_n_neighbors": "knn_n_neighbors", + "scvi_early_stopping": "scvi_early_stopping", + "scvi_early_stopping_monitor": "scvi_early_stopping_monitor", + "scvi_early_stopping_patience": "scvi_early_stopping_patience", + "scvi_early_stopping_min_delta": "scvi_early_stopping_min_delta", + "scvi_max_epochs": "scvi_max_epochs", + "scvi_reduce_lr_on_plateau": "scvi_reduce_lr_on_plateau", + "scvi_lr_factor": "scvi_lr_factor", + "scvi_lr_patience": "scvi_lr_patience", + "celltypist_majority_voting": "celltypist_majority_voting", + "celltypist_feature_selection": "celltypist_feature_selection", + "celltypist_C": "celltypist_C", + "celltypist_max_iter": "celltypist_max_iter", + "celltypist_use_SGD": "celltypist_use_SGD", + "celltypist_min_prop": "celltypist_min_prop", + "harmony_knn_n_hvg": "harmony_knn_n_hvg", + "harmony_knn_pca_num_components": "harmony_knn_pca_num_components", + "harmony_knn_theta": "harmony_knn_theta", + "scvi_knn_n_hvg": "scvi_knn_n_hvg", + "singler_input_obs_clusters": "singler_input_obs_clusters", + "singler_de_method": "singler_de_method", + "singler_de_n_genes": "singler_de_n_genes", + "singler_quantile": "singler_quantile", + "singler_fine_tune": "singler_fine_tune", + "singler_fine_tuning_threshold": "singler_fine_tuning_threshold", + "singler_prune": "singler_prune", + "celltypist_obs_predictions": "celltypist_obs_predictions", + "celltypist_obs_probability": "celltypist_obs_probability", + "harmony_knn_obs_predictions": "harmony_knn_obs_predictions", + "harmony_knn_obs_probability": "harmony_knn_obs_probability", + "harmony_knn_obsm_integrated": "harmony_knn_obsm_integrated", + "scanvi_scarches_obs_predictions": "scanvi_scarches_obs_predictions", + "scanvi_scarches_obs_probability": "scanvi_scarches_obs_probability", + "scanvi_scarches_obsm_integrated": "scanvi_scarches_obsm_integrated", + "scvi_knn_obs_predictions": "scvi_knn_obs_predictions", + "scvi_knn_obs_probability": "scvi_knn_obs_probability", + "scvi_knn_obsm_integrated": "scvi_knn_obsm_integrated", + "singler_obs_predictions": "singler_obs_predictions", + "singler_obs_probability": "singler_obs_probability", + "singler_obs_delta_next": "singler_obs_delta_next", + "singler_obs_pruned_predictions": "singler_obs_pruned_predictions", + "singler_obsm_scores": "singler_obsm_scores", + "run_consensus": "run_consensus", + "consensus_obs_predictions": "consensus_obs_predictions", + "consensus_obs_score": "consensus_obs_score", + "consensus_use_probabilities": "consensus_use_probabilities", + "consensus_tie_label": "consensus_tie_label", + "output_compression": "output_compression" + ], + toState: ["input": "output"] + ) + + // Cell types absent from the split reference are passed through unannotated + // when --allow_missing_reference_cell_type is set. + passthrough_cell_types_ch = per_cell_type_ch + | filter { id, state -> !state.annotate } + + // Concatenate the per-major-type subsets back into one unimodal file. + subtyped_ch = annotated_ch.mix(passthrough_cell_types_ch) + | map { id, state -> [state.subtype_group_id, id, state] } + | groupTuple(by: 0, sort: "hash") + | map { group_id, cell_type_ids, states -> + // Every non per-cell-type key is identical across the subsets (they all + // originate from the same modality file), so carry the first state. + def merged_state = states[0] + [ + "input": states.collect { it.input }, + "input_id": cell_type_ids + ] + [group_id, merged_state] + } + | concatenate_h5mu.run( + fromState: [ + "input": "input", + "input_id": "input_id", + "modality": "modality", + "output_compression": "output_compression" + ], + toState: ["input": "output"] + ) + + // === Merge modalities back together === + // Group the subtyped target modality with the untouched passthrough + // modalities per original sample, then merge them into one output h5mu. + output_ch = subtyped_ch.mix(passthrough_ch) + | map { id, state -> [state.original_id, state] } + | groupTuple(by: 0, sort: "hash") + | map { orig_id, states -> + def merged_state = states[0] + [ + "input": states.collect { it.input } + ] + [orig_id, merged_state] + } + | merge.run( + fromState: [ + "input": "input", + "output_compression": "output_compression" + ], + toState: ["output": "output"] + ) + | setState(["output", "_meta"]) + + emit: + output_ch +} diff --git a/src/single_cell/parallel_subtyping/nextflow.config b/src/single_cell/parallel_subtyping/nextflow.config new file mode 100644 index 0000000..21e68ab --- /dev/null +++ b/src/single_cell/parallel_subtyping/nextflow.config @@ -0,0 +1,10 @@ +manifest { + nextflowVersion = '!>=20.12.1-edge' +} + +params { + rootDir = java.nio.file.Paths.get("$projectDir/../../../").toAbsolutePath().normalize().toString() +} + +// include common settings +includeConfig("${params.rootDir}/src/configs/labels.config") diff --git a/src/single_cell/parallel_subtyping/test.nf b/src/single_cell/parallel_subtyping/test.nf new file mode 100644 index 0000000..3708dd7 --- /dev/null +++ b/src/single_cell/parallel_subtyping/test.nf @@ -0,0 +1,163 @@ +nextflow.enable.dsl=2 + +include { parallel_subtyping } from params.rootDir + "/target/nextflow/single_cell/parallel_subtyping/main.nf" +// add_id is an openpipeline dependency; it builds into the version-pinned +// dependency cache rather than target/nextflow (tag set in _viash.yaml). +include { add_id } from params.rootDir + "/target/dependencies/vsh/vsh/openpipeline/v4.2.0/nextflow/metadata/add_id/main.nf" +include { assert_h5mu_slots } from params.rootDir + "/target/_test/nextflow/test_workflows/assert_h5mu_slots/main.nf" +params.resources_test = "s3://openpipelines-bio/openpipeline_composed/resources_test/" + +workflow test_wf { + // Subtype the rna modality of a multimodal (rna + prot) sample. The query + // carries no cell-type annotation, so add_id first stamps every cell with a + // "immune" major cell type matching the reference's `compartment` column. + // parallel_subtyping then splits both query and reference on that column and + // annotates the single "immune" group with a CellTypist model trained on the + // matching reference subset, before merging everything back. + resources_test = file(params.resources_test) + + output_ch = Channel.fromList( + [ + [ + id: "subtyping_test", + input: resources_test.resolve("pbmc_1k_protein_v3/pbmc_1k_protein_v3_mms.h5mu"), + modality: "rna", + obs_major_cell_type: "compartment", + reference: resources_test.resolve("annotation_test_data/TS_Blood_filtered.h5mu"), + reference_var_gene_names: "ensemblid", + reference_obs_target: "cell_type", + reference_obs_major_cell_type: "compartment", + reference_var_input: "highly_variable", + annotation_methods: "celltypist" + ] + ]) + | view { "State at start: $it" } + | map { state -> [state.id, state] } + // Stamp the query cells with a major cell type matching the reference. + | add_id.run( + fromState: ["input": "input"], + args: ["input_id": "immune", "obs_output": "compartment", "make_observation_keys_unique": false], + toState: ["input": "output"] + ) + | parallel_subtyping + | view { "After parallel_subtyping: $it" } + | view { output -> + assert output.size() == 2 : "Outputs should contain two elements; [id, state]" + + def id = output[0] + assert id == "subtyping_test" : "Output ID should be the original sample id. Found: ${id}" + + def state = output[1] + assert state instanceof Map : "State should be a map. Found: ${state}" + assert state.containsKey("output") : "Output should contain key 'output'." + assert state.output.isFile() : "'output' should be a file." + assert state.output.toString().endsWith(".h5mu") : "Output file should end with '.h5mu'. Found: ${state.output}" + + "Output: $output" + } + // The subtyped rna modality should carry the celltypist prediction columns, + // suffixed with the reference target they were transferred from. + | assert_h5mu_slots.run( + key: "assert_rna_slots", + fromState: { id, state -> + ["input": state.output, "modality": "rna", "obs": ["celltypist_pred_cell_type", "celltypist_probability_cell_type"]] + }, + // assert_h5mu_slots emits no output; keep the state for the next assert. + toState: { id, output, state -> state } + ) + // The untouched prot modality should still be present after the merge. + | assert_h5mu_slots.run( + key: "assert_prot_modality", + fromState: { id, state -> + ["input": state.output, "modality": "prot"] + } + ) + + // Verify the test case actually ran through the assert steps. + output_ch + | toSortedList { a, b -> a[0] <=> b[0] } + | map { output_list -> + assert output_list.size() == 1 : + "output channel should contain 1 event, found ${output_list.size()}" + assert output_list.collect { ev -> ev[0] } == ["subtyping_test"] : + "expected the test case to complete; found ${output_list.collect { ev -> ev[0] }}" + } +} + +workflow test_wf_2 { + // Exercise --allow_missing_reference_cell_type + resources_test = file(params.resources_test) + + output_ch = Channel.fromList( + [ + [ + id: "allow_missing_test", + input: resources_test.resolve("pbmc_1k_protein_v3/pbmc_1k_protein_v3_mms.h5mu"), + modality: "rna", + obs_major_cell_type: "compartment", + reference: resources_test.resolve("annotation_test_data/TS_Blood_filtered.h5mu"), + reference_var_gene_names: "ensemblid", + reference_obs_target: "cell_type", + reference_obs_major_cell_type: "compartment", + reference_var_input: "highly_variable", + annotation_methods: "celltypist", + allow_missing_reference_cell_type: true + ] + ]) + | view { "State at start: $it" } + | map { state -> [state.id, state] } + // Stamp the query cells with a major cell type absent from the reference. + | add_id.run( + fromState: ["input": "input"], + args: ["input_id": "missing_cell_type", "obs_output": "compartment", "make_observation_keys_unique": false], + toState: ["input": "output"] + ) + | parallel_subtyping + | view { "After parallel_subtyping: $it" } + | view { output -> + assert output.size() == 2 : "Outputs should contain two elements; [id, state]" + + def id = output[0] + assert id == "allow_missing_test" : "Output ID should be the original sample id. Found: ${id}" + + def state = output[1] + assert state instanceof Map : "State should be a map. Found: ${state}" + assert state.containsKey("output") : "Output should contain key 'output'." + assert state.output.isFile() : "'output' should be a file." + assert state.output.toString().endsWith(".h5mu") : "Output file should end with '.h5mu'. Found: ${state.output}" + + "Output: $output" + } + // The passed-through cells are merged back unannotated: the rna modality + // retains the major cell-type column but must carry no celltypist + // predictions (the annotation step is skipped for missing cell types). + | assert_h5mu_slots.run( + key: "assert_passthrough_rna_slots", + fromState: { id, state -> + [ + "input": state.output, + "modality": "rna", + "obs": ["compartment"], + "obs_absent": ["celltypist_pred_cell_type", "celltypist_probability_cell_type"] + ] + }, + toState: { id, output, state -> state } + ) + // The untouched prot modality should still be present after the merge. + | assert_h5mu_slots.run( + key: "assert_passthrough_prot_modality", + fromState: { id, state -> + ["input": state.output, "modality": "prot"] + } + ) + + // Verify the test case actually ran through the assert steps. + output_ch + | toSortedList { a, b -> a[0] <=> b[0] } + | map { output_list -> + assert output_list.size() == 1 : + "output channel should contain 1 event, found ${output_list.size()}" + assert output_list.collect { ev -> ev[0] } == ["allow_missing_test"] : + "expected the test case to complete; found ${output_list.collect { ev -> ev[0] }}" + } +} diff --git a/src/test_workflows/assert_h5mu_slots/config.vsh.yaml b/src/test_workflows/assert_h5mu_slots/config.vsh.yaml index 4faa7c3..83342e5 100644 --- a/src/test_workflows/assert_h5mu_slots/config.vsh.yaml +++ b/src/test_workflows/assert_h5mu_slots/config.vsh.yaml @@ -3,8 +3,9 @@ namespace: "test_workflows" scope: "test" description: | Test helper. Reads a modality from an output h5mu and asserts that the expected - .obs columns, .obsm keys, .obsp keys and .uns keys are all present. Raises an - error listing any missing slots. + .obs columns, .obsm keys, .obsp keys and .uns keys are all present, and that any + slots listed as absent are not present. Raises an error listing any missing slots + (expected but absent) or unexpected slots (listed as absent but present). authors: - __merge__: /src/authors/jakub_majercik.yaml @@ -45,6 +46,28 @@ argument_groups: multiple: true description: Expected .uns keys. required: false + - name: Unexpected slots + arguments: + - name: "--obs_absent" + type: string + multiple: true + description: .obs column names that must not be present. + required: false + - name: "--obsm_absent" + type: string + multiple: true + description: .obsm keys that must not be present. + required: false + - name: "--obsp_absent" + type: string + multiple: true + description: .obsp keys that must not be present. + required: false + - name: "--uns_absent" + type: string + multiple: true + description: .uns keys that must not be present. + required: false resources: - type: python_script diff --git a/src/test_workflows/assert_h5mu_slots/script.py b/src/test_workflows/assert_h5mu_slots/script.py index 92403bb..613217f 100644 --- a/src/test_workflows/assert_h5mu_slots/script.py +++ b/src/test_workflows/assert_h5mu_slots/script.py @@ -9,6 +9,10 @@ "obsm": [], "obsp": [], "uns": [], + "obs_absent": [], + "obsm_absent": [], + "obsp_absent": [], + "uns_absent": [], } meta = {"resources_dir": "src/utils/"} ## VIASH END @@ -43,12 +47,23 @@ if absent: missing[slot] = (absent, sorted(present_keys)) -if missing: - lines = [f"Output h5mu is missing expected slots in modality '{par['modality']}':"] +# Slots listed as absent must not be present in the modality. +unexpected = {} +for slot, present_keys in present.items(): + forbidden = [key for key in (par.get(f"{slot}_absent") or []) if key] + found = [key for key in forbidden if key in present_keys] + if found: + unexpected[slot] = found + +if missing or unexpected: + lines = [f"Output h5mu has unexpected slot state in modality '{par['modality']}':"] for slot, (absent, present_keys) in missing.items(): lines.append(f" .{slot}: missing {absent}; present {present_keys}") + for slot, found in unexpected.items(): + lines.append(f" .{slot}: should be absent but present {found}") raise AssertionError("\n".join(lines)) logger.info( - "All expected slots are present in modality '%s'.", par["modality"] + "All expected slots are present and all absent slots are missing in modality '%s'.", + par["modality"], )