From 767b73bee18ff15cc9f4f450ce82d7445bcc4c7f Mon Sep 17 00:00:00 2001 From: FerriolCalvet Date: Thu, 27 Aug 2026 10:48:36 +0200 Subject: [PATCH 1/5] add plotting for dndscv -not tested --- bin/plot_selectionsideplots.py | 65 ++++++++++++++++++++++++++++++++- bin/utils_plot.py | 1 + subworkflows/local/dnds/main.nf | 1 + workflows/deepcsa.nf | 2 +- 4 files changed, 66 insertions(+), 3 deletions(-) diff --git a/bin/plot_selectionsideplots.py b/bin/plot_selectionsideplots.py index db01a109..72e6b577 100755 --- a/bin/plot_selectionsideplots.py +++ b/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) @@ -427,7 +428,44 @@ def plot_all_positive_selection(omega_truncating, 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='--') + 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 @@ -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,27 @@ 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}.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"]] + dndscv_mis["impact"] = 'missense' + dndscv_trunc = dndscv_data[["gene_name", "sample", "wnon_cv", "qtrunc_cv"]] + dndscv_trunc["impact"] = 'truncating' + dndscv_ind = dndscv_data[["gene_name", "sample", "wind_cv", "qind_cv"]] + dndscv_ind["impact"] = 'indel' + dndscv_df = pd.concat([dndscv_trunc, dndscv_mis, dndscv_ind], axis=1) + dndscv_df.columns = ["GENE", "SAMPLE", "dndscv", "pvalue"] + + 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 +758,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, diff --git a/bin/utils_plot.py b/bin/utils_plot.py index ec78f6fe..699090c4 100644 --- a/bin/utils_plot.py +++ b/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", diff --git a/subworkflows/local/dnds/main.nf b/subworkflows/local/dnds/main.nf index 0f137f5c..e1bdd892 100644 --- a/subworkflows/local/dnds/main.nf +++ b/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 diff --git a/workflows/deepcsa.nf b/workflows/deepcsa.nf index 6d12580a..5f406245 100644 --- a/workflows/deepcsa.nf +++ b/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){ From 8b562eb1daabb48b5348f775ddd4aa0d987fb69c Mon Sep 17 00:00:00 2001 From: FerriolCalvet Date: Thu, 27 Aug 2026 12:22:19 +0200 Subject: [PATCH 2/5] add dndscv to plots --- bin/plot_selectionsideplots.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/bin/plot_selectionsideplots.py b/bin/plot_selectionsideplots.py index 72e6b577..5defe1d4 100755 --- a/bin/plot_selectionsideplots.py +++ b/bin/plot_selectionsideplots.py @@ -387,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 @@ -422,7 +422,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 @@ -459,7 +459,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 @@ -591,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 ): @@ -725,18 +725,23 @@ def get_all_data(sample, outdir, print(f"Warning: Indels file {indels_file} not found. Skipping indels track.") # Check and load oncodrivefml data - dndscv_file = f"{sample}.cv.tsv" + 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"]] + 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"]] + + 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"]] + + 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=1) - dndscv_df.columns = ["GENE", "SAMPLE", "dndscv", "pvalue"] + + dndscv_df = pd.concat([dndscv_trunc, dndscv_mis, dndscv_ind], axis=0) available_tracks.append("dndscv") print(f"Loaded dNdScv data from {dndscv_file}") @@ -779,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) From 1ab7965c459d0d137d9c839f5d8bac74762b137d Mon Sep 17 00:00:00 2001 From: FerriolCalvet Date: Fri, 28 Aug 2026 08:54:17 +0200 Subject: [PATCH 3/5] add expanded plot needles -not tested --- bin/plot_needles.py | 49 ++++++++++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 16 deletions(-) diff --git a/bin/plot_needles.py b/bin/plot_needles.py index 81e30b70..3241fed9 100755 --- a/bin/plot_needles.py +++ b/bin/plot_needles.py @@ -21,7 +21,7 @@ -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') @@ -29,14 +29,22 @@ def get_counts_per_position_n_consequence(somatic_maf_file): & (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, @@ -50,7 +58,7 @@ def plot_count_track(count_df, 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']: @@ -91,7 +99,8 @@ def plot_count_track(count_df, ax_right.set_ylabel("Proportion of mutations") def plot_stacked_bar_track_binned(count_df, - gene_len, + gene_start, + gene_end, axes, colors_dict, ax=0, @@ -105,7 +114,8 @@ def plot_stacked_bar_track_binned(count_df, 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 @@ -114,6 +124,9 @@ def plot_stacked_bar_track_binned(count_df, bin_size: size of non-overlapping bins tick_every: show x-axis ticks every N bins """ + 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 +151,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, gene_end + 1, bin_size)) binned_df = binned_df.reindex(all_bins, fill_value=0) # Plot stacked bars @@ -188,12 +201,12 @@ def plot_stacked_bar_track_binned(count_df, 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 - bin_size, gene_end + bin_size) -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 +221,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 +243,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 +269,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...") + click.echo("Generating simple needle plots...") manager(sample_name, mut_file, o3d_seq_file, outdir) + click.echo("Generating expanded needle plots...") + manager(sample_name, mut_file, o3d_seq_file, outdir, expand_muts=True) + if __name__ == '__main__': main() From 1c3e78440480fd0f255b4c475d68cb6893eb4e4e Mon Sep 17 00:00:00 2001 From: FerriolCalvet Date: Fri, 28 Aug 2026 16:47:12 +0200 Subject: [PATCH 4/5] fix bugs in function --- bin/plot_needles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/plot_needles.py b/bin/plot_needles.py index 3241fed9..09470760 100755 --- a/bin/plot_needles.py +++ b/bin/plot_needles.py @@ -270,10 +270,10 @@ def manager(mutations_file, outdir, expand_muts = False): @click.option('--outdir', type=click.Path(), help='Output path for plots') def main(sample_name, mut_file, o3d_seq_file, outdir): click.echo("Generating simple needle plots...") - manager(sample_name, mut_file, o3d_seq_file, outdir) + manager(mut_file, outdir) click.echo("Generating expanded needle plots...") - manager(sample_name, mut_file, o3d_seq_file, outdir, expand_muts=True) + manager(mut_file, outdir, expand_muts=True) if __name__ == '__main__': main() From 64ec68679c062a6fda8e23d84d34bb1176f2cc2e Mon Sep 17 00:00:00 2001 From: FerriolCalvet Date: Tue, 1 Sep 2026 19:38:58 +0200 Subject: [PATCH 5/5] minor fix in stacked plotting --- bin/plot_needles.py | 41 ++++++++++++++++++++++++++--------------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/bin/plot_needles.py b/bin/plot_needles.py index 09470760..a7c45ad4 100755 --- a/bin/plot_needles.py +++ b/bin/plot_needles.py @@ -52,7 +52,24 @@ def plot_count_track(count_df, 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) @@ -98,6 +115,8 @@ 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_start, gene_end, @@ -107,8 +126,7 @@ def plot_stacked_bar_track_binned(count_df, 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. @@ -121,8 +139,8 @@ def plot_stacked_bar_track_binned(count_df, 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" @@ -151,7 +169,7 @@ def plot_stacked_bar_track_binned(count_df, ) # Ensure all bins are represented - all_bins = list(range(gene_start, gene_end + 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 @@ -194,14 +212,7 @@ 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(gene_start - bin_size, gene_end + bin_size) + axes[ax].set_xlim(gene_start - gene_len * 0.02, gene_end + gene_len * 0.02) def manager(mutations_file, outdir, expand_muts = False): @@ -227,7 +238,7 @@ def manager(mutations_file, outdir, expand_muts = False): colors_dict=metrics_colors_dictionary, indel=False, alpha=0.7 ) ax.set_title(f"{gene}") - plt.savefig(f"{outdir}/{gene}{".expanded" if expand_muts else ''}.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() @@ -251,7 +262,7 @@ def manager(mutations_file, outdir, expand_muts = False): indel=False ) ax.set_title(f"{gene}") - plt.savefig(f"{outdir}/{gene}{".expanded" if expand_muts else ''}.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()