Another fix to omega globalloc

assets/input_double_bam.csv

@@ -1,3 +1,3 @@
 sample,vcf,bam,pileup_bam,pileup_ind
-K_5_1_A_1,/workspace/datasets/transfer/ferriol_deepcsa/K_5_1_A_1.high.filtered.vcf,/workspace/datasets/transfer/ferriol_deepcsa/K_5_1_A_1.sorted.bam,/workspace/datasets/transfer/ferriol_deepcsa/allmolecules/K_5_1_A_1.sorted.bam,/workspace/datasets/transfer/ferriol_deepcsa/allmolecules/K_5_1_A_1.sorted.bam.csi
-K_6_1_A_1,/workspace/datasets/transfer/ferriol_deepcsa/K_6_1_A_1.high.filtered.vcf,/workspace/datasets/transfer/ferriol_deepcsa/K_6_1_A_1.sorted.bam,/workspace/datasets/transfer/ferriol_deepcsa/allmolecules/K_6_1_A_1.sorted.bam,/workspace/datasets/transfer/ferriol_deepcsa/allmolecules/K_6_1_A_1.sorted.bam.csi
+K_5_1_A_1,/data/bbg/datasets/transfer/ferriol_deepcsa/K_5_1_A_1.high.filtered.vcf,/data/bbg/datasets/transfer/ferriol_deepcsa/K_5_1_A_1.sorted.bam,/data/bbg/datasets/transfer/ferriol_deepcsa/allmolecules/K_5_1_A_1.sorted.bam,/data/bbg/datasets/transfer/ferriol_deepcsa/allmolecules/K_5_1_A_1.sorted.bam.csi
+K_6_1_A_1,/data/bbg/datasets/transfer/ferriol_deepcsa/K_6_1_A_1.high.filtered.vcf,/data/bbg/datasets/transfer/ferriol_deepcsa/K_6_1_A_1.sorted.bam,/data/bbg/datasets/transfer/ferriol_deepcsa/allmolecules/K_6_1_A_1.sorted.bam,/data/bbg/datasets/transfer/ferriol_deepcsa/allmolecules/K_6_1_A_1.sorted.bam.csi

bin/omega_select_mutrate.py

@@ -0,0 +1,49 @@
+#!/usr/local/bin/python
+
+
+import click
+import pandas as pd
+from read_utils import custom_na_values
+
+
+def select_syn_mutrate(mutrate_file, output_file, mode):
+    """
+    INFO
+    """
+
+    mutrate_df = pd.read_csv(mutrate_file, sep = "\t", header = 0, na_values = custom_na_values)
+
+    synonymous_mutrates_all_samples = mutrate_df[(mutrate_df["MUTTYPES"] == "SNV") &
+                                                    (mutrate_df["GENE"] != "ALL_GENES")].reset_index(drop = True)
+
+    if mode == 'mutations':
+        synonymous_mutrates_genes = synonymous_mutrates_all_samples[['GENE', 'MUTRATE_MB_ADJUSTED']]
+    elif mode == 'mutated_reads':
+        synonymous_mutrates_genes = synonymous_mutrates_all_samples[['GENE', 'MUTREADSRATE_MB_ADJUSTED']]
+
+    ## FIXME not sure if this is really needed since when called through main()
+    # the input would have already been forced to be either of the two options
+    # it might still be useful in case this was not called from main()
+    else:
+        print('unknown mode, please enter either mutations or mutated_reads')
+        exit(1) ## FIXME not sure if this is the right code to exit with
+
+    synonymous_mutrates_genes.columns = ["GENE", "MUTRATE"]
+    synonymous_mutrates_genes.to_csv(f"{output_file}",
+                                        header=True,
+                                        index=False,
+                                        sep="\t")
+
+
+@click.command()
+@click.option('--mutrates', type=click.Path(exists=True), help='Input mutation rate file')
+@click.option('--output', type=click.Path(), help='Output file')
+@click.option('--mode', type=click.Choice(['mutations', 'mutated_reads']), default='mutations')
+
+def main(mutrates, output, mode):
+    click.echo("Selecting the gene synonymous mutation rates...")
+    select_syn_mutrate(mutrates, output, mode)
+
+if __name__ == '__main__':
+    main()
+

conf/modules.config

@@ -295,7 +295,7 @@ process {
                                 ['"VAF_AM" : "le ', "${params.germline_threshold}", '"'].join('').trim(),
                                 ['"vd_VAF" : "le ', "${params.germline_threshold}", '"'].join('').trim(),
                                 ['"DEPTH" : "ge ', "${params.mutation_depth_threshold}", '"'].join('').trim(),
-                                '"FILTER" : ["notcontains NM20", "notcontains n_rich", "notcontains cohort_n_rich", "notcontains cohort_n_rich_uni",  "notcontains no_pileup_support", "notcontains other_sample_SNP", "notcontains low_mappability" ]',
+                                '"FILTER" : ["notcontains NM20", "notcontains n_rich", "notcontains cohort_n_rich", "notcontains cohort_n_rich_uni",  "notcontains no_pileup_support", "notcontains low_mappability" ]',
                                 '"VAF_distorted_expanded" : false',
                                 ].join(',\t').trim()
                             }
@@ -315,7 +315,7 @@ process {
             ext.filters     = { [['"VAF" : "le ', "${params.germline_threshold}", '"'].join('').trim(),
                                 ['"vd_VAF" : "le ', "${params.germline_threshold}", '"'].join('').trim(),
                                 ['"DEPTH" : "ge ', "${params.mutation_depth_threshold}", '"'].join('').trim(),
-                                '"FILTER" : ["notcontains NM20", "notcontains n_rich", "notcontains cohort_n_rich", "notcontains cohort_n_rich_uni",  "notcontains no_pileup_support", "notcontains other_sample_SNP", "notcontains low_mappability" ]',
+                                '"FILTER" : ["notcontains NM20", "notcontains n_rich", "notcontains cohort_n_rich", "notcontains cohort_n_rich_uni",  "notcontains no_pileup_support",  "notcontains low_mappability" ]',
                                 '"VAF_distorted_expanded" : false',
                                 ].join(',\t').trim()
                             }
@@ -489,6 +489,13 @@ process {
             ]
         }
 
+        withName: 'SYNMUTRATE' {
+            ext.mode     = 'mutations'
+        }
+        withName: 'SYNMUTREADSRATE' {
+            ext.mode     = 'mutated_reads'
+        }
+
     }
 
 

conf/test.config

@@ -27,7 +27,7 @@ params {
 
 
     omega_hotspots              = false
-    omega_hotspots_bedfile     = "/workspace/datasets/transfer/ferriol_deepcsa/kidney_panel.hotspots.bed4.bed"
+    omega_hotspots_bedfile     = "/data/bbg/datasets/transfer/ferriol_deepcsa/kidney_panel.hotspots.bed4.bed"
     hotspot_expansion           = 30
 
     oncodrivefml               = false

modules/local/bbgtools/omega/estimator/main.nf

@@ -5,7 +5,7 @@ process OMEGA_ESTIMATOR {
     label 'process_high_memory'
 
 
-    container 'docker.io/ferriolcalvet/omega:latest'
+    container 'docker.io/ferriolcalvet/omega:20241217'
 
     input:
     tuple val(meta) , path(mutabilities_table), path(mutations_table), path(depths)

modules/local/bbgtools/omega/preprocess/main.nf

@@ -5,12 +5,13 @@ process OMEGA_PREPROCESS {
     label 'process_high_memory'
 
 
-    container 'docker.io/ferriolcalvet/omega:latest'
+    container 'docker.io/ferriolcalvet/omega:20241217'
 
     input:
     tuple val(meta) , path(mutations), path(depths), path(mutation_profile)
     tuple val(meta2), path (annotated_panel)
     tuple val(meta3), path (syn_muts_global)
+    tuple val(meta4), path (mut_profile_global, stageAs: 'global_mutprofile.tsv')
 
 
     output:
@@ -26,7 +27,7 @@ process OMEGA_PREPROCESS {
     def prefix = task.ext.prefix ?: "${meta.id}"
     // TODO revise this fix
     def sample_name = prefix.tokenize('.')[0]
-    def global_loc = task.ext.global_loc ? "--absent-synonymous infer_global_custom  --relative-synonymous-muts-file ${syn_muts_global}" : "--absent-synonymous ignore"
+    def global_loc = task.ext.global_loc ? "--absent-synonymous infer_global_custom  --mutational-profile-global-file global_mutprofile.tsv --synonymous-mutrates-file ${syn_muts_global}" : "--absent-synonymous ignore"
     prefix = task.ext.global_loc ? "${prefix}.gLoc" : "${prefix}"
     """
     omega preprocessing --preprocessing-mode compute_mutabilities \\
@@ -36,7 +37,7 @@ process OMEGA_PREPROCESS {
                         --table-observed-muts mutations_per_sample_gene_impact_context.${prefix}.tsv \\
                         --mutabilities-table mutability_per_sample_gene_context.${prefix}.tsv \\
                         --synonymous-muts-table syn_muts.${prefix}.tsv \\
-                        --mutational-profile ${mutation_profile} \\
+                        --mutational-profile-file ${mutation_profile} \\
                         --single-sample ${sample_name} \\
                         ${global_loc}
     # $args -c $task.cpus
@@ -59,3 +60,15 @@ process OMEGA_PREPROCESS {
     """
 }
 
+// omega preprocessing --preprocessing-mode compute_mutabilities
+//                      --depths-file all_samples.subset_depths.tsv.gz
+//                      --mutations-file all_samples.mutations.tsv
+//                      --input-vep-postprocessed-file consensus.exons_splice_sites.tsv
+//                      --table-observed-muts mutations_per_sample_gene_impact_context.all_samples2.global_loc.gLoc.tsv
+//                      --mutabilities-table mutability_per_sample_gene_context.all_samples2.global_loc.gLoc.tsv
+//                      --synonymous-muts-table syn_muts.all_samples2.global_loc.gLoc.tsv
+//                      --mutational-profile-file all_samples.all.profile.tsv
+//                      --mutational-profile-global-file P19_0033_BTR_01.all.profile.tsv
+//                      --single-sample all_samples
+//                      --absent-synonymous infer_global_custom
+//                      --synonymous-mutrates-file mutrates_per_gene.tsv

modules/local/relative_mutrate/main.nf → modules/local/select_mutrate/main.nf

@@ -1,4 +1,4 @@
-process RELATIVE_MUTRATE {
+process SELECT_MUTRATES {
     tag "$meta.id"
     label 'process_single'
 
@@ -12,20 +12,21 @@ process RELATIVE_MUTRATE {
     tuple val(meta), path(mutation_rates)
 
     output:
-    tuple val(meta), path("*.rel_mutrates.tsv") , emit: mutrate
-    path  "versions.yml"                        , emit: versions
+    tuple val(meta), path("*.gene_mutrates.tsv") , emit: mutrate
+    path  "versions.yml"                         , emit: versions
 
 
     when:
     task.ext.when == null || task.ext.when
 
     script:
-    def args = task.ext.args ?: ''
     def prefix = task.ext.prefix ?: "${meta.id}"
+    def mode = task.ext.mode ?: "mutations"
     """
-    omega_compute_relative_mutrate.py \\
+    omega_select_mutrate.py \\
                 --mutrates ${mutation_rates} \\
-                --output ${prefix}.rel_mutrates.tsv;
+                --output ${prefix}.gene_mutrates.tsv \\
+                --mode ${mode};
 
     cat <<-END_VERSIONS > versions.yml
     "${task.process}":
@@ -35,9 +36,8 @@ process RELATIVE_MUTRATE {
 
     stub:
     def prefix = task.ext.prefix ?: "all_samples"
-    def panel_version = task.ext.panel_version ?: "${meta2.id}"
     """
-    touch ${prefix}.rel_mutrates.tsv
+    touch ${prefix}.gene_mutrates.tsv
 
     cat <<-END_VERSIONS > versions.yml
     "${task.process}":

subworkflows/local/omega/main.nf

@@ -57,6 +57,10 @@ workflow OMEGA_ANALYSIS{
     .join( profile )
     .set{ muts_n_depths_n_profile }
 
+    Channel.of([ [ id: "all_samples" ] ])
+    .join( profile ).first()
+    .set{ all_samples_mut_profile }
+
 
     if (params.omega_hotspots){
         EXPANDREGIONS(panel, hotspots_file)
@@ -68,8 +72,11 @@ workflow OMEGA_ANALYSIS{
         json_hotspots = bedfile
     }
 
-    // FIXME: here I am using bedfile as a dummy value channel
-    PREPROCESSING( muts_n_depths_n_profile, expanded_panel, bedfile)
+    // FIXME here I am using bedfile as a dummy value channel
+    PREPROCESSING( muts_n_depths_n_profile,
+                    expanded_panel,
+                    bedfile,
+                    all_samples_mut_profile)
     ch_versions = ch_versions.mix(PREPROCESSING.out.versions)
 
     PREPROCESSING.out.mutabs_n_mutations_tsv
@@ -98,14 +105,19 @@ workflow OMEGA_ANALYSIS{
 
     if (params.omega_globalloc) {
 
-        PREPROCESSINGGLOBALLOC( muts_n_depths_n_profile, expanded_panel, relative_mutationrates.first() )
+        PREPROCESSINGGLOBALLOC(muts_n_depths_n_profile,
+                                expanded_panel,
+                                relative_mutationrates.first(),
+                                all_samples_mut_profile)
         ch_versions = ch_versions.mix(PREPROCESSINGGLOBALLOC.out.versions)
 
         PREPROCESSINGGLOBALLOC.out.mutabs_n_mutations_tsv
         .join( depth )
         .set{ preprocess_globalloc_n_depths }
 
-        ESTIMATORGLOBALLOC( preprocess_globalloc_n_depths, expanded_panel, GROUPGENES.out.json_genes.first())
+        ESTIMATORGLOBALLOC(preprocess_globalloc_n_depths,
+                            expanded_panel,
+                            GROUPGENES.out.json_genes.first())
         ch_versions = ch_versions.mix(ESTIMATORGLOBALLOC.out.versions)
 
         global_loc_results = ESTIMATORGLOBALLOC.out.results

workflows/deepcsa.nf

@@ -117,7 +117,8 @@ include { TABIX_BGZIPTABIX_QUERY    as DEPTHSNONPROTCONS    } from '../modules/n
 include { TABIX_BGZIPTABIX_QUERY    as DEPTHSINTRONSCONS    } from '../modules/nf-core/tabix/bgziptabixquery/main'
 include { TABIX_BGZIPTABIX_QUERY    as DEPTHSSYNONYMOUSCONS } from '../modules/nf-core/tabix/bgziptabixquery/main'
 
-include { RELATIVE_MUTRATE          as RELMUTRATE           } from '../modules/local/relative_mutrate/main'
+include { SELECT_MUTRATES           as SYNMUTRATE           } from '../modules/local/select_mutrate/main'
+include { SELECT_MUTRATES           as SYNMUTREADSRATE      } from '../modules/local/select_mutrate/main'
 include { DNA_2_PROTEIN_MAPPING     as DNA2PROTEINMAPPING   } from '../modules/local/dna2protein/main'
 
 
@@ -257,8 +258,11 @@ workflow DEEPCSA{
         .join( MUTRATESYNONYMOUS.out.mutrates )
         .set{ all_samples_syn_mutrate }
 
-        RELMUTRATE(all_samples_syn_mutrate)
-        ch_versions = ch_versions.mix(RELMUTRATE.out.versions)
+        SYNMUTRATE(all_samples_syn_mutrate)
+        ch_versions = ch_versions.mix(SYNMUTRATE.out.versions)
+
+        SYNMUTREADSRATE(all_samples_syn_mutrate)
+        ch_versions = ch_versions.mix(SYNMUTREADSRATE.out.versions)
 
 
         // Concatenate all outputs into a single file
@@ -449,7 +453,7 @@ workflow DEEPCSA{
                     CREATEPANELS.out.exons_consensus_panel,
                     custom_groups_table,
                     hotspots_bed_file,
-                    RELMUTRATE.out.mutrate
+                    SYNMUTRATE.out.mutrate
                     )
             positive_selection_results = positive_selection_results.join(OMEGA.out.results, remainder: true)
             positive_selection_results = positive_selection_results.join(OMEGA.out.results_global, remainder: true)
@@ -463,7 +467,7 @@ workflow DEEPCSA{
                         CREATEPANELS.out.exons_consensus_panel,
                         custom_groups_table,
                         hotspots_bed_file,
-                        RELMUTRATE.out.mutrate
+                        SYNMUTREADSRATE.out.mutrate
                         )
             positive_selection_results = positive_selection_results.join(OMEGAMULTI.out.results, remainder: true)
             positive_selection_results = positive_selection_results.join(OMEGAMULTI.out.results_global, remainder: true)
@@ -477,7 +481,7 @@ workflow DEEPCSA{
                             CREATEPANELS.out.exons_consensus_panel,
                             custom_groups_table,
                             hotspots_bed_file,
-                            RELMUTRATE.out.mutrate
+                            SYNMUTRATE.out.mutrate
                             )
             ch_versions = ch_versions.mix(OMEGANONPROT.out.versions)
 
@@ -488,7 +492,7 @@ workflow DEEPCSA{
                                 CREATEPANELS.out.exons_consensus_panel,
                                 custom_groups_table,
                                 hotspots_bed_file,
-                                RELMUTRATE.out.mutrate
+                                SYNMUTREADSRATE.out.mutrate
                                 )
             ch_versions = ch_versions.mix(OMEGANONPROTMULTI.out.versions)
         }