Merge dev features to main

configs/modules.config

@@ -259,7 +259,7 @@ process {
     withName: 'MINIMAP2_ALIGN_READS' {
         tag        = { meta.assembly.build }
         ext.args   = "-x map-hifi"
-        ext.prefix = { "${meta.id}_" + (reads.getName().endsWith(".gz") ? reads.getBaseName(2) : reads.getBaseName()) }
+        ext.prefix = { "${meta.id}_${reads.getBaseName( reads.getName().endsWith(".gz") ? 2 : 1 )}" }
     }
     withName: 'MINIMAP2_ALIGN_ASSEMBLY_.*' {
         tag        = { meta.assembly.build }
@@ -307,13 +307,66 @@ process {
         ]
     }
 
+    // POLISH
+    withName: 'DVPOLISH_PBMM2_ALIGN' {
+        ext.args = '-N 1 -l 4000'
+        cpus = 20
+    }
+    // the samtools view process was modified and the
+    // qname is commented out in the process, but is used to link
+    // a bed file for -L filtering
+    withName: 'SAMTOOLS_VIEW' {
+        ext.prefix = {"${input.baseName}_${index.baseName}"}
+        cpus = 10
+        ext.args = {"-bh -F 2308 -M -L ${qname}"}
+    }
+    withName: 'SAMTOOLS_MERGE' {
+        ext.prefix = {"${meta.id}_${meta.mergeID}"}
+        cpus = 4
+    }
+    withName: 'DEEPVARIANT' {
+        cpus = 20
+        ext.args = {"--model_type=PACBIO"}
+        ext.prefix = {"${meta.id}_${meta.mergeID}"}
+    }
+    withName: 'BCFTOOLS_VIEW' {
+        cpus = 6
+        ext.args = {"-f 'PASS' -i 'GT=\"1/1\"' --no-version -Oz"}
+        ext.prefix = {"${meta.id}_${meta.mergeID}_filt"}
+    }
+    withName: 'BCFTOOLS_MERGE' {
+        cpus = 4
+        ext.args = {"--force-samples -Oz"}
+    }
+    withName: 'BCFTOOLS_CONSENSUS' {
+        ext.prefix = {"${meta.id}_consensus"}
+    }
+    withName: 'DVPOLISH_CREATE_FINALASM' {
+        ext.prefix = {"${meta.id}_dvpolish"}
+    }
+    withName: 'DVPOLISH_CHUNKFA' {
+        // split assembly file into smaller pieces of size: chunk_size
+        ext.args = 'chunk_size = 100.MB'
+    }
+    withName: 'DVPOLISH_CHUNKFA' {
+        // split assembly file into smaller pieces of size: chunk_size
+        ext.args = 'chunk_size = 100.MB'
+    }
+    withName: 'DVPOLISH_.*' {
+        tag        = { meta.assembly.build }
+        publishDir = [
+            path: { "$params.outdir/$stage.polish/dvpolish/$meta.assembly.build" },
+            mode: params.publish_mode,
+            saveAs: { filename -> filename.equals('versions.yml') ? null : filename }
+        ]
+    }
+
     // SCAFFOLD
     withName: 'BWAMEM2_MEM_SCAFFOLD' {
         ext.prefix = { "${meta.id}_${meta.assembly.build}_${reads.head().getBaseName(reads.head().name.endsWith(".gz") ? 2 : 1)}" }
         ext.args   = '-SP -T0'
         ext.args2  = { sort_bam ? "--write-index" : "" }
     }
-
     withName: 'PAIRTOOLS' {
         ext.prefix = { "${meta.id}_${meta.assembly.build}" }
         ext.args   = { "--min-mapq 20 --walks-policy 5unique --max-inter-align-gap 30" }
@@ -326,7 +379,6 @@ process {
             pattern: "*.stat"
         ]
     }
-
     withName: 'YAHS' {
         ext.prefix = { "${meta.id}_${meta.assembly.build}" }
         ext.args   = { params.hic_type == "arima_v2" ? "-e GATC,GANTC,CTNAG,TTAA" : params.hic_type == "arima_v1" ? "-e GATC,GANTC" : "" }
@@ -352,16 +404,13 @@ process {
         ext.args   = { "${params.hic_map_qv}" }  // min hic map quality
         ext.args2  = '--write-index -l1'
     }
-
     withName: 'SCAFFOLD_CURATION:SAMTOOLS_MERGE_(HIFI|HIC)' {
         ext.prefix = { "${meta.assembly.build}_merged" }
         ext.args   = '--write-index'
     }
-
     withName: 'SCAFFOLD_CURATION:BIOBAMBAM_BAMMARKDUPLICATES2' {
         ext.prefix = { "${bam.getBaseName(1)}_dupMarked" }
     }
-
     withName: 'SCAFFOLD_CURATION:BAM2BED_SORT' {
         ext.prefix = { "${bam.head().getBaseName(1)}" }
         ext.args   = '-u -F0x400 -h'            // arguments for samtools view

main.nf

@@ -22,6 +22,8 @@ include { ASSEMBLE_ORGANELLES                               } from "./subworkflo
 include { DECONTAMINATE                                     } from "./subworkflows/local/04_decontaminate/main"
 // Purge duplicates
 include { PURGE_DUPLICATES                                  } from "./subworkflows/local/05_purge_dups/main"
+// Polish
+include { DVPOLISH                                          } from "./subworkflows/local/06_polish/main"
 // Scaffold
 include { SCAFFOLD                                          } from "./subworkflows/local/07_scaffold/main"
 // Curation
@@ -184,7 +186,14 @@ workflow {
     ).dump(tag: 'Assemblies: to polish', pretty: true)
     if ( 'polish' in workflow_steps ) {
         // Run polishers
-        ch_polished_assemblies = ch_to_polish
+        DVPOLISH (
+            ch_to_polish,
+            ch_hifi,
+            BUILD_MERYL_HIFI_DATABASE.out.uniondb
+        )
+        ch_evaluate_assemblies = ch_evaluate_assemblies.mix( DVPOLISH.out.assemblies )
+        ch_polished_assemblies = DVPOLISH.out.assemblies
+        //TODO: logs, versions
     } else {
         ch_polished_assemblies = ch_to_polish
     }

modules/local/dvpolish/chunkfa.nf

@@ -0,0 +1,68 @@
+process DVPOLISH_CHUNKFA {
+    tag "$meta.id"
+    label 'process_single'
+
+    conda "${moduleDir}/environment.yml"
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/ubuntu:20.04' :
+        'nf-core/ubuntu:20.04' }"
+
+
+    input:
+    tuple val(meta), path(fai)
+
+    output:
+    tuple val(meta), path ("*.bed", arity: '1..*')        , emit: bed
+    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}"
+
+    // Split the string by spaces, then iterate to create a map
+    // assumption all arguments are provided like key1 = value1 key2 = value2 ...
+    def args_map = [:]
+    args.replaceAll(' = ', '=').split().each {
+        def (key, value) = it.split('=')
+        args_map[key.trim()] = value.trim()
+    }
+    def chunk_size = MemoryUnit.of(args_map['chunk_size'] ?: '90MB').toBytes()
+
+    """
+    # convert chunk size into base pairs
+    awk -v chunk_size_inBases=${chunk_size} -v prefix=$prefix 'BEGIN {
+        block=1
+        cum_basecount=0
+    }{
+        output_file = sprintf("%s_chunk_%d.bed", prefix, block)
+        printf("%s\\t0\\t%s\\n", \$1, \$2) > output_file
+        cum_basecount+=\$2
+
+        if (cum_basecount >= chunk_size_inBases)
+        {
+            cum_basecount=0
+            block+=1
+        }
+    }' $fai
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        dvpolish: \$(awk --version |& sed '1!d')
+    END_VERSIONS
+    """
+
+    stub:
+    def args = task.ext.args ?: ''
+    def prefix = task.ext.prefix ?: "${meta.id}"
+    """
+    touch ${prefix}_chunk_1.bed
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        dvpolish: \$(awk --version |& sed '1!d')
+    END_VERSIONS
+    """
+}

modules/local/dvpolish/createFinalAsm.nf

@@ -0,0 +1,78 @@
+process DVPOLISH_CREATE_FINALASM {
+    tag "$meta.id"
+    label 'process_single'
+
+    conda 'bioconda::seqkit=2.8.2'
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/seqkit:2.8.2--h9ee0642_0' :
+        'nf-core/seqkit:2.8.2--h9ee0642_0' }"
+
+    input:
+    tuple val(meta), path(unpol_fasta), path(unpol_merqury_csv)  // meta map, unpolished assembly, corresponding merqury qv file
+    tuple val(meta2), path(pol_fasta), path(pol_merqury_csv)     // meta map, polished assembly, corresponding merqury qv file
+
+    output:
+    tuple val(meta), path('*.fasta.gz')  , emit: fasta_gz
+    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 unpol_name = unpol_fasta.getBaseName()
+    def unpol_ext  = unpol_fasta.getExtension()
+    def pol_name   = pol_fasta.getBaseName()
+    def pol_ext    = pol_fasta.getExtension()
+
+    """
+    nl_unpol_ASM=\$(wc -l < ${unpol_merqury_csv})
+    nl_pol_ASM=\$(wc -l < ${pol_merqury_csv})
+
+    if [[ \${nl_unpol_ASM} -ne \${nl_pol_ASM} ]]
+    then 
+        >&2 echo "[ERROR] DVPOLISH_CREATE_FINALASM: merqury files have different lines: ${unpol_merqury_csv}: \${nl_unpol_ASM} != ${pol_merqury_csv}: \${nl_pol_ASM}"
+        exit 1
+    fi
+
+    #split unpolished assembly by sequence ID
+    seqkit split -i -O unpolished_asm ${unpol_fasta}
+
+    #split polished assembly by sequence ID
+    seqkit split -i -O polished_asm ${pol_fasta}
+
+    l=1 
+    while [[ \$l -le \${nl_pol_ASM} ]]
+    do 
+        # read qv string in line p from unpolished and polished merqury files into bash array l_uasm and l_pasm respectively
+        IFS='\t' read -r -a l_uasm <<< "\$(sed -n \${l}p ${unpol_merqury_csv})"
+        IFS='\t' read -r -a l_pasm <<< "\$(sed -n \${l}p ${pol_merqury_csv})"
+
+        # check if the contig names (column 1) are the same
+        if [[ "\${l_uasm[0]}" != "\${l_pasm[0]}" ]]
+        then 
+            >&2 echo "[ERROR] DVPOLISH_CREATE_FINALASM: merqury files are not in the same order!"
+            >&2 echo "[ERROR]        file:  ${unpol_fasta} line \$l: \${l_uasm[*]}"
+            >&2 echo "[ERROR]        file:  ${pol_fasta}   line \$l: \${l_pasm[*]}"
+            exit 2
+        fi 
+
+        # compare number of errorneous kmers (column 2)
+        if [[ \${l_uasm[1]} -le \${l_pasm[1]} ]]    # unpolished assembly has fewer errors, or no difference -> go with the unpolished assembly 
+        then 
+            cat unpolished_asm/${unpol_name}.part_\${l_uasm[0]}.${unpol_ext}
+            >&2 echo "[WARNING] DVPOLISH_CREATE_FINALASM: unpolished contig \${l_uasm[0]} has better or equal QV: \${l_uasm[3]} vs \${l_pasm[3]}"
+        else                                      # polished assembly has fewer errors
+            cat polished_asm/${pol_name}.part_\${l_pasm[0]}.${pol_ext}
+        fi
+
+        l=\$((l+1))
+    done | gzip -c > ${prefix}.fasta.gz
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        dvpolish: \$(seqkit version)
+    END_VERSIONS
+    """
+}

modules/local/dvpolish/environment.yml

@@ -0,0 +1,7 @@
+name: cat_fastq
+channels:
+  - conda-forge
+  - bioconda
+  - defaults
+dependencies:
+  - conda-forge::sed=4.7

modules/local/dvpolish/pbmm2_align.nf

@@ -0,0 +1,42 @@
+process DVPOLISH_PBMM2_ALIGN {
+    tag "$meta.id"
+    label 'process_medium'
+
+    // Note: the versions here need to match the versions used in pbmm2/index
+    conda 'bioconda::pbmm2=1.13.1'
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/pbmm2:1.13.1--h9ee0642_0' :
+        'biocontainers/pbmm2:1.13.1--h9ee0642_0' }"
+
+    input:
+    tuple val(meta), path(reads)
+    tuple val(meta2), path(reference)
+
+    output:
+    tuple val(meta), path("*.bam"), path("*.bai"), emit: bam_bai
+    path "versions.yml"           , emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def args  = task.ext.args ?: ''
+    def args2 = task.ext.args2 ?: ''
+
+    def out_name_part1 = reference.name.endsWith(".gz") ? reference.getBaseName(2) : reference.baseName
+    def out_name_part2 = reads.name.endsWith(".gz") ?  reads.getBaseName(2) : reads.baseName
+
+    """
+    pbmm2 align --sort \\
+        $args \\
+        -j $task.cpus \\
+        "$reference" \\
+        "$reads" \\
+        ${out_name_part1}_${out_name_part2}.bam
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        \$(pbmm2 --version 2>&1 | head -n 1)
+    END_VERSIONS
+    """
+}

modules/local/dvpolish/pbmm2_index.nf

@@ -0,0 +1,34 @@
+process DVPOLISH_PBMM2_INDEX {
+    label 'process_medium'
+
+    // Note: the versions here need to match the versions used in pbmm2/align
+    conda 'bioconda::pbmm2=1.13.1'
+    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
+        'https://depot.galaxyproject.org/singularity/pbmm2:1.13.1--h9ee0642_0' :
+        'biocontainers/pbmm2:1.13.1--h9ee0642_0' }"
+
+    input:
+    tuple val(meta), path(fasta)
+
+    output:
+    tuple val(meta), path("*.mmi"), emit: index
+    path "versions.yml"           , emit: versions
+
+    when:
+    task.ext.when == null || task.ext.when
+
+    script:
+    def args = task.ext.args ?: ''
+    """
+    pbmm2 index\\
+        -j $task.cpus \\
+        $args \\
+        $fasta \\
+        ${fasta.baseName}.mmi
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        \$(pbmm2 --version 2>&1 | head -n 1)
+    END_VERSIONS
+    """
+}