feat: Add genome subsampling (#63)

* feat: Add module to extract distances

* feat: Incorporate subsampling into main workflows

* test: Add test for subsampling subworkflow

* chore: Remove commented out code

* test: Add subsetting to poppunk test

* test: Add subsetting to test profile

* fix: Change way of setting distance output

* refactor: Extract distances should be part of subsetting

* test: Change test to use popdb

* test: Fix poppunk trace size

* refactor: Update default threshold

* feat: Write removed genomes to a file

* test: Update similarity

* feat: Add poppunk entry

* fix: Change variable declaration

* refactor: Change calculus

* test: Update subsampling test

conf/modules.config

@@ -252,6 +252,14 @@ process {
         ]
     }
 
+    withName: POPPUNK_EXTRACT_DISTANCES {
+        publishDir = [
+            path: { "${params.outdir}/poppunk_results/distances/" },
+            mode: params.publish_dir_mode,
+            saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+        ]
+    }
+
     //  pangenomics
 
     withName: PANAROO_RUN {

conf/test.config

@@ -19,6 +19,7 @@ params {
     db_cache = false
     use_full_alignment = false
     use_fasttree = true
+    enable_subsetting = true
     use_prokka = true
     skip_kraken = true
     skip_poppunk = true

main.nf

@@ -35,6 +35,7 @@ include { ARETE } from './workflows/arete'
 include { ASSEMBLY } from './workflows/arete'
 include { ANNOTATION } from './workflows/arete'
 include { QUALITYCHECK } from './workflows/arete'
+include { POPPUNK } from './workflows/arete'
 
 
 //
@@ -57,6 +58,11 @@ workflow annotation {
 workflow assembly_qc {
     QUALITYCHECK()
 }
+
+workflow poppunk {
+    POPPUNK()
+}
+
 /*
 ========================================================================================
     RUN ALL WORKFLOWS

modules/local/poppunk/extractdistances/main.nf

@@ -0,0 +1,32 @@
+process POPPUNK_EXTRACT_DISTANCES {
+    label 'process_high'
+
+    conda (params.enable_conda ? "bioconda::poppunk=2.6.0" : null)
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/poppunk:2.6.0--py39h9b916c0_0':
+        'quay.io/biocontainers/poppunk:2.6.0--py39h9b916c0_0' }"
+
+    input:
+
+    path poppunk_db
+
+    output:
+
+    path "poppunk_db_distances.tsv", emit: poppunk_distances
+    path "versions.yml", emit: versions
+
+    script:
+    def args = task.ext.args ?: ''
+
+    """
+    poppunk_extract_distances.py \\
+        --distances $poppunk_db/${poppunk_db}.dists \\
+        --output poppunk_db_distances.tsv \\
+        $args
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        poppunk: \$(echo \$(poppunk --version 2>&1) | sed 's/^.*poppunk //;')
+    END_VERSIONS
+    """
+}

modules/local/poppunk/extractdistances/meta.yml

@@ -0,0 +1,33 @@
+name: poppunk_extract_distances
+description: Extract pairwise distances between genomes
+keywords:
+  - genomes
+  - poppunk
+  - database
+  - distances
+tools:
+  - poppunk:
+    description: Population partitioning using nucleotide k-mers
+    homepage: https://poppunk.net/
+    documentation: https://poppunk.readthedocs.io/en/latest/
+    tool_dev_url: https://github.com/bacpop/PopPUNK
+    doi: "doi:10.1101/gr.241455.118"
+    licence: ["Apache-2.0"]
+
+input:
+  - poppunk_db:
+      type: directory
+      description: Directory containing the PopPunk database
+
+output:
+  - poppunk_distances:
+      type: file
+      description: TSV file containing pairwise distances
+      pattern: "poppunk_db_distances.tsv"
+  - versions:
+      type: file
+      description: File containing software versions
+      pattern: "versions.yml"
+
+authors:
+  - "@jvfe"

nextflow.config

@@ -20,6 +20,11 @@ params {
     poppunk_model              = null
     run_poppunk_qc             = false
 
+    // Subsampling parameters
+    enable_subsetting          = false
+    core_similarity            = 99.9
+    accessory_similarity       = 99
+
     // Phylogenomics parameters
     use_full_alignment         = false
     use_fasttree               = true

nextflow_schema.json

@@ -158,6 +158,23 @@
                     "type": "boolean",
                     "description": "Whether to run the QC step for PopPunk",
                     "fa_icon": "fas fa-industry"
+                },
+                "enable_subsetting": {
+                    "type": "boolean",
+                    "fa_icon": "fas fa-tasks",
+                    "description": "Enable subsetting workflow based on genome similarity"
+                },
+                "core_similarity": {
+                    "type": "number",
+                    "default": 99.99,
+                    "fa_icon": "fas fa-clone",
+                    "description": "Similarity threshold for core genomes"
+                },
+                "accessory_similarity": {
+                    "type": "number",
+                    "default": 99,
+                    "fa_icon": "far fa-clone",
+                    "description": "Similarity threshold for accessory genes"
                 }
             }
         },

subworkflows/local/poppunk.nf

@@ -51,4 +51,5 @@ workflow RUN_POPPUNK {
         poppunk_version = POPPUNK_FITMODEL.out.versions.ifEmpty(null)
         poppunk_results = poppunk_results
         poppunk_visualisations = POPPUNK_VISUALISE.out.poppunk_visualizations
+        poppunk_db = poppunk_db
 }

subworkflows/local/subsample.nf

@@ -0,0 +1,36 @@
+include { POPPUNK_EXTRACT_DISTANCES } from '../../modules/local/poppunk/extractdistances/main'
+
+workflow SUBSET_GENOMES {
+
+    take:
+        genome_assemblies
+        poppunk_db
+
+    main:
+
+        POPPUNK_EXTRACT_DISTANCES(poppunk_db)
+        POPPUNK_EXTRACT_DISTANCES.out.poppunk_distances.set{ poppunk_distances }
+
+        def core_threshold = 100.0 - params.core_similarity
+        def accessory_threshold = 100.0 - params.accessory_similarity
+
+        poppunk_distances
+            .splitCsv(header: true, sep: '\t')
+            .filter { row -> (row.Core.toFloat() * 100) < core_threshold  && (row.Accessory.toFloat() * 100) < accessory_threshold }
+            .map { row -> row.Query }
+            .set { genomes_to_remove }
+
+        genomes_to_remove
+            .unique()
+            .collectFile(newLine: true)
+            .collectFile(name: 'removed_genomes.txt', storeDir: "${params.outdir}/poppunk_results")
+
+        genome_assemblies
+            .combine (genomes_to_remove.collect().map { [it] })
+            .filter { meta, path, to_remove -> !(meta.id in to_remove) }
+            .map { it[0, 1] }
+            .set { filtered_genomes }
+
+    emit:
+        filtered_genomes = filtered_genomes
+}

tests/subworkflows/local/subsample.nf.test

@@ -0,0 +1,44 @@
+nextflow_workflow {
+
+    name "Test Workflow SUBSET_GENOMES"
+    script "subworkflows/local/subsample.nf"
+    workflow "SUBSET_GENOMES"
+
+    test("Subsampling subworkflow runs without failures") {
+
+        when {
+            params {
+                // define parameters here. Example:
+                outdir = "$outputDir"
+                core_similarity = 99
+                accessory_similarity = 95
+            }
+            workflow {
+                """
+                // define inputs of the workflow here. Example:
+                input[0] = Channel.of(
+                    [[id:'SRR14022735'], "$baseDir/test/SRR14022735_T1.scaffolds.fa"],
+                    [[id:'SRR14022737'], "$baseDir/test/SRR14022737_T1.scaffolds.fa"],
+                    [[id:'SRR14022754'], "$baseDir/test/SRR14022754_T1.scaffolds.fa"],
+                    [[id:'SRR14022764'], "$baseDir/test/SRR14022764_T1.scaffolds.fa"],
+                )
+                input[1] = Channel.fromPath("$baseDir/test/popdb")
+                """
+            }
+        }
+
+        then {
+            assert workflow.success
+            assert workflow.out.filtered_genomes.size() == 3
+            assert workflow.out.filtered_genomes.toList() == [
+                [['id':'SRR14022735'], "$baseDir/test/SRR14022735_T1.scaffolds.fa"],
+                [['id':'SRR14022754'], "$baseDir/test/SRR14022754_T1.scaffolds.fa"],
+                [['id':'SRR14022764'], "$baseDir/test/SRR14022764_T1.scaffolds.fa"]
+            ]
+
+
+        }
+
+    }
+
+}