diff --git a/bin/omega_vs_dndscv_qc.py b/bin/omega_vs_dndscv_qc.py new file mode 100755 index 00000000..a473687c --- /dev/null +++ b/bin/omega_vs_dndscv_qc.py @@ -0,0 +1,161 @@ +#!/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 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(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_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(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_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.shape) + + # Process dNdScv table + dndscv_cv_df = process_dndscv_table(dndscv_df) + print('Processed dNdScv table:') + 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.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) + +if __name__ == '__main__': + main() \ No newline at end of file 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 a0a77334..82356105 100644 --- a/conf/results_outputs.config +++ b/conf/results_outputs.config @@ -118,6 +118,22 @@ process { ] } + withName: PLOTOMEGAVSDNDSCV { + publishDir = [ + 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}', + ] + } + withName: APPLYOMEGAQC { publishDir = [ [ @@ -126,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/annotate_omega/main.nf b/modules/local/plot/qc/annotate_omega/main.nf index c1aa8932..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("*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 new file mode 100644 index 00000000..d11ac355 --- /dev/null +++ b/modules/local/plot/qc/omega_vs_dndscv/main.nf @@ -0,0 +1,42 @@ +process PLOT_OMEGA_VS_DNDSCV { + + tag "all_samples" + label 'process_low' + label 'deepcsa_core' + + input: + path (all_omegas) + path (dndscv_cv) + path (compiled_flagged) + + output: + path("**.pdf") , optional: true , emit: plots + path("**.tsv") , optional: true , emit: tables + path "versions.yml" , topic: versions + + script: + """ + omega_vs_dndscv_qc.py \\ + --input-omega-file ${all_omegas} \\ + --input-dndscv-file ${dndscv_cv} \\ + --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/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..26579a4e --- /dev/null +++ b/modules/local/plot/qc/omega_vs_omegaglobal/main.nf @@ -0,0 +1,42 @@ +process PLOT_OMEGA_VS_OMEGAGLOBAL { + + 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 ab183920..9466023b 100644 --- a/subworkflows/local/plotting_qc/main.nf +++ b/subworkflows/local/plotting_qc/main.nf @@ -1,8 +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_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 { @@ -17,6 +20,7 @@ workflow PLOTTING_QC { panel groups_definition group_name + dndscv_cv // all_samples_depth // all_groups // full_panel_rich @@ -49,12 +53,29 @@ workflow PLOTTING_QC { all_omegas_globalloc ) + APPLYOMEGAQC(all_omegas, PLOTMUTDENSITYQC.out.compiled_flagged.collect()) + + // 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, APPLYOMEGAQC.out.flagged_synonymous_cases) + omega_vs_dndscv_plots = PLOTOMEGAVSDNDSCV.out.plots + } + + // 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 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 = omega_vs_dndscv_plots + omega_vs_omegaglobal = omega_vs_omegaglobal_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,