add consensus masking

rmcolq · Sep 22, 2023 · 3e77c2d · 3e77c2d
1 parent af60d54
commit 3e77c2d
Show file tree

Hide file tree

Showing 5 changed files with 82 additions and 8 deletions.
diff --git a/conf/params.config b/conf/params.config
@@ -51,6 +51,7 @@ params {
 
     reference_segment_sep = "NNNNNNNNNNNNNNNNNNNNNNNNNNNNNN"
     reference_min_count = 10
+    mask_depth = 20
 
     assembler = null
 
@@ -73,6 +74,6 @@ params {
       ]
       agent = null
       container = "rmcolq/pantheon"
-      container_sha = "sha256:7b547c61219e9b6c5821cf9a9b50ca184613608f172fc7ae156de9b5a142f19a"
+      container_sha = "sha256:6f08c4c419929666b62e8b0699c588d23c22ebe18803fdcf207f24d6e4c68631"
     }
 }
diff --git a/environment.yml b/environment.yml
@@ -14,4 +14,9 @@ dependencies:
   - tabix
   - mako
   - entrez-direct
+  - minimap2
+  - samtools=1.17
+  - bedtools
+  - pip:
+      - git+https://github.com/Desperate-Dan/maskara.git
 
diff --git a/modules/assemble_taxa.nf b/modules/assemble_taxa.nf
@@ -69,3 +69,74 @@ process medaka_consensus {
         mv medaka_output/consensus.fasta "assembly_${taxon_id}.fa"
         """
 }
+
+process map_to_consensus {
+
+    label 'process_low'
+
+    conda 'bioconda::minimmap2=v2.26 bioconda::samtools=v1.17'
+    container "${params.wf.container}@${params.wf.container_sha}"
+
+    input:
+        tuple val(taxon_id), val(barcode_id), path(reads), path(assembly)
+    output:
+        tuple val(taxon_id), val(barcode_id), path(reads), path(assembly), path("out.bam"), path("out.bam.csi")
+    script:
+        """
+        minimap2 -x map-ont -a ${assembly} ${reads} | samtools sort -o out.bam --write-index -
+        """
+}
+
+process define_mask {
+
+    label 'process_low'
+
+    container "${params.wf.container}@${params.wf.container_sha}"
+
+    input:
+        tuple val(taxon_id), val(barcode_id), path(reads), path(assembly), path(bam), path(index)
+    output:
+        tuple val(taxon_id), val(barcode_id), path(reads), path(assembly), path("mask.bed.tsv")
+    script:
+        """
+        str=\$(head -n1 ${assembly})
+        ref="\${str: 1}"
+        maskara -d ${params.mask_depth} -r \$ref -o mask.bed ${bam}
+        """
+}
+
+process apply_mask {
+
+    label 'process_low'
+
+    publishDir path: "${params.outdir}/${barcode_id}/assemblies", mode: 'copy', pattern: 'assembly_*.fa'
+
+    conda 'bioconda::bedtools=v2.31.0'
+    container "${params.wf.container}@${params.wf.container_sha}"
+
+    input:
+        tuple val(taxon_id), val(barcode_id), path(reads), path(assembly), path(mask)
+    output:
+        tuple val(taxon_id), val(barcode_id), path(reads), path("assembly_${taxon_id}.masked.fa")
+    script:
+        """
+        bedtools maskfasta -fi ${assembly} -bed ${mask} -fo "assembly_${taxon_id}.masked.fa"
+        """
+}
+
+workflow assemble {
+    take:
+        unique_id
+        ch_assembly
+    main:
+        subset_references(ch_assembly)
+        check_subset(subset_references.out.all)
+        medaka_consensus(check_subset.out)
+        map_to_consensus(medaka_consensus.out)
+        define_mask(map_to_consensus.out)
+        apply_mask(define_mask.out)
+        subset_references.out.summary.collectFile(name: "${params.outdir}/${unique_id}/reference_by_barcode.csv", keepHeader: true, skip:1).set{ assembly_summary }
+    emit:
+        completed = apply_mask.out
+        summary = assembly_summary
+    }
diff --git a/nextflow.config b/nextflow.config
@@ -32,7 +32,7 @@ manifest {
     description     = 'Investigate and assemble metagenomic sequence data from human respiratory infections.'
     mainScript      = 'main.nf'
     nextflowVersion = '>=20.10.0'
-    version         = 'v0.5.0'
+    version         = 'v0.5.1'
 }
 
 profiles {

diff --git a/subworkflows/process_run.nf b/subworkflows/process_run.nf
@@ -1,6 +1,6 @@
 // workflow to extract reads and create assemblies from a wf-metagenomic or scylla run
 include { move_or_compress; unpack_taxonomy; extract_reads; check_reads; taxid_from_name; download_references_by_taxid; filter_references } from '../modules/preprocess'
-include { subset_references; check_subset; medaka_consensus } from '../modules/assemble_taxa'
+include { assemble } from '../modules/assemble_taxa'
 include { generate_assembly_report } from '../modules/generate_report'
 
 EXTENSIONS = ["fastq", "fastq.gz", "fq", "fq.gz"]
@@ -69,11 +69,8 @@ workflow process_run {
               .set{ ch_reads }
 
             ch_assembly = ch_reads.combine(ch_references, by: 0)
-            subset_references(ch_assembly)
-            check_subset(subset_references.out.all)
-            medaka_consensus(check_subset.out)
-            subset_references.out.summary.collectFile(name: "${params.outdir}/${unique_id}/reference_by_barcode.csv", keepHeader: true, skip:1).set{ assembly_summary }
+            assemble(unique_id, ch_assembly)
 
-            generate_assembly_report(unique_id, input_summary, reference_summary, extract_summary, assembly_summary)
+            generate_assembly_report(unique_id, input_summary, reference_summary, extract_summary, assemble.out.summary)
         }
 }