From e6a0f5059964b277ec64e87162e2c738db21b1d3 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Thu, 2 Jul 2026 12:40:05 +0200 Subject: [PATCH 01/14] add draft for omega vs dndscv plot qc --- bin/omega_vs_dndscv_qc.py | 155 ++++++++++++++++++ modules/local/plot/qc/omega_vs_dndscv/main.nf | 44 +++++ subworkflows/local/plotting_qc/main.nf | 4 + 3 files changed, 203 insertions(+) create mode 100644 bin/omega_vs_dndscv_qc.py create mode 100644 modules/local/plot/qc/omega_vs_dndscv/main.nf diff --git a/bin/omega_vs_dndscv_qc.py b/bin/omega_vs_dndscv_qc.py new file mode 100644 index 00000000..6f9d5f50 --- /dev/null +++ b/bin/omega_vs_dndscv_qc.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python + +# Needed basic packages +import click +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np +import seaborn as sns +import os +import time +import json +from scipy import stats +from scipy.stats import pearsonr + + +# Functions +def process_dndscv_table(dndscv_df): + + # Separate missense and truncating variants for dNdScv table in two separated tables + missense_df = dndscv_df[['sample', 'gene_name', 'n_mis', 'wmis_cv', 'qmis_cv']].copy() + missense_df = missense_df.rename(columns={'n_mis': 'mutations', 'wmis_cv': 'dnds', 'qmis_cv': 'pvalue', 'gene_name': 'gene'}) + missense_df['impact'] = 'missense' + + truncating_df = dndscv_df[['sample', 'gene_name', 'n_non', 'n_spl', 'wnon_cv', 'qtrunc_cv']].copy() + truncating_df['n_trunc'] = truncating_df['n_non'] + truncating_df['n_spl'] + truncating_df = truncating_df.rename(columns={'n_trunc': 'mutations', 'wnon_cv': 'dnds', 'qtrunc_cv': 'pvalue', 'gene_name': 'gene'}) + truncating_df['impact'] = 'truncating' + + # Then concat the two tables + dndscv_cv_df = pd.concat([missense_df, truncating_df.drop(columns=['n_non', 'n_spl'])], axis=0) + + return dndscv_cv_df + + +def filter_flagged_genes_per_group(df, flagged_cases_df): + # Optimized: Vectorized filtering using a multi-index instead of .apply row-by-row + flagged_index = pd.MultiIndex.from_frame(flagged_cases_df[['cohort', 'gene']].rename(columns={'cohort': 'sample'})) + df_index = pd.MultiIndex.from_frame(df[['sample', 'gene']]) + + filtered_df = df[~df_index.isin(flagged_index)].copy() + return filtered_df + + +def apply_correlation_and_plotting(df, samples_group, output_dir): + + os.makedirs(output_dir, exist_ok=True) # Ensure output dir exists + + for group in sorted(samples_group): + subset = df[df['sample'] == group].copy() + subset[['dnds_omega', 'dnds_dndscv']] = subset[['dnds_omega', 'dnds_dndscv']].fillna(0) + print(subset.head()) + + if subset.empty: + print(f'No data available for sample: {group}') + continue + + fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharex=False, sharey=False) + + impacts = ['missense', 'truncating'] + + for i, impact in enumerate(impacts): + ax = axes[i] + impact_subset = subset[subset['impact'] == impact] + impact_subset_reg = impact_subset[(impact_subset['dnds_omega'] > 0) & (impact_subset['dnds_dndscv'] > 0)] + + if impact_subset.empty: + ax.set_title(f'No {impact} data for {group}') + continue + + # Scatter plot + sns.scatterplot(data=impact_subset, x='dnds_omega', y='dnds_dndscv', hue='gene', ax=ax) + + # Regression line + sns.regplot(data=impact_subset_reg, x='dnds_omega', y='dnds_dndscv', scatter=False, color='red', label='Regression Line', ax=ax) + + # X=Y line + all_vals = pd.concat([impact_subset['dnds_omega'], impact_subset['dnds_dndscv']]) + min_val, max_val = all_vals.min(), all_vals.max() + ax.plot([min_val, max_val], [min_val, max_val], color='black', linestyle='--', label='x=y') + + if impact_subset_reg.shape[0] > 2: + # Compute correlation + corr, pval = pearsonr(impact_subset_reg['dnds_omega'], impact_subset_reg['dnds_dndscv']) + + ax.set_title(f'{impact.capitalize()} mutations\nPearson R: {corr:.3f} (p={pval:.3e})') + else: + ax.set_title(f'{impact.capitalize()} mutations\nNot enough data for correlation') + + ax.set_xlabel('dN/dS (Omega)') + ax.set_ylabel('dN/dS (dNdScv)') + ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + ax.legend().remove() + + fig.suptitle(f'Comparison for {group}') + fig.tight_layout() + + output_path = os.path.join(output_dir, f"{group}_omega_vs_dndscv_qc_summary_plot.pdf") + with PdfPages(output_path) as pdf: + pdf.savefig(fig) + + plt.close(fig) + + return + + +@click.command() +@click.option("--input-omega-file", required=True, type=click.Path(exists=True), + help="Directory containing selection/omega/all_omegas.tsv") + +@click.option("--input-dndscv-file", required=True, type=click.Path(exists=True), + help="Directory containing selection/dndscv/cv/all_dNdScv.cv.tsv") + +@click.option("--output-dir", required=True, type=click.Path(writable=True), + help="Directory where output files will be written") + +@click.option("--flagged-genes-omega", required=True, type=click.Path(writable=True), + help="Directory where flagged_omega genes were stored from /qc/omega_flagged/debug.syn_flagged_gene.tsv") + + +# Main function +def main(input_omega_file, input_dndscv_file, output_dir, flagged_genes_omega): + + # Read tables + omega_df = pd.read_table(input_omega_file, sep='\t') + dndscv_df = pd.read_table(input_dndscv_file, sep='\t') + flagged_genes_df = pd.read_table(flagged_genes_omega, sep='\t') + + # Extract sample groups + sample_groups = flagged_genes_df['cohort'].unique().tolist() + + # Apply basic filter for omega table + omega_filtered_df = omega_df[~(omega_df['gene'].str.contains("--")) & # discards exons if they are included in the analysis + (omega_df['impact'].isin(['missense', 'truncating']))].copy() + print('Filtered omega table:') + print(omega_filtered_df.head()) + + # Process dNdScv table + dndscv_cv_df = process_dndscv_table(dndscv_df) + print('Processed dNdScv table:') + print(dndscv_cv_df.head()) + + # Merge omega and dNdScv tables + dndscv_n_omegas_df = omega_filtered_df.merge(dndscv_cv_df, on=['sample', 'gene', 'impact'], how='outer', suffixes=('_omega', '_dndscv')) + print('Merged omega and dNdScv table:') + print(dndscv_n_omegas_df.head()) + + # Filter flagged omega genes within groups + filtered_dndscv_n_omegas_df = filter_flagged_genes_per_group(dndscv_n_omegas_df, flagged_genes_df) + + # Apply function to plot omega values vs dndscv per group and compute correlation + apply_correlation_and_plotting(filtered_dndscv_n_omegas_df, sample_groups, output_dir) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/modules/local/plot/qc/omega_vs_dndscv/main.nf b/modules/local/plot/qc/omega_vs_dndscv/main.nf new file mode 100644 index 00000000..fccf4b78 --- /dev/null +++ b/modules/local/plot/qc/omega_vs_dndscv/main.nf @@ -0,0 +1,44 @@ +process PLOT_OMEGA_VS_DNDSCV { + + tag "${group_name}" + label 'process_low' + + label 'deepcsa_core' + + input: + path (all_omegas) + path (groups_json) + val (group_name) + path (compiled_flagged) + + output: + path("**.pdf") , optional: true , emit: plots + path "versions.yml" , topic: versions + + script: + """ + mkdir ${group_name}.plots + omega_vs_dndscv_qc.py \\ + --input-omega-file ${all_mutdensities} \\ + --input-dndscv-file ${} \\ + --output-dir ${group_name}.plots \\ + --flagged-genes-omega ${compiled_flagged} + + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "all_samples" + """ + touch ${prefix}.pdf + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ +} diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index ab183920..1f3eb9a5 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -3,6 +3,7 @@ include { PLOT_MUTDENSITY_QC as PLOTMUTDENSITYQC } from '../../../mod include { PLOT_METRICS_VS_DEPTH_QC as PLOTMETRICSVSDEPTHQC } from '../../../modules/local/plot/qc/metrics_vs_depth/main' include { ANNOTATE_OMEGA_QC as APPLYOMEGAQC } from '../../../modules/local/plot/qc/annotate_omega/main' include { PLOT_MUTATION_SPECIFIC as PLOTMUTATIONSPECIFIC } from '../../../modules/local/plot/qc/mutation_specific/main' +include { PLOT_OMEGA_VS_DNDSCV as PLOTOMEGAVSDNDSCV } from '../../../modules/local/plot/qc/omega_vs_dndscv/main' workflow PLOTTING_QC { @@ -49,6 +50,8 @@ workflow PLOTTING_QC { all_omegas_globalloc ) + PLOTOMEGAVSDNDSCV(all_omegas, groups_definition, group_name) + APPLYOMEGAQC(all_omegas, PLOTMUTDENSITYQC.out.compiled_flagged.collect()) emit: @@ -56,5 +59,6 @@ workflow PLOTTING_QC { metrics_vs_depth_plots = PLOTMETRICSVSDEPTHQC.out.plots metrics_vs_depth_tables = PLOTMETRICSVSDEPTHQC.out.tables flagged_omegas = APPLYOMEGAQC.out.all_omegas_annotated + omega_vs_dndscv_qc = PLOTOMEGAVSDNDSCV.out.plots } From 3287ffe13508f7bba5396dbf98e2d9cfaa9cbd7a Mon Sep 17 00:00:00 2001 From: efiguerola Date: Thu, 2 Jul 2026 12:56:15 +0200 Subject: [PATCH 02/14] fixed and commented --- modules/local/plot/qc/omega_vs_dndscv/main.nf | 5 +++-- subworkflows/local/plotting_qc/main.nf | 4 +++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/modules/local/plot/qc/omega_vs_dndscv/main.nf b/modules/local/plot/qc/omega_vs_dndscv/main.nf index fccf4b78..8561df95 100644 --- a/modules/local/plot/qc/omega_vs_dndscv/main.nf +++ b/modules/local/plot/qc/omega_vs_dndscv/main.nf @@ -7,6 +7,7 @@ process PLOT_OMEGA_VS_DNDSCV { input: path (all_omegas) + path (dndscv_cv) path (groups_json) val (group_name) path (compiled_flagged) @@ -19,8 +20,8 @@ process PLOT_OMEGA_VS_DNDSCV { """ mkdir ${group_name}.plots omega_vs_dndscv_qc.py \\ - --input-omega-file ${all_mutdensities} \\ - --input-dndscv-file ${} \\ + --input-omega-file ${all_omegas} \\ + --input-dndscv-file ${dndscv_cv} \\ --output-dir ${group_name}.plots \\ --flagged-genes-omega ${compiled_flagged} diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index 1f3eb9a5..629eeb63 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -18,6 +18,8 @@ workflow PLOTTING_QC { panel groups_definition group_name + dndscv_cv //needs to be defined, where? + compiled_flagged //needs to be extracted from a process from the same module, how? // all_samples_depth // all_groups // full_panel_rich @@ -50,7 +52,7 @@ workflow PLOTTING_QC { all_omegas_globalloc ) - PLOTOMEGAVSDNDSCV(all_omegas, groups_definition, group_name) + PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, compiled_flagged) APPLYOMEGAQC(all_omegas, PLOTMUTDENSITYQC.out.compiled_flagged.collect()) From ab5375967df35841731eeb92bd081299075c484a Mon Sep 17 00:00:00 2001 From: efiguerola Date: Thu, 2 Jul 2026 15:51:54 +0200 Subject: [PATCH 03/14] fixed input calling variables --- subworkflows/local/plotting_qc/main.nf | 7 ++++--- workflows/deepcsa.nf | 5 ++++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index 629eeb63..955d9eda 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -18,8 +18,7 @@ workflow PLOTTING_QC { panel groups_definition group_name - dndscv_cv //needs to be defined, where? - compiled_flagged //needs to be extracted from a process from the same module, how? + dndscv_cv // all_samples_depth // all_groups // full_panel_rich @@ -52,9 +51,11 @@ workflow PLOTTING_QC { all_omegas_globalloc ) - PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, compiled_flagged) APPLYOMEGAQC(all_omegas, PLOTMUTDENSITYQC.out.compiled_flagged.collect()) + + PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, APPLYOMEGAQC.out.flagged_cases) + emit: mutdensity_plots = PLOTMUTDENSITYQC.out.plots diff --git a/workflows/deepcsa.nf b/workflows/deepcsa.nf index 87930718..6231b6d9 100644 --- a/workflows/deepcsa.nf +++ b/workflows/deepcsa.nf @@ -162,6 +162,7 @@ workflow DEEPCSA { all_mutdensities_file = channel.empty() all_adjusted_mutdensities_file = channel.value(file("${projectDir}/assets/placeholder_no_file.tsv", checkIfExists: true)) all_compiled_stabilities = channel.empty() + dndscv_table = channel.empty() // if the user wants to use custom gene groups, import the gene groups table // otherwise I am using the input csv as a dummy value channel @@ -451,6 +452,7 @@ workflow DEEPCSA { CREATEPANELS.out.exons_consensus_panel, params.fasta ) + dndscv_table = DNDS.out.all_dndscv_results } if (params.omega){ @@ -613,7 +615,8 @@ workflow DEEPCSA { // TABLE2GROUP.out.json_allgroups, CREATEPANELS.out.exons_consensus_panel, TABLE2GROUP.out.json_allgroups.first(), - group_keys_ch + group_keys_ch, + dndscv_table // CREATEPANELS.out.panel_annotated_rich, // seqinfo_df, // CREATEPANELS.out.domains_in_panel, From c78dbb72279f722b76232a857ca22ee47520365f Mon Sep 17 00:00:00 2001 From: efiguerola Date: Thu, 2 Jul 2026 16:11:43 +0200 Subject: [PATCH 04/14] added condition to only run PLOT_OMEGA_VS_DNDSCV whenever there DNDS step was enabled in the pipeline --- subworkflows/local/plotting_qc/main.nf | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index 955d9eda..7e7ba26c 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -54,7 +54,12 @@ workflow PLOTTING_QC { APPLYOMEGAQC(all_omegas, PLOTMUTDENSITYQC.out.compiled_flagged.collect()) - PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, APPLYOMEGAQC.out.flagged_cases) + // Only run omega vs dndscv QC plot when dndscv results are available + omega_vs_dndscv_plots = channel.empty() + if (params.dnds) { + PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, APPLYOMEGAQC.out.flagged_cases) + omega_vs_dndscv_plots = PLOTOMEGAVSDNDSCV.out.plots + } emit: @@ -62,6 +67,6 @@ workflow PLOTTING_QC { metrics_vs_depth_plots = PLOTMETRICSVSDEPTHQC.out.plots metrics_vs_depth_tables = PLOTMETRICSVSDEPTHQC.out.tables flagged_omegas = APPLYOMEGAQC.out.all_omegas_annotated - omega_vs_dndscv_qc = PLOTOMEGAVSDNDSCV.out.plots + omega_vs_dndscv_qc = omega_vs_dndscv_plots } From 9e2e42ec8a850893b430035898cef2c3b6d5dcaa Mon Sep 17 00:00:00 2001 From: efiguerola Date: Thu, 2 Jul 2026 16:37:07 +0200 Subject: [PATCH 05/14] fixed script execution and minor corrections --- bin/omega_vs_dndscv_qc.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) mode change 100644 => 100755 bin/omega_vs_dndscv_qc.py diff --git a/bin/omega_vs_dndscv_qc.py b/bin/omega_vs_dndscv_qc.py old mode 100644 new mode 100755 index 6f9d5f50..7f297ef5 --- a/bin/omega_vs_dndscv_qc.py +++ b/bin/omega_vs_dndscv_qc.py @@ -7,8 +7,6 @@ from matplotlib.backends.backend_pdf import PdfPages import numpy as np import seaborn as sns -import os -import time import json from scipy import stats from scipy.stats import pearsonr @@ -114,7 +112,7 @@ def apply_correlation_and_plotting(df, samples_group, output_dir): @click.option("--output-dir", required=True, type=click.Path(writable=True), help="Directory where output files will be written") -@click.option("--flagged-genes-omega", required=True, type=click.Path(writable=True), +@click.option("--flagged-genes-omega", required=True, type=click.Path(exists=True), help="Directory where flagged_omega genes were stored from /qc/omega_flagged/debug.syn_flagged_gene.tsv") From c08dbb7391b6ffde6688fce905ea0366f5c8bbb6 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Fri, 3 Jul 2026 15:44:50 +0200 Subject: [PATCH 06/14] added debuging on flagged_cases --- modules/local/plot/qc/annotate_omega/main.nf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/local/plot/qc/annotate_omega/main.nf b/modules/local/plot/qc/annotate_omega/main.nf index c1aa8932..6e5a706d 100644 --- a/modules/local/plot/qc/annotate_omega/main.nf +++ b/modules/local/plot/qc/annotate_omega/main.nf @@ -11,7 +11,7 @@ process ANNOTATE_OMEGA_QC { output: path("*flagged_annotated.tsv") , emit: all_omegas_annotated - path("*flagged.tsv") , optional: true, emit: flagged_cases + path("debug_*flagged*.tsv") , optional: true, emit: flagged_cases path("*.tsv") , optional: true, emit: files path("*.png") , optional: true, emit: plots path "versions.yml" , topic: versions From a7570bc0ad8610684fc0d954dd0c3ef9a15bfff6 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Fri, 3 Jul 2026 16:49:06 +0200 Subject: [PATCH 07/14] fixed debug typo --- modules/local/plot/qc/annotate_omega/main.nf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/local/plot/qc/annotate_omega/main.nf b/modules/local/plot/qc/annotate_omega/main.nf index 6e5a706d..5a9a11f1 100644 --- a/modules/local/plot/qc/annotate_omega/main.nf +++ b/modules/local/plot/qc/annotate_omega/main.nf @@ -11,7 +11,7 @@ process ANNOTATE_OMEGA_QC { output: path("*flagged_annotated.tsv") , emit: all_omegas_annotated - path("debug_*flagged*.tsv") , optional: true, emit: flagged_cases + path("debug.*flagged*.tsv") , optional: true, emit: flagged_cases path("*.tsv") , optional: true, emit: files path("*.png") , optional: true, emit: plots path "versions.yml" , topic: versions From 0c39d3f1cf4c39568c4b039c38da51e0fe3184b3 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Mon, 6 Jul 2026 15:21:05 +0200 Subject: [PATCH 08/14] fix error on taking flagged genes so it takes only the synonymous flagged file --- subworkflows/local/plotting_qc/main.nf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index 7e7ba26c..7f5de5fd 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -57,7 +57,7 @@ workflow PLOTTING_QC { // Only run omega vs dndscv QC plot when dndscv results are available omega_vs_dndscv_plots = channel.empty() if (params.dnds) { - PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, APPLYOMEGAQC.out.flagged_cases) + PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, APPLYOMEGAQC.out.flagged_cases.flatten().filter { it.name == 'debug.syn_flagged_gene.tsv' }) omega_vs_dndscv_plots = PLOTOMEGAVSDNDSCV.out.plots } From 8f69ce5d1615f6f191c93013d4d40e21e2c25d87 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Mon, 6 Jul 2026 15:27:55 +0200 Subject: [PATCH 09/14] fix error on missing import os --- bin/omega_vs_dndscv_qc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/bin/omega_vs_dndscv_qc.py b/bin/omega_vs_dndscv_qc.py index 7f297ef5..24814b3b 100755 --- a/bin/omega_vs_dndscv_qc.py +++ b/bin/omega_vs_dndscv_qc.py @@ -2,6 +2,7 @@ # Needed basic packages import click +import os import pandas as pd import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages From e7763d362b0390fa5674c7d336bf837601da85e4 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Mon, 6 Jul 2026 15:52:23 +0200 Subject: [PATCH 10/14] fixed output storage, added prints in the script and also stored the omega-dndscv filtered table --- bin/omega_vs_dndscv_qc.py | 15 +++++++++++---- conf/results_outputs.config | 8 ++++++++ modules/local/plot/qc/omega_vs_dndscv/main.nf | 1 + 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/bin/omega_vs_dndscv_qc.py b/bin/omega_vs_dndscv_qc.py index 24814b3b..a473687c 100755 --- a/bin/omega_vs_dndscv_qc.py +++ b/bin/omega_vs_dndscv_qc.py @@ -48,7 +48,7 @@ def apply_correlation_and_plotting(df, samples_group, output_dir): for group in sorted(samples_group): subset = df[df['sample'] == group].copy() subset[['dnds_omega', 'dnds_dndscv']] = subset[['dnds_omega', 'dnds_dndscv']].fillna(0) - print(subset.head()) + print(f'Subset shape for {group}: {subset.shape[0]}') if subset.empty: print(f'No data available for sample: {group}') @@ -132,21 +132,28 @@ def main(input_omega_file, input_dndscv_file, output_dir, flagged_genes_omega): omega_filtered_df = omega_df[~(omega_df['gene'].str.contains("--")) & # discards exons if they are included in the analysis (omega_df['impact'].isin(['missense', 'truncating']))].copy() print('Filtered omega table:') - print(omega_filtered_df.head()) + print(omega_filtered_df.shape) # Process dNdScv table dndscv_cv_df = process_dndscv_table(dndscv_df) print('Processed dNdScv table:') - print(dndscv_cv_df.head()) + print(dndscv_cv_df.shape) # Merge omega and dNdScv tables dndscv_n_omegas_df = omega_filtered_df.merge(dndscv_cv_df, on=['sample', 'gene', 'impact'], how='outer', suffixes=('_omega', '_dndscv')) print('Merged omega and dNdScv table:') - print(dndscv_n_omegas_df.head()) + print(dndscv_n_omegas_df.shape) # Filter flagged omega genes within groups filtered_dndscv_n_omegas_df = filter_flagged_genes_per_group(dndscv_n_omegas_df, flagged_genes_df) + # Export filtered table + print('Filtered table from flagged genes, shape:') + print(filtered_dndscv_n_omegas_df.shape) + + filtered_output_path = os.path.join(output_dir, "filtered_omega_dndscv_table.tsv") + filtered_dndscv_n_omegas_df.to_csv(filtered_output_path, sep='\t', index=False) + # Apply function to plot omega values vs dndscv per group and compute correlation apply_correlation_and_plotting(filtered_dndscv_n_omegas_df, sample_groups, output_dir) diff --git a/conf/results_outputs.config b/conf/results_outputs.config index a0a77334..87fab20b 100644 --- a/conf/results_outputs.config +++ b/conf/results_outputs.config @@ -118,6 +118,14 @@ process { ] } + withName: PLOTOMEGAVSDNDSCV { + publishDir = [ + path: { "${params.outdir}/qc/omega_vs_dndscv" }, + mode: params.publish_dir_mode, + pattern: '**{tsv,pdf,png}', + ] + } + withName: APPLYOMEGAQC { publishDir = [ [ diff --git a/modules/local/plot/qc/omega_vs_dndscv/main.nf b/modules/local/plot/qc/omega_vs_dndscv/main.nf index 8561df95..7e36a11a 100644 --- a/modules/local/plot/qc/omega_vs_dndscv/main.nf +++ b/modules/local/plot/qc/omega_vs_dndscv/main.nf @@ -14,6 +14,7 @@ process PLOT_OMEGA_VS_DNDSCV { output: path("**.pdf") , optional: true , emit: plots + path("**.tsv") , optional: true , emit: tables path "versions.yml" , topic: versions script: From 8a1885f6d2c8758405197aa77e548d7eb9ffadab Mon Sep 17 00:00:00 2001 From: efiguerola Date: Tue, 7 Jul 2026 11:00:24 +0200 Subject: [PATCH 11/14] fixed output name --- modules/local/plot/qc/annotate_omega/main.nf | 11 ++++++----- modules/local/plot/qc/omega_vs_dndscv/main.nf | 7 ++----- subworkflows/local/plotting_qc/main.nf | 2 +- 3 files changed, 9 insertions(+), 11 deletions(-) diff --git a/modules/local/plot/qc/annotate_omega/main.nf b/modules/local/plot/qc/annotate_omega/main.nf index 5a9a11f1..bc91915f 100644 --- a/modules/local/plot/qc/annotate_omega/main.nf +++ b/modules/local/plot/qc/annotate_omega/main.nf @@ -10,11 +10,12 @@ process ANNOTATE_OMEGA_QC { path (compiled_flagged_cases) output: - path("*flagged_annotated.tsv") , emit: all_omegas_annotated - path("debug.*flagged*.tsv") , optional: true, emit: flagged_cases - path("*.tsv") , optional: true, emit: files - path("*.png") , optional: true, emit: plots - path "versions.yml" , topic: versions + path("*flagged_annotated.tsv") , emit: all_omegas_annotated + path("debug.*flagged*.tsv") , optional: true, emit: flagged_cases + path("debug.syn_flagged_gene.tsv") , optional: true, emit: flagged_synonymous_cases + path("*.tsv") , optional: true, emit: files + path("*.png") , optional: true, emit: plots + path "versions.yml" , topic: versions script: """ diff --git a/modules/local/plot/qc/omega_vs_dndscv/main.nf b/modules/local/plot/qc/omega_vs_dndscv/main.nf index 7e36a11a..ee5fd98c 100644 --- a/modules/local/plot/qc/omega_vs_dndscv/main.nf +++ b/modules/local/plot/qc/omega_vs_dndscv/main.nf @@ -1,6 +1,6 @@ process PLOT_OMEGA_VS_DNDSCV { - tag "${group_name}" + tag "all_samples" label 'process_low' label 'deepcsa_core' @@ -8,8 +8,6 @@ process PLOT_OMEGA_VS_DNDSCV { input: path (all_omegas) path (dndscv_cv) - path (groups_json) - val (group_name) path (compiled_flagged) output: @@ -19,11 +17,10 @@ process PLOT_OMEGA_VS_DNDSCV { script: """ - mkdir ${group_name}.plots omega_vs_dndscv_qc.py \\ --input-omega-file ${all_omegas} \\ --input-dndscv-file ${dndscv_cv} \\ - --output-dir ${group_name}.plots \\ + --output-dir . \\ --flagged-genes-omega ${compiled_flagged} diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index 7f5de5fd..fe532552 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -57,7 +57,7 @@ workflow PLOTTING_QC { // Only run omega vs dndscv QC plot when dndscv results are available omega_vs_dndscv_plots = channel.empty() if (params.dnds) { - PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, groups_definition, group_name, APPLYOMEGAQC.out.flagged_cases.flatten().filter { it.name == 'debug.syn_flagged_gene.tsv' }) + PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, APPLYOMEGAQC.out.flagged_synonymous_cases) omega_vs_dndscv_plots = PLOTOMEGAVSDNDSCV.out.plots } From dfbc7f02d928a881077a9f9e7b7a10bd3aecc4df Mon Sep 17 00:00:00 2001 From: efiguerola Date: Tue, 7 Jul 2026 15:01:48 +0200 Subject: [PATCH 12/14] added omega omegaglobal comparison code --- bin/omega_vs_omegaglobal_qc.py | 145 ++++++++++++++++++ bin/plot_explore_variability.py | 1 + conf/results_outputs.config | 12 +- modules/local/plot/qc/omega_vs_dndscv/main.nf | 3 +- .../plot/qc/omega_vs_omegaglobal/main.nf | 42 +++++ subworkflows/local/plotting_qc/main.nf | 14 +- 6 files changed, 208 insertions(+), 9 deletions(-) create mode 100755 bin/omega_vs_omegaglobal_qc.py create mode 100644 modules/local/plot/qc/omega_vs_omegaglobal/main.nf diff --git a/bin/omega_vs_omegaglobal_qc.py b/bin/omega_vs_omegaglobal_qc.py new file mode 100755 index 00000000..6bfb71af --- /dev/null +++ b/bin/omega_vs_omegaglobal_qc.py @@ -0,0 +1,145 @@ +#!/usr/bin/env python + +# Needed basic packages +import click +import os +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np +import seaborn as sns +import json +from scipy import stats +from scipy.stats import pearsonr + + +# Functions +def filter_flagged_genes_per_group(df, flagged_cases_df): + # Optimized: Vectorized filtering using a multi-index instead of .apply row-by-row + flagged_index = pd.MultiIndex.from_frame(flagged_cases_df[['cohort', 'gene']].rename(columns={'cohort': 'sample'})) + df_index = pd.MultiIndex.from_frame(df[['sample', 'gene']]) + + filtered_df = df[~df_index.isin(flagged_index)].copy() + return filtered_df + + +def apply_correlation_and_plotting(df, samples_group, output_dir): + + os.makedirs(output_dir, exist_ok=True) # Ensure output dir exists + + for group in sorted(samples_group): + subset = df[df['sample'] == group].copy() + subset[['dnds_omega', 'dnds_omegaglobal']] = subset[['dnds_omega', 'dnds_omegaglobal']].fillna(0) + print(f'Subset shape for {group}: {subset.shape[0]}') + + if subset.empty: + print(f'No data available for sample: {group}') + continue + + fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharex=False, sharey=False) + + impacts = ['missense', 'truncating'] + + for i, impact in enumerate(impacts): + ax = axes[i] + impact_subset = subset[subset['impact'] == impact] + impact_subset_reg = impact_subset[(impact_subset['dnds_omega'] > 0) & (impact_subset['dnds_omegaglobal'] > 0)] + + if impact_subset.empty: + ax.set_title(f'No {impact} data for {group}') + continue + + # Scatter plot + sns.scatterplot(data=impact_subset, x='dnds_omega', y='dnds_omegaglobal', hue='gene', ax=ax) + + # Regression line + sns.regplot(data=impact_subset_reg, x='dnds_omega', y='dnds_omegaglobal', scatter=False, color='red', label='Regression Line', ax=ax) + + # X=Y line + all_vals = pd.concat([impact_subset['dnds_omega'], impact_subset['dnds_omegaglobal']]) + min_val, max_val = all_vals.min(), all_vals.max() + ax.plot([min_val, max_val], [min_val, max_val], color='black', linestyle='--', label='x=y') + + if impact_subset_reg.shape[0] > 2: + # Compute correlation + corr, pval = pearsonr(impact_subset_reg['dnds_omega'], impact_subset_reg['dnds_omegaglobal']) + + ax.set_title(f'{impact.capitalize()} mutations\nPearson R: {corr:.3f} (p={pval:.3e})') + else: + ax.set_title(f'{impact.capitalize()} mutations\nNot enough data for correlation') + + ax.set_xlabel('dN/dS (Omega)') + ax.set_ylabel('dN/dS (Omegaglobal)') + ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + ax.legend().remove() + + fig.suptitle(f'Comparison for {group}') + fig.tight_layout() + + output_path = os.path.join(output_dir, f"{group}_omega_vs_omegaglobal_qc_summary_plot.pdf") + with PdfPages(output_path) as pdf: + pdf.savefig(fig) + + plt.close(fig) + + return + + +@click.command() +@click.option("--input-omega-file", required=True, type=click.Path(exists=True), + help="Directory containing selection/omega/all_omegas.tsv") + +@click.option("--input-omegaglobal-file", required=True, type=click.Path(exists=True), + help="Directory containing selection/omegaglobal/all_omegaglobal.tsv") + +@click.option("--output-dir", required=True, type=click.Path(writable=True), + help="Directory where output files will be written") + +@click.option("--flagged-genes-omega", required=True, type=click.Path(exists=True), + help="Directory where flagged_omega genes were stored from /qc/omega_flagged/debug.syn_flagged_gene.tsv") + + +# Main function +def main(input_omega_file, input_omegaglobal_file, output_dir, flagged_genes_omega): + + # Read tables + omega_df = pd.read_table(input_omega_file, sep='\t') + omegaglobal_df = pd.read_table(input_omegaglobal_file, sep='\t') + flagged_genes_df = pd.read_table(flagged_genes_omega, sep='\t') + + # Extract sample groups + sample_groups = flagged_genes_df['cohort'].unique().tolist() + + # Apply basic filter for omega and omegaglobal table + omega_filtered_df = omega_df[~(omega_df['gene'].str.contains("--")) & # discards exons if they are included in the analysis + (omega_df['impact'].isin(['missense', 'truncating']))].copy() + print('Filtered omega table:') + print(omega_filtered_df.shape) + + # Apply basic filter for omega and omegaglobal table + omegaglobal_filtered_df = omegaglobal_df[~(omegaglobal_df['gene'].str.contains("--")) & # discards exons if they are included in the analysis + (omegaglobal_df['impact'].isin(['missense', 'truncating']))].copy() + print('Filtered omegaglobal table:') + print(omegaglobal_filtered_df.shape) + + + # Merge omega and omegaglobal tables + omega_vs_omegaglobal_df = omega_filtered_df.merge(omegaglobal_filtered_df, on=['sample', 'gene', 'impact'], how='outer', suffixes=('_omega', '_omegaglobal')) + print('Merged omega and omegaglobal table:') + print(omega_vs_omegaglobal_df.shape) + + # Filter flagged omega genes within groups + filtered_omega_vs_omegaglobal_df = filter_flagged_genes_per_group(omega_vs_omegaglobal_df, flagged_genes_df) + + # Export filtered table + print('Filtered table from flagged genes, shape:') + print(filtered_omega_vs_omegaglobal_df.shape) + + filtered_output_path = os.path.join(output_dir, "filtered_omega_omegaglobal_table.tsv") + filtered_omega_vs_omegaglobal_df.to_csv(filtered_output_path, sep='\t', index=False) + + # Apply function to plot omega values vs omegaglobal per group and compute correlation + apply_correlation_and_plotting(filtered_omega_vs_omegaglobal_df, sample_groups, output_dir) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/bin/plot_explore_variability.py b/bin/plot_explore_variability.py index 0ba22504..466d2ae2 100755 --- a/bin/plot_explore_variability.py +++ b/bin/plot_explore_variability.py @@ -252,6 +252,7 @@ def main(outdir, panel_regions, samples_json, all_groups_json, mutdensities, adj except Exception as e: print("Error in the process", e) + # pca plot for mutdensity and mutdensity adjusted if __name__ == '__main__': main() diff --git a/conf/results_outputs.config b/conf/results_outputs.config index 87fab20b..82356105 100644 --- a/conf/results_outputs.config +++ b/conf/results_outputs.config @@ -120,7 +120,15 @@ process { withName: PLOTOMEGAVSDNDSCV { publishDir = [ - path: { "${params.outdir}/qc/omega_vs_dndscv" }, + path: { "${params.outdir}/qc/omegaqc/omega_vs_dndscv" }, + mode: params.publish_dir_mode, + pattern: '**{tsv,pdf,png}', + ] + } + + withName: PLOTOMEGAVSOMEGAGLOBAL { + publishDir = [ + path: { "${params.outdir}/qc/omegaqc/omega_vs_omegaglobal" }, mode: params.publish_dir_mode, pattern: '**{tsv,pdf,png}', ] @@ -134,7 +142,7 @@ process { pattern: 'omega.flagged_annotated.tsv' ], [ - path: { "${params.outdir}/qc/omega_flagged" }, + path: { "${params.outdir}/qc/omegaqc/omega_flagged" }, mode: params.publish_dir_mode, pattern: '*.{tsv,pdf,png}', saveAs: { filename -> diff --git a/modules/local/plot/qc/omega_vs_dndscv/main.nf b/modules/local/plot/qc/omega_vs_dndscv/main.nf index ee5fd98c..d11ac355 100644 --- a/modules/local/plot/qc/omega_vs_dndscv/main.nf +++ b/modules/local/plot/qc/omega_vs_dndscv/main.nf @@ -2,7 +2,6 @@ process PLOT_OMEGA_VS_DNDSCV { tag "all_samples" label 'process_low' - label 'deepcsa_core' input: @@ -12,7 +11,7 @@ process PLOT_OMEGA_VS_DNDSCV { output: path("**.pdf") , optional: true , emit: plots - path("**.tsv") , optional: true , emit: tables + path("**.tsv") , optional: true , emit: tables path "versions.yml" , topic: versions script: diff --git a/modules/local/plot/qc/omega_vs_omegaglobal/main.nf b/modules/local/plot/qc/omega_vs_omegaglobal/main.nf new file mode 100644 index 00000000..7012b8d4 --- /dev/null +++ b/modules/local/plot/qc/omega_vs_omegaglobal/main.nf @@ -0,0 +1,42 @@ +process PLOTOMEGAVSOMEGAGLOBAL { + + tag "all_samples" + label 'process_low' + label 'deepcsa_core' + + input: + path (all_omegas) + path (all_omegas_globalloc) + path (compiled_flagged) + + output: + path("**.pdf") , optional: true , emit: plots + path("**.tsv") , optional: true , emit: tables + path "versions.yml" , topic: versions + + script: + """ + omega_vs_omegaglobal_qc.py \\ + --input-omega-file ${all_omegas} \\ + --input-omegaglobal-file ${all_omegas_globalloc} \\ + --output-dir . \\ + --flagged-genes-omega ${compiled_flagged} + + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "all_samples" + """ + touch ${prefix}.pdf + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ +} diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index fe532552..b96cc3ae 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -1,9 +1,11 @@ -include { PLOT_MUTDENSITY_QC as PLOTMUTDENSITYQC } from '../../../modules/local/plot/qc/mutation_densities/main' -include { PLOT_METRICS_VS_DEPTH_QC as PLOTMETRICSVSDEPTHQC } from '../../../modules/local/plot/qc/metrics_vs_depth/main' -include { ANNOTATE_OMEGA_QC as APPLYOMEGAQC } from '../../../modules/local/plot/qc/annotate_omega/main' -include { PLOT_MUTATION_SPECIFIC as PLOTMUTATIONSPECIFIC } from '../../../modules/local/plot/qc/mutation_specific/main' -include { PLOT_OMEGA_VS_DNDSCV as PLOTOMEGAVSDNDSCV } from '../../../modules/local/plot/qc/omega_vs_dndscv/main' +include { PLOT_MUTDENSITY_QC as PLOTMUTDENSITYQC } from '../../../modules/local/plot/qc/mutation_densities/main' +include { PLOT_METRICS_VS_DEPTH_QC as PLOTMETRICSVSDEPTHQC } from '../../../modules/local/plot/qc/metrics_vs_depth/main' +include { ANNOTATE_OMEGA_QC as APPLYOMEGAQC } from '../../../modules/local/plot/qc/annotate_omega/main' +include { PLOT_MUTATION_SPECIFIC as PLOTMUTATIONSPECIFIC } from '../../../modules/local/plot/qc/mutation_specific/main' +include { PLOT_OMEGA_VS_DNDSCV as PLOTOMEGAVSDNDSCV } from '../../../modules/local/plot/qc/omega_vs_dndscv/main' +include { PLOT_OMEGA_VS_OMEGAGLOBAL as PLOTOMEGAVSOMEGAGLOBAL} from '../../../modules/local/plot/qc/omega_vs_omegaglobal/main' + workflow PLOTTING_QC { @@ -61,6 +63,7 @@ workflow PLOTTING_QC { omega_vs_dndscv_plots = PLOTOMEGAVSDNDSCV.out.plots } + PLOTOMEGAVSOMEGAGLOBAL(all_omegas, all_omegas_globalloc, APPLYOMEGAQC.out.flagged_synonymous_cases) emit: mutdensity_plots = PLOTMUTDENSITYQC.out.plots @@ -68,5 +71,6 @@ workflow PLOTTING_QC { metrics_vs_depth_tables = PLOTMETRICSVSDEPTHQC.out.tables flagged_omegas = APPLYOMEGAQC.out.all_omegas_annotated omega_vs_dndscv_qc = omega_vs_dndscv_plots + omega_vs_omegaglobal = PLOTOMEGAVSOMEGAGLOBAL.out.plots } From dd9e556a644ac3eef221cdfe837ab13ca7c21765 Mon Sep 17 00:00:00 2001 From: efiguerola Date: Tue, 7 Jul 2026 15:13:06 +0200 Subject: [PATCH 13/14] fix module name --- modules/local/plot/qc/omega_vs_omegaglobal/main.nf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/local/plot/qc/omega_vs_omegaglobal/main.nf b/modules/local/plot/qc/omega_vs_omegaglobal/main.nf index 7012b8d4..26579a4e 100644 --- a/modules/local/plot/qc/omega_vs_omegaglobal/main.nf +++ b/modules/local/plot/qc/omega_vs_omegaglobal/main.nf @@ -1,4 +1,4 @@ -process PLOTOMEGAVSOMEGAGLOBAL { +process PLOT_OMEGA_VS_OMEGAGLOBAL { tag "all_samples" label 'process_low' From f3fc0387fe0f6715bda74e3977fc3f85d3845def Mon Sep 17 00:00:00 2001 From: Bet Figuerola <154223352+efigb@users.noreply.github.com> Date: Tue, 14 Jul 2026 12:19:40 +0200 Subject: [PATCH 14/14] Potential fix for pull request finding Make this plot conditional (similar to the dNdScv plot) and emit an empty channel when not applicable. Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- subworkflows/local/plotting_qc/main.nf | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/subworkflows/local/plotting_qc/main.nf b/subworkflows/local/plotting_qc/main.nf index b96cc3ae..9466023b 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -62,8 +62,13 @@ workflow PLOTTING_QC { PLOTOMEGAVSDNDSCV(all_omegas, dndscv_cv, APPLYOMEGAQC.out.flagged_synonymous_cases) omega_vs_dndscv_plots = PLOTOMEGAVSDNDSCV.out.plots } - - PLOTOMEGAVSOMEGAGLOBAL(all_omegas, all_omegas_globalloc, APPLYOMEGAQC.out.flagged_synonymous_cases) + + // Only run omega vs omegaglobal QC plot when omegaglobal results are available + omega_vs_omegaglobal_plots = channel.empty() + if (params.omega_globalloc) { + PLOTOMEGAVSOMEGAGLOBAL(all_omegas, all_omegas_globalloc, APPLYOMEGAQC.out.flagged_synonymous_cases) + omega_vs_omegaglobal_plots = PLOTOMEGAVSOMEGAGLOBAL.out.plots + } emit: mutdensity_plots = PLOTMUTDENSITYQC.out.plots @@ -71,6 +76,6 @@ workflow PLOTTING_QC { metrics_vs_depth_tables = PLOTMETRICSVSDEPTHQC.out.tables flagged_omegas = APPLYOMEGAQC.out.all_omegas_annotated omega_vs_dndscv_qc = omega_vs_dndscv_plots - omega_vs_omegaglobal = PLOTOMEGAVSOMEGAGLOBAL.out.plots + omega_vs_omegaglobal = omega_vs_omegaglobal_plots }