Add SNV variant scoring by MIVMIR, GICAM models

assets/rank_model_mivmir.ini

@@ -0,0 +1,121 @@
+[Version]
+  version = 1.39
+  name = rank_model_mivmir
+
+[Categories]
+
+ [[inheritance_models]]
+   category_aggregation = min
+
+ [[variant_call_quality_filter]]
+   category_aggregation = sum
+
+[model_score]
+  category = variant_call_quality_filter
+  data_type = integer
+  description = Inheritance model score
+  field = INFO
+  info_key = ModelScore
+  record_rule = min
+  separators = ',',':',
+
+  [[not_reported]]
+    score = 0
+
+  [[low_qual]]
+    score = -5
+    lower = 0
+    upper = 10
+
+  [[medium_qual]]
+    score = -2
+    lower = 10
+    upper = 20
+
+  [[high_qual]]
+    score = 0
+    lower = 20
+    upper = 300
+
+[genetic_models]
+  category = inheritance_models
+  data_type = string
+  description = Inheritance models followed for the variant
+  field = INFO
+  info_key = GeneticModels
+  record_rule = max
+  separators = ',', ':', '|',
+
+ [[ad]]
+   priority = 1
+   score = 1
+   string = 'AD'
+
+ [[ad_dn]]
+   score = 1
+   priority = 1
+   string = 'AD_dn'
+
+ [[ar]]
+   score = 1
+   priority = 1
+   string = 'AR_hom'
+
+ [[ar_dn]]
+   score = 1
+   priority = 1
+   string = 'AR_hom_dn'
+
+ [[ar_comp]]
+   score = 1
+   priority = 1
+   string = 'AR_comp'
+
+ [[ar_comp_dn]]
+   score = 1
+   priority = 1
+   string = 'AR_comp_dn'
+
+ [[xr]]
+   score = 1
+   priority = 1
+   string = 'XR'
+
+ [[xr_dn]]
+   score = 1
+   priority = 1
+   string = 'XR_dn'
+
+ [[xd]]
+   score = 1
+   priority = 1
+   string = 'XD'
+
+ [[xd_dn]]
+   score = 1
+   priority = 1
+   string = 'XD_dn'
+
+ [[not_reported]]
+   score = -12
+
+[filter]
+  category = variant_call_quality_filter
+  data_type = string
+  description = The filters for the variant
+  field = FILTER
+  record_rule = min
+  separators = ';',
+
+  [[not_reported]]
+    score = 0
+
+  [[pass]]
+    score = 3
+    priority = 1
+    string = 'PASS'
+
+  [[dot]]
+    score = 3
+    priority = 2
+    string = '.'

conf/modules/rank_variants.config

@@ -107,6 +107,24 @@ process {
         ]
     }
 
+    withName: '.*RANK_VARIANTS_SNV:TABIX_BGZIPTABIX_GICAM' {
+        ext.prefix = { "${meta.id}_snv_ranked_gicam_${meta.set}" }
+    }
+
+    withName: '.*RANK_VARIANTS_SNV:TABIX_BGZIP_GENMOD_GICAM' {
+        ext.prefix = { "${meta.id}_snv_ranked_${meta.set}" }
+        publishDir = [
+            path: { "${params.outdir}/rank_and_filter" },
+            mode: params.publish_dir_mode,
+            saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+        ]
+    }
+
+    withName: '.*RANK_VARIANTS_SNV:BCFTOOLS_MERGE_GENMOD_GICAM' {
+        ext.args   = { "--columns VrsModelPrediction,VrsModelExplanation,GICAM" }
+        ext.prefix = { "${meta.id}_snv_ranked_${meta.set}" }
+    }
+
     withName: '.*RANK_VARIANTS_SNV:TABIX_TABIX' {
         publishDir = [
             path: { "${params.outdir}/rank_and_filter" },

conf/test.config

@@ -52,6 +52,7 @@ params {
     reduced_penetrance              = params.pipelines_testdata_base_path + '/reference/reduced_penetrance.tsv'
     score_config_mt                 = params.pipelines_testdata_base_path + '/reference/rank_model_snv.ini'
     score_config_snv                = params.pipelines_testdata_base_path + '/reference/rank_model_snv.ini'
+    score_config_genmod_gicam_snv   = 'assets/rank_model_mivmir.ini'
     score_config_sv                 = params.pipelines_testdata_base_path + '/reference/rank_model_sv.ini'
     svdb_query_dbs                  = params.pipelines_testdata_base_path + '/reference/svdb_querydb_files.csv'
     target_bed                      = params.pipelines_testdata_base_path + '/reference/target.bed'

conf/test_full.config

@@ -51,6 +51,7 @@ params {
     reduced_penetrance              = params.pipelines_testdata_base_path + 'reference/reduced_penetrance.tsv'
     score_config_mt                 = params.pipelines_testdata_base_path + 'reference/rank_model_snv.ini'
     score_config_snv                = params.pipelines_testdata_base_path + 'reference/rank_model_snv.ini'
+    score_config_genmod_gicam_snv   = 'assets/rank_model_mivmir.ini'
     score_config_sv                 = params.pipelines_testdata_base_path + 'reference/rank_model_sv.ini'
     svdb_query_dbs                  = params.pipelines_testdata_base_path + 'reference/svdb_querydb_files.csv'
     target_bed                      = params.pipelines_testdata_base_path + 'reference/target.bed'

modules/local/gicam/main.nf

@@ -0,0 +1,40 @@
+process GICAM_INFER {
+    // https://github.com/Clinical-Genomics/rdds/tree/master/src/rdds/gicam
+
+    tag "${meta.id}"
+    label 'process_single'
+
+    container "docker.io/clinicalgenomics/rdds_mivmir:v1.12.0-rc3"
+
+    beforeScript "mkdir ${task.workDir}/rdds-tmp"
+    afterScript "rm -r ${task.workDir}/rdds-tmp"
+    containerOptions {[
+        workflow.containerEngine.equals("singularity") ? "--bind ${task.workDir}/rdds-tmp:/rdds/tmp" : "",
+        workflow.containerEngine.equals("docker") ? "--tmpfs /rdds/tmp": "",
+        ""
+    ].minus("").join(" ")}
+
+    input:
+    tuple val(meta), path(input_vcf)
+
+    output:
+    tuple val(meta), path('*-predictions.vcf'), emit: vcf
+    path "versions.yml",                        emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def VERSION = 'v1.12.0-rc3'
+    """
+    . /opt/pyenv/bin/activate
+    export PYTHONPATH=/rdds/src
+    python3 -m rdds.gicam infer-vcf --cpu_cores ${task.cpus} ${input_vcf}
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        gicam: ${VERSION}
+        python: \$(python --version | sed 's/Python //g')
+    END_VERSIONS
+    """
+}

modules/local/mivmir/main.nf

@@ -0,0 +1,40 @@
+process MIVMIR_INFER {
+    // https://github.com/Clinical-Genomics/rdds/tree/master/src/rdds/variant_rank_score
+
+    tag "${meta.id}"
+    label 'process_single'
+
+    container "docker.io/clinicalgenomics/rdds_mivmir:v1.12.0-rc3"
+
+    beforeScript "mkdir ${task.workDir}/rdds-tmp"
+    afterScript "rm -r ${task.workDir}/rdds-tmp"
+    containerOptions {[
+        workflow.containerEngine.equals("singularity") ? "--bind ${task.workDir}/rdds-tmp:/rdds/tmp" : "",
+        workflow.containerEngine.equals("docker") ? "--tmpfs /rdds/tmp": "",
+        ""
+    ].minus("").join(" ")}
+
+    input:
+    tuple val(meta), path(input_vcf)
+
+    output:
+    tuple val(meta), path('*-predictions.vcf'), emit: vcf
+    path "versions.yml",                        emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def VERSION = 'v1.12.0-rc3'
+    """
+    . /opt/pyenv/bin/activate
+    export PYTHONPATH=/rdds/src
+    python3 -m rdds.variant_rank_score predict-on-vcf  --cpu_cores ${task.cpus} ${input_vcf}
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        mivmir: ${VERSION}
+        python: \$(python --version | sed 's/Python //g')
+    END_VERSIONS
+    """
+}

nextflow.config

@@ -68,6 +68,7 @@ params {
     sequence_dictionary             = null
     score_config_mt                 = null
     score_config_snv                = null
+    score_config_genmod_gicam_snv   = 'assets/rank_model_mivmir.ini'
     score_config_sv                 = null
     sdf                             = null
     svdb_query_bedpedbs             = null

nextflow_schema.json

@@ -350,6 +350,13 @@
                     "fa_icon": "fas fa-file",
                     "description": "SV rank model config file for genmod."
                 },
+              "score_config_genmod_gicam_snv": {
+                "type": "string",
+                "exists": true,
+                "format": "path",
+                "fa_icon": "fas fa-file",
+                "description": "Rank model config file for genmod in conjunction with MIVMIR-GICAM scoring."
+              },
                 "sdf": {
                     "type": "string",
                     "exists": true,

subworkflows/local/rank_variants.nf

@@ -2,21 +2,29 @@
 // A subworkflow to score and rank variants.
 //
 
-include { GENMOD_ANNOTATE  } from '../../modules/nf-core/genmod/annotate/main'
-include { GENMOD_MODELS    } from '../../modules/nf-core/genmod/models/main'
-include { GENMOD_SCORE     } from '../../modules/nf-core/genmod/score/main'
-include { GENMOD_COMPOUND  } from '../../modules/nf-core/genmod/compound/main'
-include { BCFTOOLS_SORT    } from '../../modules/nf-core/bcftools/sort/main'
-include { TABIX_BGZIP      } from '../../modules/nf-core/tabix/bgzip/main'
-include { TABIX_TABIX      } from '../../modules/nf-core/tabix/tabix/main'
+include { GENMOD_ANNOTATE                                       } from '../../modules/nf-core/genmod/annotate/main'
+include { GENMOD_MODELS                                         } from '../../modules/nf-core/genmod/models/main'
+include { GENMOD_SCORE                                          } from '../../modules/nf-core/genmod/score/main'
+include { GENMOD_SCORE as GENMOD_SCORE_FOR_GICAM                } from '../../modules/nf-core/genmod/score/main'
+include { GENMOD_COMPOUND                                       } from '../../modules/nf-core/genmod/compound/main'
+include { MIVMIR_INFER                                          } from '../../modules/local/mivmir/main'
+include { GICAM_INFER                                           } from '../../modules/local/gicam/main'
+include { BCFTOOLS_SORT                                         } from '../../modules/nf-core/bcftools/sort/main'
+include { TABIX_BGZIP                                           } from '../../modules/nf-core/tabix/bgzip/main'
+include { TABIX_BGZIP as TABIX_BGZIP_GENMOD_GICAM               } from '../../modules/nf-core/tabix/bgzip/main'
+include { TABIX_BGZIPTABIX as TABIX_BGZIPTABIX_GICAM            } from '../../modules/nf-core/tabix/bgziptabix/main'
+include { TABIX_TABIX                                           } from '../../modules/nf-core/tabix/tabix/main'
+include { BCFTOOLS_ANNOTATE as BCFTOOLS_MERGE_GENMOD_GICAM      } from '../../modules/nf-core/bcftools/annotate/main'
 
 workflow RANK_VARIANTS {
 
     take:
-        ch_vcf                // channel: [mandatory] [ val(meta), path(vcf) ]
-        ch_pedfile            // channel: [mandatory] [ path(ped) ]
-        ch_reduced_penetrance // channel: [mandatory] [ path(pentrance) ]
-        ch_score_config       // channel: [mandatory] [ path(ini) ]
+        ch_vcf                          // channel: [mandatory] [ val(meta), path(vcf) ]
+        ch_pedfile                      // channel: [mandatory] [ path(ped) ]
+        ch_reduced_penetrance           // channel: [mandatory] [ path(pentrance) ]
+        ch_score_config                 // channel: [mandatory] [ path(ini) ]
+        ch_genmod_gicam_score_config    // channel: [mandatory] [ path(ini) ]
+        rank_with_mivmir_gicam          // value
 
     main:
         ch_versions = Channel.empty()
@@ -31,13 +39,34 @@ workflow RANK_VARIANTS {
 
         GENMOD_SCORE(ch_score_in, ch_score_config)
 
-        GENMOD_COMPOUND(GENMOD_SCORE.out.vcf)
-
-        BCFTOOLS_SORT(GENMOD_COMPOUND.out.vcf) // SV file needs to be sorted before indexing
+        // Run MIVMIR - GICAM scoring (not supported for MT SNVs and SVs)
+        if (rank_with_mivmir_gicam) {
+            GENMOD_SCORE_FOR_GICAM(ch_score_in, ch_genmod_gicam_score_config)
+            MIVMIR_INFER(GENMOD_SCORE_FOR_GICAM.out.vcf)
+            GICAM_INFER(MIVMIR_INFER.out.vcf)
+            TABIX_BGZIPTABIX_GICAM(GICAM_INFER.out.vcf)
+        }
 
+        GENMOD_COMPOUND(GENMOD_SCORE.out.vcf)
         TABIX_BGZIP(GENMOD_COMPOUND.out.vcf) //run only for SNVs
 
-        ch_vcf = TABIX_BGZIP.out.output.mix(BCFTOOLS_SORT.out.vcf)
+        // Merge Genmod and MIVMIR-GICAM scores
+        if (rank_with_mivmir_gicam) {
+            TABIX_BGZIP.out.output
+                .join(TABIX_BGZIPTABIX_GICAM.out.gz_tbi, failOnMismatch: true)
+                .map {meta, vcf_genmod, vcf_gicam, vcf_index_gicam -> return [ meta, vcf_genmod, [], vcf_gicam, vcf_index_gicam ]}
+                .set {ch_merge_genmod_gicam}
+            BCFTOOLS_MERGE_GENMOD_GICAM(ch_merge_genmod_gicam, [])
+            TABIX_BGZIP_GENMOD_GICAM(BCFTOOLS_MERGE_GENMOD_GICAM.out.vcf)
+        }
+
+        BCFTOOLS_SORT(GENMOD_COMPOUND.out.vcf) // SV file needs to be sorted before indexing
+        // Mix SNVs and SVs
+        if (rank_with_mivmir_gicam) {
+            ch_vcf = TABIX_BGZIP_GENMOD_GICAM.out.output.mix(BCFTOOLS_SORT.out.vcf)
+        } else {
+            ch_vcf = TABIX_BGZIP.out.output.mix(BCFTOOLS_SORT.out.vcf)
+        }
 
         TABIX_TABIX (ch_vcf)
 

workflows/raredisease.nf

@@ -203,6 +203,8 @@ workflow RAREDISEASE {
                                                                             : Channel.value([])
     ch_score_config_snv         = params.score_config_snv                   ? Channel.fromPath(params.score_config_snv).collect()
                                                                             : Channel.value([])
+    ch_genmod_gicam_score_config = params.score_config_genmod_gicam_snv     ? Channel.fromPath(params.score_config_genmod_gicam_snv).collect()
+                                                                            : Channel.value([])
     ch_score_config_sv          = params.score_config_sv                    ? Channel.fromPath(params.score_config_sv).collect()
                                                                             : Channel.value([])
     ch_sdf                      = params.sdf                                ? Channel.fromPath(params.sdf).map{it -> [[id:it.simpleName],it]}.collect()
@@ -547,7 +549,9 @@ workflow RAREDISEASE {
                 ch_ranksnv_nuclear_in,
                 ch_pedfile,
                 ch_reduced_penetrance,
-                ch_score_config_snv
+                ch_score_config_snv,
+                ch_genmod_gicam_score_config,
+                true
             )
             ch_versions = ch_versions.mix(RANK_VARIANTS_SNV.out.versions)
         }
@@ -615,7 +619,9 @@ workflow RAREDISEASE {
                 ch_ranksnv_mt_in,
                 ch_pedfile,
                 ch_reduced_penetrance,
-                ch_score_config_mt
+                ch_score_config_mt,
+                ch_genmod_gicam_score_config,
+                false
             )
             ch_versions = ch_versions.mix(RANK_VARIANTS_MT.out.versions)
         }
@@ -706,7 +712,9 @@ workflow RAREDISEASE {
                 ch_ranksnv_sv_in,
                 ch_pedfile,
                 ch_reduced_penetrance,
-                ch_score_config_sv
+                ch_score_config_sv,
+                ch_genmod_gicam_score_config,
+                false
             )
             ch_versions = ch_versions.mix(RANK_VARIANTS_SV.out.versions)
         }