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..187a4b5 100644
--- a/modules/local/compare/gliph2.nf
+++ b/modules/local/compare/gliph2.nf
@@ -8,16 +8,22 @@ 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'
     
     script:
     """
+    mkdir -p ${patient}
+
     Rscript - <<EOF
     #!/usr/bin/env Rscript
 
@@ -30,27 +36,80 @@ process GLIPH2_TURBOGLIPH {
     colnames(df)[colnames(df) == "duplicate_count"] <- "counts"
     df[,'patient'] <- df[,'sample']
 
-    result <- turboGliph::gliph2(
-        cdr3_sequences = df,
-        result_folder = "./${patient}",
-        lcminp = ${params.local_min_pvalue},
-        sim_depth = ${params.simulation_depth},
-        kmer_mindepth = ${params.kmer_min_depth},
-        lcminove = ${params.local_min_OVE},
-        all_aa_interchangeable = FALSE,
-        n_cores = ${task.cpus}
-    )
-
-    df3 <- read.csv('${patient}/cluster_member_details.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
-    df3[,'sample'] <- df3[,'patient']
-    df3 <- merge(df3, df[, c("CDR3b", "TRBV", "sample", 'counts')], by = c("CDR3b", "TRBV", "sample", 'counts'), all.x = TRUE)
-    df3 <- df3[, c('CDR3b', 'TRBV', 'TRBJ', 'counts', 'sample', 'tag', 'seq_ID', 'ultCDR3b')]
+    # 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,
+            result_folder = "./${patient}",
+            lcminp = ${params.local_min_pvalue},
+            sim_depth = ${params.simulation_depth},
+            kmer_mindepth = ${params.kmer_min_depth},
+            lcminove = ${params.local_min_OVE},
+            all_aa_interchangeable = FALSE,
+            n_cores = ${task.cpus}
+        )
+        TRUE
+    }, error = function(e) {
+        message("turboGliph::gliph2() failed for ${patient}, treating as no clusters found: ", conditionMessage(e))
+        FALSE
+    })
+
+    # 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']
+        tmp <- merge(tmp, df[, c("CDR3b", "TRBV", "sample", 'counts')], by = c("CDR3b", "TRBV", "sample", 'counts'), all.x = TRUE)
+        tmp[, c('CDR3b', 'TRBV', 'TRBJ', 'counts', 'sample', 'tag', 'seq_ID', 'ultCDR3b')]
+    }, error = function(e) {
+        data.frame(CDR3b=character(), TRBV=character(), TRBJ=character(), counts=integer(),
+                   sample=character(), tag=character(), seq_ID=character(), ultCDR3b=character())
+    })
     write.table(df3, "${patient}/cluster_member_details.txt", sep = "\t", row.names = FALSE, quote = FALSE)
+
+    # 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()),
+                    '${patient}/all_motifs.txt', sep = "\t", row.names = FALSE, quote = FALSE)
+    }
+    if (!file.exists('${patient}/clone_network.txt')) {
+        file.create('${patient}/clone_network.txt')
+    }
+    if (!file.exists('${patient}/global_similarities.txt')) {
+        write.table(data.frame(cluster_tag=character(), cluster_size=integer(), unique_CDR3b=integer(),
+                                num_in_ref=integer(), fisher.score=double(), aa_at_position=character(),
+                                TRBV=character(), CDR3b=character()),
+                    '${patient}/global_similarities.txt', sep = "\t", row.names = FALSE, quote = FALSE)
+    }
+    if (!file.exists('${patient}/convergence_groups.txt')) {
+        file.create('${patient}/convergence_groups.txt')
+    }
+    if (length(Sys.glob('${patient}/local_similarities_*.txt')) == 0) {
+        write.table(data.frame(motif=character(), num_in_sample=integer(), num_in_ref=integer(),
+                                fisher.score=double(), num_fold=double(), start=integer(), stop=integer(),
+                                members=character()),
+                    '${patient}/local_similarities_none.txt', sep = "\t", row.names = FALSE, quote = FALSE)
+    }
+    if (!file.exists('${patient}/parameter.txt')) {
+        file.create('${patient}/parameter.txt')
+    }
     EOF
 
     # 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..1c2cbf4 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -75,6 +75,13 @@ 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"
+    template_patient_clustering_on = "${projectDir}/notebooks/template_patient_clustering_on.qmd"
+    template_patient_clustering_off = "${projectDir}/notebooks/template_patient_clustering_off.qmd"
 }
 
 includeConfig 'conf/base.config'
diff --git a/nextflow_schema.json b/nextflow_schema.json
index 26d3c20..4c3ad79 100644
--- a/nextflow_schema.json
+++ b/nextflow_schema.json
@@ -120,6 +120,34 @@
         "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."
+        },
+        "template_patient_clustering_on": {
+          "type": "string",
+          "description": "Path to patient-level clustering (GIANA/GLIPH2) notebook template, used when patient workflow_level is run."
+        },
+        "template_patient_clustering_off": {
+          "type": "string",
+          "description": "Path to placeholder notebook template used when patient workflow_level is not run."
         }
       }
     },
diff --git a/notebooks/template_details_part1.qmd b/notebooks/template_details_part1.qmd
index d615125..0aef09d 100644
--- a/notebooks/template_details_part1.qmd
+++ b/notebooks/template_details_part1.qmd
@@ -6,7 +6,7 @@ format:
     toc: true
     toc_depth: 3
     code-fold: true
-    embed-resources: false
+    embed-resources: true
     number-sections: true
     smooth-scroll: true
     grid:
diff --git a/notebooks/template_details_part2.qmd b/notebooks/template_details_part2.qmd
index e66a05a..c02e088 100644
--- a/notebooks/template_details_part2.qmd
+++ b/notebooks/template_details_part2.qmd
@@ -6,7 +6,7 @@ format:
     toc: true
     toc_depth: 3
     code-fold: true
-    embed-resources: false
+    embed-resources: true
     number-sections: true
     smooth-scroll: true
     grid:
@@ -103,6 +103,4 @@ This pipeline can be used to analyze both **single-cell and bulk TCR data**. Ple
 
 {{< include ./template_sharing.qmd >}}
 
-{{< include ./template_giana.qmd >}}
-
-{{< include ./template_gliph.qmd >}}
+{{< include ./template_patient_clustering.qmd >}}
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. "<sample>_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..ec4894b 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]):
@@ -172,8 +143,13 @@ 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 "##" 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:
+            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]]
@@ -190,7 +166,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()
 
@@ -248,8 +224,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]]
@@ -261,7 +242,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 +291,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
 
@@ -324,18 +303,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')
@@ -350,7 +336,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 +346,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 +355,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 +404,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 +412,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 +423,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 +551,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..d77ea33 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,27 +43,32 @@ 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"
+
+    # 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
 
-    if os.path.exists(cd_csv):
+    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
@@ -102,89 +100,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 "<br>".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 "<br>".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"<b>CDR3:</b> {g.vs[i]['name']}<br><b>Found in:</b><br>{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"<b>CDR3:</b> {g.vs[i]['name']}<br><b>Found in:</b><br>{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 +256,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 +338,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_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 >}}
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..342f79e 100644
--- a/notebooks/template_sharing.qmd
+++ b/notebooks/template_sharing.qmd
@@ -48,8 +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
 
 
 # - Loading data
@@ -81,16 +79,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 +194,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 +213,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..e02a2d4 100644
--- a/subworkflows/local/patient.nf
+++ b/subworkflows/local/patient.nf
@@ -45,9 +45,31 @@ 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
         )
     }
+
+    // .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([])
+        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/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/<subdir>/<file> 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/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 f962d0b..bf9742f 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([])
+
+    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 reads its inputs from a project_dir/<project_name>/
+    // <subdir>/<file> 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).
+    //
+    // .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)
+        .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)
+
+    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)
+
+    // 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)
+
+    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
+    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"),
+        patient_clustering_notebook
+    ]
+    def part2_staged_layout_extra = [
+        ["template_patient_clustering.qmd", patient_clustering_notebook.name]
+    ]
+
+    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] })
+            .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") 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],
+                    ["${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] } +
+                    part2_staged_layout_extra
+                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]
+                    ] + part2_staged_layout_extra
+                )
+            }
+    }
+    ch_reports = ch_reports.mix(ch_details_part2_report)
 
     REPORT( ch_reports )