From 6bcad7428389b647e6643fcae84bd2e0ebd98f86 Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Tue, 21 Jul 2026 17:23:19 -0400 Subject: [PATCH 1/4] Wire report notebooks into pipeline and fix real render bugs Extends RENDER_NOTEBOOK/REPORT wiring to template_discovery_brief.qmd, template_details_part1.qmd, and template_details_part2.qmd, and fixes bugs surfaced by actually executing all 7 underlying notebooks (template_sample, template_overlap, template_sharing, template_giana, template_gliph, template_pheno_sc, template_pheno_bulk) against real pipeline-generated data, including a full end-to-end Nextflow run, rather than assuming correctness from static review. Pipeline wiring: - Add named emits for sample/patient/compare subworkflow outputs needed by the new reports (v_family, j_family, tcrdist/olga/vdjdb files, convergence files, tcrpheno files, shared_cdr3, giana/gliph2 outputs). - Stage {{< include >}} sibling .qmd files alongside their parent notebook in RENDER_NOTEBOOK, since Nextflow only auto-stages the named notebook file. - Give GLIPH2_TURBOGLIPH outputs patient-prefixed filenames (_all_motifs.txt, etc.) to avoid Nextflow's fatal input-file-name collision when staging multiple patients' outputs together, matching the existing GIANA_CALC convention. Notebooks read these prefixed names directly instead of requiring a rename step. - Wire CONVERT.out.sample_map_converted into the discovery brief's staged files so the VDJdb section reads AIRR-converted files instead of falling back to raw adaptive TSVs. - Make template_discovery_brief.qmd's phenotype section conditional: a single generic {{< include ./template_pheno.qmd >}} is spliced in, and the pipeline stages whichever real notebook applies under that name - template_pheno_sc.qmd for cellranger input with sobject_gex supplied, template_pheno_bulk.qmd otherwise - since Quarto's include shortcode has no native runtime conditional. - Register template_pheno_sc/template_pheno_bulk as pipeline params. Notebook fixes (found via end-to-end render testing): - Fix systemic legacy-vs-AIRR column naming bug (CDR3b/counts/TRBV/TRBJ vs junction_aa/duplicate_count/v_call/j_call) across template_discovery_brief, template_sample, template_overlap, and template_sharing, wherever the data originates from AIRR/concat outputs. Left legacy names intact where they genuinely belong (GIANA/GLIPH2 module outputs, tcrdist3/sample_stats schemas). - Fix convergence file merge silently producing zero rows by handling the filename suffix when parsing sample IDs. - Fix 3 VDJdb plotting functions and template_pheno_bulk.qmd's phenotype-composition functions silently using hardcoded column names (subject_id/alias/timepoint_order) instead of the notebook's configured subject_col/alias_col/timepoint_order_col. - Fix VDJdb converted-file lookup to match by sample name (matching how CONVERT_ADAPTIVE actually names its outputs) instead of reconstructing a filename from the raw input basename. - Fix template_pheno_bulk.qmd's weighted-phenotype merge, which joined TCRPHENO scores against clone counts on a CDR3/junction_aa column that TCRPHENO's real R-script output never contains - confirmed by running the actual tcrpheno R package - and rewrite it to join on sequence_id, which both files share. - Point template_pheno_sc.qmd's phenotype clonotype source at annotate/concatenated_cdr3_sorted.tsv (matching the rest of the pipeline) instead of a compare_phenotype/ path no module produces, and remove dead samplesheet_phenotype.csv loading code. - Remove genuinely dead imports (logomaker and others) across all notebooks, verified individually via pyflakes plus manual chain tracing to avoid removing imports still needed by earlier-executing cells in the Quarto include splice. Testing: - Update render_notebook.nf.test for the new [dest, src] staged_layout tuple format, with cases for nested paths and renamed basenames. --- .gitignore | 5 +- modules/local/compare/giana.nf | 4 +- modules/local/compare/gliph2.nf | 19 +- modules/local/report/render_notebook.nf | 21 +- modules/local/sample/tcrspecificity.nf | 6 +- nextflow.config | 5 + nextflow_schema.json | 20 ++ notebooks/template_discovery_brief.qmd | 121 ++++++----- notebooks/template_giana.qmd | 51 +---- notebooks/template_gliph.qmd | 180 +++++++++-------- notebooks/template_overlap.qmd | 32 ++- notebooks/template_pheno_bulk.qmd | 22 +- notebooks/template_pheno_sc.qmd | 18 +- notebooks/template_sample.qmd | 60 ++---- notebooks/template_sharing.qmd | 15 +- subworkflows/local/compare.nf | 3 + subworkflows/local/patient.nf | 20 +- subworkflows/local/sample.nf | 15 +- .../local/report/render_notebook.nf.test | 78 ++++++- workflows/tcrtoolkit.nf | 191 +++++++++++++++++- 20 files changed, 584 insertions(+), 302 deletions(-) diff --git a/.gitignore b/.gitignore index 9c9b98a..e9ff6c1 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ runs/* results* work out-* +out/* ## nf-test .nf-test/ @@ -34,4 +35,6 @@ notebooks/* tmp ## vscode -.vscode/* \ No newline at end of file +.vscode/* + +.e2e_test_tmp/ diff --git a/modules/local/compare/giana.nf b/modules/local/compare/giana.nf index 0537f8a..fdaea2c 100644 --- a/modules/local/compare/giana.nf +++ b/modules/local/compare/giana.nf @@ -9,8 +9,8 @@ process GIANA_CALC { val threshold_vgene output: - path "${patient}_VgeneScores.txt" - path "${patient}_giana.txt" + path "${patient}_VgeneScores.txt", emit: 'vgene_scores' + path "${patient}_giana.txt", emit: 'giana_output' // path "giana_EncodingMatrix.txt" script: diff --git a/modules/local/compare/gliph2.nf b/modules/local/compare/gliph2.nf index 01b76e7..ec71b91 100644 --- a/modules/local/compare/gliph2.nf +++ b/modules/local/compare/gliph2.nf @@ -8,11 +8,15 @@ process GLIPH2_TURBOGLIPH { tuple val(patient), path(concat_cdr3) output: - path "${patient}/all_motifs.txt", emit: 'all_motifs' - path "${patient}/clone_network.txt", emit: 'clone_network' - path "${patient}/cluster_member_details.txt", emit: 'cluster_member_details' + // all_motifs/clone_network/cluster_member_details/global_similarities are + // copied to patient-prefixed top-level names because multiple patients' + // outputs share the same basename otherwise (e.g. "all_motifs.txt"), which + // collides when several patients' files are staged together downstream. + tuple val(patient), path("${patient}_all_motifs.txt"), emit: 'all_motifs' + tuple val(patient), path("${patient}_clone_network.txt"), emit: 'clone_network' + tuple val(patient), path("${patient}_cluster_member_details.txt"), emit: 'cluster_member_details' path "${patient}/convergence_groups.txt", emit: 'convergence_groups' - path "${patient}/global_similarities.txt", emit: 'global_similarities' + tuple val(patient), path("${patient}_global_similarities.txt"), emit: 'global_similarities' path "${patient}/local_similarities.txt", emit: 'local_similarities' path "${patient}/parameter.txt", emit: 'gliph2_parameters' @@ -51,6 +55,13 @@ process GLIPH2_TURBOGLIPH { # Rename local_similarities file to standardize output name input_file="${patient}/local_similarities_*.txt" cat \$input_file > ${patient}/local_similarities.txt + + # Copy to patient-prefixed top-level names to avoid basename collisions + # when multiple patients' outputs are staged together downstream. + cp ${patient}/all_motifs.txt ${patient}_all_motifs.txt + cp ${patient}/clone_network.txt ${patient}_clone_network.txt + cp ${patient}/cluster_member_details.txt ${patient}_cluster_member_details.txt + cp ${patient}/global_similarities.txt ${patient}_global_similarities.txt """ } diff --git a/modules/local/report/render_notebook.nf b/modules/local/report/render_notebook.nf index a7f8136..78b5e69 100644 --- a/modules/local/report/render_notebook.nf +++ b/modules/local/report/render_notebook.nf @@ -4,7 +4,16 @@ process RENDER_NOTEBOOK { label 'process_single' input: - tuple path(notebook), path(files) // path(files) just stages files in root dir + // path(files) stages files flat in the root dir; staged_layout optionally + // symlinks them into a project_dir-style subdirectory tree for notebooks that + // read from a nested project_dir/subdir/file layout instead of bare filenames. + // It's a list of [dest_path, source_basename] pairs (e.g. + // ["sample/sample_stats.csv", "sample_stats.csv"]) rather than just a dest + // path, because source and dest basenames can legitimately differ - e.g. + // per-patient files staged under a shared basename (like GLIPH2's + // "patientA_all_motifs.txt") get symlinked to their unprefixed, per-patient + // subdirectory destination ("gliph2/patientA/all_motifs.txt"). + tuple path(notebook), path(files), val(staged_layout) val project_name val workflow_cmd @@ -12,12 +21,22 @@ process RENDER_NOTEBOOK { path "${notebook.getBaseName()}.html", emit: report_html script: + def stage_cmds = staged_layout.collect { dest, src -> + "mkdir -p \"\$(dirname '${dest}')\"; ln -sf \"\$PWD/${src}\" '${dest}'" + }.join('\n ') + def project_dir_arg = staged_layout ? "-P project_dir:'.'" : '' """ + ${stage_cmds} ## render qmd report to html quarto render ${notebook} \\ -P project_name:${project_name} \\ -P workflow_cmd:'${workflow_cmd}' \\ -P sample_table:${file(params.samplesheet)} \\ + -P subject_col:'${params.subject_col}' \\ + -P timepoint_col:'${params.timepoint_col}' \\ + -P timepoint_order_col:'${params.timepoint_order_col}' \\ + -P alias_col:'${params.alias_col}' \\ + ${project_dir_arg} \\ --to html """ diff --git a/modules/local/sample/tcrspecificity.nf b/modules/local/sample/tcrspecificity.nf index 30ce184..fd8c514 100644 --- a/modules/local/sample/tcrspecificity.nf +++ b/modules/local/sample/tcrspecificity.nf @@ -20,9 +20,9 @@ process VDJDB_VDJMATCH { path(ref_db) output: - path("${sample_meta.sample}.vdjmatch.txt") - path("${sample_meta.sample}.annot.summary.txt") - path "logs/${sample_meta.sample}.vdjmatch.log" + path("${sample_meta.sample}.vdjmatch.txt"), emit: 'vdjmatch_txt' + path("${sample_meta.sample}.annot.summary.txt"), emit: 'annot_summary' + path "logs/${sample_meta.sample}.vdjmatch.log", emit: 'log' script: def memGb = (task.memory.toMega() * 0.8 / 1024).intValue() diff --git a/nextflow.config b/nextflow.config index ccc4534..fac5172 100644 --- a/nextflow.config +++ b/nextflow.config @@ -75,6 +75,11 @@ params { //reports template_qc = "${projectDir}/notebooks/template_qc.qmd" + template_discovery_brief = "${projectDir}/notebooks/template_discovery_brief.qmd" + template_details_part1 = "${projectDir}/notebooks/template_details_part1.qmd" + template_details_part2 = "${projectDir}/notebooks/template_details_part2.qmd" + template_pheno_sc = "${projectDir}/notebooks/template_pheno_sc.qmd" + template_pheno_bulk = "${projectDir}/notebooks/template_pheno_bulk.qmd" } includeConfig 'conf/base.config' diff --git a/nextflow_schema.json b/nextflow_schema.json index 26d3c20..a40abbf 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -120,6 +120,26 @@ "template_qc": { "type": "string", "description": "Path to QC notebook template." + }, + "template_discovery_brief": { + "type": "string", + "description": "Path to discovery brief notebook template." + }, + "template_details_part1": { + "type": "string", + "description": "Path to details (part 1) notebook template." + }, + "template_details_part2": { + "type": "string", + "description": "Path to details (part 2) notebook template." + }, + "template_pheno_sc": { + "type": "string", + "description": "Path to single-cell phenotype notebook template." + }, + "template_pheno_bulk": { + "type": "string", + "description": "Path to bulk (TCRpheno) phenotype notebook template." } } }, diff --git a/notebooks/template_discovery_brief.qmd b/notebooks/template_discovery_brief.qmd index 7956ad9..ad1c3ae 100644 --- a/notebooks/template_discovery_brief.qmd +++ b/notebooks/template_discovery_brief.qmd @@ -475,20 +475,20 @@ def calculate_overlaps(df, meta): # Ensure patient mapping is available sample_to_patient = meta.set_index('sample')[subject_col].to_dict() df[subject_col] = df['sample'].map(sample_to_patient) - df = df.dropna(subset=[subject_col, 'CDR3b', 'counts']) + df = df.dropna(subset=[subject_col, 'junction_aa', 'duplicate_count']) samples = sorted(df['sample'].unique()) patients = sorted(df[subject_col].unique()) - clones = sorted(df['CDR3b'].unique()) - + clones = sorted(df['junction_aa'].unique()) + clone_map = {c: i for i, c in enumerate(clones)} n_clones = len(clones) # --- 2. MORISITA (Sample-level, Frequency-based) --- sample_map = {s: i for i, s in enumerate(samples)} row_idx_s = df['sample'].map(sample_map).values - col_idx_s = df['CDR3b'].map(clone_map).values - values_s = df['counts'].values + col_idx_s = df['junction_aa'].map(clone_map).values + values_s = df['duplicate_count'].values mat_s = sparse.coo_matrix((values_s, (row_idx_s, col_idx_s)), shape=(len(samples), n_clones)).tocsr() @@ -501,12 +501,12 @@ def calculate_overlaps(df, meta): morisita_df = pd.DataFrame(morisita, index=samples, columns=samples) # --- 3. JACCARD (Patient-level, Binary-based) --- - # Aggregate counts by patient/CDR3b to get unique clones per patient - patient_df = df.groupby([subject_col, 'CDR3b']).size().reset_index() - + # Aggregate counts by patient/junction_aa to get unique clones per patient + patient_df = df.groupby([subject_col, 'junction_aa']).size().reset_index() + pat_map = {p: i for i, p in enumerate(patients)} row_idx_p = patient_df[subject_col].map(pat_map).values - col_idx_p = patient_df['CDR3b'].map(clone_map).values + col_idx_p = patient_df['junction_aa'].map(clone_map).values # Binary matrix: 1 if patient has clone, else 0 mat_p_bin = sparse.coo_matrix((np.ones(len(patient_df)), (row_idx_p, col_idx_p)), @@ -644,8 +644,8 @@ def create_upset_tabs(df, patient_list): patient_df = df[df[subject_col] == patient_id] timepoints = sorted(patient_df[timepoint_col].unique()) - # Extracting clones (CDR3b) by timepoint - clones_by_timepoint = {tp: set(patient_df[patient_df[timepoint_col] == tp]['CDR3b']) for tp in timepoints} + # Extracting clones (junction_aa) by timepoint + clones_by_timepoint = {tp: set(patient_df[patient_df[timepoint_col] == tp]['junction_aa']) for tp in timepoints} upset_data = upsetplot.from_contents(clones_by_timepoint) @@ -740,11 +740,11 @@ for patient_id in upset_patients: patient_df = concat_df[concat_df[subject_col] == patient_id] # 1. Aggregate counts for each clone at each timepoint - # (Using 'counts' as the abundance metric based on your previous structures) - clone_tp_data = patient_df.groupby(['CDR3b', timepoint_col])['counts'].sum().reset_index() + # (Using 'duplicate_count' as the abundance metric based on your previous structures) + clone_tp_data = patient_df.groupby(['junction_aa', timepoint_col])['duplicate_count'].sum().reset_index() # 2. Calculate the number of timepoints and list them for each clone - clone_stats = clone_tp_data.groupby('CDR3b').agg( + clone_stats = clone_tp_data.groupby('junction_aa').agg( n_timepoints=(timepoint_col, 'nunique'), timepoints_present=(timepoint_col, lambda x: ', '.join(sorted(x.astype(str)))) ) @@ -757,12 +757,12 @@ for patient_id in upset_patients: # 4. Pivot the data to get frequencies (counts) as separate columns per timepoint persistent_clones = persistent_stats.index - persistent_counts_df = clone_tp_data[clone_tp_data['CDR3b'].isin(persistent_clones)] + persistent_counts_df = clone_tp_data[clone_tp_data['junction_aa'].isin(persistent_clones)] wide_counts = persistent_counts_df.pivot( - index='CDR3b', + index='junction_aa', columns=timepoint_col, - values='counts' + values='duplicate_count' ).fillna(0) # Fill missing timepoints with 0 counts # Prefix the timepoint columns so they are clearly identifiable as counts @@ -874,17 +874,17 @@ from statsmodels.stats.multitest import multipletests clonotypes_df = concat_df.copy() clonotypes_df = clonotypes_df.dropna(subset=[timepoint_col, 'origin']) -clone_cols = ['CDR3b', 'TRBV', 'TRBJ'] +clone_cols = ['junction_aa', 'v_call', 'j_call'] merge_cols = clone_cols + [subject_col, 'origin'] clonotypes_df['timepoint_rank'] = clonotypes_df[timepoint_order_col] clonotypes_df = clonotypes_df.sort_values('timepoint_rank') # Grouping by subject, origin, and timepoint for accurate total counts -sample_total_counts = clonotypes_df.groupby([subject_col, 'origin', timepoint_col])['counts'].sum().reset_index() -sample_total_counts.rename(columns={'counts': 'total_counts'}, inplace=True) +sample_total_counts = clonotypes_df.groupby([subject_col, 'origin', timepoint_col])['duplicate_count'].sum().reset_index() +sample_total_counts.rename(columns={'duplicate_count': 'total_counts'}, inplace=True) clonotypes_df = pd.merge(clonotypes_df, sample_total_counts, on=[subject_col, 'origin', timepoint_col]) -clonotypes_df['frequency'] = clonotypes_df['counts'] / clonotypes_df['total_counts'] +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): @@ -911,8 +911,8 @@ for (subject, origin), subject_df in clonotypes_df.groupby([subject_col, 'origin continue freq_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='frequency', fill_value=0) - counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='counts', fill_value=0) - subject_total_counts = subject_df.groupby(timepoint_col)['counts'].sum() + counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='duplicate_count', fill_value=0) + subject_total_counts = subject_df.groupby(timepoint_col)['duplicate_count'].sum() for t_pre, t_post in itertools.combinations(timepoints, 2): # 1. Vectorized DataFrame creation instead of iterrows @@ -1036,8 +1036,8 @@ for (subject, origin), subject_df in clonotypes_df.groupby([subject_col, 'origin continue freq_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='frequency', fill_value=0) - counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='counts', fill_value=0) - subject_total_counts = subject_df.groupby(timepoint_col)['counts'].sum() + counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='duplicate_count', fill_value=0) + subject_total_counts = subject_df.groupby(timepoint_col)['duplicate_count'].sum() for t_pre, t_post in itertools.combinations(timepoints, 2): # 1. Vectorize @@ -1149,7 +1149,7 @@ import pandas as pd # --- 1. Identify "New" vs "Pre-existing" via Cumulative History --- -present_clones = clonotypes_df[clonotypes_df['counts'] > 0].copy() +present_clones = clonotypes_df[clonotypes_df['duplicate_count'] > 0].copy() present_clones = present_clones.sort_values([subject_col, timepoint_order_col]) dynamic_results = [] @@ -1161,7 +1161,7 @@ if 'detailed_comparisons_table' in locals() and not detailed_comparisons_table.e (detailed_comparisons_table['fold_change'] > fc_thold_up) & (pd.to_numeric(detailed_comparisons_table['q_value'], errors='coerce') < signif_thold) ] - expanded_set = set(zip(sig_exp[subject_col], sig_exp['origin'], sig_exp['t_post'], sig_exp['CDR3b'])) + expanded_set = set(zip(sig_exp[subject_col], sig_exp['origin'], sig_exp['t_post'], sig_exp['junction_aa'])) contracted_set = set() if 'detailed_comparisons_table_cont' in locals() and not detailed_comparisons_table_cont.empty: @@ -1169,7 +1169,7 @@ if 'detailed_comparisons_table_cont' in locals() and not detailed_comparisons_ta (detailed_comparisons_table_cont['fold_change'] < fc_thold_down) & (pd.to_numeric(detailed_comparisons_table_cont['q_value'], errors='coerce') < signif_thold) ] - contracted_set = set(zip(sig_cont[subject_col], sig_cont['origin'], sig_cont['t_post'], sig_cont['CDR3b'])) + contracted_set = set(zip(sig_cont[subject_col], sig_cont['origin'], sig_cont['t_post'], sig_cont['junction_aa'])) # Iterate over both subject AND origin @@ -1179,7 +1179,7 @@ for (subject, origin), subj_df in present_clones.groupby([subject_col, 'origin'] for rank in timepoints: tp_name = subj_df[subj_df[timepoint_order_col] == rank][timepoint_col].iloc[0] - current_clones = subj_df[subj_df[timepoint_order_col] == rank]['CDR3b'].unique() + current_clones = subj_df[subj_df[timepoint_order_col] == rank]['junction_aa'].unique() for clone in current_clones: status = "Stable" # Default @@ -1202,7 +1202,7 @@ for (subject, origin), subj_df in present_clones.groupby([subject_col, 'origin'] 'origin': origin, # Ensure origin is appended here! timepoint_col: tp_name, timepoint_order_col: rank, - 'CDR3b': clone, + 'junction_aa': clone, 'Status': status }) @@ -1356,19 +1356,19 @@ In this visualization, the total unique repertoire at each timepoint is broken d # --- Save Detailed Clone Categories with Info for further analysis --- # Merge the categorical statuses back with the original data to keep counts, frequencies, and gene usage -export_cols = [subject_col, timepoint_col, 'CDR3b', 'TRBV', 'TRBJ', 'counts', 'frequency'] +export_cols = [subject_col, timepoint_col, 'junction_aa', 'v_call', 'j_call', 'duplicate_count', 'frequency'] export_cols = [c for c in export_cols if c in present_clones.columns] # Safely keep only existing columns detailed_export_df = pd.merge( plot_df, # Contains subject, timepoint, timepoint_order, CDR3b, Status present_clones[export_cols], - on=[subject_col, timepoint_col, 'CDR3b'], + on=[subject_col, timepoint_col, 'junction_aa'], how='left' ) # Sort it logically: By Patient -> Chronological Timepoint -> Status -> Largest Clones first detailed_export_df = detailed_export_df.sort_values( - by=[subject_col, timepoint_order_col, 'Status', 'counts'], + by=[subject_col, timepoint_order_col, 'Status', 'duplicate_count'], ascending=[True, True, True, False] ) ``` @@ -1411,10 +1411,10 @@ if 'detailed_comparisons_table_cont' in locals() and not detailed_comparisons_ta # Combine them if sig_stats_list: all_sig_stats = pd.concat(sig_stats_list, ignore_index=True) - stats_to_merge = all_sig_stats[[subject_col, 'origin', 't_post', 'CDR3b', 't_pre', 'freq_pre', 'freq_post', 'p_value']] + stats_to_merge = all_sig_stats[[subject_col, 'origin', 't_post', 'junction_aa', 't_pre', 'freq_pre', 'freq_post', 'p_value']] else: # Fallback if there are absolutely zero significant clones in the entire dataset - stats_to_merge = pd.DataFrame(columns=[subject_col, 'origin', 't_post', 'CDR3b', 't_pre', 'freq_pre', 'freq_post', 'p_value']) + stats_to_merge = pd.DataFrame(columns=[subject_col, 'origin', 't_post', 'junction_aa', 't_pre', 'freq_pre', 'freq_post', 'p_value']) # 4. Merge the plotting classifications with the exact statistical comparisons # Note: A clone might match multiple 't_pre' comparisons. This merge will mathematically branch @@ -1422,16 +1422,16 @@ else: merged_df = pd.merge( dynamic_clones_df, stats_to_merge, - on=[subject_col, 'origin', 't_post', 'CDR3b'], + on=[subject_col, 'origin', 't_post', 'junction_aa'], how='left' ) # 5. Fetch actual frequencies for "New" clones (since they don't have a Fisher test row) -current_freqs = clonotypes_df[[subject_col, 'origin', timepoint_col, 'CDR3b', 'frequency', 'TRBV', 'TRBJ']].rename( +current_freqs = clonotypes_df[[subject_col, 'origin', timepoint_col, 'junction_aa', 'frequency', 'v_call', 'j_call']].rename( columns={timepoint_col: 't_post', 'frequency': 'current_frequency'} ).drop_duplicates() -merged_df = pd.merge(merged_df, current_freqs, on=[subject_col, 'origin', 't_post', 'CDR3b'], how='left') +merged_df = pd.merge(merged_df, current_freqs, on=[subject_col, 'origin', 't_post', 'junction_aa'], how='left') # Fill in the blanks for the "New" clones merged_df['freq_post'] = merged_df['freq_post'].fillna(merged_df['current_frequency']) @@ -1441,7 +1441,7 @@ merged_df['t_pre'] = merged_df['t_pre'].fillna("Not Present Previously") # 6. Clean up and Rename Columns for the final export final_cols = [ - subject_col, 'origin', 'CDR3b', 'TRBV', 'TRBJ', + subject_col, 'origin', 'junction_aa', 'v_call', 'j_call', 'Status', 't_post', 't_pre', 'freq_post', 'freq_pre', 'p_value' ] export_df = merged_df[[c for c in final_cols if c in merged_df.columns]].copy() @@ -1556,7 +1556,7 @@ if 'detailed_comparisons_table' in locals() and not detailed_comparisons_table.e (pd.to_numeric(detailed_comparisons_table['q_value'], errors='coerce') < signif_thold) ] # Store tuple of (subject, origin, t_pre, t_post, CDR3b) - sig_exp_lookup = set(zip(sig_exp[subject_col], sig_exp['origin'], sig_exp['t_pre'], sig_exp['t_post'], sig_exp['CDR3b'])) + sig_exp_lookup = set(zip(sig_exp[subject_col], sig_exp['origin'], sig_exp['t_pre'], sig_exp['t_post'], sig_exp['junction_aa'])) sig_cont_lookup = set() if 'detailed_comparisons_table_cont' in locals() and not detailed_comparisons_table_cont.empty: @@ -1564,7 +1564,7 @@ if 'detailed_comparisons_table_cont' in locals() and not detailed_comparisons_ta (detailed_comparisons_table_cont['fold_change'] < fc_thold_down) & (pd.to_numeric(detailed_comparisons_table_cont['q_value'], errors='coerce') < signif_thold) ] - sig_cont_lookup = set(zip(sig_cont[subject_col], sig_cont['origin'], sig_cont['t_pre'], sig_cont['t_post'], sig_cont['CDR3b'])) + sig_cont_lookup = set(zip(sig_cont[subject_col], sig_cont['origin'], sig_cont['t_pre'], sig_cont['t_post'], sig_cont['junction_aa'])) # Define our color palette color_discrete_map = { @@ -1600,7 +1600,7 @@ for subject, subj_df in clonotypes_df.groupby(subject_col): continue pivot = pair_df.pivot_table( - index=['CDR3b', 'TRBV', 'TRBJ'], + index=['junction_aa', 'v_call', 'j_call'], columns=timepoint_col, values='frequency', fill_value=0 @@ -1629,7 +1629,7 @@ for subject, subj_df in clonotypes_df.groupby(subject_col): # Assign highlight status using our sets def assign_status(row): - key = (subject, origin, t1, t2, row['CDR3b']) + key = (subject, origin, t1, t2, row['junction_aa']) if key in sig_exp_lookup: return "Expanded" elif key in sig_cont_lookup: @@ -1651,8 +1651,8 @@ for subject, subj_df in clonotypes_df.groupby(subject_col): color='Status', color_discrete_map=color_discrete_map, hover_data={ - 'CDR3b': True, - 'TRBV': True, + 'junction_aa': True, + 'v_call': True, t1: ':.6f', # Show the true frequency on hover (including actual 0s) t2: ':.6f', 'freq_t1_plot': False, @@ -1908,7 +1908,7 @@ else: (detailed_export_df[subject_col] == patient) & (detailed_export_df['Status'] == 'Expanded') ] - expanded_tcr_set = set(patient_expanded_df['CDR3b'].tolist()) + expanded_tcr_set = set(patient_expanded_df['junction_aa'].tolist()) # 2. Iterate through connected components and keep only those with at least one expanded clone keep_vids = [] @@ -1988,12 +1988,9 @@ for _, row in meta.iterrows(): filepath = row['file'] if os.path.exists(convert_dir) and os.listdir(convert_dir): - base_name = os.path.basename(row['file']) - name_part, ext_part = os.path.splitext(base_name) - new_file_name = f"{name_part}_airr{ext_part}" - - # Update the path to point to convert_dir - target_file = os.path.join(convert_dir, new_file_name) + # Converted files are named by sample (e.g. "_airr.tsv"), not by + # the original raw input filename. + target_file = os.path.join(convert_dir, f"{row['sample']}_airr.tsv") else: target_file = row['file'] @@ -2200,7 +2197,7 @@ import numpy as np exp_df = detailed_export_df[detailed_export_df['Status'] == 'Expanded'].copy() rename_dict = {} -if 'CDR3b' in exp_df.columns: rename_dict['CDR3b'] = 'cdr3aa' +if 'junction_aa' in exp_df.columns: rename_dict['junction_aa'] = 'cdr3aa' if subject_col in exp_df.columns and subject_col != subject_col: rename_dict[subject_col] = subject_col exp_df = exp_df.rename(columns=rename_dict) @@ -2349,18 +2346,18 @@ import pandas as pd def get_public_clones_df(concat_df, pgen_df): """Calculates and formats the full table of public (shared) clones.""" # Create the "Wide" Matrix (One Column per Patient) - patient_map = concat_df.groupby(['CDR3b', subject_col])[alias_col].apply( + patient_map = concat_df.groupby(['junction_aa', subject_col])[alias_col].apply( lambda x: ', '.join(sorted(x.unique().astype(str))) ).unstack(fill_value='-') # Calculate Summary Stats (Total Patients) - stats_df = concat_df.groupby('CDR3b').agg( + stats_df = concat_df.groupby('junction_aa').agg( total_patients=(subject_col, 'nunique') ).reset_index() # Merge Everything Together - df = pd.merge(stats_df, pgen_df[['CDR3b', 'pgen']], on='CDR3b', how='inner') - df = pd.merge(df, patient_map, on='CDR3b', how='inner') + df = pd.merge(stats_df, pgen_df[['junction_aa', 'pgen']], on='junction_aa', how='inner') + df = pd.merge(df, patient_map, on='junction_aa', how='inner') # Filter for public clones (>1 patient) and Sort df = df[df['total_patients'] > 1].copy() @@ -2369,7 +2366,7 @@ def get_public_clones_df(concat_df, pgen_df): # Formatting for Display/Export df['pgen'] = df['pgen'].apply(lambda x: f"{x:.2e}") df.rename(columns={ - 'CDR3b': 'CDR3b Sequence', + 'junction_aa': 'CDR3b Sequence', 'total_patients': 'No_Shared_Individuals', # Updated column name 'pgen': 'Pgen' }, inplace=True) @@ -2454,10 +2451,10 @@ def plot_public_clones_upset(concat_df, subject_col, min_shared_patients=2): # --- DEFENSIVE FIX 1: Drop NaNs --- # Missing sequences cause duplicate MultiIndex entries, breaking Matplotlib text rendering - df_clean = concat_df.dropna(subset=[subject_col, 'CDR3b']).copy() + df_clean = concat_df.dropna(subset=[subject_col, 'junction_aa']).copy() # --- 1. Extract Unique Clones per Patient --- - patient_clones = df_clean.groupby(subject_col)['CDR3b'].unique().to_dict() + 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()} # --- 2. Convert to UpSet Multi-Index Format --- @@ -2526,6 +2523,4 @@ plot_public_clones_upset(concat_df, subject_col=subject_col, min_shared_patients ``` **Figure 8. UpSet plot of public TCR clones**. This plot illustrates the intersection of unique CDR3b sequences shared among multiple individuals. The horizontal bars on the left indicate the total number of unique clones identified in each individual patient. In the central matrix, connected black dots denote the specific combination of patients sharing a subset of clones. The vertical bars at the top represent the intersection size, indicating the exact number of public clones shared exclusively by the patients marked in the corresponding column below. -{{< include ./template_pheno_sc.qmd >}} - -{{< include ./template_pheno_bulk.qmd >}} \ No newline at end of file +{{< include ./template_pheno.qmd >}} \ No newline at end of file diff --git a/notebooks/template_giana.qmd b/notebooks/template_giana.qmd index e8f1b7e..b850832 100644 --- a/notebooks/template_giana.qmd +++ b/notebooks/template_giana.qmd @@ -18,41 +18,13 @@ notebook_results_dir = f"{project_dir}/notebook-analysis/" #| code-fold: true # 1. Load Packages -from IPython.display import Image -import os -import datetime -import sys import pandas as pd -import math -import matplotlib.pyplot as plt -import seaborn as sns -from matplotlib.colors import LinearSegmentedColormap import plotly.express as px import plotly.graph_objects as go -import glob import itertools -import h5py -import igraph as ig -import matplotlib.ticker as ticker import numpy as np import scipy.cluster.hierarchy as sch from IPython.display import display, Markdown -from scipy.sparse import csr_matrix -from scipy.stats import gaussian_kde -from scipy.stats import entropy -from scipy.stats import skew, wasserstein_distance -from scipy.stats import pearsonr -from sklearn.decomposition import PCA -from sklearn.preprocessing import StandardScaler -from scipy import sparse -from sklearn.preprocessing import normalize -from scipy.stats import wilcoxon -from scipy.stats import mannwhitneyu -from itertools import combinations -from scipy.spatial.distance import pdist -from scipy.cluster.hierarchy import linkage, leaves_list -import plotly.figure_factory as ff -import re import warnings warnings.filterwarnings( @@ -89,7 +61,6 @@ GIANA stands out from other TCR analysis tools for several key reasons: import pandas as pd import plotly.graph_objects as go -import os from IPython.display import display, Markdown def plot_giana_tables_per_patient(meta_df, giana_dir, subject_col, timepoint_col, min_unique_thresholds=[0, 1]): @@ -190,7 +161,7 @@ def plot_giana_tables_per_patient(meta_df, giana_dir, subject_col, timepoint_col # Aggregate cluster stats across all samples for this patient cluster_stats = df.groupby('cluster').agg( - Total_Counts=('counts', 'sum'), + Total_Counts=('duplicate_count', 'sum'), Unique_TCRs=('CDR3b', 'nunique') ).reset_index() @@ -261,7 +232,7 @@ for patient in meta[subject_col].unique(): # 1. Aggregate cluster stats (All clusters, no min_unique threshold) cluster_stats = df.groupby('cluster').agg( - Total_Counts=('counts', 'sum'), + Total_Counts=('duplicate_count', 'sum'), Unique_TCRs=('CDR3b', 'nunique') ).reset_index() @@ -310,8 +281,6 @@ However, simply identifying the existence of these top clusters is only the firs import os import pandas as pd import plotly.graph_objects as go -import matplotlib.pyplot as plt -from upsetplot import from_contents, plot as upset_plot import warnings import itertools @@ -350,7 +319,7 @@ giana_df = giana_df.dropna(subset=[subject_col]) # Level 2: Cluster Level (Innermost Ring) df_cluster = giana_df.groupby([alias_col, 'cluster']).agg( - total_counts=('counts', 'sum'), + total_counts=('duplicate_count', 'sum'), n_unique_tcrs=('CDR3b', 'nunique') ).reset_index() @@ -360,7 +329,7 @@ df_cluster['label'] = df_cluster['cluster'].astype(str) # Level 1: Alias Level (Middle Ring) df_sample = giana_df.groupby(alias_col).agg( - total_counts=('counts', 'sum'), + total_counts=('duplicate_count', 'sum'), n_unique_tcrs=('CDR3b', 'nunique') ).reset_index() @@ -369,7 +338,7 @@ df_sample['parent'] = 'All Samples' df_sample['label'] = df_sample[alias_col] # Level 0: Root Level (Center) -total_counts_root = giana_df['counts'].sum() +total_counts_root = giana_df['duplicate_count'].sum() unique_tcrs_root = giana_df['CDR3b'].nunique() df_root = pd.DataFrame({ @@ -418,7 +387,7 @@ fig.show() # Level 2: Alias Level (Outer Ring) df_alias_level = (giana_df.groupby(['cluster', alias_col]) - .agg(total_counts=('counts', 'sum'), n_unique_tcrs=('CDR3b', 'nunique')) + .agg(total_counts=('duplicate_count', 'sum'), n_unique_tcrs=('CDR3b', 'nunique')) .reset_index() .assign(id=lambda x: x['cluster'].astype(str) + '_' + x[alias_col], parent=lambda x: x['cluster'].astype(str), @@ -426,7 +395,7 @@ df_alias_level = (giana_df.groupby(['cluster', alias_col]) # Level 1: Cluster Level (Middle Ring) df_cluster_level = (giana_df.groupby('cluster') - .agg(total_counts=('counts', 'sum'), n_unique_tcrs=('CDR3b', 'nunique')) + .agg(total_counts=('duplicate_count', 'sum'), n_unique_tcrs=('CDR3b', 'nunique')) .reset_index() .assign(id=lambda x: x['cluster'].astype(str), parent='All Clusters', @@ -437,7 +406,7 @@ df_root = pd.DataFrame([{ 'id': 'All Clusters', 'parent': '', 'label': 'All Clusters', - 'total_counts': giana_df['counts'].sum(), + 'total_counts': giana_df['duplicate_count'].sum(), 'n_unique_tcrs': giana_df['CDR3b'].nunique() }]) @@ -565,9 +534,9 @@ def create_patient_tabs_abundance_scatter(df): print(f"## {patient_id}\n") patient_df = df[df[subject_col] == patient_id] - sample_depths = patient_df.groupby(timepoint_col)['counts'].sum() + sample_depths = patient_df.groupby(timepoint_col)['duplicate_count'].sum() patient_pivot = patient_df.pivot_table( - index='cluster', columns=timepoint_col, values='counts', aggfunc='sum' + index='cluster', columns=timepoint_col, values='duplicate_count', aggfunc='sum' ).fillna(0) patient_cpm = patient_pivot.div(sample_depths, axis=1) * 1e6 diff --git a/notebooks/template_gliph.qmd b/notebooks/template_gliph.qmd index 267a9cd..349a4ac 100644 --- a/notebooks/template_gliph.qmd +++ b/notebooks/template_gliph.qmd @@ -22,18 +22,11 @@ gliph2_dir = f"{project_dir}/gliph2/" ```{python} #| include: false -from IPython.display import Image, HTML, display import os -import datetime import pandas as pd import matplotlib.pyplot as plt import plotly.express as px import plotly.graph_objects as go -import igraph as ig -import logomaker -import io -import base64 -import json import warnings # - Loading data @@ -50,10 +43,10 @@ global_similarities_list = [] for patient in meta[subject_col].unique(): p_dir = f"{gliph2_dir}/{patient}" - cd_csv = f"{p_dir}/cluster_member_details.txt" - cn_csv = f"{p_dir}/clone_network.txt" - am_csv = f"{p_dir}/all_motifs.txt" - gs_csv = f"{p_dir}/global_similarities.txt" + cd_csv = f"{p_dir}/{patient}_cluster_member_details.txt" + cn_csv = f"{p_dir}/{patient}_clone_network.txt" + am_csv = f"{p_dir}/{patient}_all_motifs.txt" + gs_csv = f"{p_dir}/{patient}_global_similarities.txt" if os.path.exists(cd_csv): df = pd.read_csv(cd_csv, sep='\t') @@ -102,89 +95,92 @@ import pandas as pd import igraph as ig # --- 1. Prepare Data & Hover Info --- -merged_meta = pd.merge(cluster_details, meta[['sample', alias_col]], on='sample', how='left').dropna(subset=[subject_col]) - -# A. Categorize Nodes (Public vs Private) -cdr3_subject_counts = merged_meta.groupby('CDR3b')[subject_col].nunique() -public_cdr3s = set(cdr3_subject_counts[cdr3_subject_counts > 1].index) - -# B. Create Detailed Hover Text -def create_hover_string(group): - patient_summaries = [] - for subj, subgroup in group.groupby(subject_col): - aliases = ", ".join(sorted(subgroup[alias_col].unique().astype(str))) - patient_summaries.append(f"{subj} ({aliases})") - return "
".join(patient_summaries) - -hover_map = merged_meta.groupby('CDR3b').apply(create_hover_string, include_groups=False).to_dict() - -# C. Create Category Map -private_df = merged_meta[~merged_meta['CDR3b'].isin(public_cdr3s)] -category_map = dict(zip(private_df['CDR3b'], private_df[subject_col])) - -def get_node_category(cdr3): - if cdr3 in public_cdr3s: return "Shared" - return category_map.get(cdr3, "Unknown") - -# --- 2. Build Graph --- -# Isolate just the 4 core GLIPH columns in case extra were appended during concatenation -clone_network = clone_network.iloc[:, :4] -clone_network.columns = ["source", "target", "type", "group"] -clone_network = clone_network[clone_network['type'] != 'singleton'] - -edges = clone_network[["source", "target"]].values.tolist() -g = ig.Graph.TupleList(edges, directed=False) +if cluster_details.empty or clone_network.empty: + print("No GLIPH2 data available (use_gliph2 disabled, or no clusters found for any patient) - skipping network plot.") +else: + merged_meta = pd.merge(cluster_details, meta[['sample', alias_col]], on='sample', how='left').dropna(subset=[subject_col]) + + # A. Categorize Nodes (Public vs Private) + cdr3_subject_counts = merged_meta.groupby('CDR3b')[subject_col].nunique() + public_cdr3s = set(cdr3_subject_counts[cdr3_subject_counts > 1].index) + + # B. Create Detailed Hover Text + def create_hover_string(group): + patient_summaries = [] + for subj, subgroup in group.groupby(subject_col): + aliases = ", ".join(sorted(subgroup[alias_col].unique().astype(str))) + patient_summaries.append(f"{subj} ({aliases})") + return "
".join(patient_summaries) + + hover_map = merged_meta.groupby('CDR3b').apply(create_hover_string, include_groups=False).to_dict() + + # C. Create Category Map + private_df = merged_meta[~merged_meta['CDR3b'].isin(public_cdr3s)] + category_map = dict(zip(private_df['CDR3b'], private_df[subject_col])) + + def get_node_category(cdr3): + if cdr3 in public_cdr3s: return "Shared" + return category_map.get(cdr3, "Unknown") + + # --- 2. Build Graph --- + # Isolate just the 4 core GLIPH columns in case extra were appended during concatenation + clone_network = clone_network.iloc[:, :4] + clone_network.columns = ["source", "target", "type", "group"] + clone_network = clone_network[clone_network['type'] != 'singleton'] + + edges = clone_network[["source", "target"]].values.tolist() + g = ig.Graph.TupleList(edges, directed=False) + + g.vs['category'] = [get_node_category(v['name']) for v in g.vs] + g.vs['hover_detail'] = [hover_map.get(v['name'], "No Data") for v in g.vs] + + layout = g.layout("fr") + coords = list(map(tuple, layout)) + + # --- 3. Plotting --- + fig = go.Figure() + + edge_x, edge_y = [], [] + for edge in g.es: + src, tgt = edge.tuple + x0, y0 = coords[src] + x1, y1 = coords[tgt] + edge_x.extend([x0, x1, None]) + edge_y.extend([y0, y1, None]) -g.vs['category'] = [get_node_category(v['name']) for v in g.vs] -g.vs['hover_detail'] = [hover_map.get(v['name'], "No Data") for v in g.vs] - -layout = g.layout("fr") -coords = list(map(tuple, layout)) - -# --- 3. Plotting --- -fig = go.Figure() - -edge_x, edge_y = [], [] -for edge in g.es: - src, tgt = edge.tuple - x0, y0 = coords[src] - x1, y1 = coords[tgt] - edge_x.extend([x0, x1, None]) - edge_y.extend([y0, y1, None]) - -fig.add_trace(go.Scatter( - x=edge_x, y=edge_y, line=dict(width=0.5, color='#cccccc'), - hoverinfo='none', mode='lines', showlegend=False -)) + fig.add_trace(go.Scatter( + x=edge_x, y=edge_y, line=dict(width=0.5, color='#cccccc'), + hoverinfo='none', mode='lines', showlegend=False + )) -categories = sorted(list(set(g.vs['category']))) -base_colors = px.colors.qualitative.Bold -color_map = {cat: base_colors[i % len(base_colors)] for i, cat in enumerate(categories)} -color_map["Shared"] = "black" + categories = sorted(list(set(g.vs['category']))) + base_colors = px.colors.qualitative.Bold + color_map = {cat: base_colors[i % len(base_colors)] for i, cat in enumerate(categories)} + color_map["Shared"] = "black" -for cat in categories: - indices = [i for i, v in enumerate(g.vs) if v['category'] == cat] - if not indices: continue + for cat in categories: + indices = [i for i, v in enumerate(g.vs) if v['category'] == cat] + if not indices: continue - node_x = [coords[i][0] for i in indices] - node_y = [coords[i][1] for i in indices] - hover_texts = [f"CDR3: {g.vs[i]['name']}
Found in:
{g.vs[i]['hover_detail']}" for i in indices] + node_x = [coords[i][0] for i in indices] + node_y = [coords[i][1] for i in indices] + hover_texts = [f"CDR3: {g.vs[i]['name']}
Found in:
{g.vs[i]['hover_detail']}" for i in indices] - node_size = 18 if cat == "Shared" else 12 - line_width = 2 if cat == "Shared" else 1 + node_size = 18 if cat == "Shared" else 12 + line_width = 2 if cat == "Shared" else 1 - fig.add_trace(go.Scatter( - x=node_x, y=node_y, mode='markers', name=str(cat), text=hover_texts, hoverinfo='text', - marker=dict(color=color_map[cat], size=node_size, line=dict(width=line_width, color='white')) - )) + fig.add_trace(go.Scatter( + x=node_x, y=node_y, mode='markers', name=str(cat), text=hover_texts, hoverinfo='text', + marker=dict(color=color_map[cat], size=node_size, line=dict(width=line_width, color='white')) + )) -fig.update_layout( - title='TCRB Network', plot_bgcolor='white', width=1000, height=800, - xaxis=dict(showgrid=False, zeroline=False, showticklabels=False), - yaxis=dict(showgrid=False, zeroline=False, showticklabels=False), - margin=dict(t=40, b=20, l=10, r=10), legend_title_text='Category' -) -fig.show() + fig.update_layout( + title='TCRB Network', plot_bgcolor='white', width=1000, height=800, + xaxis=dict(showgrid=False, zeroline=False, showticklabels=False), + yaxis=dict(showgrid=False, zeroline=False, showticklabels=False), + margin=dict(t=40, b=20, l=10, r=10), legend_title_text='Category' + ) + fig.show() ``` @@ -255,7 +251,10 @@ def create_patient_histogram_tabs(all_motifs, subject_col): print(":::\n") # Run -create_patient_histogram_tabs(all_motifs, subject_col) +if all_motifs.empty: + print("No GLIPH2 motif data available - skipping.") +else: + create_patient_histogram_tabs(all_motifs, subject_col) ``` @@ -334,7 +333,10 @@ def create_nested_motif_tabs(all_motifs, subject_col): print(":::\n") # Run -create_nested_motif_tabs(all_motifs, subject_col) +if all_motifs.empty: + print("No GLIPH2 motif data available - skipping.") +else: + create_nested_motif_tabs(all_motifs, subject_col) ``` diff --git a/notebooks/template_overlap.qmd b/notebooks/template_overlap.qmd index 2b7961a..2a81cd7 100644 --- a/notebooks/template_overlap.qmd +++ b/notebooks/template_overlap.qmd @@ -22,7 +22,6 @@ concat_csv = f"{project_dir}/annotate/concatenated_cdr3_sorted.tsv" # Load Packages from IPython.display import Image, display, Markdown, HTML from matplotlib.colors import LinearSegmentedColormap -from io import StringIO from scipy.sparse import csr_matrix from scipy.stats import gaussian_kde, fisher_exact from statsmodels.stats.multitest import multipletests @@ -30,7 +29,6 @@ from scipy.cluster.hierarchy import linkage, leaves_list from scipy.spatial.distance import pdist import igraph as ig -import logomaker import base64 import datetime @@ -44,16 +42,12 @@ import matplotlib.ticker as ticker import numpy as np import os import pandas as pd -import pathlib import plotly.express as px import plotly.graph_objects as go import scipy.cluster.hierarchy as sch import seaborn as sns -import shutil import sys -import upsetplot import warnings -from pathlib import Path import re ## Reading sample metadata @@ -142,17 +136,17 @@ from statsmodels.stats.multitest import multipletests clonotypes_df = concat_df.copy() clonotypes_df = clonotypes_df.dropna(subset=[timepoint_col, 'origin']) -clone_cols = ['CDR3b', 'TRBV', 'TRBJ'] +clone_cols = ['junction_aa', 'v_call', 'j_call'] merge_cols = clone_cols + [subject_col, 'origin'] clonotypes_df['timepoint_rank'] = clonotypes_df[timepoint_order_col] clonotypes_df = clonotypes_df.sort_values('timepoint_rank') # Grouping by subject, origin, and timepoint for accurate total counts -sample_total_counts = clonotypes_df.groupby([subject_col, 'origin', timepoint_col])['counts'].sum().reset_index() -sample_total_counts.rename(columns={'counts': 'total_counts'}, inplace=True) +sample_total_counts = clonotypes_df.groupby([subject_col, 'origin', timepoint_col])['duplicate_count'].sum().reset_index() +sample_total_counts.rename(columns={'duplicate_count': 'total_counts'}, inplace=True) clonotypes_df = pd.merge(clonotypes_df, sample_total_counts, on=[subject_col, 'origin', timepoint_col]) -clonotypes_df['frequency'] = clonotypes_df['counts'] / clonotypes_df['total_counts'] +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): @@ -179,8 +173,8 @@ for (subject, origin), subject_df in clonotypes_df.groupby([subject_col, 'origin continue freq_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='frequency', fill_value=0) - counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='counts', fill_value=0) - subject_total_counts = subject_df.groupby(timepoint_col)['counts'].sum() + counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='duplicate_count', fill_value=0) + subject_total_counts = subject_df.groupby(timepoint_col)['duplicate_count'].sum() for t_pre, t_post in itertools.combinations(timepoints, 2): # 1. Vectorized DataFrame creation instead of iterrows @@ -304,8 +298,8 @@ for (subject, origin), subject_df in clonotypes_df.groupby([subject_col, 'origin continue freq_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='frequency', fill_value=0) - counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='counts', fill_value=0) - subject_total_counts = subject_df.groupby(timepoint_col)['counts'].sum() + counts_pivot = subject_df.pivot_table(index=clone_cols, columns=timepoint_col, values='duplicate_count', fill_value=0) + subject_total_counts = subject_df.groupby(timepoint_col)['duplicate_count'].sum() for t_pre, t_post in itertools.combinations(timepoints, 2): # 1. Vectorize @@ -500,7 +494,7 @@ else: sorted_tps += [t for t in valid_tps if t not in sorted_tps] # Set Index to CDR3b - matrix_df = orig_df.set_index('CDR3b')[sorted_tps] + matrix_df = orig_df.set_index('junction_aa')[sorted_tps] # Drop completely empty rows matrix_df = matrix_df.loc[~(matrix_df == 0).all(axis=1)] @@ -708,7 +702,7 @@ def get_top_ids(df, direction): subset = group[group['fold_change'] < 1] # Sort by P-value (lowest first) -> Take top 15 Unique CDR3s - top_clones = subset.sort_values('p_value', ascending=True)['CDR3b'].drop_duplicates().head(15).tolist() + top_clones = subset.sort_values('p_value', ascending=True)['junction_aa'].drop_duplicates().head(15).tolist() id_map[name] = set(top_clones) return id_map @@ -718,7 +712,7 @@ contracted_map = get_top_ids(stats_df, 'down') # --- 2. Prepare Master Frequency Data --- # Aggregate by CDR3 including 'origin' so we can split traces later -agg_cols = [subject_col, 'origin', 'CDR3b', timepoint_col] if 'origin' in clonotypes_df.columns else [subject_col, 'CDR3b', timepoint_col] +agg_cols = [subject_col, 'origin', 'junction_aa', timepoint_col] if 'origin' in clonotypes_df.columns else [subject_col, 'junction_aa', timepoint_col] cdr3_df = clonotypes_df.groupby(agg_cols)['frequency'].sum().reset_index() # Log Transform (with pseudocount) @@ -748,7 +742,7 @@ def generate_plot_quarto_tabs(id_map, title_prefix, line_color): def get_status(row): subj = row[subject_col] orig = row['origin'] if 'origin' in row else None - seq = row['CDR3b'] + seq = row['junction_aa'] # Check map using the correct key format key = (subj, orig) if orig is not None else subj @@ -801,7 +795,7 @@ def generate_plot_quarto_tabs(id_map, title_prefix, line_color): # Generate temporary figure to extract traces temp_fig = px.line( orig_df, x=timepoint_col, y="log_freq", color="status", - line_group="CDR3b", hover_name="CDR3b", + line_group="junction_aa", hover_name="junction_aa", color_discrete_map=color_map, category_orders={timepoint_col: orig_tps} ) diff --git a/notebooks/template_pheno_bulk.qmd b/notebooks/template_pheno_bulk.qmd index 4970300..a17d883 100644 --- a/notebooks/template_pheno_bulk.qmd +++ b/notebooks/template_pheno_bulk.qmd @@ -148,9 +148,9 @@ def create_phenotype_tabs(meta_df): # --- Merge Metadata --- meta_row = meta_df[meta_df['sample'] == sample_id] if not meta_row.empty: - composition['subject_id'] = meta_row.iloc[0]['subject_id'] - composition['alias'] = meta_row.iloc[0]['alias'] - composition['sort_order'] = meta_row.iloc[0]['timepoint_order'] if 'timepoint_order' in meta_df.columns else 0 + composition['subject_id'] = meta_row.iloc[0][subject_col] + composition['alias'] = meta_row.iloc[0][alias_col] + composition['sort_order'] = meta_row.iloc[0][timepoint_order_col] if timepoint_order_col in meta_df.columns else 0 else: composition['subject_id'] = np.nan composition['alias'] = sample_id @@ -321,10 +321,10 @@ def create_weighted_phenotype_plot(meta_df): # Skip samples not in metadata, or define defaults continue - subject_id = meta_row.iloc[0]['subject_id'] - alias = meta_row.iloc[0]['alias'] + subject_id = meta_row.iloc[0][subject_col] + alias = meta_row.iloc[0][alias_col] # Default to 0 if order column is missing - tp_order = meta_row.iloc[0]['timepoint_order'] if 'timepoint_order' in meta_df.columns else 0 + tp_order = meta_row.iloc[0][timepoint_order_col] if timepoint_order_col in meta_df.columns else 0 # --- Load Data --- pheno_df = pd.read_csv(pheno_file_path, sep='\t') @@ -339,10 +339,12 @@ def create_weighted_phenotype_plot(meta_df): continue # --- Merge & Weight --- - # Match phenotype (junction_aa) with counts (CDR3b) - merged_df = pd.merge(pheno_df, sample_counts_df[['CDR3b', 'counts']], - left_on='junction_aa', right_on='CDR3b', how='inner') - + # TCRPHENO output only carries 'sequence_id' (no CDR3/junction_aa column), + # which matches the same "{sample}|sequence{N}" IDs in concat_csv. + merged_df = pd.merge(pheno_df, sample_counts_df[['sequence_id', 'duplicate_count']], + on='sequence_id', how='inner') + merged_df.rename(columns={'duplicate_count': 'counts'}, inplace=True) + merged_df = merged_df.dropna(subset=phenotype_cols + ['counts']) if merged_df.empty: diff --git a/notebooks/template_pheno_sc.qmd b/notebooks/template_pheno_sc.qmd index eaf8f6c..2f18455 100644 --- a/notebooks/template_pheno_sc.qmd +++ b/notebooks/template_pheno_sc.qmd @@ -5,11 +5,10 @@ project_dir=f"{project_dir}" -# Define files +# Define files concat_csv = f"{project_dir}/annotate/concatenated_cdr3_sorted.tsv" -concat_csv_pheno=f"{project_dir}/compare_phenotype/concatenated_cdr3.tsv" -sample_table_pheno=f"{project_dir}/pipeline_info/samplesheet_phenotype.csv" +concat_csv_pheno = f"{project_dir}/annotate/concatenated_cdr3_sorted.tsv" # Define dirs tcrpheno_dir = f"{project_dir}/tcrpheno/" @@ -46,7 +45,6 @@ from scipy import stats from itertools import combinations from scipy.stats import mannwhitneyu from plotly.subplots import make_subplots -from pathlib import Path import re meta = pd.read_csv(sample_table, sep=',') @@ -58,14 +56,8 @@ concat_df = concat_df.merge(meta[['sample', subject_col, alias_col, timepoint_co # Phenotype files -# Samplesheet -metadata_pheno_df = pd.read_csv(sample_table_pheno, sep=',', header=0, index_col="file") -metadata_pheno_df.index = metadata_pheno_df.index.map(lambda p: Path(p).name) -metadata_pheno_df.index = metadata_pheno_df.index.str.replace('_pseudobulk', '', regex=False) -metadata_pheno_df.index = metadata_pheno_df.index.str.replace('_phenotype.tsv', '', regex=False) - # Importing sample metadata (pheno-ps) -clonotypes_pheno_df = pd.read_csv(concat_csv_pheno, sep='\t', header=0, index_col=0).reset_index() +clonotypes_pheno_df = pd.read_csv(concat_csv_pheno, sep='\t', header=0) clonotypes_pheno_df['sample_phenotype'] = clonotypes_pheno_df['sample'] # Add 'sample_phenotype' and 'phenotype' samples = sorted(meta['sample'].astype(str).unique(), key=len, reverse=True) @@ -437,13 +429,13 @@ import sys warnings.filterwarnings("ignore") def create_patient_upset_figure(df, patient_id): - df_clean = df.dropna(subset=['phenotype', 'CDR3b']).copy() + df_clean = df.dropna(subset=['phenotype', 'junction_aa']).copy() phenotypes = sorted(df_clean['phenotype'].unique()) if len(phenotypes) < 2: return None clones_by_pheno = { - p: list(df_clean[df_clean['phenotype'] == p]['CDR3b'].unique()) + p: list(df_clean[df_clean['phenotype'] == p]['junction_aa'].unique()) for p in phenotypes } upset_data = from_contents(clones_by_pheno) diff --git a/notebooks/template_sample.qmd b/notebooks/template_sample.qmd index aadde37..e49837f 100644 --- a/notebooks/template_sample.qmd +++ b/notebooks/template_sample.qmd @@ -27,27 +27,14 @@ convergence_dir = f"{project_dir}/convergence/" ```{python} #| code-fold: true -import datetime import glob -import itertools -import math -import os -import sys import h5py -import igraph as ig -import matplotlib.pyplot as plt -import matplotlib.ticker as ticker +import os import numpy as np import pandas as pd import plotly.express as px import plotly.graph_objects as go -import scipy.cluster.hierarchy as sch -import seaborn as sns -from IPython.display import Image -from matplotlib.colors import LinearSegmentedColormap from scipy.sparse import csr_matrix -from scipy.stats import gaussian_kde -import plotly.io as pio from scipy import stats from itertools import combinations from scipy.stats import mannwhitneyu @@ -104,7 +91,6 @@ $$f_i = \frac{\text{Read count for clone } i}{\text{Total reads for all producti #| echo: false import pandas as pd -import plotly.graph_objects as go def create_expansion_stacked_barplot_per_individual(df, subject_col, origin_col, timepoint_order_col): """ @@ -112,11 +98,11 @@ def create_expansion_stacked_barplot_per_individual(df, subject_col, origin_col, Generates one Quarto tab per individual with dropdowns for origin. Includes a horizontal threshold line for sequencing depth (3000 counts). """ - # --- Ensure 'counts' column is numeric --- + # --- Ensure 'duplicate_count' column is numeric --- try: - df['counts'] = pd.to_numeric(df['counts']) + df['duplicate_count'] = pd.to_numeric(df['duplicate_count']) except ValueError: - print("Error: 'counts' column could not be converted to a numeric type.") + print("Error: 'duplicate_count' column could not be converted to a numeric type.") return subjects = [s for s in df[subject_col].unique() if pd.notna(s)] @@ -138,11 +124,11 @@ def create_expansion_stacked_barplot_per_individual(df, subject_col, origin_col, for sample_name in unique_samples: sample_df = origin_df[origin_df['sample'] == sample_name].copy() - total_reads = sample_df['counts'].sum() + total_reads = sample_df['duplicate_count'].sum() if total_reads == 0: continue - sample_df['frequency'] = (sample_df['counts'] / total_reads) * 100 + sample_df['frequency'] = (sample_df['duplicate_count'] / total_reads) * 100 highly_expanded = sample_df[sample_df['frequency'] > 1] expanded = sample_df[(sample_df['frequency'] >= 0.1) & (sample_df['frequency'] <= 1)] @@ -360,7 +346,6 @@ This table is useful because it reveals the identities of the dominant clones dr #| output: asis #| echo: false -import plotly.graph_objects as go def create_top_clones_tables_by_patient(df, n_clones=15): """ @@ -393,21 +378,21 @@ def create_top_clones_tables_by_patient(df, n_clones=15): # Build a table trace for each timepoint for i, tp in enumerate(unique_timepoints): tp_df = pat_df[pat_df[timepoint_col] == tp].copy() - tp_df = tp_df[tp_df['counts'] > 0] + tp_df = tp_df[tp_df['duplicate_count'] > 0] if tp_df.empty: continue # Grab the top clones - top_n_clones = tp_df.nlargest(n_clones, 'counts') + top_n_clones = tp_df.nlargest(n_clones, 'duplicate_count') # Safely extract columns (avoids KeyErrors if TRBV/TRBJ are missing in some runs) - cols_to_show = ['CDR3b', 'TRBV', 'TRBJ', 'counts'] + cols_to_show = ['junction_aa', 'v_call', 'j_call', 'duplicate_count'] available_cols = [c for c in cols_to_show if c in top_n_clones.columns] display_df = top_n_clones[available_cols].copy() # Clean up headers for display - display_df.columns = [c.capitalize() if c == 'counts' else c for c in display_df.columns] + display_df.columns = [c.capitalize() if c == 'duplicate_count' else c for c in display_df.columns] fig.add_trace( go.Table( @@ -489,8 +474,6 @@ create_top_clones_tables_by_patient(concat_df, n_clones=15) import pandas as pd import numpy as np import plotly.graph_objects as go -from plotly.subplots import make_subplots -from scipy.stats import mannwhitneyu def create_violin_with_heatmap_inset(df, meta_df, origin_col='origin'): """ @@ -503,16 +486,16 @@ def create_violin_with_heatmap_inset(df, meta_df, origin_col='origin'): # --- 1. Data Prep --- try: - df['counts'] = pd.to_numeric(df['counts']) + df['duplicate_count'] = pd.to_numeric(df['duplicate_count']) except ValueError: - print("Error: 'counts' column could not be converted to a numeric type.") + print("Error: 'duplicate_count' column could not be converted to a numeric type.") return # Filter for valid counts - plot_df = df[df['counts'] > 0].copy() + plot_df = df[df['duplicate_count'] > 0].copy() if plot_df.empty: return - plot_df['log10_counts'] = np.log10(plot_df['counts']) + plot_df['log10_counts'] = np.log10(plot_df['duplicate_count']) # --- 2. Merge with Metadata --- # Assuming df already has the metadata merged based on your previous code structure @@ -958,9 +941,11 @@ conv_data = [] for f in files: filename = os.path.basename(f) sample_id = filename.replace('_tcr_convergence.tsv', '') - if sample_id.endswith('_pseudobulk'): + if sample_id.endswith('_pseudobulk'): sample_id = sample_id[:-11] - if sample_id.endswith('_airr'): + if sample_id.endswith('_airr'): + sample_id = sample_id[:-5] + if sample_id.endswith('_cdr3'): sample_id = sample_id[:-5] try: @@ -1098,7 +1083,6 @@ print(":::\n") import plotly.express as px import plotly.graph_objects as go -from scipy.stats import mannwhitneyu from itertools import combinations import pandas as pd @@ -1574,8 +1558,6 @@ import glob import pandas as pd import numpy as np import plotly.graph_objects as go -from plotly.subplots import make_subplots -from scipy.stats import mannwhitneyu def create_pgen_tabs(meta_df, origin_col='origin'): """ @@ -2331,7 +2313,7 @@ def plot_vdjdb_species_specificity(VDJdb_dir, meta_df, subject_col='subject_id') print(":::\n") # Run the function -plot_vdjdb_species_specificity(VDJdb_dir, meta) +plot_vdjdb_species_specificity(VDJdb_dir, meta, subject_col=subject_col) ``` **Figure 14: Top 15 recognized epitopes identified in bulk-TCR data.** The number of unique TCR clonotypes in the sample that match a known epitope in the VDJdb database is shown on the x-axis. The plot highlights the most abundant antigen specificities inferred from the T-cell repertoire. @@ -2515,7 +2497,7 @@ def plot_vdjdb_sunburst_hierarchy(VDJdb_dir, meta_df, subject_col='subject_id'): print(":::\n") # Run the function -plot_vdjdb_sunburst_hierarchy(VDJdb_dir, meta) +plot_vdjdb_sunburst_hierarchy(VDJdb_dir, meta, subject_col=subject_col) ``` **Figure 15: TCR Repertoire Mapping to specific gene epitopes:** Sunburst chart visualizing the relationships between T-cell receptor (TCR) sequences and their known antigen specificities from the VDJdb curated database. The central ring represents the species, branching out to their specific genes, and finally to the individual TCR sequences from your dataset that match those specificities. The size of each segment corresponds to the number of TCRs associated with that particular specificity. @@ -2691,7 +2673,7 @@ def plot_vdjdb_sunburst_reversed(VDJdb_dir, meta_df, subject_col='subject_id'): print(":::\n") # Run the function -plot_vdjdb_sunburst_reversed(VDJdb_dir, meta) +plot_vdjdb_sunburst_reversed(VDJdb_dir, meta, subject_col=subject_col) ``` **Figure 16: Gene Epitopes Mapping to TCR sequences:** Sunburst chart visualizing the relationships antigens and T-cell receptor (TCR) sequences. The central ring represents TCRs sequences, branching out to their specific genes, and finally to the species the antigen belongs to. Only TCRs with a Freq>2 are shown, for visualization purposes. diff --git a/notebooks/template_sharing.qmd b/notebooks/template_sharing.qmd index 0bd0e6b..11dd0ef 100644 --- a/notebooks/template_sharing.qmd +++ b/notebooks/template_sharing.qmd @@ -48,7 +48,6 @@ import base64 import json from IPython.display import HTML, display import warnings -import matplotlib.pyplot as plt from upsetplot import from_contents, plot as upset_plot @@ -81,16 +80,16 @@ This plot is generated by first calculating the frequency of each unique TCR clo #| fig-cap: "**TCR Sharing by Maximum Clonal Expansion.** Number of samples where a TCR has an exact match on the aminoacid level. Color represents clonal expansion category using the highest frequency across all samples." # --- 1. Load and Process Actual Data --- -# We assume 'concat_df' has columns ['CDR3b', 'counts', 'sample'] +# We assume 'concat_df' has columns ['junction_aa', 'duplicate_count', 'sample'] # We assume 'meta' has columns ['sample', subject_col] raw_df = concat_df.copy() # Rename columns to standard internal names -raw_df.rename(columns={'CDR3b': 'cdr3_sequence', 'sample': 'sample_id'}, inplace=True) +raw_df.rename(columns={'junction_aa': 'cdr3_sequence', 'sample': 'sample_id'}, inplace=True) # Calculate frequency per SAMPLE first (Expansion is a property of a specific physical sample) -sample_total_counts = raw_df.groupby('sample_id')['counts'].transform('sum') -raw_df['frequency'] = raw_df['counts'] / sample_total_counts +sample_total_counts = raw_df.groupby('sample_id')['duplicate_count'].transform('sum') +raw_df['frequency'] = raw_df['duplicate_count'] / sample_total_counts # --- CRITICAL STEP: Merge with Metadata to get Patient IDs --- # We merge on 'sample_id' (which matches 'sample' in meta) @@ -196,14 +195,14 @@ import plotly.express as px merged_df = pd.merge(concat_df, meta[['sample', subject_col, timepoint_col]], on='sample', how='left') # Group by TCR (CDR3b) and aggregate the exact data we need for plotting and hovering -agg_df = merged_df.groupby('CDR3b').agg( +agg_df = merged_df.groupby('junction_aa').agg( total_patients=(subject_col, 'nunique'), individuals=(subject_col, lambda x: ', '.join(sorted(set(x.dropna().astype(str))))), timepoints=(timepoint_col, lambda x: ', '.join(sorted(set(x.dropna().astype(str))))) ).reset_index() # Merge this with your existing generation probability dataframe -plot_df = pd.merge(prob_generation_df[['CDR3b', 'pgen']], agg_df, on='CDR3b', how='inner') +plot_df = pd.merge(prob_generation_df[['junction_aa', 'pgen']], agg_df, on='junction_aa', how='inner') # Calculate the log10 of the generation probability plot_df['log10_pgen'] = np.log10(plot_df['pgen']) @@ -215,7 +214,7 @@ fig = px.scatter( plot_df, x='log10_pgen', y='total_patients', - hover_name='CDR3b', # Puts the sequence at the top of the tooltip in bold + hover_name='junction_aa', # Puts the sequence at the top of the tooltip in bold hover_data={ 'log10_pgen': ':.2f', # Format to 2 decimal places 'total_patients': False, # Hide this because it's already obvious from the Y-axis diff --git a/subworkflows/local/compare.nf b/subworkflows/local/compare.nf index 9e2453e..0265bf6 100644 --- a/subworkflows/local/compare.nf +++ b/subworkflows/local/compare.nf @@ -33,4 +33,7 @@ workflow COMPARE { TCRSHARING_SCATTERPLOT( TCRSHARING_CALC.out.shared_cdr3 ) + + emit: + shared_cdr3 = TCRSHARING_CALC.out.shared_cdr3 } \ No newline at end of file diff --git a/subworkflows/local/patient.nf b/subworkflows/local/patient.nf index aa93022..423a920 100644 --- a/subworkflows/local/patient.nf +++ b/subworkflows/local/patient.nf @@ -45,9 +45,27 @@ workflow PATIENT { params.threshold_vgene ) - if(params.use_gliph2) { + // Each gliph2_* emit is a list of [patient, file] pairs - kept separate per + // output type (rather than mixed together) so downstream staging can map + // each pair back to its known target leaf-name (e.g. "all_motifs.txt") + // without having to parse it back out of the patient-prefixed filename. + // + // NOTE: this Nextflow version's strict-syntax workflow output collection + // requires emit values to be direct .out expressions - referencing a + // pre-computed local `def` variable in `emit:` fails at definition time + // with "Missing workflow output parameter", even outside any conditional. + // So the params.use_gliph2 ternary has to live in the emit line itself, + // referencing GLIPH2_TURBOGLIPH.out directly. + if (params.use_gliph2) { GLIPH2_TURBOGLIPH( PATIENT_CONCATENATE.out.patient_cdr3 ) } + + emit: + giana_files = GIANA_CALC.out.giana_output.collect() + gliph2_all_motifs = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.all_motifs.collect() : channel.value([]) + gliph2_clone_network = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.clone_network.collect() : channel.value([]) + gliph2_cluster_member_details = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.cluster_member_details.collect() : channel.value([]) + gliph2_global_similarities = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.global_similarities.collect() : channel.value([]) } \ No newline at end of file diff --git a/subworkflows/local/sample.nf b/subworkflows/local/sample.nf index 28f7b67..a3f53fe 100644 --- a/subworkflows/local/sample.nf +++ b/subworkflows/local/sample.nf @@ -111,5 +111,18 @@ workflow SAMPLE { VDJDB_VDJMATCH (processed_samples, VDJDB_GET.out.ref_db) emit: - sample_csv = SAMPLE_CALC.out.sample_csv + sample_csv = SAMPLE_CALC.out.sample_csv + v_family = SAMPLE_CALC_PIVOT.out.v_family_wide + j_family = SAMPLE_CALC_PIVOT.out.j_family_wide + tcrdist_files = TCRDIST3_MATRIX.out.tcrdist_output.map { _sample_meta, dist -> dist } + .mix( TCRDIST3_MATRIX.out.clone_df ) + .flatten() + .collect() + olga_files = OLGA_SAMPLE_MERGE.out.olga_pgen.map { _sample_meta, pgen -> pgen } + .collect() + vdjdb_files = VDJDB_VDJMATCH.out.vdjmatch_txt + .mix( VDJDB_VDJMATCH.out.annot_summary ) + .collect() + convergence_files = CONVERGENCE.out.convergence_output.collect() + tcrpheno_files = TCRPHENO.out.tcrpheno_output.map { _sample_meta, f -> f }.collect() } \ No newline at end of file diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index eba24f6..13b4931 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -21,7 +21,83 @@ nextflow_process { """ input[0] = [ file("${projectDir}/tests/fixtures/test.qmd"), - file("${projectDir}/tests/fixtures") + file("${projectDir}/tests/fixtures"), + [] + ] + input[1] = "TCRtoolkit" + input[2] = "nextflow run main.nf" + """ + } + } + + then { + assert process.success + with(process.out.report_html) { + def html_file = path(get(0)) + assert html_file.exists() + assert html_file.getFileName().toString().endsWith('.html') + } + } + } + + test("Should stage files into a nested project_dir layout") { + // Verifies the staged_layout mechanism used by notebooks that read from a + // project_dir// tree (e.g. template_discovery_brief.qmd) + // instead of bare relative filenames. + + tag "notebook" + tag "staged_layout" + + when { + params { + samplesheet = "${projectDir}/tests/fixtures/valid_samplesheet.csv" + container = System.getenv("CONTAINER") ?: "ghcr.io/karchinlab/tcrtoolkit:main" + } + + process { + """ + input[0] = [ + file("${projectDir}/tests/fixtures/test.qmd"), + [file("${projectDir}/tests/fixtures/valid_samplesheet.csv")], + [["myproj/sample/valid_samplesheet.csv", "valid_samplesheet.csv"]] + ] + input[1] = "TCRtoolkit" + input[2] = "nextflow run main.nf" + """ + } + } + + then { + assert process.success + with(process.out.report_html) { + def html_file = path(get(0)) + assert html_file.exists() + assert html_file.getFileName().toString().endsWith('.html') + } + } + } + + test("Should stage a file under a renamed dest basename") { + // Verifies source and dest basenames can differ - e.g. per-patient files + // staged under a shared/prefixed basename (like GLIPH2's + // "patientA_all_motifs.txt") get symlinked to an unprefixed dest name + // ("gliph2/patientA/all_motifs.txt"). + + tag "notebook" + tag "staged_layout" + + when { + params { + samplesheet = "${projectDir}/tests/fixtures/valid_samplesheet.csv" + container = System.getenv("CONTAINER") ?: "ghcr.io/karchinlab/tcrtoolkit:main" + } + + process { + """ + input[0] = [ + file("${projectDir}/tests/fixtures/test.qmd"), + [file("${projectDir}/tests/fixtures/valid_samplesheet.csv")], + [["myproj/gliph2/patientA/all_motifs.txt", "valid_samplesheet.csv"]] ] input[1] = "TCRtoolkit" input[2] = "nextflow run main.nf" diff --git a/workflows/tcrtoolkit.nf b/workflows/tcrtoolkit.nf index f962d0b..32d4843 100644 --- a/workflows/tcrtoolkit.nf +++ b/workflows/tcrtoolkit.nf @@ -107,24 +107,203 @@ workflow TCRTOOLKIT { ) } - // Report - works on channel of tuples [report name, [report files]] + // Report - works on channel of tuples [notebook template, files to stage, staged directory layout] ch_reports = channel.empty() - // QC report requires sample-level aggregate outputs. - ch_qc_report = SAMPLE.out.sample_csv + def sample_stats_agg = SAMPLE.out.sample_csv .collectFile(name: "sample_stats.csv", keepHeader: true, skip: 1, sort: true) + + // template_discovery_brief.qmd's VDJdb section falls back to reading raw + // (unconverted) input files when convert_dir is empty, which only has AIRR- + // standard frequency columns for input_format 'airr' - for 'adaptive'/ + // 'cellranger' it needs the AIRR-converted files, so those are only staged + // when CONVERT actually ran (mirrors the same input_format branch above). + def convert_files = (input_format == 'adaptive' || input_format == 'cellranger') + ? CONVERT.out.sample_map_converted.map { _meta, f -> f }.collect() + : channel.value([]) + + // template_discovery_brief.qmd includes a single generic template_pheno.qmd - + // RENDER_NOTEBOOK stages whichever real notebook applies under that shared + // name (see staged_layout below), since Quarto's {{< include >}} shortcode is + // a static textual splice with no native runtime if/else. pheno_sc needs the + // per-cell/per-phenotype pseudobulk files, which only exist for cellranger + // input with sobject_gex supplied; everything else gets pheno_bulk, which + // only needs TCRPHENO output (produced for every sample regardless of format). + def use_pheno_sc = (input_format == 'cellranger' && params.sobject_gex) + def pheno_notebook = use_pheno_sc ? file(params.template_pheno_sc) : file(params.template_pheno_bulk) + // .collect() on a channel that never emits (channel.empty(), which is what + // CONVERT.out.pseudobulk_phenotype_files is outside the cellranger+sobject_gex + // case) never emits either - not even an empty list - which would silently + // starve every downstream .combine() in ch_discovery_report and make + // RENDER_NOTEBOOK(template_discovery_brief) never get scheduled at all (no + // error, just silently skipped - confirmed directly with a minimal repro). + // Guard with the same channel.value([]) fallback already used for convert_files. + def pseudobulk_pheno_files = use_pheno_sc + ? CONVERT.out.pseudobulk_phenotype_files.map { _meta, files -> files }.flatten().collect() + : channel.value([]) + + // QC report requires sample-level aggregate outputs. + ch_qc_report = sample_stats_agg .combine(ANNOTATE.out.concat_cdr3_sorted) .map { sample_stats_csv, concat_cdr3_sorted -> tuple( file(params.template_qc), [sample_stats_csv, - concat_cdr3_sorted] + concat_cdr3_sorted], + [] ) } ch_reports = ch_reports.mix(ch_qc_report) - // Another report - // ch_reports = ch_reports.mix( ... ) + // Discovery brief requires sample-level aggregate outputs, TCRdist3/VDJdb + // per-sample outputs, and cross-sample TCR sharing - it reads them from a + // project_dir/// layout, so staged_layout tells + // RENDER_NOTEBOOK where to symlink each staged file within the task work dir. + // Each staged_layout entry is a [dest_path, source_basename] pair - source and + // dest basenames match here, but don't always (see gliph2 below). + // + // NOTE: .collect()-produced list channels are wrapped in an extra list via + // `.map { l -> [l] }` before every .combine() below - .combine() otherwise + // flattens a list-valued item's contents into the tuple instead of keeping + // it as one element (confirmed directly: combining a scalar channel with an + // unwrapped .collect() channel of 3 items produces a 4-element flat tuple, + // not a 2-element [scalar, list] tuple). + ch_discovery_report = sample_stats_agg + .combine(ANNOTATE.out.concat_cdr3_sorted) + .combine(COMPARE.out.shared_cdr3) + .combine(SAMPLE.out.tcrdist_files.map { l -> [l] }) + .combine(SAMPLE.out.vdjdb_files.map { l -> [l] }) + .combine(convert_files.map { l -> [l] }) + .combine(SAMPLE.out.tcrpheno_files.map { l -> [l] }) + .combine(pseudobulk_pheno_files.map { l -> [l] }) + .map { sample_stats_csv, concat_cdr3_sorted, shared_cdr3, tcrdist_files, vdjdb_files, convert_files_l, tcrpheno_files, pseudobulk_files -> + def report_files = [sample_stats_csv, concat_cdr3_sorted, shared_cdr3, pheno_notebook] + + tcrdist_files + vdjdb_files + convert_files_l + tcrpheno_files + pseudobulk_files + def staged_layout = [ + ["${params.project_name}/sample/${sample_stats_csv.name}", sample_stats_csv.name], + ["${params.project_name}/annotate/${concat_cdr3_sorted.name}", concat_cdr3_sorted.name], + ["${params.project_name}/tcrsharing/${shared_cdr3.name}", shared_cdr3.name], + ["template_pheno.qmd", pheno_notebook.name] + ] + tcrdist_files.collect { f -> ["${params.project_name}/tcrdist3/${f.name}", f.name] } + + vdjdb_files.collect { f -> ["${params.project_name}/vdjdb/${f.name}", f.name] } + + convert_files_l.collect { f -> ["${params.project_name}/convert/${f.name}", f.name] } + + tcrpheno_files.collect { f -> ["${params.project_name}/tcrpheno/${f.name}", f.name] } + + pseudobulk_files.collect { f -> ["${params.project_name}/pseudobulk/${f.name}", f.name] } + tuple( + file(params.template_discovery_brief), + report_files, + staged_layout + ) + } + ch_reports = ch_reports.mix(ch_discovery_report) + + // Details report (part 1) includes template_sample.qmd, which reads sample-level + // aggregate outputs, the V/J gene-family pivot tables, and per-sample TCRdist3, + // OLGA, VDJdb, and convergence outputs from the same project_dir// + // / layout as the discovery brief. + ch_details_part1_report = sample_stats_agg + .combine(ANNOTATE.out.concat_cdr3_sorted) + .combine(SAMPLE.out.v_family) + .combine(SAMPLE.out.j_family) + .combine(SAMPLE.out.tcrdist_files.map { l -> [l] }) + .combine(SAMPLE.out.olga_files.map { l -> [l] }) + .combine(SAMPLE.out.vdjdb_files.map { l -> [l] }) + .combine(SAMPLE.out.convergence_files.map { l -> [l] }) + .map { sample_stats_csv, concat_cdr3_sorted, v_family, j_family, tcrdist_files, olga_files, vdjdb_files, convergence_files -> + // template_details_part1.qmd pulls in template_sample.qmd via a + // {{< include >}} shortcode resolved relative to its own directory, + // so that sibling file has to be staged alongside it too - Nextflow + // only stages the single notebook file named in the tuple otherwise. + def include_files = [file("${file(params.template_details_part1).parent}/template_sample.qmd")] + def report_files = [sample_stats_csv, concat_cdr3_sorted, v_family, j_family] + + tcrdist_files + olga_files + vdjdb_files + convergence_files + include_files + def staged_layout = [ + ["${params.project_name}/sample/${sample_stats_csv.name}", sample_stats_csv.name], + ["${params.project_name}/annotate/${concat_cdr3_sorted.name}", concat_cdr3_sorted.name], + ["${params.project_name}/sample/${v_family.name}", v_family.name], + ["${params.project_name}/sample/${j_family.name}", j_family.name] + ] + tcrdist_files.collect { f -> ["${params.project_name}/tcrdist3/${f.name}", f.name] } + + olga_files.collect { f -> ["${params.project_name}/olga/${f.name}", f.name] } + + vdjdb_files.collect { f -> ["${params.project_name}/vdjdb/${f.name}", f.name] } + + convergence_files.collect { f -> ["${params.project_name}/convergence/${f.name}", f.name] } + tuple( + file(params.template_details_part1), + report_files, + staged_layout + ) + } + ch_reports = ch_reports.mix(ch_details_part1_report) + + // Details report (part 2) includes template_overlap.qmd, template_sharing.qmd, + // template_giana.qmd, and template_gliph.qmd. The latter two need patient-level + // clonotype clustering outputs, which only exist when 'patient' is in + // workflow_level; their own os.path.exists guards degrade gracefully when those + // files are absent, so those files are only staged when PATIENT actually ran. + def details_part2_base = sample_stats_agg + .combine(ANNOTATE.out.concat_cdr3_sorted) + .combine(COMPARE.out.shared_cdr3) + + // template_details_part2.qmd pulls in these 4 files via {{< include >}} + // shortcodes resolved relative to its own directory, so they have to be + // staged alongside it too - see the equivalent part1 comment above. + def part2_notebooks_dir = file(params.template_details_part2).parent + def part2_include_files = [ + file("${part2_notebooks_dir}/template_overlap.qmd"), + file("${part2_notebooks_dir}/template_sharing.qmd"), + file("${part2_notebooks_dir}/template_giana.qmd"), + file("${part2_notebooks_dir}/template_gliph.qmd") + ] + + if (levels.contains('patient')) { + ch_details_part2_report = details_part2_base + .combine(PATIENT.out.giana_files.map { l -> [l] }) + .combine(PATIENT.out.gliph2_all_motifs.map { l -> [l] }) + .combine(PATIENT.out.gliph2_clone_network.map { l -> [l] }) + .combine(PATIENT.out.gliph2_cluster_member_details.map { l -> [l] }) + .combine(PATIENT.out.gliph2_global_similarities.map { l -> [l] }) + .map { sample_stats_csv, concat_cdr3_sorted, shared_cdr3, giana_files, + all_motifs_pairs, clone_network_pairs, cluster_member_pairs, global_sim_pairs -> + def report_files = [sample_stats_csv, concat_cdr3_sorted, shared_cdr3] + giana_files + + all_motifs_pairs.collect { p -> p[1] } + + clone_network_pairs.collect { p -> p[1] } + + cluster_member_pairs.collect { p -> p[1] } + + global_sim_pairs.collect { p -> p[1] } + + part2_include_files + // gliph2 files keep their patient-prefixed basename (e.g. + // "patientA_all_motifs.txt", set in GLIPH2_TURBOGLIPH to avoid + // collisions across patients when staged) all the way through - + // template_gliph.qmd reads that same prefixed name from within + // each patient's subdir, so source and dest basenames match here. + def staged_layout = [ + ["${params.project_name}/sample/${sample_stats_csv.name}", sample_stats_csv.name], + ["${params.project_name}/annotate/${concat_cdr3_sorted.name}", concat_cdr3_sorted.name], + ["${params.project_name}/tcrsharing/${shared_cdr3.name}", shared_cdr3.name] + ] + giana_files.collect { f -> ["${params.project_name}/giana/${f.name}", f.name] } + + all_motifs_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + + clone_network_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + + cluster_member_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + + global_sim_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + tuple( + file(params.template_details_part2), + report_files, + staged_layout + ) + } + } else { + ch_details_part2_report = details_part2_base + .map { sample_stats_csv, concat_cdr3_sorted, shared_cdr3 -> + tuple( + file(params.template_details_part2), + [sample_stats_csv, concat_cdr3_sorted, shared_cdr3] + part2_include_files, + [ + ["${params.project_name}/sample/${sample_stats_csv.name}", sample_stats_csv.name], + ["${params.project_name}/annotate/${concat_cdr3_sorted.name}", concat_cdr3_sorted.name], + ["${params.project_name}/tcrsharing/${shared_cdr3.name}", shared_cdr3.name] + ] + ) + } + } + ch_reports = ch_reports.mix(ch_details_part2_report) REPORT( ch_reports ) From 3b8621e24fb1a68aff671675e2c115216330d23c Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Wed, 22 Jul 2026 16:07:07 -0400 Subject: [PATCH 2/4] - Add timepoint_order column to the minimal-example samplesheet (1=Base, 2=Post), required by timepoint_order_col in the report notebooks. Update the matching header assertion in samplesheet_check.nf.test. - Make GIANA/GLIPH2 inclusion in template_details_part2.qmd conditional on the patient workflow_level: splice in a single generic template_patient_clustering.qmd, which workflows/tcrtoolkit.nf resolves to either template_patient_clustering_on.qmd (includes both notebooks) or a placeholder off.qmd, mirroring the existing template_pheno.qmd mechanism. Avoids a hard crash in template_giana.qmd when patient clustering never ran, and a wall of "no data" warnings otherwise. - Fix subworkflows/local/patient.nf's gliph2_* emits: plain .collect() on tuple(val(patient), path) channels flattens the tuples by default, corrupting downstream [patient, file] pair indexing (surfaced as a cryptic "No such variable: name" error). Fixed with .collect(flat: false). - Handle GLIPH2/GIANA finding no significant results, which isn't always a clean "no results" case but can also leave output files missing or malformed: - gliph2.nf: wrap turboGliph::gliph2() in tryCatch (it can throw an uncaught internal dplyr type-mismatch error) and backfill empty placeholders for any of its output files left missing. - template_giana.qmd: catch EmptyDataError when GIANA writes only metadata comments and no data for a patient, at all 3 read sites, and guard the final concat against an all-empty patient list. - template_gliph.qmd: skip genuinely 0-byte output files (not just missing ones), which also crash pd.read_csv. - Set embed-resources: true in template_details_part1/part2.qmd, matching the other reports - RENDER_NOTEBOOK only publishes the bare .html, so the separate CSS/JS assets Quarto writes under embed-resources: false were never shipped, leaving these two reports unstyled. --- modules/local/compare/gliph2.nf | 90 +++++++++++++++---- nextflow.config | 2 + nextflow_schema.json | 8 ++ notebooks/template_details_part1.qmd | 2 +- notebooks/template_details_part2.qmd | 6 +- notebooks/template_giana.qmd | 38 ++++++-- notebooks/template_gliph.qmd | 17 ++-- notebooks/template_sharing.qmd | 1 - subworkflows/local/patient.nf | 14 ++- .../samplesheet/samplesheet_check.nf.test | 2 +- .../test_data/minimal-example/samplesheet.csv | 72 +++++++-------- workflows/tcrtoolkit.nf | 32 +++++-- 12 files changed, 199 insertions(+), 85 deletions(-) diff --git a/modules/local/compare/gliph2.nf b/modules/local/compare/gliph2.nf index ec71b91..b0ea684 100644 --- a/modules/local/compare/gliph2.nf +++ b/modules/local/compare/gliph2.nf @@ -22,6 +22,8 @@ process GLIPH2_TURBOGLIPH { script: """ + mkdir -p ${patient} + Rscript - <}} -{{< include ./template_giana.qmd >}} - -{{< include ./template_gliph.qmd >}} +{{< include ./template_patient_clustering.qmd >}} diff --git a/notebooks/template_giana.qmd b/notebooks/template_giana.qmd index b850832..2bf0893 100644 --- a/notebooks/template_giana.qmd +++ b/notebooks/template_giana.qmd @@ -143,8 +143,16 @@ def plot_giana_tables_per_patient(meta_df, giana_dir, subject_col, timepoint_col display(Markdown(f"*[WARNING] No GIANA file found for patient '{patient}' (expected: {giana_file}). Skipping.*\n\n")) continue - # Load this patient's GIANA output - giana_df = pd.read_csv(giana_file, comment='#', sep='\t') + # GIANA writes only its "##"-prefixed metadata/header comment lines + # (no data rows at all) when it finds zero clusters for a patient's + # samples - after comment='#' strips those, pandas has nothing left + # to parse and raises EmptyDataError. Skip gracefully like any other + # "no data" case instead of crashing the whole report. + try: + giana_df = pd.read_csv(giana_file, comment='#', sep='\t') + except pd.errors.EmptyDataError: + display(Markdown(f"*[INFO] No GIANA clusters found for patient '{patient}'. Skipping.*\n\n")) + continue # Merge with all samples belonging to this patient patient_meta = meta_df[meta_df[subject_col] == patient][['sample', subject_col, timepoint_col]] @@ -219,8 +227,13 @@ for patient in meta[subject_col].unique(): if not os.path.exists(giana_file): continue - # Load this patient's GIANA output - giana_df = pd.read_csv(giana_file, comment='#', sep='\t') + # See the comment on the equivalent read above: GIANA writes only + # comment lines with zero clusters, which leaves nothing for pandas to + # parse once comment='#' strips them. + try: + giana_df = pd.read_csv(giana_file, comment='#', sep='\t') + except pd.errors.EmptyDataError: + continue # Merge with all samples belonging to this patient patient_meta = meta[meta[subject_col] == patient][['sample', subject_col, timepoint_col]] @@ -293,18 +306,25 @@ for patient in meta[subject_col].unique(): giana_file = os.path.join(giana_dir, f"{patient}_giana.txt") if os.path.exists(giana_file): - patient_df = pd.read_csv(giana_file, comment='#', sep='\t') - - # CRITICAL FIX: Make cluster IDs unique to the patient + # See the comment on the first read in this notebook: GIANA writes + # only comment lines with zero clusters, which leaves nothing for + # pandas to parse once comment='#' strips them. + try: + patient_df = pd.read_csv(giana_file, comment='#', sep='\t') + except pd.errors.EmptyDataError: + print(f"*(No GIANA clusters found for: {giana_file})*") + continue + + # CRITICAL FIX: Make cluster IDs unique to the patient # to prevent false merging across different GIANA runs patient_df['cluster'] = str(patient) + "_c" + patient_df['cluster'].astype(str) - + all_giana_data.append(patient_df) else: print(f"*(No GIANA file found for: {giana_file})*") # Combine into a single master dataframe -giana_df = pd.concat(all_giana_data, ignore_index=True) +giana_df = pd.concat(all_giana_data, ignore_index=True) if all_giana_data else pd.DataFrame() # Merge metadata giana_df = giana_df.merge(meta[['sample', subject_col, timepoint_col, alias_col]], on='sample', how='left') diff --git a/notebooks/template_gliph.qmd b/notebooks/template_gliph.qmd index 349a4ac..d4cb1bb 100644 --- a/notebooks/template_gliph.qmd +++ b/notebooks/template_gliph.qmd @@ -48,22 +48,29 @@ for patient in meta[subject_col].unique(): am_csv = f"{p_dir}/{patient}_all_motifs.txt" gs_csv = f"{p_dir}/{patient}_global_similarities.txt" - if os.path.exists(cd_csv): + # GLIPH2 writes these as genuinely empty (0-byte) files, not just + # header-only, when it finds no clusters for a patient (small/sparse + # samples) - pd.read_csv raises EmptyDataError on those, so an + # os.path.exists check alone isn't enough; skip zero-size files too. + def _has_content(path): + return os.path.exists(path) and os.path.getsize(path) > 0 + + if _has_content(cd_csv): df = pd.read_csv(cd_csv, sep='\t') df[subject_col] = patient cluster_details_list.append(df) - if os.path.exists(cn_csv): + if _has_content(cn_csv): # clone_network has no header - df = pd.read_csv(cn_csv, sep='\t', header=None) + df = pd.read_csv(cn_csv, sep='\t', header=None) clone_network_list.append(df) - if os.path.exists(am_csv): + if _has_content(am_csv): df = pd.read_csv(am_csv, sep='\t') df[subject_col] = patient all_motifs_list.append(df) - if os.path.exists(gs_csv): + if _has_content(gs_csv): df = pd.read_csv(gs_csv, sep='\t') df[subject_col] = patient # CRITICAL: Prefix cluster tags to prevent cross-patient collisions diff --git a/notebooks/template_sharing.qmd b/notebooks/template_sharing.qmd index 11dd0ef..342f79e 100644 --- a/notebooks/template_sharing.qmd +++ b/notebooks/template_sharing.qmd @@ -48,7 +48,6 @@ import base64 import json from IPython.display import HTML, display import warnings -from upsetplot import from_contents, plot as upset_plot # - Loading data diff --git a/subworkflows/local/patient.nf b/subworkflows/local/patient.nf index 423a920..1434ddf 100644 --- a/subworkflows/local/patient.nf +++ b/subworkflows/local/patient.nf @@ -62,10 +62,16 @@ workflow PATIENT { ) } + // channel .collect() flattens tuple(val, path) emissions by default (e.g. + // [patientA, fileA, patientB, fileB] instead of [[patientA, fileA], + // [patientB, fileB]]) - confirmed via a minimal repro. workflows/tcrtoolkit.nf + // indexes these as [patient, file] pairs (p[0]/p[1]) downstream when staging + // per-patient gliph2 outputs, so flat: false is required to preserve that + // shape instead of silently interleaving patients and files. emit: giana_files = GIANA_CALC.out.giana_output.collect() - gliph2_all_motifs = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.all_motifs.collect() : channel.value([]) - gliph2_clone_network = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.clone_network.collect() : channel.value([]) - gliph2_cluster_member_details = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.cluster_member_details.collect() : channel.value([]) - gliph2_global_similarities = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.global_similarities.collect() : channel.value([]) + gliph2_all_motifs = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.all_motifs.collect(flat: false) : channel.value([]) + gliph2_clone_network = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.clone_network.collect(flat: false) : channel.value([]) + gliph2_cluster_member_details = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.cluster_member_details.collect(flat: false) : channel.value([]) + gliph2_global_similarities = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.global_similarities.collect(flat: false) : channel.value([]) } \ No newline at end of file diff --git a/tests/modules/local/samplesheet/samplesheet_check.nf.test b/tests/modules/local/samplesheet/samplesheet_check.nf.test index 92328c9..575d059 100644 --- a/tests/modules/local/samplesheet/samplesheet_check.nf.test +++ b/tests/modules/local/samplesheet/samplesheet_check.nf.test @@ -55,7 +55,7 @@ nextflow_process { def lines = path(samplesheet_utf8.get(0)).readLines() assert lines.size() == 36 // Header + 35 samples - assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,origin,file" + assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,timepoint_order,origin,file" } } } diff --git a/tests/test_data/minimal-example/samplesheet.csv b/tests/test_data/minimal-example/samplesheet.csv index 5f89b10..08403f7 100644 --- a/tests/test_data/minimal-example/samplesheet.csv +++ b/tests/test_data/minimal-example/samplesheet.csv @@ -1,36 +1,36 @@ -sample,alias,patient,timepoint,origin,file -Patient01_Base,Patient01_Base,Patient01,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient01_Base.tsv -Patient02_Base,Patient02_Base,Patient02,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient02_Base.tsv -Patient03_Base,Patient03_Base,Patient03,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient03_Base.tsv -Patient04_Base,Patient04_Base,Patient04,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient04_Base.tsv -Patient05_Base,Patient05_Base,Patient05,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient05_Base.tsv -Patient06_Base,Patient06_Base,Patient06,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient06_Base.tsv -Patient07_Base,Patient07_Base,Patient07,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient07_Base.tsv -Patient08_Base,Patient08_Base,Patient08,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient08_Base.tsv -Patient09_Base,Patient09_Base,Patient09,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient09_Base.tsv -Patient10_Base,Patient10_Base,Patient10,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient10_Base.tsv -Patient11_Base,Patient11_Base,Patient11,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient11_Base.tsv -Patient12_Base,Patient12_Base,Patient12,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient12_Base.tsv -Patient13_Base,Patient13_Base,Patient13,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient13_Base.tsv -Patient14_Base,Patient14_Base,Patient14,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient14_Base.tsv -Patient15_Base,Patient15_Base,Patient15,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient15_Base.tsv -Patient16_Base,Patient16_Base,Patient16,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Base.tsv -Patient16_Post,Patient16_Post,Patient16,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Post.tsv -Patient17_Base,Patient17_Base,Patient17,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Base.tsv -Patient17_Post,Patient17_Post,Patient17,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Post.tsv -Patient18_Base,Patient18_Base,Patient18,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Base.tsv -Patient18_Post,Patient18_Post,Patient18,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Post.tsv -Patient19_Base,Patient19_Base,Patient19,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Base.tsv -Patient19_Post,Patient19_Post,Patient19,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Post.tsv -Patient20_Base,Patient20_Base,Patient20,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Base.tsv -Patient20_Post,Patient20_Post,Patient20,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Post.tsv -Patient21_Base,Patient21_Base,Patient21,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Base.tsv -Patient21_Post,Patient21_Post,Patient21,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Post.tsv -Patient22_Base,Patient22_Base,Patient22,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Base.tsv -Patient22_Post,Patient22_Post,Patient22,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Post.tsv -Patient23_Base,Patient23_Base,Patient23,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Base.tsv -Patient23_Post,Patient23_Post,Patient23,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Post.tsv -Patient24_Base,Patient24_Base,Patient24,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Base.tsv -Patient24_Post,Patient24_Post,Patient24,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Post.tsv -Patient25_Base,Patient25_Base,Patient25,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Base.tsv -Patient25_Post,Patient25_Post,Patient25,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Post.tsv +sample,alias,patient,timepoint,timepoint_order,origin,file +Patient01_Base,Patient01_Base,Patient01,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient01_Base.tsv +Patient02_Base,Patient02_Base,Patient02,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient02_Base.tsv +Patient03_Base,Patient03_Base,Patient03,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient03_Base.tsv +Patient04_Base,Patient04_Base,Patient04,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient04_Base.tsv +Patient05_Base,Patient05_Base,Patient05,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient05_Base.tsv +Patient06_Base,Patient06_Base,Patient06,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient06_Base.tsv +Patient07_Base,Patient07_Base,Patient07,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient07_Base.tsv +Patient08_Base,Patient08_Base,Patient08,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient08_Base.tsv +Patient09_Base,Patient09_Base,Patient09,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient09_Base.tsv +Patient10_Base,Patient10_Base,Patient10,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient10_Base.tsv +Patient11_Base,Patient11_Base,Patient11,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient11_Base.tsv +Patient12_Base,Patient12_Base,Patient12,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient12_Base.tsv +Patient13_Base,Patient13_Base,Patient13,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient13_Base.tsv +Patient14_Base,Patient14_Base,Patient14,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient14_Base.tsv +Patient15_Base,Patient15_Base,Patient15,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient15_Base.tsv +Patient16_Base,Patient16_Base,Patient16,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Base.tsv +Patient16_Post,Patient16_Post,Patient16,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Post.tsv +Patient17_Base,Patient17_Base,Patient17,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Base.tsv +Patient17_Post,Patient17_Post,Patient17,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Post.tsv +Patient18_Base,Patient18_Base,Patient18,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Base.tsv +Patient18_Post,Patient18_Post,Patient18,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Post.tsv +Patient19_Base,Patient19_Base,Patient19,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Base.tsv +Patient19_Post,Patient19_Post,Patient19,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Post.tsv +Patient20_Base,Patient20_Base,Patient20,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Base.tsv +Patient20_Post,Patient20_Post,Patient20,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Post.tsv +Patient21_Base,Patient21_Base,Patient21,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Base.tsv +Patient21_Post,Patient21_Post,Patient21,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Post.tsv +Patient22_Base,Patient22_Base,Patient22,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Base.tsv +Patient22_Post,Patient22_Post,Patient22,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Post.tsv +Patient23_Base,Patient23_Base,Patient23,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Base.tsv +Patient23_Post,Patient23_Post,Patient23,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Post.tsv +Patient24_Base,Patient24_Base,Patient24,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Base.tsv +Patient24_Post,Patient24_Post,Patient24,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Post.tsv +Patient25_Base,Patient25_Base,Patient25,Base,1,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Base.tsv +Patient25_Post,Patient25_Post,Patient25,Post,2,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Post.tsv diff --git a/workflows/tcrtoolkit.nf b/workflows/tcrtoolkit.nf index 32d4843..4539921 100644 --- a/workflows/tcrtoolkit.nf +++ b/workflows/tcrtoolkit.nf @@ -235,15 +235,24 @@ workflow TCRTOOLKIT { ch_reports = ch_reports.mix(ch_details_part1_report) // Details report (part 2) includes template_overlap.qmd, template_sharing.qmd, - // template_giana.qmd, and template_gliph.qmd. The latter two need patient-level - // clonotype clustering outputs, which only exist when 'patient' is in - // workflow_level; their own os.path.exists guards degrade gracefully when those - // files are absent, so those files are only staged when PATIENT actually ran. + // and a single generic template_patient_clustering.qmd - mirroring the + // template_pheno.qmd trick above, since template_giana.qmd/template_gliph.qmd + // need patient-level clonotype clustering outputs that only exist when + // 'patient' is in workflow_level. Rather than always splicing them in and + // relying on their internal os.path.exists guards (GIANA's own concat call + // crashes outright when no per-patient file exists at all, i.e. whenever + // PATIENT didn't run), Nextflow stages the "on" wrapper (which includes both) + // or the "off" placeholder under that shared name depending on workflow_level. def details_part2_base = sample_stats_agg .combine(ANNOTATE.out.concat_cdr3_sorted) .combine(COMPARE.out.shared_cdr3) - // template_details_part2.qmd pulls in these 4 files via {{< include >}} + def run_patient_clustering = levels.contains('patient') + def patient_clustering_notebook = run_patient_clustering + ? file(params.template_patient_clustering_on) + : file(params.template_patient_clustering_off) + + // template_details_part2.qmd pulls in these files via {{< include >}} // shortcodes resolved relative to its own directory, so they have to be // staged alongside it too - see the equivalent part1 comment above. def part2_notebooks_dir = file(params.template_details_part2).parent @@ -251,10 +260,14 @@ workflow TCRTOOLKIT { file("${part2_notebooks_dir}/template_overlap.qmd"), file("${part2_notebooks_dir}/template_sharing.qmd"), file("${part2_notebooks_dir}/template_giana.qmd"), - file("${part2_notebooks_dir}/template_gliph.qmd") + file("${part2_notebooks_dir}/template_gliph.qmd"), + patient_clustering_notebook + ] + def part2_staged_layout_extra = [ + ["template_patient_clustering.qmd", patient_clustering_notebook.name] ] - if (levels.contains('patient')) { + if (run_patient_clustering) { ch_details_part2_report = details_part2_base .combine(PATIENT.out.giana_files.map { l -> [l] }) .combine(PATIENT.out.gliph2_all_motifs.map { l -> [l] }) @@ -282,7 +295,8 @@ workflow TCRTOOLKIT { all_motifs_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + clone_network_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + cluster_member_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + - global_sim_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + global_sim_pairs.collect { p -> ["${params.project_name}/gliph2/${p[0]}/${p[1].name}", p[1].name] } + + part2_staged_layout_extra tuple( file(params.template_details_part2), report_files, @@ -299,7 +313,7 @@ workflow TCRTOOLKIT { ["${params.project_name}/sample/${sample_stats_csv.name}", sample_stats_csv.name], ["${params.project_name}/annotate/${concat_cdr3_sorted.name}", concat_cdr3_sorted.name], ["${params.project_name}/tcrsharing/${shared_cdr3.name}", shared_cdr3.name] - ] + ] + part2_staged_layout_extra ) } } From 6f9410deb00ea279cb1423f385afd8b2ed6bcc67 Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Wed, 22 Jul 2026 16:53:02 -0400 Subject: [PATCH 3/4] Removed unnecessary commented code --- modules/local/compare/gliph2.nf | 28 ++++++-------------- notebooks/template_giana.qmd | 7 ++--- notebooks/template_gliph.qmd | 6 ++--- subworkflows/local/patient.nf | 10 +++---- workflows/tcrtoolkit.nf | 46 ++++++++++++--------------------- 5 files changed, 32 insertions(+), 65 deletions(-) diff --git a/modules/local/compare/gliph2.nf b/modules/local/compare/gliph2.nf index b0ea684..187a4b5 100644 --- a/modules/local/compare/gliph2.nf +++ b/modules/local/compare/gliph2.nf @@ -36,15 +36,9 @@ process GLIPH2_TURBOGLIPH { colnames(df)[colnames(df) == "duplicate_count"] <- "counts" df[,'patient'] <- df[,'sample'] - # turboGliph::gliph2() can throw an uncaught internal error - not just - # print "no clusters/similarities found" - on edge-case-shaped input, e.g. - # a dplyr::left_join() type mismatch between the tag column of sample_stats - # (factor) and the tag column of ref_stats (double) when very few - # candidate motifs pass its internal filtering (seen on real Patient02 - # data: crashes partway through "Part 2: Searching for global - # similarities", after only all_motifs.txt and local_similarities_*.txt - # were written). Catch this and fall back to empty results for the whole - # patient rather than failing the task. + # gliph2() can throw an uncaught internal error (e.g. a dplyr::left_join() + # type mismatch) rather than just reporting "no results found" - catch it + # and fall back to empty results for the whole patient. gliph2_ok <- tryCatch({ turboGliph::gliph2( cdr3_sequences = df, @@ -62,11 +56,9 @@ process GLIPH2_TURBOGLIPH { FALSE }) - # gliph2() also doesn't error on the narrower "no significant clusters" - # case - it just writes cluster_member_details.txt as a single blank line - # with no header, which read.csv rejects outright ("no lines available in - # input"). The same fallback covers that case, a totally missing file (if - # gliph2() failed before reaching this point), and the tryCatch above. + # Also covers the narrower "no significant clusters" case, where gliph2() + # writes cluster_member_details.txt as a single blank line with no header, + # which read.csv rejects ("no lines available in input"). df3 <- tryCatch({ tmp <- read.csv('${patient}/cluster_member_details.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE) tmp[,'sample'] <- tmp[,'patient'] @@ -78,12 +70,8 @@ process GLIPH2_TURBOGLIPH { }) write.table(df3, "${patient}/cluster_member_details.txt", sep = "\t", row.names = FALSE, quote = FALSE) - # If gliph2() failed partway through (or before writing anything), some of - # its other declared output files may be entirely missing - write minimal - # placeholders for whichever ones aren't already there so Nextflow's own - # output declarations (which require every path to exist) don't turn this - # into a hard pipeline failure on top of the already-handled analysis - # failure. + # Backfill placeholders for any other declared output file gliph2() left + # missing, so Nextflow's output declarations don't fail the task. if (!file.exists('${patient}/all_motifs.txt')) { write.table(data.frame(motif=character(), num_in_sample=integer(), num_in_ref=integer(), fisher.score=double(), num_fold=double()), diff --git a/notebooks/template_giana.qmd b/notebooks/template_giana.qmd index 2bf0893..ec4894b 100644 --- a/notebooks/template_giana.qmd +++ b/notebooks/template_giana.qmd @@ -143,11 +143,8 @@ def plot_giana_tables_per_patient(meta_df, giana_dir, subject_col, timepoint_col display(Markdown(f"*[WARNING] No GIANA file found for patient '{patient}' (expected: {giana_file}). Skipping.*\n\n")) continue - # GIANA writes only its "##"-prefixed metadata/header comment lines - # (no data rows at all) when it finds zero clusters for a patient's - # samples - after comment='#' strips those, pandas has nothing left - # to parse and raises EmptyDataError. Skip gracefully like any other - # "no data" case instead of crashing the whole report. + # GIANA writes only "##" comment lines when it finds zero clusters, + # which leaves pandas nothing to parse after comment='#' strips them. try: giana_df = pd.read_csv(giana_file, comment='#', sep='\t') except pd.errors.EmptyDataError: diff --git a/notebooks/template_gliph.qmd b/notebooks/template_gliph.qmd index d4cb1bb..d77ea33 100644 --- a/notebooks/template_gliph.qmd +++ b/notebooks/template_gliph.qmd @@ -48,10 +48,8 @@ for patient in meta[subject_col].unique(): am_csv = f"{p_dir}/{patient}_all_motifs.txt" gs_csv = f"{p_dir}/{patient}_global_similarities.txt" - # GLIPH2 writes these as genuinely empty (0-byte) files, not just - # header-only, when it finds no clusters for a patient (small/sparse - # samples) - pd.read_csv raises EmptyDataError on those, so an - # os.path.exists check alone isn't enough; skip zero-size files too. + # Some GLIPH2 outputs are genuinely 0-byte files when it finds no + # clusters, which crashes pd.read_csv - skip empty files too, not just missing. def _has_content(path): return os.path.exists(path) and os.path.getsize(path) > 0 diff --git a/subworkflows/local/patient.nf b/subworkflows/local/patient.nf index 1434ddf..e02a2d4 100644 --- a/subworkflows/local/patient.nf +++ b/subworkflows/local/patient.nf @@ -62,12 +62,10 @@ workflow PATIENT { ) } - // channel .collect() flattens tuple(val, path) emissions by default (e.g. - // [patientA, fileA, patientB, fileB] instead of [[patientA, fileA], - // [patientB, fileB]]) - confirmed via a minimal repro. workflows/tcrtoolkit.nf - // indexes these as [patient, file] pairs (p[0]/p[1]) downstream when staging - // per-patient gliph2 outputs, so flat: false is required to preserve that - // shape instead of silently interleaving patients and files. + // .collect() flattens tuple(val, path) emissions by default (e.g. + // [patientA, fileA, patientB, fileB] instead of [[patientA, fileA], ...]), + // which corrupts the [patient, file] pair indexing used downstream in + // workflows/tcrtoolkit.nf - flat: false preserves the pair shape. emit: giana_files = GIANA_CALC.out.giana_output.collect() gliph2_all_motifs = params.use_gliph2 ? GLIPH2_TURBOGLIPH.out.all_motifs.collect(flat: false) : channel.value([]) diff --git a/workflows/tcrtoolkit.nf b/workflows/tcrtoolkit.nf index 4539921..bf9742f 100644 --- a/workflows/tcrtoolkit.nf +++ b/workflows/tcrtoolkit.nf @@ -142,7 +142,6 @@ workflow TCRTOOLKIT { ? CONVERT.out.pseudobulk_phenotype_files.map { _meta, files -> files }.flatten().collect() : channel.value([]) - // QC report requires sample-level aggregate outputs. ch_qc_report = sample_stats_agg .combine(ANNOTATE.out.concat_cdr3_sorted) .map { sample_stats_csv, concat_cdr3_sorted -> @@ -155,19 +154,15 @@ workflow TCRTOOLKIT { } ch_reports = ch_reports.mix(ch_qc_report) - // Discovery brief requires sample-level aggregate outputs, TCRdist3/VDJdb - // per-sample outputs, and cross-sample TCR sharing - it reads them from a - // project_dir/// layout, so staged_layout tells - // RENDER_NOTEBOOK where to symlink each staged file within the task work dir. - // Each staged_layout entry is a [dest_path, source_basename] pair - source and - // dest basenames match here, but don't always (see gliph2 below). + // Discovery brief reads its inputs from a project_dir// + // / layout, so staged_layout tells RENDER_NOTEBOOK where to + // symlink each staged file - each entry is a [dest_path, source_basename] + // pair (source and dest basenames usually match, but not always - see + // gliph2 below). // - // NOTE: .collect()-produced list channels are wrapped in an extra list via - // `.map { l -> [l] }` before every .combine() below - .combine() otherwise - // flattens a list-valued item's contents into the tuple instead of keeping - // it as one element (confirmed directly: combining a scalar channel with an - // unwrapped .collect() channel of 3 items produces a 4-element flat tuple, - // not a 2-element [scalar, list] tuple). + // .collect()-produced list channels are wrapped via `.map { l -> [l] }` + // before every .combine() below - otherwise .combine() flattens the + // list's contents into the tuple instead of keeping it as one element. ch_discovery_report = sample_stats_agg .combine(ANNOTATE.out.concat_cdr3_sorted) .combine(COMPARE.out.shared_cdr3) @@ -197,10 +192,6 @@ workflow TCRTOOLKIT { } ch_reports = ch_reports.mix(ch_discovery_report) - // Details report (part 1) includes template_sample.qmd, which reads sample-level - // aggregate outputs, the V/J gene-family pivot tables, and per-sample TCRdist3, - // OLGA, VDJdb, and convergence outputs from the same project_dir// - // / layout as the discovery brief. ch_details_part1_report = sample_stats_agg .combine(ANNOTATE.out.concat_cdr3_sorted) .combine(SAMPLE.out.v_family) @@ -234,15 +225,11 @@ workflow TCRTOOLKIT { } ch_reports = ch_reports.mix(ch_details_part1_report) - // Details report (part 2) includes template_overlap.qmd, template_sharing.qmd, - // and a single generic template_patient_clustering.qmd - mirroring the - // template_pheno.qmd trick above, since template_giana.qmd/template_gliph.qmd - // need patient-level clonotype clustering outputs that only exist when - // 'patient' is in workflow_level. Rather than always splicing them in and - // relying on their internal os.path.exists guards (GIANA's own concat call - // crashes outright when no per-patient file exists at all, i.e. whenever - // PATIENT didn't run), Nextflow stages the "on" wrapper (which includes both) - // or the "off" placeholder under that shared name depending on workflow_level. + // template_giana.qmd/template_gliph.qmd need patient-level clustering + // outputs that only exist when 'patient' is in workflow_level (GIANA's + // concat crashes outright otherwise), so template_patient_clustering.qmd + // is resolved to an "on" (includes both) or "off" (placeholder) wrapper, + // mirroring the template_pheno.qmd trick above. def details_part2_base = sample_stats_agg .combine(ANNOTATE.out.concat_cdr3_sorted) .combine(COMPARE.out.shared_cdr3) @@ -283,10 +270,9 @@ workflow TCRTOOLKIT { global_sim_pairs.collect { p -> p[1] } + part2_include_files // gliph2 files keep their patient-prefixed basename (e.g. - // "patientA_all_motifs.txt", set in GLIPH2_TURBOGLIPH to avoid - // collisions across patients when staged) all the way through - - // template_gliph.qmd reads that same prefixed name from within - // each patient's subdir, so source and dest basenames match here. + // "patientA_all_motifs.txt") all the way through, since + // template_gliph.qmd reads that same name from each patient's + // subdir - source and dest basenames match here. def staged_layout = [ ["${params.project_name}/sample/${sample_stats_csv.name}", sample_stats_csv.name], ["${params.project_name}/annotate/${concat_cdr3_sorted.name}", concat_cdr3_sorted.name], From 78bb82fe17360b0bf83b8b599dde0972773d10d7 Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Wed, 22 Jul 2026 19:17:51 -0400 Subject: [PATCH 4/4] Add missing template_patient_clustering_on/off.qmd notebooks These were created in 3b8621e but silently blocked from git add by the stale `notebooks/*` .gitignore rule, so they only ever existed on disk locally. CI checks out from git, so the RENDER_NOTEBOOK symlink for template_patient_clustering.qmd pointed at a nonexistent file, and Quarto's {{< include >}} failed with "could not find file". Co-Authored-By: Claude Sonnet 5 --- notebooks/template_patient_clustering_off.qmd | 5 +++++ notebooks/template_patient_clustering_on.qmd | 3 +++ 2 files changed, 8 insertions(+) create mode 100644 notebooks/template_patient_clustering_off.qmd create mode 100644 notebooks/template_patient_clustering_on.qmd diff --git a/notebooks/template_patient_clustering_off.qmd b/notebooks/template_patient_clustering_off.qmd new file mode 100644 index 0000000..505afce --- /dev/null +++ b/notebooks/template_patient_clustering_off.qmd @@ -0,0 +1,5 @@ +## Patient-level TCR clustering (GIANA/GLIPH2) + +::: {.callout-note} +This section requires patient-level clonotype clustering (GIANA and GLIPH2), which was not run for this project. Include `patient` in `--workflow_level` to generate it. +::: diff --git a/notebooks/template_patient_clustering_on.qmd b/notebooks/template_patient_clustering_on.qmd new file mode 100644 index 0000000..f56d020 --- /dev/null +++ b/notebooks/template_patient_clustering_on.qmd @@ -0,0 +1,3 @@ +{{< include ./template_giana.qmd >}} + +{{< include ./template_gliph.qmd >}}