Add sample TCR convergence

Author: dltamayo

modules/local/convergence.nf

@@ -0,0 +1,60 @@
+process CONVERGENCE {
+    tag "${sample_meta[0]}"
+    label 'process_low'
+    container "ghcr.io/break-through-cancer/bulktcr:latest"
+
+    input:
+    tuple val(sample_meta), path(count_table)
+
+    output:
+    path "${count_table.baseName}_tcr_convergence.tsv", emit: "convergence_output"
+    path "${count_table.baseName}_tcr_convergence_histogram.png"
+
+    script:
+    """
+    # Extract vector of cdr3 aa, dropping null values
+    python - <<EOF
+    import pandas as pd
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    # Load your TCR data (make sure the file has 'cdr3_aa' and 'cdr3_nt' columns)
+    df = pd.read_csv("${count_table}", sep="\t", usecols=["aminoAcid", "nucleotide"])
+    df = df.dropna(subset=["aminoAcid"])
+
+    # Group by amino acid sequence and count unique nucleotide sequences (convergence)
+    convergence_df = (
+        df.groupby("aminoAcid")["nucleotide"]
+        .nunique()
+        .reset_index(name="convergence")
+    )
+
+    # Sort by convergence count, descending
+    convergence_df = convergence_df.sort_values(by="convergence", ascending=False)
+
+    # Export
+    convergence_df.to_csv("${count_table.baseName}_tcr_convergence.tsv", sep="\t", index=False)
+    
+    # Plot histogram
+    convergence = convergence_df['convergence'].values
+    average_convergence = convergence_df["convergence"].mean()
+
+    # Create integer bin edges from 0 to max(data)
+    bins = np.arange(min(convergence), max(convergence) + 2)  # +2 to include the last value as a bin edge
+
+    plt.figure(figsize=(8, 5))
+    plt.hist(convergence, bins=bins, edgecolor='black', align='left')
+    plt.xticks(bins[:-1])  # whole number positions and labels
+    plt.yscale('log')
+
+    plt.xlabel('TCR Convergence Number')
+    plt.ylabel('TCR Convergence Frequency (log scale)')
+    plt.title(f'${count_table.baseName} TCR Convergence Histogram, Average: {average_convergence:.2f}')
+
+    # Save to file
+    plt.savefig("${count_table.baseName}_tcr_convergence_histogram.png", dpi=300, bbox_inches="tight")
+    plt.close()
+
+    EOF
+    """
+}

subworkflows/local/sample.nf

@@ -10,6 +10,7 @@ include { SAMPLE_PLOT } from '../../modules/local/sample_plot'
 include { TCRDIST3_MATRIX } from '../../modules/local/tcrdist3_matrix'
 include { TCRDIST3_PLOT } from '../../modules/local/tcrdist3_plot'
 include { OLGA } from '../../modules/local/olga'
+include { CONVERGENCE } from '../../modules/local/convergence'
 
 /*
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -69,6 +70,8 @@ workflow SAMPLE {
 
     OLGA ( sample_map )
 
+    CONVERGENCE ( sample_map )
+    
     // emit:
     // sample_stats_csv
     // v_family_csv