Add support for variant scoring by MIVMIR, GICAM models

Signed-off-by: Tor Björgen <tor.bjorgen@scilifelab.se>

Author: torbjorgen

CHANGELOG.md

@@ -30,6 +30,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Env variable NXF_SINGULARITY_NEW_PID_NAMESPACE = false to accommodate hisat2 running with latest Nextflow and Singularity [#775](https://github.com/nf-core/raredisease/pull/775)
 - Parameter `exclude_alt` to filter alignments to alt/unplaced contigs after alignment using samtools view, retaining only primary chromosomes (GRCh37: 1-22,X,Y,MT / GRCh38: chr1-chr22,chrX,chrY,chrM). Note that enabling this will restrict variant calling to these chromosomes [#803](https://github.com/nf-core/raredisease/pull/803)]
 - Parameters `save_all_mapped_as_cram` and `save_noalt_mapped_as_cram` to replace `save_mapped_as_cram`, allowing independent control over publishing unfiltered and alt-filtered alignment files as CRAM [#807](https://github.com/nf-core/raredisease/pull/807)
+- Add SNV scoring by MIVMIR, GICAM models [#816](https://github.com/nf-core/raredisease/pull/812)
 
 ### `Changed`
 

README.md

@@ -175,7 +175,7 @@ For more details about the output files and reports, please refer to the
 
 nf-core/raredisease was written in a collaboration between the Clinical Genomics nodes in Sweden, with major contributions from [Ramprasad Neethiraj](https://github.com/ramprasadn), [Anders Jemt](https://github.com/jemten), [Lucia Pena Perez](https://github.com/Lucpen), and [Mei Wu](https://github.com/projectoriented) at Clinical Genomics Stockholm.
 
-Additional contributors were [Sima Rahimi](https://github.com/sima-r), [Gwenna Breton](https://github.com/Gwennid) and [Emma Västerviga](https://github.com/EmmaCAndersson) (Clinical Genomics Gothenburg); [Halfdan Rydbeck](https://github.com/hrydbeck) and [Lauri Mesilaakso](https://github.com/ljmesi) (Clinical Genomics Linköping); [Subazini Thankaswamy Kosalai](https://github.com/sysbiocoder) (Clinical Genomics Örebro); [Annick Renevey](https://github.com/rannick), [Peter Pruisscher](https://github.com/peterpru) and [Eva Caceres](https://github.com/fevac) (Clinical Genomics Stockholm); [Ryan Kennedy](https://github.com/ryanjameskennedy) (Clinical Genomics Lund); [Anders Sune Pedersen](https://github.com/asp8200) (Danish National Genome Center) and [Lucas Taniguti](https://github.com/lmtani).
+Additional contributors were [Sima Rahimi](https://github.com/sima-r), [Gwenna Breton](https://github.com/Gwennid) and [Emma Västerviga](https://github.com/EmmaCAndersson) (Clinical Genomics Gothenburg); [Halfdan Rydbeck](https://github.com/hrydbeck) and [Lauri Mesilaakso](https://github.com/ljmesi) (Clinical Genomics Linköping); [Subazini Thankaswamy Kosalai](https://github.com/sysbiocoder) (Clinical Genomics Örebro); [Annick Renevey](https://github.com/rannick), [Peter Pruisscher](https://github.com/peterpru), [Eva Caceres](https://github.com/fevac) and [Tor Björgen](https://github.com/torbjorgen) (Clinical Genomics Stockholm); [Ryan Kennedy](https://github.com/ryanjameskennedy) (Clinical Genomics Lund); [Anders Sune Pedersen](https://github.com/asp8200) (Danish National Genome Center) and [Lucas Taniguti](https://github.com/lmtani).
 
 We thank the nf-core community for their extensive assistance in the development of this pipeline.
 

conf/modules/rank_variants.config

@@ -83,6 +83,31 @@ process {
         ext.prefix = { "${meta.id}_snv_ranked_${meta.set}" }
     }
 
+    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 MivmirScore,MivmirExplanation,GicamScore" }
+        ext.prefix = { "${meta.id}_snv_ranked_${meta.set}" }
+    }
+
+    withName: '.*RANK_VARIANTS_SNV:TABIX_TABIX' {
+        publishDir = [
+            path: { "${params.outdir}/rank_and_filter" },
+            mode: params.publish_dir_mode,
+            saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+        ]
+    }
 }
 
 //

docs/output.md

@@ -69,6 +69,7 @@ The pipeline is built using [Nextflow](https://www.nextflow.io/) and processes d
     - [Filtering and ranking](#filtering-and-ranking)
       - [Filter_vep](#filter_vep)
       - [GENMOD](#genmod)
+      - [MIVMIR, GICAM](#mivmir-gicam)
     - [Mobile element analysis](#mobile-element-analysis)
       - [Calling mobile elements](#calling-mobile-elements)
       - [Annotating mobile elements](#annotating-mobile-elements)
@@ -527,6 +528,14 @@ We recommend using vcfanno to annotate SNVs with precomputed CADD scores (files
 
 </details>
 
+#### MIVMIR, GICAM
+[MIVMIR](../modules/local/mivmir/meta.yml) and [GICAM](../modules/local/gicam/meta.yml) are two machine learning models used to
+infer a pathogenicity score for SNVs. In essence, MIVMIR infer SNV pathogenicity and GICAM improves precision for
+duo, trio, ... analysis. MIVMIR, GICAM can be enabled by setting the `--rank_with_mivmir_gicam` feature flag and
+adds annotations `INFO/MivmirScore`, `INFO/MivmirExplanation`, `INFO/GicamScore`.
+Only `<case_id>_snv_ranked_<research|clinical>.vcf.gz` contains the above annotations.
+Refer to the module documentation `.yml` for more information on required inputs and output formats.
+
 ### Mobile element analysis
 
 #### Calling mobile elements

docs/usage.md

@@ -287,7 +287,7 @@ The mandatory and optional parameters for each category are tabulated below.
 | vcfanno_resources<sup>2</sup>        | vcfanno_lua                                    |
 | vcfanno_toml<sup>3</sup>             | vep_filters/vep_filters_scout_fmt<sup>10</sup> |
 | vep_cache_version                    | cadd_resources<sup>11</sup>                    |
-| vep_cache<sup>4</sup>                |                                                |
+| vep_cache<sup>4</sup>                | rank_with_mivmir_gicam<sup>12</sup>            |
 | gnomad_af<sup>5</sup>                |                                                |
 | score_config_snv<sup>6</sup>         |                                                |
 | variant_consequences_snv<sup>7</sup> |                                                |
@@ -307,6 +307,7 @@ no header and the following columns: `CHROM POS REF_ALLELE,ALT_ALLELE AF`. Sampl
 <sup>9</sup>Used by GENMOD while modeling the variants. Contains a list of loci that show [reduced penetrance](https://medlineplus.gov/genetics/understanding/inheritance/penetranceexpressivity/) in people. Sample file [here](https://github.com/nf-core/test-datasets/blob/raredisease/reference/reduced_penetrance.tsv).<br />
 <sup>10</sup> This file contains a list of candidate genes (with [HGNC](https://www.genenames.org/) IDs) that is used to split the variants into candidate variants and research variants. Research variants contain all the variants, while candidate variants are a subset of research variants and are associated with candidate genes. Sample file [here](https://github.com/nf-core/test-datasets/blob/raredisease/reference/hgnc.txt). Not required if `--skip_subworkflows generate_clinical_set` is set. To skip this splitting entirely, add `generate_clinical_set` to `--skip_subworkflows`.<br />
 <sup>11</sup>Path to a folder containing cadd annotations. Equivalent of the data/annotations/ folder described [here](https://github.com/kircherlab/CADD-scripts/#manual-installation), and it is used to calculate CADD scores for small indels. <br />
+<sup>12</sup> Enable variant scoring using MIVMIR, GICAM machine learning models.
 
 :::note
 We use CADD only to annotate small indels. To annotate SNVs with precomputed CADD scores, pass the file containing CADD scores as a resource to vcfanno instead. Files containing the precomputed CADD scores for SNVs can be downloaded from [here](https://cadd.gs.washington.edu/download) (download files listed under the description: "All possible SNVs of GRCh3<7/8>/hg3<7/8>")

main.nf

@@ -105,6 +105,7 @@ workflow NFCORE_RAREDISEASE {
     val_readcount_intervals
     val_reduced_penetrance
     val_rtg_truthvcfs
+    val_rank_with_mivmir_gicam
     val_run_mt_for_wes
     val_run_rtgvcfeval
     val_sambamba_regions
@@ -508,6 +509,7 @@ workflow NFCORE_RAREDISEASE {
         val_mt_subsample_rd,
         val_mt_subsample_seed,
         val_platform,
+        val_rank_with_mivmir_gicam,
         val_run_mt_for_wes,
         val_run_rtgvcfeval,
         val_sample_id_map,
@@ -619,6 +621,7 @@ workflow {
         params.readcount_intervals,
         params.reduced_penetrance,
         params.rtg_truthvcfs,
+        params.rank_with_mivmir_gicam,
         params.run_mt_for_wes,
         params.run_rtgvcfeval,
         params.sambamba_regions,

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_high'
+
+    container "docker.io/clinicalgenomics/rdds_mivmir:v1.12.0-rc6"
+
+    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-rc6'
+    """
+    . /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/gicam/meta.yml

@@ -0,0 +1,51 @@
+name: gicam
+description: Machine learning tool to improve precision for duo, trio, ... analysis
+keywords:
+  - score
+  - ranking
+  - gicam
+tools:
+  - mivmir:
+      description: Model for SNV ranking in conjunction with MIVMIR.
+        This model improves precision for duo, trio, ... analysis situations by reducing MIVMIR scores for variants that's
+        not following the appropriate GENMOD genetic inheritance model.
+        Applied to MIVMIR scores as a post-processing step.
+
+        VCF key inputs to the model are
+        - MivmirScore, (0, 1)
+        - RankScoreNormalized, (0, 1) as produced by GENMOD using rank_model_genmod_gicam.ini custom scoring config
+
+        The tool adds one key to the VCF
+          - GicamScore (0, 1) where 1.0 inferred pathogenic.
+      homepage: https://github.com/clinicalgenomics/rdds
+      documentation: https://github.com/Clinical-Genomics/rdds/tree/master/src/rdds/gicam
+      doi: ""
+      licence: ["MIT"]
+      identifier: ""
+input:
+  - - meta:
+        type: map
+        description: |
+          Groovy Map containing sample information
+          e.g. [ id:'test', single_end:false ]
+    - input_vcf:
+        type: file
+        description: vcf file
+        pattern: "*.{vcf}"
+        ontologies: []
+output:
+  vcf:
+    - - meta:
+          type: map
+          description: |
+            Groovy Map containing sample information
+            e.g. [ id:'test', single_end:false ]
+      - "*-predictions.vcf":
+          type: file
+          description: Scored output VCF file
+          pattern: "*.{vcf}"
+          ontologies: [ ]
+authors:
+  - "@torbjorgen"
+maintainers:
+  - "@torbjorgen"

modules/local/gicam/rank_model_genmod_gicam.ini

@@ -0,0 +1,121 @@
+[Version]
+  version = 1
+  name = rank_model_for_gicam
+
+[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 = '.'

modules/local/gicam/tests/main.nf.test

@@ -0,0 +1,31 @@
+nextflow_process {
+
+    name "Test Process GICAM_INFER"
+    script "modules/local/gicam/main.nf"
+    process "GICAM_INFER"
+
+    test("Test GICAM inference on annotated VCF") {
+
+        when {
+            params {
+            }
+            process {
+                """
+                input[0] = [
+                    [ id:'test', single_end:false ], // meta map
+                    // TODO@torbjorgen: Update data path to main repo
+                    // VCF as annotated by rank_variants subworkflow (including mivmir and genmod custom inheritance config), prior to gicam inference call
+                    file('https://raw.githubusercontent.com/torbjorgen/nfcore-test-datasets/mivmir-gicam-test-data/testdata/justhusky_snv_gicam.vcf', checkIfExists: true)
+                ]
+                """
+            }
+        }
+
+        then {
+            assert process.success
+            assert snapshot(process.out).match()
+        }
+
+    }
+
+}