Merge pull request #30 from KarchinLab/dltamayo-dev

Implement TCR convergence and plots, code cleanup

Author: favorov

.cirro/preprocess.py

@@ -17,26 +17,19 @@
 ds.logger.info("Checking samplesheet parameter")
 ds.logger.info(ds.samplesheet)
 samplesheet = ds.samplesheet
-
-ds.logger.info("Dropping incorrect file path & Merging ds.files w samplesheets")
-samplesheet = samplesheet.drop(columns=['file'])
-samplesheet2 = samplesheet.merge(ds.files, on='sample', how='left')
-
-samplesheet2.to_csv('samplesheet.csv', index=None)
+samplesheet.to_csv('samplesheet.csv', index=None)
 ds.add_param("samplesheet", "samplesheet.csv")
 
 
 # 3. Set workflow_level value based on form input
 ds.logger.info("Setting workflow_level")
-if ds.params['sample_lvl'] == ds.params['compare_lvl'] == ds.params['cluster_lvl'] == True:
+if ds.params['sample_lvl'] == ds.params['compare_lvl'] == True:
     workflow_level = ['complete']
 else:
-    workflow_lvls = ['sample', 'compare', 'cluster']
-    chosen_lvls = [ds.params['sample_lvl'], ds.params['compare_lvl'], ds.params['cluster_lvl']]
+    workflow_lvls = ['sample', 'compare']
+    chosen_lvls = [ds.params['sample_lvl'], ds.params['compare_lvl']]
     workflow_level = [i for i, j in zip(workflow_lvls, chosen_lvls) if j]
 
 ds.add_param('workflow_level', ','.join(workflow_level))
 
 ds.logger.info(ds.params)
-
-## 

modules/local/compare_clonal_publicity.nf

@@ -6,12 +6,99 @@ process COMPARE_CLONAL_PUBLICITY {
     path concat_cdr3
 
     output:
-    path "cdr3_sharing.tsv", emit: "shared_cdr3"
+    path "cdr3_sharing_pgen.tsv", emit: "shared_cdr3"
     path "sample_mapping.tsv", emit: "sample_mapping"
+    path "sharing_histogram.png"
+    path "sharing_pgen_scatterplot.png"
 
     script:
     """
-    # Concatenate input Adaptive files and process metadata
-    compare_clonal_publicity.py $concat_cdr3
+    python - <<EOF
+    import pandas as pd
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    # Load data
+    df = pd.read_csv("${concat_cdr3}", sep="\t")
+
+    # Step 1: Map samples to integers
+    sample_mapping = {sample: i + 1 for i, sample in enumerate(df['sample'].unique())}
+    df['sample_id'] = df['sample'].map(sample_mapping)
+
+    # Step 2: Group by CDR3b and aggregate sample_ids
+    grouped = (
+        df.groupby('CDR3b')['sample_id']
+        .apply(lambda x: sorted(set(x)))  # remove duplicates if any
+        .reset_index()
+    )
+
+    # Step 3: Add comma-separated list and total count
+    grouped['samples_present'] = grouped['sample_id'].apply(lambda x: ",".join(map(str, x)))
+    grouped['total_samples'] = grouped['sample_id'].apply(len)
+
+    # Step 4: Final output — drop raw list
+    final_df = grouped[['CDR3b', 'total_samples', 'samples_present']]
+    final_df = final_df.sort_values(by='total_samples', axis=0, ascending=False)
+
+    # Step 5: Export both outputs
+    final_df.to_csv("cdr3_sharing.tsv", sep="\t", index=False)
+
+    # Also export the sample mapping
+    sample_map_df = pd.DataFrame.from_dict(sample_mapping, orient='index', columns=['sample_id']).reset_index()
+    sample_map_df.columns = ['patient', 'sample_id']
+    sample_map_df.to_csv("sample_mapping.tsv", sep="\t", index=False)
+
+
+    # Plot histogram
+    sharing = final_df['total_samples'].values
+
+    # Create integer bin edges from 0 to max(data)
+    bins = np.arange(min(sharing), max(sharing) + 2)  # +2 to include the last value as a bin edge
+
+    plt.figure(figsize=(8, 5))
+    plt.hist(sharing, bins=bins, edgecolor='black', align='left')
+    plt.xticks(bins[:-1])  # whole number positions and labels
+    plt.yscale('log')
+
+    plt.xlabel('Number of Shared Samples')
+    plt.ylabel('TCR Sequence Frequency (log scale)')
+    plt.title('TCR Sharing Histogram')
+
+    # Save to file
+    plt.savefig("sharing_histogram.png", dpi=300, bbox_inches="tight")
+    plt.close()
+    EOF
+
+    olga-compute_pgen --humanTRB -i cdr3_sharing.tsv -o pgen_sharing.tsv
+
+    python - <<EOF
+    import pandas as pd
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from matplotlib.ticker import MaxNLocator
+
+    # Load TSVs for shared cdr3s and corresponding pgen values
+    left_df = pd.read_csv('pgen_sharing.tsv', sep='\t', header=None, usecols=[0, 1], names=['CDR3b', 'pgen'])
+    right_df = pd.read_csv('cdr3_sharing.tsv', sep='\t')
+
+    # Drop rows where pgen == 0 and merge
+    left_df = left_df[left_df['pgen'] != 0]
+    merged_df = pd.merge(left_df, right_df, on='CDR3b', how='left')
+    merged_df.to_csv('cdr3_sharing_pgen.tsv', sep='\t', index=False)
+
+    # Create scatter plot with log-transform pgen
+    merged_df["log10_pgen"] = np.log10(merged_df["pgen"])
+    plt.figure(figsize=(8, 6))
+    plt.grid(True)
+    plt.scatter(merged_df["log10_pgen"], merged_df["total_samples"], c='blue', alpha=0.7)
+    plt.gca().yaxis.set_major_locator(MaxNLocator(integer=True))
+
+    plt.xlabel("log10(Probability)")
+    plt.ylabel("Number of Shared Samples")
+    plt.title("Scatterplot of Shared TCRs vs log10(Generation Probability)")
+    plt.tight_layout()
+    plt.savefig("sharing_pgen_scatterplot.png", dpi=300, bbox_inches="tight")
+    plt.close()
+    EOF
     """
 }

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
+    """
+}

modules/local/olga.nf

@@ -2,30 +2,55 @@ process OLGA {
     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_generation_probabilities.tsv", emit: "olga_output"
-    
+    path "${count_table.baseName}_tcr_pgen.tsv", emit: "olga_output"
+    path "${count_table.baseName}_tcr_pgen_histogram.png"
+
     script:
     """
     # Extract vector of cdr3 aa, dropping null values
-    
     cat > dropAA.py <<EOF
-    
     import pandas as pd
-    
+
     df = pd.read_csv("${count_table}", sep="\t")
     df = df.dropna(subset=["aminoAcid"])
     df = df["aminoAcid"]
     df.to_csv("output.tsv", sep="\t", index=False, header=False)
-    
     EOF
-    
+
     python dropAA.py
-    
-    olga-compute_pgen --humanTRB -i output.tsv -o "${count_table.baseName}_tcr_generation_probabilities.tsv"
+
+    olga-compute_pgen --humanTRB -i output.tsv -o "${count_table.baseName}_tcr_pgen.tsv"
+
+    python - <<EOF
+    import pandas as pd
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    # Load TSV with no header
+    df = pd.read_csv('${count_table.baseName}_tcr_pgen.tsv', sep='\t', header=None, usecols=[0, 1], names=['CDR3b', 'probability'])
+    
+    # Drop rows where pgen is 0
+    df = df[df['probability'] != 0]
+    log_probs = np.log10(df['probability'])
+
+    # Plot histogram
+    plt.figure(figsize=(8, 5))
+    plt.hist(log_probs, bins=30, density=True, edgecolor='black')
+
+    # Label with LaTeX formatting
+    plt.xlabel('log_10 Generation Probability')
+    plt.ylabel('Probability Density')
+    plt.title(f'${count_table.baseName} TCR Generation Probability Histogram')
+    # plt.grid(True)
+
+    # Save to file
+    plt.savefig("${count_table.baseName}_tcr_pgen_histogram.png", dpi=300, bbox_inches="tight")
+    plt.close()
+    EOF
     """
-}
+}

modules/local/tcrdist3_matrix.nf

@@ -9,7 +9,7 @@ process TCRDIST3_MATRIX {
     path ref_db
 
     output:
-    path "${count_table.baseName}_distance_matrix.csv", emit: 'distance_matrix'
+    tuple val(sample_meta), path("${count_table.baseName}_distance_matrix.csv"), emit: 'tcr_output'
     path "${count_table.baseName}_clone_df.csv", emit: 'clone_df'
 
     script:

modules/local/tcrdist3_plot.nf

@@ -0,0 +1,33 @@
+process TCRDIST3_PLOT {
+    tag "${sample_meta[0]}"
+    label 'process_high'
+    label 'process_high_memory'
+    container "ghcr.io/break-through-cancer/bulktcr:latest"
+
+    input:
+    tuple val(sample_meta), path(distance_matrix)
+
+    output:
+    path "${sample_meta[0]}_pairwise_distance_distribution.png", emit: 'beta_histogram'
+
+    script:
+    """
+    python - <<EOF
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    distances = np.loadtxt("${distance_matrix}", delimiter=',')
+    lower_triangle = distances[np.tril_indices(distances.shape[0], k=-1)]
+    counts, bin_edges = np.histogram(lower_triangle, bins=100)
+
+    plt.figure(figsize=(8, 5))
+    plt.bar((bin_edges[:-1] + bin_edges[1:]) / 2, counts, width=np.diff(bin_edges), edgecolor='black')
+    plt.xlabel("Pairwise Distance")
+    plt.ylabel("Frequency (log scale)")
+    plt.yscale("log")
+    plt.title("Distribution of Beta Chain Pairwise Distances - ${sample_meta[0]}")
+    plt.savefig("${sample_meta[0]}_pairwise_distance_distribution.png", dpi=300, bbox_inches="tight")
+    plt.close()
+    EOF
+    """
+}

subworkflows/local/sample.nf

@@ -8,7 +8,9 @@
 include { SAMPLE_CALC } from '../../modules/local/sample_calc'
 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'
 
 /*
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -53,6 +55,10 @@ workflow SAMPLE {
         file(params.db_path)
     )
 
+    TCRDIST3_PLOT(
+        TCRDIST3_MATRIX.out.tcr_output
+    )
+
     /////// =================== PLOT SAMPLE ===================  ///////
 
     SAMPLE_PLOT (
@@ -64,6 +70,8 @@ workflow SAMPLE {
 
     OLGA ( sample_map )
 
+    CONVERGENCE ( sample_map )
+    
     // emit:
     // sample_stats_csv
     // v_family_csv