Add dndscv plots and expanded needle plots

bin/plot_needles.py

@@ -21,36 +21,61 @@
 
 
 
-def get_counts_per_position_n_consequence(somatic_maf_file):
+def get_counts_per_position_n_consequence(somatic_maf_file, expanded = False):
     somatic_maf = pd.read_table(somatic_maf_file, na_values = custom_na_values)
 
     somatic_maf_clean = somatic_maf[(somatic_maf["TYPE"] == 'SNV')
                                     & (~somatic_maf["FILTER.not_in_exons"])
                                     & (somatic_maf['canonical_Protein_position'] != '-')
                                     ].reset_index(drop = True)
     somatic_maf_clean['canonical_Protein_position'] = somatic_maf_clean['canonical_Protein_position'].astype(int)
-    counts_per_position = somatic_maf_clean.groupby(by = ["SAMPLE_ID", "canonical_SYMBOL", 'canonical_Consequence_broader', 'canonical_Protein_position'])['ALT_DEPTH'].size().to_frame('Count').reset_index()
+    if expanded:
+        somatic_maf_clean.loc[somatic_maf_clean["ALT_DEPTH"] > 5, "ALT_DEPTH"] = 5
+        counts_per_position = somatic_maf_clean.groupby(by = ["SAMPLE_ID", "canonical_SYMBOL", 'canonical_Consequence_broader', 'canonical_Protein_position'])['ALT_DEPTH'].sum().to_frame('Count').reset_index()
+        click.echo("Using expanded counts")
+
+    else:
+        counts_per_position = somatic_maf_clean.groupby(by = ["SAMPLE_ID", "canonical_SYMBOL", 'canonical_Consequence_broader', 'canonical_Protein_position'])['ALT_DEPTH'].size().to_frame('Count').reset_index()
+
     counts_per_position.columns = ["SAMPLE_ID", 'Gene', 'Consequence', 'Pos', 'Count']
 
     return counts_per_position
 
 
 def plot_count_track(count_df,
-                        gene_len,
+                     gene_start,
+                     gene_end,
                         axes,
                         colors_dict,
                         ax=0,
                         alpha=1,
                         indel=False,
                         n_batches = 10
                     ):
+    """
+    Plots needle plot of mutation counts by position.
 
+    Parameters:
+        count_df: DataFrame with ['Pos', 'Consequence', 'Count'] columns
+        gene_start: Start position of the gene
+        gene_end: End position of the gene
+        axes: matplotlib axes array
+        colors_dict: dictionary mapping consequence -> color
+        ax: index of subplot
+        alpha: transparency
+        indel: whether to include 'indel' consequence
+        n_batches: number of batches to split the data into for plotting
+
+    """
+    assert gene_start < gene_end, "gene_start must be less than gene_end"
+
+    gene_len = gene_end - gene_start + 1
     # Shuffle the data and split into batches
     shuffled_df = count_df.sample(frac=1, random_state=42).reset_index(drop=True)
     batches = np.array_split(shuffled_df, n_batches)
 
     legend_list = []
-    pos_df = pd.DataFrame({"Pos" : range(1, gene_len+1)})
+    pos_df = pd.DataFrame({"Pos" : range(gene_start, gene_end + 1)})
 
     for batch_idx, batch in enumerate(batches):
         for cnsq in ['nonsense', 'missense', 'synonymous']:
@@ -90,30 +115,36 @@ def plot_count_track(count_df,
         ax_right.set_yticklabels(["0.00"] * len(axes[ax].get_yticks()))
         ax_right.set_ylabel("Proportion of mutations")
 
+    axes[ax].set_xlim(gene_start - gene_len * 0.02, gene_end + gene_len * 0.02)
+
 def plot_stacked_bar_track_binned(count_df,
-                                    gene_len,
+                                    gene_start,
+                                    gene_end,
                                     axes,
                                     colors_dict,
                                     ax=0,
                                     alpha=1,
                                     indel=False,
                                     min_bin_size=3,
-                                    num_bins = 100,
-                                    num_ticks=5):
+                                    num_bins = 100):
     """
     Plots stacked barplot of mutation counts binned by position.
 
     Parameters:
         count_df: DataFrame with ['Pos', 'Consequence', 'Count'] columns
-        gene_len: Length of the protein sequence
+        gene_start: Start position of the gene
+        gene_end: End position of the gene
         axes: matplotlib axes array
         colors_dict: dictionary mapping consequence -> color
         ax: index of subplot
         alpha: transparency
         indel: whether to include 'indel' consequence
-        bin_size: size of non-overlapping bins
-        tick_every: show x-axis ticks every N bins
+        min_bin_size: minimum size of non-overlapping bins
+        num_bins: number of bins to create
     """
+    assert gene_start < gene_end, "gene_start must be less than gene_end"
+
+    gene_len = gene_end - gene_start + 1
 
     # Compute bin_size or result to default
     candidate_bin_size = max(1, gene_len // num_bins)
@@ -138,7 +169,7 @@ def plot_stacked_bar_track_binned(count_df,
     )
 
     # Ensure all bins are represented
-    all_bins = list(range(1, gene_len + 1, bin_size))
+    all_bins = list(range((gene_start-1)//bin_size * bin_size + 1, gene_end + 1, bin_size))
     binned_df = binned_df.reindex(all_bins, fill_value=0)
 
     # Plot stacked bars
@@ -181,19 +212,12 @@ def plot_stacked_bar_track_binned(count_df,
         ax_right.set_yticklabels(["0.00"] * len(axes[ax].get_yticks()))
         ax_right.set_ylabel("Proportion of mutations")
 
-    # Sparse x-ticks
-    tick_every = len(all_bins) // num_ticks
-    sparse_ticks = all_bins[::tick_every]
-    sparse_ticks = [x-1 for x in sparse_ticks]
-
-    axes[ax].set_xticks(sparse_ticks)
-    axes[ax].set_xticklabels(sparse_ticks)
-    axes[ax].set_xlim(0, gene_len + bin_size)
+    axes[ax].set_xlim(gene_start - gene_len * 0.02, gene_end + gene_len * 0.02)
 
 
-def manager(sample_name, mutations_file, o3d_seq_file, outdir):
+def manager(mutations_file, outdir, expand_muts = False):
 
-    counts_per_position = get_counts_per_position_n_consequence(mutations_file)
+    counts_per_position = get_counts_per_position_n_consequence(mutations_file, expanded=expand_muts)
 
     gene_order = sorted(pd.unique(counts_per_position["Gene"]))
 
@@ -208,13 +232,13 @@ def manager(sample_name, mutations_file, o3d_seq_file, outdir):
             fig, ax = plt.subplots(1, 1, figsize=(5, 1.2))
             plot_count_track(
                 mut_count_df,
-                gene_len=mut_count_df["Pos"].max(),  # FIXME: this is not ideal, the max position is the biggest position with mutation
-
+                gene_start=mut_count_df["Pos"].min(),
+                gene_end=mut_count_df["Pos"].max(),
                 axes=[ax], ax=0,
                 colors_dict=metrics_colors_dictionary, indel=False, alpha=0.7
             )
             ax.set_title(f"{gene}")
-            plt.savefig(f"{outdir}/{gene}.needle.pdf", bbox_inches='tight', dpi=100)
+            plt.savefig(f"{outdir}/{gene}{'_expanded' if expand_muts else ''}.needle.pdf", bbox_inches='tight', dpi=100)
             plt.show()
             plt.close()
 
@@ -230,14 +254,15 @@ def manager(sample_name, mutations_file, o3d_seq_file, outdir):
             fig, ax = plt.subplots(1, 1, figsize=(5, 1.2))
             plot_stacked_bar_track_binned(
                 count_df=mut_count_df,
-                gene_len=mut_count_df["Pos"].max(),  # FIXME: this is not ideal, the max position is the biggest position with mutation
+                gene_start=mut_count_df["Pos"].min(),
+                gene_end=mut_count_df["Pos"].max(),
                 axes=[ax], ax=0,
                 colors_dict=metrics_colors_dictionary,
                 alpha=1,
                 indel=False
             )
             ax.set_title(f"{gene}")
-            plt.savefig(f"{outdir}/{gene}.stacked.pdf", bbox_inches='tight', dpi=100)
+            plt.savefig(f"{outdir}/{gene}{'_expanded' if expand_muts else ''}.stacked.pdf", bbox_inches='tight', dpi=100)
             plt.show()
             plt.close()
 
@@ -255,8 +280,11 @@ def manager(sample_name, mutations_file, o3d_seq_file, outdir):
 @click.option('--o3d_seq_file', type=click.Path(exists=True), help='Input Oncodrive3D sequence df file')
 @click.option('--outdir', type=click.Path(), help='Output path for plots')
 def main(sample_name, mut_file, o3d_seq_file, outdir):
-    click.echo("Plotting omega results...")
-    manager(sample_name, mut_file, o3d_seq_file, outdir)
+    click.echo("Generating simple needle plots...")
+    manager(mut_file, outdir)
+
+    click.echo("Generating expanded needle plots...")
+    manager(mut_file, outdir, expand_muts=True)
 
 if __name__ == '__main__':
     main()

bin/plot_selectionsideplots.py

@@ -333,16 +333,17 @@ def plot_all_positive_selection(omega_truncating,
                                 indels_panel_df,
                                 oncodrive3d_data_scores,
                                 oncodrivefml_data,
+                                dndscv_data,
                                 gene_order,
                                 title = None,
                                 pvalue_thres = 0.05,
                                 linewidth_def = 0.6,
-                                tracks = ("omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml")
+                                tracks = ("omega_trunc", "omega_mis", "dndscv", "oncodrive3d", "oncodrivefml")
                                 ):
 
     num_genes = len(gene_order)
     # Determine which tracks to plot and their order
-    all_tracks = ["omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml", "indels"]
+    all_tracks = ["omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml", "dndscv", "indels"]
     plot_tracks = [t for t in all_tracks if t in tracks]
     n_tracks = len(plot_tracks)
 
@@ -386,7 +387,7 @@ def plot_all_positive_selection(omega_truncating,
         else:
             ax.set_xticks(range(len(gene_order)))
             ax.set_xticklabels([])
-        ax.axhline(1, color='black', linestyle='--')
+        ax.axhline(1, color='black', linestyle='--', linewidth=0.5)
         ax.spines['top'].set_visible(False)
         ax.spines['right'].set_visible(False)
         ax_idx += 1
@@ -421,13 +422,50 @@ def plot_all_positive_selection(omega_truncating,
         else:
             ax.set_xticks(range(len(gene_order)))
             ax.set_xticklabels([])
-        ax.axhline(1, color='black', linestyle='--')
+        ax.axhline(1, color='black', linestyle='--', linewidth=0.5)
         ax.spines['top'].set_visible(False)
         ax.spines['right'].set_visible(False)
         ax_idx += 1
 
 
+    if "dndscv" in plot_tracks and dndscv_data is not None:
+        ax = axes[ax_idx]
+        dndscv_sig = dndscv_data[dndscv_data["pvalue"] <= pvalue_thres].reset_index(drop = True)
+        dndscv_notsig = dndscv_data[dndscv_data["pvalue"] > pvalue_thres].reset_index(drop = True)
+        sns.barplot(data=dndscv_notsig, x='GENE', y='dndscv',
+                    ax=ax, alpha=1,
+                    fill = False,
+                    legend = False,
+                    linewidth = linewidth_def,
+                    order = gene_order,
+                    hue = 'impact',
+                    hue_order = ["missense", "truncating", "indel"],
+                    palette = metrics_colors_dictionary)
+        sns.barplot(data=dndscv_sig, x='GENE', y='dndscv',
+                    ax=ax, alpha=1,
+                    legend = False,
+                    linewidth = linewidth_def,
+                    order = gene_order,
+                    hue = 'impact',
+                    hue_order = ["missense", "truncating", "indel"],
+                    palette = metrics_colors_dictionary
+                    )
+        ax.set_xlabel('')
+        ax.set_ylabel('dNdScv', rotation = 0, labelpad=17, verticalalignment = 'center')
+        # Only set xticklabels on last axis
+        if ax_idx == n_tracks - 1:
+            ax.set_xticks(range(num_genes))
+            ax.set_xticklabels(gene_order, rotation=90)
+        else:
+            ax.set_xticks(range(len(gene_order)))
+            ax.set_xticklabels([])
+        ax.axhline(1, color='black', linestyle='--', linewidth=0.5)
+        ax.spines['top'].set_visible(False)
+        ax.spines['right'].set_visible(False)
+        ax_idx += 1
+
 
+
     if "oncodrive3d" in plot_tracks and oncodrive3d_data_scores is not None:
         ax = axes[ax_idx]
         df = oncodrive3d_data_scores
@@ -553,7 +591,7 @@ def plot_all_positive_selection(omega_truncating,
 
 def get_all_data(sample, outdir,
                  pvaluee = 0.05,
-                 tracks=("omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml"),
+                 tracks=("omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml", "dndscv"),
                  gene_order = None
                  ):
 
@@ -563,6 +601,7 @@ def get_all_data(sample, outdir,
     omega_missense = None
     oncodrive3d_data_scores = None
     indels_panel_df = None
+    dndscv_df = None
     global_omega_decreasing = []
     available_tracks = []
 
@@ -685,6 +724,32 @@ def get_all_data(sample, outdir,
     else:
         print(f"Warning: Indels file {indels_file} not found. Skipping indels track.")
 
+    # Check and load oncodrivefml data
+    dndscv_file = f"{sample}.dNdScv.cv.tsv"
+    if os.path.exists(dndscv_file) and "dndscv" in tracks:
+        try:
+            dndscv_data = pd.read_table(dndscv_file)
+            dndscv_mis = dndscv_data[["gene_name", "sample", "wmis_cv", "qmis_cv"]].copy()
+            dndscv_mis.columns = ["GENE", "SAMPLE", "dndscv", "pvalue"]
+            dndscv_mis["impact"] = 'missense'
+
+            dndscv_trunc = dndscv_data[["gene_name", "sample", "wnon_cv", "qtrunc_cv"]].copy()
+            dndscv_trunc.columns = ["GENE", "SAMPLE", "dndscv", "pvalue"]
+            dndscv_trunc["impact"] = 'truncating'
+
+            dndscv_ind = dndscv_data[["gene_name", "sample", "wind_cv", "qind_cv"]].copy()
+            dndscv_ind.columns = ["GENE", "SAMPLE", "dndscv", "pvalue"]
+            dndscv_ind["impact"] = 'indel'
+
+            dndscv_df = pd.concat([dndscv_trunc, dndscv_mis, dndscv_ind], axis=0)
+
+            available_tracks.append("dndscv")
+            print(f"Loaded dNdScv data from {dndscv_file}")
+        except Exception as e:
+            print(f"Warning: Failed to load dNdScv data: {e}")
+    else:
+        print(f"Warning: dNdScv file {dndscv_file} not found. Skipping dNdScv track.")
+
     # Check if we have any data to plot
     if not available_tracks:
         print("Warning: No data files found for any of the requested tracks. Skipping plot generation.")
@@ -698,6 +763,7 @@ def get_all_data(sample, outdir,
                                             indels_panel_df,
                                             oncodrive3d_data_scores,
                                             oncodrivefml_data,
+                                            dndscv_df,
                                             global_omega_decreasing if gene_order is None else gene_order,
                                             title = sample,
                                             pvalue_thres = pvaluee,
@@ -718,7 +784,7 @@ def main(sample_name, outdir, include_indels):
     try:
         generate_all_side_figures(sample_name, f"{outdir}/side_figures")
         # By default, exclude indels unless --include_indels is set
-        tracks = ("omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml")
+        tracks = ("omega_trunc", "omega_mis", "oncodrive3d", "oncodrivefml", "dndscv")
         if include_indels:
             tracks = tracks + ("indels",)
         get_all_data(sample_name, outdir, tracks=tracks)

bin/utils_plot.py

@@ -143,6 +143,7 @@ def plot_profile(frequencies, title='Mutational profile',
                                 "omega_miss_tert" : "#f5840c",
                                 "omega_synon_tert": "#378c12",
                                 "nonsense" : "#FA5E32",
+                                "truncating" : "#FA5E32",                                
                                 "synonymous" : "#89E4A2",
                                 "missense"  : "#FABE4A",
                                 "indel"       : "#ECC4F7",

subworkflows/local/dnds/main.nf

@@ -45,6 +45,7 @@ workflow DNDS {
     local_results.collectFile(name: "all_dNdScv.local.tsv", storeDir:"${params.outdir}/selection/dndscv/local", skip: 1, keepHeader: true).set{ all_dndscv_local_results }
 
     emit:
+    dnds_cv_per_sample = DNDSRUN.out.results_cv
     all_dndscv_results
     all_dndscv_global_results
     all_dndscv_local_results

workflows/deepcsa.nf

@@ -467,7 +467,6 @@ workflow DEEPCSA {
         positive_selection_results = positive_selection_results.join(ONCODRIVE3D.out.results_pos, remainder: true)
     }
 
-    // if (params.expected_mutated_cells & params.dnds){
     if (params.dnds){
         DNDS(mutations_in_exons,
                     DEPTHSEXONSCONS.out.subset,
@@ -476,6 +475,7 @@ workflow DEEPCSA {
                     params.fasta
                     )
         dndscv_table = DNDS.out.all_dndscv_results
+        positive_selection_results = positive_selection_results.join(DNDS.out.dnds_cv_per_sample, remainder: true)
     }
 
     if (params.omega){