Fix bulk Cirro pipeline config and OOM/resource issues

.cirro/bulk_analysis/preprocess.py

@@ -25,14 +25,4 @@
 samplesheet.to_csv('samplesheet.csv', index=None)
 ds.add_param("samplesheet", "samplesheet.csv")
 
-
-# 3. Set workflow_level value based on form input
-ds.logger.info("Setting workflow_level")
-
-levels = ['convert', 'sample', 'patient', 'compare']
-flags = [ds.params['convert_lvl'], ds.params['sample_lvl'], ds.params['patient_lvl'], ds.params['compare_lvl']]
-workflow_level = [lvl for lvl, flag in zip(levels, flags) if flag]
-
-ds.add_param('workflow_level', ','.join(workflow_level))
-
 ds.logger.info(ds.params)

.cirro/bulk_analysis/process-input.json

@@ -1,9 +1,8 @@
 {
     "input_format": "airr",
-    "convert_lvl": false,
-    "sample_lvl": "$.params.dataset.paramJson.sample_lvl",
-    "patient_lvl": "$.params.dataset.paramJson.patient_lvl",
-    "compare_lvl": "$.params.dataset.paramJson.compare_lvl",
+    "run_sample": "$.params.dataset.paramJson.sample_lvl",
+    "run_patient": "$.params.dataset.paramJson.patient_lvl",
+    "run_compare": "$.params.dataset.paramJson.compare_lvl",
     "olga_chunk_length": "$.params.dataset.paramJson.olga_chunk_length",
     "matrix_sparsity": "sparse",
     "distance_metric": "$.params.dataset.paramJson.distance_metric",

.cirro/convert_adaptive/preprocess.py

@@ -25,14 +25,4 @@
 samplesheet.to_csv('samplesheet.csv', index=None)
 ds.add_param("samplesheet", "samplesheet.csv")
 
-
-# 3. Set workflow_level value based on form input
-ds.logger.info("Setting workflow_level")
-
-levels = ['convert', 'sample', 'compare']
-flags = [ds.params['convert_lvl'], ds.params['sample_lvl'], ds.params['compare_lvl']]
-workflow_level = [lvl for lvl, flag in zip(levels, flags) if flag]
-
-ds.add_param('workflow_level', ','.join(workflow_level))
-
 ds.logger.info(ds.params)

.cirro/convert_adaptive/process-input.json

@@ -1,7 +1,7 @@
 {
     "input_format": "adaptive",
-    "convert_lvl": true,
-    "sample_lvl": false,
-    "compare_lvl": false,
+    "run_convert": true,
+    "run_sample": false,
+    "run_compare": false,
     "outdir": "$.params.dataset.s3|/data/"
 }

bin/samplesheet.py

@@ -9,7 +9,7 @@ def samplesheet(samplesheet):
     ss.to_csv('samplesheet_utf8.csv', index=False, encoding='utf-8-sig')
 
     stats = ss.describe()
-    stats.to_csv('samplesheet_stats.csv', index=False, encoding='utf-8-sig')
+    stats.to_csv('samplesheet_stats.csv', index=True, encoding='utf-8-sig')
 
     print(ss.head())
 

conf/modules.config

@@ -29,6 +29,10 @@ process {
             mode: params.publish_dir_mode,
             saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
         ]
+        // template_details_sample builds dense per-sample TCRdist matrices (OOM'd at 16GB); compare/discovery_brief process full VDJdb match files.
+        cpus   = { ['template_details_sample', 'template_details_compare', 'template_discovery_brief'].contains(notebook.baseName) ? 16 * task.attempt : 4 * task.attempt }
+        memory = { ['template_details_sample', 'template_details_compare', 'template_discovery_brief'].contains(notebook.baseName) ? 64.GB * task.attempt : 16.GB * task.attempt }
+        time   = { ['template_details_sample', 'template_details_compare', 'template_discovery_brief'].contains(notebook.baseName) ? 16.h : 4.h }
     }
 
 }

modules/local/report/render_notebook.nf

@@ -1,7 +1,6 @@
 // Generic process to render a Quarto notebook to HTML
 process RENDER_NOTEBOOK {
     tag "${notebook.getBaseName()}"
-    label 'process_single'
 
     input:
     // path(files) stages files flat in the root dir; staged_layout optionally
@@ -16,6 +15,8 @@ process RENDER_NOTEBOOK {
     tuple path(notebook), path(files), val(staged_layout)
     val project_name
     val workflow_cmd
+    // Fixed name avoids colliding with a same-named file in `files`.
+    path samplesheet, stageAs: 'render_notebook_samplesheet.csv'
 
     output:
     path "${notebook.getBaseName()}.html", emit: report_html
@@ -34,7 +35,7 @@ process RENDER_NOTEBOOK {
     quarto render ${notebook} \\
         -P project_name:${project_name} \\
         -P workflow_cmd:'${workflow_cmd}' \\
-        -P sample_table:${file(params.samplesheet)} \\
+        -P sample_table:${samplesheet} \\
         -P subject_col:'${params.subject_col}' \\
         -P timepoint_col:'${params.timepoint_col}' \\
         -P timepoint_order_col:'${params.timepoint_order_col}' \\

nextflow_schema.json

@@ -214,19 +214,19 @@
           "default": true
         },
         "local_min_pvalue": {
-          "type": "string",
-          "default": 0.001
+          "type": ["string", "number"],
+          "default": "0.001"
         },
         "simulation_depth": {
           "type": "string",
           "default": 1000
         },
         "kmer_min_depth": {
-          "type": "string",
-          "default": 3
+          "type": ["string", "integer"],
+          "default": "3"
         },
         "local_min_OVE": {
-          "type": "string",
+          "type": ["string", "integer"],
           "default": "c(1000, 100, 10)"
         },
         "all_aa_interchangeable": {

notebooks/template_discovery_brief.qmd

@@ -148,6 +148,10 @@ print('Date and time:         ' + str(datetime.datetime.now()))
 
 meta = pd.read_csv(sample_table, sep=',')
 
+# alias_col is optional; default to sample name if missing.
+if alias_col not in meta.columns:
+    meta[alias_col] = meta['sample']
+
 # timepoint_order_col isn't a samplesheet column - compute and inject it here.
 # timepoint_order (an ordered comma-separated list, e.g. "Base,Week4,EOT") ranks
 # timepoints by position in that list; any timepoint not listed sorts after the
@@ -892,7 +896,7 @@ def create_styled_table(df, title="", max_height='400px', formatter=None):
 import pandas as pd
 import numpy as np
 import itertools
-from scipy.stats import fisher_exact
+from scipy.stats import hypergeom
 from statsmodels.stats.multitest import multipletests
 
 # --- Setup & Pre-computation ---
@@ -911,11 +915,13 @@ sample_total_counts.rename(columns={'duplicate_count': 'total_counts'}, inplace=
 clonotypes_df = pd.merge(clonotypes_df, sample_total_counts, on=[subject_col, 'origin', timepoint_col])
 clonotypes_df['frequency'] = clonotypes_df['duplicate_count'] / clonotypes_df['total_counts']
 
-# Helper to run fisher exact test on rows
-def run_fisher(row, alt_hyp):
-    table = [[row['count_post'], row['total_post'] - row['count_post']],
-             [row['count_pre'], row['total_pre'] - row['count_pre']]]
-    return fisher_exact(table, alternative=alt_hyp)[1]
+# One-sided Fisher's exact test reduces to a single vectorized hypergeometric call since total_pre/total_post are fixed margins.
+def run_fisher_vectorized(count_post, count_pre, total_post, total_pre, alt_hyp):
+    M = total_pre + total_post
+    n = count_post + count_pre
+    if alt_hyp == 'greater':
+        return hypergeom.sf(count_post - 1, M, n, total_post)
+    return hypergeom.cdf(count_post, M, n, total_post)
 
 
 # ==========================================
@@ -968,7 +974,9 @@ for (subject, origin), subject_df in clonotypes_df.groupby([subject_col, 'origin
             continue
 
         # 5. Apply Fisher's Exact test only to the remaining rows
-        comp_df['p_value'] = comp_df.apply(run_fisher, axis=1, alt_hyp='greater')
+        comp_df['p_value'] = run_fisher_vectorized(
+            comp_df['count_post'], comp_df['count_pre'], comp_df['total_post'], comp_df['total_pre'], 'greater'
+        )
 
         # Add metadata
         comp_df[subject_col] = subject
@@ -1090,7 +1098,9 @@ for (subject, origin), subject_df in clonotypes_df.groupby([subject_col, 'origin
             continue
 
         # 5. Apply Fisher
-        comp_df['p_value'] = comp_df.apply(run_fisher, axis=1, alt_hyp='less')
+        comp_df['p_value'] = run_fisher_vectorized(
+            comp_df['count_post'], comp_df['count_pre'], comp_df['total_post'], comp_df['total_pre'], 'less'
+        )
 
         # Add metadata
         comp_df[subject_col] = subject
@@ -1652,17 +1662,12 @@ for subject, subj_df in clonotypes_df.groupby(subject_col):
             pivot['freq_t1_plot'] = pivot[t1].replace(0, min_freq)
             pivot['freq_t2_plot'] = pivot[t2].replace(0, min_freq)
 
-            # Assign highlight status using our sets
-            def assign_status(row):
-                key = (subject, origin, t1, t2, row['junction_aa'])
-                if key in sig_exp_lookup:
-                    return "Expanded"
-                elif key in sig_cont_lookup:
-                    return "Contracted"
-                else:
-                    return "Not Significant"
-
-            pivot['Status'] = pivot.apply(assign_status, axis=1)
+            # Avoids the per-row Series-construction overhead of .apply(axis=1).
+            keys = [(subject, origin, t1, t2, j) for j in pivot['junction_aa']]
+            pivot['Status'] = [
+                "Expanded" if key in sig_exp_lookup else ("Contracted" if key in sig_cont_lookup else "Not Significant")
+                for key in keys
+            ]
 
             # Sort so that colored points render on top of gray points
             pivot['sort_order'] = pivot['Status'].map({"Not Significant": 0, "Expanded": 1, "Contracted": 2})
@@ -1869,24 +1874,27 @@ else:
                 if len(tps) > 1: return COLOR_MULTI, 'Multiple Timepoints'
                 else: return timepoint_colors.get(tps[0], '#888888'), tps[0]
 
-            node_properties[['color', 'category']] = node_properties['timepoints'].apply(
-                lambda x: pd.Series(get_color_and_category(x))
-            )
+            colors_categories = [get_color_and_category(x) for x in node_properties['timepoints']]
+            node_properties['color'] = [c for c, _ in colors_categories]
+            node_properties['category'] = [cat for _, cat in colors_categories]
 
-            node_properties['size'] = node_properties['max_frequency'].apply(lambda x: max(10, np.sqrt(x * 100) * 15))
+            node_properties['size'] = np.maximum(10, np.sqrt(node_properties['max_frequency'] * 100) * 15)
 
             tcr_to_idx = {tcr: i for i, tcr in enumerate(node_properties['cdr3_b_aa'])}
 
             g = ig.Graph()
             g.add_vertices(len(node_properties))
             g.vs['name'] = node_properties['cdr3_b_aa']
-            g.vs['size'] = node_properties['size'] 
+            g.vs['size'] = node_properties['size']
             g.vs['color'] = node_properties['color']
             g.vs['category'] = node_properties['category']
-            g.vs['title'] = node_properties.apply(
-                lambda row: f"TCR: {row['cdr3_b_aa']}<br>Peak Freq: {(row['max_frequency']*100):.3f}%<br>Peak Count: {row['max_count']}<br>TPs: {row['timepoints']}",
-                axis=1
-            ).tolist()
+            g.vs['title'] = [
+                f"TCR: {cdr3}<br>Peak Freq: {freq*100:.3f}%<br>Peak Count: {count}<br>TPs: {tps}"
+                for cdr3, freq, count, tps in zip(
+                    node_properties['cdr3_b_aa'], node_properties['max_frequency'],
+                    node_properties['max_count'], node_properties['timepoints']
+                )
+            ]
 
             patient_samples = patient_df['sample'].unique()
             for sample_id in patient_samples:
@@ -2480,7 +2488,11 @@ def plot_public_clones_upset(concat_df, subject_col, min_shared_patients=2):
     # --- 1. Extract Unique Clones per Patient ---
     patient_clones = df_clean.groupby(subject_col)['junction_aa'].unique().to_dict()
     patient_contents = {str(patient): list(clones) for patient, clones in patient_clones.items()}
-
+
+    if len(patient_contents) < 2:
+        print("Only one patient present; skipping UpSet plot (needs ≥ 2 patients to show intersections).")
+        return
+
     # --- 2. Convert to UpSet Multi-Index Format ---
     upset_data = from_contents(patient_contents)
 

notebooks/template_giana.qmd

@@ -37,6 +37,11 @@ warnings.filterwarnings(
 # Loading data
 ## reading sample metadata
 meta = pd.read_csv(sample_table, sep=',')
+
+# alias_col is optional; default to sample name if missing.
+if alias_col not in meta.columns:
+    meta[alias_col] = meta['sample']
+
 meta_cols = meta.columns.tolist()
 
 ```

notebooks/template_gliph.qmd

@@ -52,6 +52,10 @@ if timepoint_order_col not in meta.columns:
             rank_map = {t: r for r, t in enumerate(sorted(unique_timepoints, key=str))}
     meta[timepoint_order_col] = meta[timepoint_col].map(rank_map)
 
+# alias_col is optional; default to sample name if missing.
+if alias_col not in meta.columns:
+    meta[alias_col] = meta['sample']
+
 concat_df = pd.read_csv(concat_csv, sep='\t')
 concat_df = concat_df.merge(meta[['sample', 'origin', timepoint_col, timepoint_order_col, alias_col]], on='sample', how='left')