From 58532d0b92898f8231c00383a1212a1f919a5675 Mon Sep 17 00:00:00 2001 From: Jason Stajich Date: Mon, 29 Jun 2026 09:41:30 -0700 Subject: [PATCH 1/5] =?UTF-8?q?feat(#6):=20extract=20GENOME=5FCLEAN=20+=20?= =?UTF-8?q?GENOME=5FCLEAN=5FBATCH=20=E2=86=92=20modules/local/=20+=20CLEAN?= =?UTF-8?q?=5FGENOMES=20subworkflow?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two processes moved out of the funannotate.nf monolith: modules/local/genome_clean.nf — per-genome FCS-GX purge + length-filter modules/local/genome_clean_batch.nf — batched variant (amortizes /dev/shm staging) Composed into subworkflows/local/clean_genomes.nf (CLEAN_GENOMES), which encapsulates the batching decision (clean_batch_size / skip_fcs), the jobs_to_clean pre-filter, the sentinel-gate clean_done_ch, and the clean_genome_ch reconstruction. funannotate.nf reduces to: CLEAN_GENOMES(jobs, taxondb_ch) def clean_genome_ch = CLEAN_GENOMES.out.genomes Dual-label pattern: 'genome_clean' (provisioning) + 'process_high'/'process_long' (tier defaults); inline cpus/memory/time removed (overridden by withName in conf/profile_annotate.config). withName patterns in profile_annotate.config updated from 'GENOME_CLEAN' to '.*:GENOME_CLEAN' (regex) to match the subworkflow-namespaced process names. genomeFile() is a top-level def in funannotate.nf, not accessible from included subworkflow files. clean_genomes.nf inlines the equivalent two-liner logic (prefer .fa.gz if non-empty, else .fa) rather than refactoring to lib/ now. Stub-run gate: 15 completed, 0 failed (CLEAN_GENOMES:GENOME_CLEAN_BATCH correctly namespaced; SETUP_AUGUSTUS_CONFIG skips) Co-Authored-By: Claude Sonnet 4.6 Claude-Session: https://claude.ai/code/session_01AbWwmm9v6dJn13sgAJn3Sa --- conf/profile_annotate.config | 4 +- funannotate.nf | 281 +--------------------------- modules/local/genome_clean.nf | 85 +++++++++ modules/local/genome_clean_batch.nf | 125 +++++++++++++ subworkflows/local/clean_genomes.nf | 74 ++++++++ 5 files changed, 289 insertions(+), 280 deletions(-) create mode 100644 modules/local/genome_clean.nf create mode 100644 modules/local/genome_clean_batch.nf create mode 100644 subworkflows/local/clean_genomes.nf diff --git a/conf/profile_annotate.config b/conf/profile_annotate.config index a8ed12b..c6acd23 100644 --- a/conf/profile_annotate.config +++ b/conf/profile_annotate.config @@ -126,7 +126,7 @@ process { queue = 'short'; cpus = 1; memory = '4 GB'; time = '30m' } - withName: 'GENOME_CLEAN' { + withName: '.*:GENOME_CLEAN' { // Full FCS-GX purge needs highmem (stages ~470 GB gx db into /dev/shm); // with --skip_fcs it's a light decompress + length-filter. queue = params.skip_fcs ? 'short' : 'highmem' @@ -135,7 +135,7 @@ process { time = params.skip_fcs ? '1h' : '6h' } - withName: 'GENOME_CLEAN_BATCH' { + withName: '.*:GENOME_CLEAN_BATCH' { // Stages the ~470 GB FCS-GX db into /dev/shm once, then cleans up to // clean_batch_size genomes sequentially — needs highmem and a long walltime. queue = 'highmem'; cpus = 16; memory = '500 GB'; time = '7d' diff --git a/funannotate.nf b/funannotate.nf index 94c7d87..16e65c5 100644 --- a/funannotate.nf +++ b/funannotate.nf @@ -33,229 +33,6 @@ // FUNANNOTATE_PREDICT receives: tuple val(meta), val(genome_fa) // → emits metadata: tuple val(meta) -// Download and extract NCBI taxdump once; storeDir caches it at params.taxondb so -// subsequent runs skip this entirely. -process GENOME_CLEAN { - label 'genome_clean' - tag "${meta.id}" - - // container '/rhome/jstajich/projects/AAFTF/AAFTF_v0.6.1-signed.sif' - - // Nextflow skips this task when input_clean_genomes/.fa.gz already exists. - storeDir "${launchDir}/input_clean_genomes" - - cpus 16 - memory '450 GB' - time '6h' - - input: - tuple val(meta), path(genome_gz), val(taxondb) - - output: - tuple val(meta), path("${meta.asmid}.fa.gz"), emit: genome - - script: - def out = meta.id - def asmid = meta.asmid - def taxonid = meta.taxonid - """ - if [ ! -f "${genome_gz}" ]; then - echo "ERROR: genome_gz not found at path: ${genome_gz}" >&2 - exit 1 - fi - - source /etc/profile.d/modules.sh 2>/dev/null || true - SCRATCH=\$(printf '%s' "\${SCRATCH:-.}" | tr -d '\\n\\r') - - echo "[INFO] Decompressing genome for ${asmid}..." - # Accept either a gzipped (NCBI_ASM .fna.gz) or plain (local GENOME column) FASTA. - if printf '%s' "${genome_gz}" | grep -qiE '\\.gz\$'; then - pigz -dc ${genome_gz} > \$SCRATCH/${asmid}.raw.fa - else - cat ${genome_gz} > \$SCRATCH/${asmid}.raw.fa - fi - - if [ "${params.skip_fcs}" = "true" ]; then - # --skip_fcs: bypass AAFTF FCS-GX contaminant purge (no 470 GB gxdb needed); - # just length-filter the assembly. - echo "[INFO] --skip_fcs set: skipping FCS-GX purge for ${asmid}" - ${params.clean_script} --len ${params.min_contig_len} \ - -i \$SCRATCH/${asmid}.raw.fa -o ${asmid}.fa - else - # Ensure /dev/shm/gxdb is present on this node; register for cleanup when done. - source ${params.fcs_shm_script} - TAXONKIT_DB=${taxondb} - phylum=\$(echo ${taxonid} | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{p}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | cut -f2 | uniq | head -n 1) - if [ -z "\$phylum" ]; then - phylum=\$(echo ${taxonid} | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{K}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | uniq | cut -f2 | head -n 1) - # weird we are getting 2 lines from name2taxid when input is Fungi add the uniq/head -n 1 to ensure only one line - fi - echo "[INFO] Phylum for ${asmid} (taxonid=${taxonid}): \$phylum" - echo "[INFO] FCS-GX purge + cleaning genome for ${asmid}..." - AAFTF fcs_gx_purge --db /dev/shm/gxdb/all \ - -i \$SCRATCH/${asmid}.raw.fa --cpus ${task.cpus} \ - -o \$SCRATCH/${asmid}.purge.fasta \ - -t "\$phylum" -w \$SCRATCH/${asmid}.fcs_report - mkdir -p ${launchDir}/input_clean_genomes/clean - cat \$SCRATCH/${asmid}.purge.fasta | \ - ${params.clean_script} --len ${params.min_contig_len} > ${asmid}.fa - pigz \$SCRATCH/${asmid}.purge.fasta - [ -f \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv ] && pigz \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv - mv \$SCRATCH/${asmid}.purge.fasta.gz ${launchDir}/input_clean_genomes/clean/ - [ -f \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv.gz ] && \ - mv \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv.gz ${launchDir}/input_clean_genomes/clean/ - fi - # Deliver the clean genome gzip-compressed to save space in input_clean_genomes; - # downstream tools inflate it on the fly (they cannot read a gzipped FASTA via -i). - pigz -f ${asmid}.fa - echo "[INFO] Clean genome written: ${asmid}.fa.gz (\$(du -sh ${asmid}.fa.gz | cut -f1))" - rm -f \$SCRATCH/${asmid}.raw.fa - """ - - stub: - def asmid = meta.asmid - """ - echo ">stub_${asmid}" | pigz -c > ${asmid}.fa.gz - mkdir -p ${launchDir}/input_clean_genomes/clean - touch ${launchDir}/input_clean_genomes/clean/${asmid}.purge.fasta - touch ${launchDir}/input_clean_genomes/clean/${asmid}.purge.fcs_gx-taxonomy.tsv - """ -} - -// Batched variant of GENOME_CLEAN. Receives a LIST of per-genome tuples and stages -// the FCS-GX database into /dev/shm ONCE (~30 min) via the shared label provisioning, -// then cleans every genome in the batch sequentially against that in-memory DB. This -// amortizes the expensive staging step over ~clean_batch_size genomes instead of -// paying it per genome. (Only used when params.skip_fcs is false — with FCS skipped -// there is no DB to stage and per-genome GENOME_CLEAN is used instead.) -// -// Outputs are written directly to ${launchDir}/input_clean_genomes/.fa (the same -// location GENOME_CLEAN's storeDir uses) plus a per-batch manifest listing every cleaned -// assembly. Genomes whose .fa already exists are skipped, so a killed/retried batch -// resumes without redoing finished assemblies. -// -// Uses the 'genome_clean' label, whose beforeScript loads miniconda3/conda + taxonkit + -// AAFTF, so the script stays tool-agnostic (no inline module loads), matching GENOME_CLEAN. -process GENOME_CLEAN_BATCH { - label 'genome_clean' - tag "clean_batch_${task.index}" - - cpus 16 - memory '450 GB' - time '7d' - - input: - tuple val(items), val(taxondb) - - output: - path "clean_batch_*.manifest.tsv", emit: manifest - - script: - def batch_tsv = items.collect { row -> "${row[0].asmid}\t${row[1]}\t${row[0].taxonid}" }.join('\n') - """ - set -uo pipefail - source /etc/profile.d/modules.sh 2>/dev/null || true - - SCRATCH=\$(printf '%s' "\${SCRATCH:-.}" | tr -d '\\n\\r') - TAXONKIT_DB=${taxondb} - DEST=${launchDir}/input_clean_genomes - mkdir -p \$DEST/clean - - MANIFEST=clean_batch_${task.index}.manifest.tsv - : > \$MANIFEST - - cat > batch.tsv <<'BATCH_EOF' -${batch_tsv} -BATCH_EOF - - n_total=\$(grep -c . batch.tsv || true) - echo "[INFO] batch ${task.index}: \$n_total genomes to consider" - - # Stage the FCS-GX DB into /dev/shm ONCE for the whole batch (~30 min). FCS_GX_KEEP_SHM=1 - # tells the staging script not to register its own per-shell EXIT cleanup; we remove the - # RAM copy ourselves when the batch finishes (or aborts) via the trap below. - export FCS_GX_KEEP_SHM=1 - source ${params.fcs_shm_script} - trap 'rm -rf /dev/shm/gxdb 2>/dev/null || true' EXIT - if [ ! -f /dev/shm/gxdb/all.gxi ]; then - echo "[ERROR] FCS-GX DB not staged into /dev/shm/gxdb; aborting batch" >&2 - exit 1 - fi - - i=0 - while IFS=\$'\\t' read -r asmid gz taxonid; do - [ -z "\$asmid" ] && continue - i=\$((i+1)) - target=\$DEST/\${asmid}.fa.gz - if [ -s "\$target" ]; then - echo "[\$i/\$n_total][SKIP] \$asmid already cleaned" - printf '%s\\t%s\\n' "\$asmid" "\$target" >> \$MANIFEST - continue - elif [ -s "\$DEST/\${asmid}.fa" ]; then - # Back-compat: a prior run may have left an uncompressed .fa. - echo "[\$i/\$n_total][SKIP] \$asmid already cleaned (uncompressed)" - printf '%s\\t%s\\n' "\$asmid" "\$DEST/\${asmid}.fa" >> \$MANIFEST - continue - fi - if [ ! -f "\$gz" ]; then - echo "[\$i/\$n_total][WARN] missing genome for \$asmid: \$gz" >&2 - continue - fi - - phylum=\$(echo \$taxonid | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{p}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | cut -f2 | uniq | head -n 1) - if [ -z "\$phylum" ]; then - phylum=\$(echo \$taxonid | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{K}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | uniq | cut -f2 | head -n 1) - fi - echo "[\$i/\$n_total][INFO] \$asmid taxonid=\$taxonid phylum=\$phylum" - - # Accept gzipped (NCBI_ASM .fna.gz) or plain (local GENOME column) FASTA input. - if printf '%s' "\$gz" | grep -qiE '\\.gz\$'; then - pigz -dc "\$gz" > \$SCRATCH/\${asmid}.raw.fa - else - cat "\$gz" > \$SCRATCH/\${asmid}.raw.fa - fi - if AAFTF fcs_gx_purge --db /dev/shm/gxdb/all \\ - -i \$SCRATCH/\${asmid}.raw.fa --cpus ${task.cpus} \\ - -o \$SCRATCH/\${asmid}.purge.fasta \\ - -t "\$phylum" -w \$SCRATCH/\${asmid}.fcs_report ; then - cat \$SCRATCH/\${asmid}.purge.fasta | ${params.clean_script} --len ${params.min_contig_len} > \$SCRATCH/\${asmid}.clean.fa \\ - && pigz -c \$SCRATCH/\${asmid}.clean.fa > \${target}.tmp \\ - && mv \${target}.tmp \$target - rm -f \$SCRATCH/\${asmid}.clean.fa - echo "[\$i/\$n_total][OK] \$asmid -> \$target (\$(du -sh \$target | cut -f1))" - pigz -f \$SCRATCH/\${asmid}.purge.fasta - [ -f \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv ] && pigz -f \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv - mv \$SCRATCH/\${asmid}.purge.fasta.gz \$DEST/clean/ 2>/dev/null || true - [ -f \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv.gz ] && mv \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv.gz \$DEST/clean/ - printf '%s\\t%s\\n' "\$asmid" "\$target" >> \$MANIFEST - else - echo "[\$i/\$n_total][FAIL] fcs_gx_purge failed for \$asmid" >&2 - fi - rm -f \$SCRATCH/\${asmid}.raw.fa \$SCRATCH/\${asmid}.purge.fasta - done < batch.tsv - - echo "[INFO] batch ${task.index} complete: \$(grep -c . \$MANIFEST || echo 0) cleaned genomes in manifest" - """ - - stub: - def batch_tsv = items.collect { row -> "${row[0].asmid}\t${row[1]}\t${row[0].taxonid}" }.join('\n') - """ - DEST=${launchDir}/input_clean_genomes - mkdir -p \$DEST/clean - MANIFEST=clean_batch_${task.index}.manifest.tsv - : > \$MANIFEST - cat > batch.tsv <<'BATCH_EOF' -${batch_tsv} -BATCH_EOF - while IFS=\$'\\t' read -r asmid gz taxonid; do - [ -z "\$asmid" ] && continue - echo ">stub_\${asmid}" | pigz -c > \$DEST/\${asmid}.fa.gz - touch \$DEST/clean/\${asmid}.purge.fasta \$DEST/clean/\${asmid}.purge.fcs_gx-taxonomy.tsv - printf '%s\\t%s\\n' "\$asmid" "\$DEST/\${asmid}.fa.gz" >> \$MANIFEST - done < batch.tsv - """ -} - // Soft-mask each assembly using funannotate mask with tantan. // storeDir caches the masked FASTA alongside the clean genome. process MASKREPEAT_TANTAN_RUN { @@ -1638,6 +1415,7 @@ include { validateParameters; paramsSummaryLog; paramsHelp } from 'plugin/nf-sch include { ASM_STATS } from './modules/local/asm_stats' include { INPUT_CHECK } from './subworkflows/local/input_check' include { SETUP_DBS } from './subworkflows/local/setup_dbs' +include { CLEAN_GENOMES } from './subworkflows/local/clean_genomes' include { ANTISMASH_RUN } from './modules/local/antismash_run' include { INTERPROSCAN_RUN } from './modules/local/interproscan_run' include { SIGNALP_RUN } from './modules/local/signalp_run' @@ -1677,63 +1455,10 @@ workflow { SETUP_DBS() def taxondb_ch = SETUP_DBS.out.taxondb - // Only clean genomes whose cleaned .fa does not already exist. This keeps batches - // from being padded with finished genomes — a batch that is entirely cleaned is never - // scheduled, so it never pays the ~30-min /dev/shm staging cost. (GENOME_CLEAN_BATCH - // also re-checks per genome at runtime, which handles partial completion on retry.) - def jobs_to_clean = jobs.filter { meta, gz -> - !genomeFile("${launchDir}/input_clean_genomes/${meta.asmid}.fa").exists() - } - - // Genome cleaning. The FCS-GX DB staging into /dev/shm costs ~30 min per task, so by - // default we batch genomes (clean_batch_size, default 1000) into single SLURM jobs that - // stage the DB once and then clean every genome in the batch. Set clean_batch_size = 0 - // (or --skip_fcs, where there is no DB to amortize) to fall back to one SLURM job per - // genome via GENOME_CLEAN. clean_done_ch gates downstream on cleaning finishing; - // ifEmpty([]) ensures it still emits (so downstream runs) when every genome was already - // clean and nothing was scheduled. - def clean_done_ch - int clean_batch_size = params.clean_batch_size as int - if (clean_batch_size > 0 && !params.skip_fcs.toBoolean()) { - // Wrap each collated batch (a List of per-genome tuples) in a single-element list - // so .combine() appends taxondb as the 2nd tuple element instead of spreading the - // batch's rows into the tuple (which would break GENOME_CLEAN_BATCH's - // `tuple val(items), val(taxondb)` declaration). - def clean_batches = jobs_to_clean.collate(clean_batch_size).map { batch -> [ batch ] } - GENOME_CLEAN_BATCH(clean_batches.combine(taxondb_ch)) - clean_done_ch = GENOME_CLEAN_BATCH.out.manifest.collect().ifEmpty([]) - } else { - GENOME_CLEAN(jobs_to_clean.combine(taxondb_ch)) - clean_done_ch = GENOME_CLEAN.out.genome.map { meta, gz -> gz }.collect().ifEmpty([]) - } + CLEAN_GENOMES(jobs, taxondb_ch) if (!params.only_clean.toBoolean()) { - // Re-attach the cleaned genome to its full per-sample metadata. The cleaned genome - // lands at input_clean_genomes/.fa.gz (or .fa for legacy runs); genomeFile() - // resolves whichever exists. We rebuild from the jobs channel and gate on - // clean_done_ch (combine waits until all cleaning is done). genome_fa is emitted as - // an absolute-path string so downstream val(genome_fa) processes reference the file - // directly without Nextflow re-staging it. - // Note: combine on jobs (a [meta,gz] tuple) appends clean_done_ch as a third element; - // we drop gz and the sentinel together. Keeping the tuple avoids emitting a bare Map - // (which Nextflow's combine may flatten unexpectedly). - def clean_genome_ch = jobs - .combine(clean_done_ch) // gate: blocks until cleaning done - // Use it[0] (meta) rather than destructuring: combine with ifEmpty([]) produces - // a variable-length tuple — 2 elements when the sentinel is empty, 3+ when it - // carries collected paths. A fixed-arity { meta, _gz, _sentinel -> } fails on - // the 2-element case, so we index into the element directly. - .map { tuple(it[0], genomeFile("${launchDir}/input_clean_genomes/${it[0].asmid}.fa")) } - .filter { meta, g -> - if (!g.exists()) { - log.warn "No cleaned genome for ${meta.id} (asmid=${meta.asmid}) — skipping downstream" - return false - } - return true - } - .map { meta, genome_fa -> - tuple(meta, genome_fa.toAbsolutePath().toString()) - } + def clean_genome_ch = CLEAN_GENOMES.out.genomes // ── Generate assembly statistics (for earlgrey_mask.nf SELECT_REPS) ──────── // Generate asm_stats.tsv if --gen_asm_stats is true and the file doesn't exist. diff --git a/modules/local/genome_clean.nf b/modules/local/genome_clean.nf new file mode 100644 index 0000000..3f74e7d --- /dev/null +++ b/modules/local/genome_clean.nf @@ -0,0 +1,85 @@ +// FCS-GX contaminant purge + length-filter for a single genome assembly. +// Uses storeDir so the task is skipped if the cleaned .fa.gz already exists. +// Runs under the 'genome_clean' provisioning label (AAFTF + taxonkit + FCS-GX); +// resources are overridden per-profile by withName: '.*:GENOME_CLEAN' in +// conf/profile_annotate.config (500 GB highmem with FCS, 8 GB short without). +process GENOME_CLEAN { + label 'genome_clean' + label 'process_high' + tag "${meta.id}" + + storeDir "${launchDir}/input_clean_genomes" + + input: + tuple val(meta), path(genome_gz), val(taxondb) + + output: + tuple val(meta), path("${meta.asmid}.fa.gz"), emit: genome + + script: + def out = meta.id + def asmid = meta.asmid + def taxonid = meta.taxonid + """ + if [ ! -f "${genome_gz}" ]; then + echo "ERROR: genome_gz not found at path: ${genome_gz}" >&2 + exit 1 + fi + + source /etc/profile.d/modules.sh 2>/dev/null || true + SCRATCH=\$(printf '%s' "\${SCRATCH:-.}" | tr -d '\\n\\r') + + echo "[INFO] Decompressing genome for ${asmid}..." + # Accept either a gzipped (NCBI_ASM .fna.gz) or plain (local GENOME column) FASTA. + if printf '%s' "${genome_gz}" | grep -qiE '\\.gz\$'; then + pigz -dc ${genome_gz} > \$SCRATCH/${asmid}.raw.fa + else + cat ${genome_gz} > \$SCRATCH/${asmid}.raw.fa + fi + + if [ "${params.skip_fcs}" = "true" ]; then + # --skip_fcs: bypass AAFTF FCS-GX contaminant purge (no 470 GB gxdb needed); + # just length-filter the assembly. + echo "[INFO] --skip_fcs set: skipping FCS-GX purge for ${asmid}" + ${params.clean_script} --len ${params.min_contig_len} \ + -i \$SCRATCH/${asmid}.raw.fa -o ${asmid}.fa + else + # Ensure /dev/shm/gxdb is present on this node; register for cleanup when done. + source ${params.fcs_shm_script} + TAXONKIT_DB=${taxondb} + phylum=\$(echo ${taxonid} | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{p}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | cut -f2 | uniq | head -n 1) + if [ -z "\$phylum" ]; then + phylum=\$(echo ${taxonid} | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{K}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | uniq | cut -f2 | head -n 1) + # weird we are getting 2 lines from name2taxid when input is Fungi add the uniq/head -n 1 to ensure only one line + fi + echo "[INFO] Phylum for ${asmid} (taxonid=${taxonid}): \$phylum" + echo "[INFO] FCS-GX purge + cleaning genome for ${asmid}..." + AAFTF fcs_gx_purge --db /dev/shm/gxdb/all \ + -i \$SCRATCH/${asmid}.raw.fa --cpus ${task.cpus} \ + -o \$SCRATCH/${asmid}.purge.fasta \ + -t "\$phylum" -w \$SCRATCH/${asmid}.fcs_report + mkdir -p ${launchDir}/input_clean_genomes/clean + cat \$SCRATCH/${asmid}.purge.fasta | \ + ${params.clean_script} --len ${params.min_contig_len} > ${asmid}.fa + pigz \$SCRATCH/${asmid}.purge.fasta + [ -f \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv ] && pigz \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv + mv \$SCRATCH/${asmid}.purge.fasta.gz ${launchDir}/input_clean_genomes/clean/ + [ -f \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv.gz ] && \ + mv \$SCRATCH/${asmid}.purge.fcs_gx-taxonomy.tsv.gz ${launchDir}/input_clean_genomes/clean/ + fi + # Deliver the clean genome gzip-compressed to save space in input_clean_genomes; + # downstream tools inflate it on the fly (they cannot read a gzipped FASTA via -i). + pigz -f ${asmid}.fa + echo "[INFO] Clean genome written: ${asmid}.fa.gz (\$(du -sh ${asmid}.fa.gz | cut -f1))" + rm -f \$SCRATCH/${asmid}.raw.fa + """ + + stub: + def asmid = meta.asmid + """ + echo ">stub_${asmid}" | pigz -c > ${asmid}.fa.gz + mkdir -p ${launchDir}/input_clean_genomes/clean + touch ${launchDir}/input_clean_genomes/clean/${asmid}.purge.fasta + touch ${launchDir}/input_clean_genomes/clean/${asmid}.purge.fcs_gx-taxonomy.tsv + """ +} diff --git a/modules/local/genome_clean_batch.nf b/modules/local/genome_clean_batch.nf new file mode 100644 index 0000000..054a166 --- /dev/null +++ b/modules/local/genome_clean_batch.nf @@ -0,0 +1,125 @@ +// Batched FCS-GX purge for a list of genome assemblies. +// Stages the ~470 GB FCS-GX database into /dev/shm ONCE (~30 min), then +// cleans every genome in the batch sequentially against that in-memory DB, +// amortizing the staging cost over up to clean_batch_size genomes. +// Only used when params.skip_fcs is false; GENOME_CLEAN handles the per-genome +// fallback (clean_batch_size=0 or --skip_fcs). +// +// Resources overridden per-profile by withName: '.*:GENOME_CLEAN_BATCH' in +// conf/profile_annotate.config (500 GB highmem, 7d walltime). +process GENOME_CLEAN_BATCH { + label 'genome_clean' + label 'process_long' + tag "clean_batch_${task.index}" + + input: + tuple val(items), val(taxondb) + + output: + path "clean_batch_*.manifest.tsv", emit: manifest + + script: + def batch_tsv = items.collect { row -> "${row[0].asmid}\t${row[1]}\t${row[0].taxonid}" }.join('\n') + """ + set -uo pipefail + source /etc/profile.d/modules.sh 2>/dev/null || true + + SCRATCH=\$(printf '%s' "\${SCRATCH:-.}" | tr -d '\\n\\r') + TAXONKIT_DB=${taxondb} + DEST=${launchDir}/input_clean_genomes + mkdir -p \$DEST/clean + + MANIFEST=clean_batch_${task.index}.manifest.tsv + : > \$MANIFEST + + cat > batch.tsv <<'BATCH_EOF' +${batch_tsv} +BATCH_EOF + + n_total=\$(grep -c . batch.tsv || true) + echo "[INFO] batch ${task.index}: \$n_total genomes to consider" + + # Stage the FCS-GX DB into /dev/shm ONCE for the whole batch (~30 min). FCS_GX_KEEP_SHM=1 + # tells the staging script not to register its own per-shell EXIT cleanup; we remove the + # RAM copy ourselves when the batch finishes (or aborts) via the trap below. + export FCS_GX_KEEP_SHM=1 + source ${params.fcs_shm_script} + trap 'rm -rf /dev/shm/gxdb 2>/dev/null || true' EXIT + if [ ! -f /dev/shm/gxdb/all.gxi ]; then + echo "[ERROR] FCS-GX DB not staged into /dev/shm/gxdb; aborting batch" >&2 + exit 1 + fi + + i=0 + while IFS=\$'\\t' read -r asmid gz taxonid; do + [ -z "\$asmid" ] && continue + i=\$((i+1)) + target=\$DEST/\${asmid}.fa.gz + if [ -s "\$target" ]; then + echo "[\$i/\$n_total][SKIP] \$asmid already cleaned" + printf '%s\\t%s\\n' "\$asmid" "\$target" >> \$MANIFEST + continue + elif [ -s "\$DEST/\${asmid}.fa" ]; then + # Back-compat: a prior run may have left an uncompressed .fa. + echo "[\$i/\$n_total][SKIP] \$asmid already cleaned (uncompressed)" + printf '%s\\t%s\\n' "\$asmid" "\$DEST/\${asmid}.fa" >> \$MANIFEST + continue + fi + if [ ! -f "\$gz" ]; then + echo "[\$i/\$n_total][WARN] missing genome for \$asmid: \$gz" >&2 + continue + fi + + phylum=\$(echo \$taxonid | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{p}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | cut -f2 | uniq | head -n 1) + if [ -z "\$phylum" ]; then + phylum=\$(echo \$taxonid | taxonkit --data-dir \$TAXONKIT_DB lineage | taxonkit --data-dir \$TAXONKIT_DB reformat -f "{K}" --output-ambiguous-result | cut -f3 | taxonkit --data-dir \$TAXONKIT_DB name2taxid | uniq | cut -f2 | head -n 1) + fi + echo "[\$i/\$n_total][INFO] \$asmid taxonid=\$taxonid phylum=\$phylum" + + # Accept gzipped (NCBI_ASM .fna.gz) or plain (local GENOME column) FASTA input. + if printf '%s' "\$gz" | grep -qiE '\\.gz\$'; then + pigz -dc "\$gz" > \$SCRATCH/\${asmid}.raw.fa + else + cat "\$gz" > \$SCRATCH/\${asmid}.raw.fa + fi + if AAFTF fcs_gx_purge --db /dev/shm/gxdb/all \\ + -i \$SCRATCH/\${asmid}.raw.fa --cpus ${task.cpus} \\ + -o \$SCRATCH/\${asmid}.purge.fasta \\ + -t "\$phylum" -w \$SCRATCH/\${asmid}.fcs_report ; then + cat \$SCRATCH/\${asmid}.purge.fasta | ${params.clean_script} --len ${params.min_contig_len} > \$SCRATCH/\${asmid}.clean.fa \\ + && pigz -c \$SCRATCH/\${asmid}.clean.fa > \${target}.tmp \\ + && mv \${target}.tmp \$target + rm -f \$SCRATCH/\${asmid}.clean.fa + echo "[\$i/\$n_total][OK] \$asmid -> \$target (\$(du -sh \$target | cut -f1))" + pigz -f \$SCRATCH/\${asmid}.purge.fasta + [ -f \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv ] && pigz -f \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv + mv \$SCRATCH/\${asmid}.purge.fasta.gz \$DEST/clean/ 2>/dev/null || true + [ -f \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv.gz ] && mv \$SCRATCH/\${asmid}.purge.fcs_gx-taxonomy.tsv.gz \$DEST/clean/ + printf '%s\\t%s\\n' "\$asmid" "\$target" >> \$MANIFEST + else + echo "[\$i/\$n_total][FAIL] fcs_gx_purge failed for \$asmid" >&2 + fi + rm -f \$SCRATCH/\${asmid}.raw.fa \$SCRATCH/\${asmid}.purge.fasta + done < batch.tsv + + echo "[INFO] batch ${task.index} complete: \$(grep -c . \$MANIFEST || echo 0) cleaned genomes in manifest" + """ + + stub: + def batch_tsv = items.collect { row -> "${row[0].asmid}\t${row[1]}\t${row[0].taxonid}" }.join('\n') + """ + DEST=${launchDir}/input_clean_genomes + mkdir -p \$DEST/clean + MANIFEST=clean_batch_${task.index}.manifest.tsv + : > \$MANIFEST + cat > batch.tsv <<'BATCH_EOF' +${batch_tsv} +BATCH_EOF + while IFS=\$'\\t' read -r asmid gz taxonid; do + [ -z "\$asmid" ] && continue + echo ">stub_\${asmid}" | pigz -c > \$DEST/\${asmid}.fa.gz + touch \$DEST/clean/\${asmid}.purge.fasta \$DEST/clean/\${asmid}.purge.fcs_gx-taxonomy.tsv + printf '%s\\t%s\\n' "\$asmid" "\$DEST/\${asmid}.fa.gz" >> \$MANIFEST + done < batch.tsv + """ +} diff --git a/subworkflows/local/clean_genomes.nf b/subworkflows/local/clean_genomes.nf new file mode 100644 index 0000000..6dd9398 --- /dev/null +++ b/subworkflows/local/clean_genomes.nf @@ -0,0 +1,74 @@ +/* + * CLEAN_GENOMES — FCS-GX contaminant purge + length-filter + * + * Accepts the full jobs channel and the taxondb gate channel from SETUP_DBS. + * Handles the batching decision internally (params.clean_batch_size > 0 and + * !params.skip_fcs → GENOME_CLEAN_BATCH; otherwise per-genome GENOME_CLEAN). + * + * Emits: + * genomes — channel: tuple(val(meta), val(genome_fa)) + * genome_fa is an absolute-path String to the cleaned .fa/.fa.gz; + * samples with no cleaned file are silently dropped with a warning. + */ + +include { GENOME_CLEAN } from './../../modules/local/genome_clean' +include { GENOME_CLEAN_BATCH } from './../../modules/local/genome_clean_batch' + +workflow CLEAN_GENOMES { + + take: + jobs // channel: tuple(val(meta), path(gz)) + taxondb // channel: val(String taxondb_path) + + main: + // Only submit cleaning jobs for genomes not already in input_clean_genomes/. + // Mirrors genomeFile() from funannotate.nf: prefer .fa.gz if non-empty, else .fa. + // GENOME_CLEAN_BATCH also re-checks per genome at runtime for partial retry. + def jobs_to_clean = jobs.filter { meta, gz -> + def fgz = file("${launchDir}/input_clean_genomes/${meta.asmid}.fa.gz") + def fa = file("${launchDir}/input_clean_genomes/${meta.asmid}.fa") + !(fgz.exists() && fgz.size() > 0) && !fa.exists() + } + + def clean_done_ch + int clean_batch_size = params.clean_batch_size as int + if (clean_batch_size > 0 && !params.skip_fcs.toBoolean()) { + // Wrap each collated batch in a single-element list so .combine() appends + // taxondb as the 2nd tuple element rather than spreading the batch rows. + def clean_batches = jobs_to_clean.collate(clean_batch_size).map { batch -> [ batch ] } + GENOME_CLEAN_BATCH(clean_batches.combine(taxondb)) + clean_done_ch = GENOME_CLEAN_BATCH.out.manifest.collect().ifEmpty([]) + } else { + GENOME_CLEAN(jobs_to_clean.combine(taxondb)) + clean_done_ch = GENOME_CLEAN.out.genome.map { meta, gz -> gz }.collect().ifEmpty([]) + } + + // Re-attach the cleaned genome path to each sample's meta. Gate on clean_done_ch + // so this channel only opens once all cleaning tasks have finished. + // Use it[0] index access (not fixed-arity destructuring) because combine with + // ifEmpty([]) produces a variable-length tuple (2 elements when the sentinel is + // empty, 3+ when it carries collected paths). + // Mirrors genomeFile() from funannotate.nf: prefer .fa.gz if non-empty, else .fa. + def ch_genomes = jobs + .combine(clean_done_ch) + .map { row -> + def meta = row[0] + def fgz = file("${launchDir}/input_clean_genomes/${meta.asmid}.fa.gz") + def fa = file("${launchDir}/input_clean_genomes/${meta.asmid}.fa") + def g = (fgz.exists() && fgz.size() > 0) ? fgz : fa + tuple(meta, g) + } + .filter { meta, g -> + if (!g.exists()) { + log.warn "No cleaned genome for ${meta.id} (asmid=${meta.asmid}) — skipping downstream" + return false + } + return true + } + .map { meta, genome_fa -> + tuple(meta, genome_fa.toAbsolutePath().toString()) + } + + emit: + genomes = ch_genomes +} From 9e4b3d1439febe60e5e35cd5292137074aa1ab5e Mon Sep 17 00:00:00 2001 From: Jason Stajich Date: Mon, 29 Jun 2026 09:52:11 -0700 Subject: [PATCH 2/5] =?UTF-8?q?feat(#7):=20extract=20MASKREPEAT=5FTANTAN?= =?UTF-8?q?=5FRUN=20=E2=86=92=20modules/local/=20+=20MASK=5FGENOME=20subwo?= =?UTF-8?q?rkflow?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Process moved to modules/local/maskrepeat_tantan_run.nf; labels updated to dual pattern: 'funannotate' (provisioning) + 'process_medium' (tier default, overridden by withName '.*:MASKREPEAT_TANTAN_RUN' in profile_annotate.config). MASK_GENOME subworkflow (subworkflows/local/mask_genome.nf) encapsulates the run_repeatmasker if/else and the fallback-to-prior-masked-or-clean logic. funannotate.nf reduces to: MASK_GENOME(clean_genome_ch) def predict_genome_ch = MASK_GENOME.out.genomes genomeFile() unavailability in subworkflows: fallback path inlines equivalent logic (prefer .masked.fasta.gz if non-empty, else .masked.fasta). Stub-run gate: 15 completed, 0 failed Co-Authored-By: Claude Sonnet 4.6 Claude-Session: https://claude.ai/code/session_01AbWwmm9v6dJn13sgAJn3Sa --- conf/profile_annotate.config | 2 +- funannotate.nf | 70 ++------------------------ modules/local/maskrepeat_tantan_run.nf | 41 +++++++++++++++ subworkflows/local/mask_genome.nf | 47 +++++++++++++++++ 4 files changed, 94 insertions(+), 66 deletions(-) create mode 100644 modules/local/maskrepeat_tantan_run.nf create mode 100644 subworkflows/local/mask_genome.nf diff --git a/conf/profile_annotate.config b/conf/profile_annotate.config index c6acd23..e130bbd 100644 --- a/conf/profile_annotate.config +++ b/conf/profile_annotate.config @@ -141,7 +141,7 @@ process { queue = 'highmem'; cpus = 16; memory = '500 GB'; time = '7d' } - withName: 'MASKREPEAT_TANTAN_RUN' { + withName: '.*:MASKREPEAT_TANTAN_RUN' { queue = 'short'; cpus = 8; memory = '16 GB'; time = '2h' } diff --git a/funannotate.nf b/funannotate.nf index 16e65c5..040315d 100644 --- a/funannotate.nf +++ b/funannotate.nf @@ -22,8 +22,6 @@ // // GENOME_CLEAN receives: tuple val(meta), path(genome_gz), val(taxondb) // → emits: tuple val(meta), path(genome_fa) [storeDir writes input_clean_genomes/.fa.gz] -// MASKREPEAT_TANTAN_RUN receives: tuple val(meta), val(genome_fa) -// → emits: tuple val(meta), path(masked_fa) [storeDir caches input_clean_genomes/.masked.fasta.gz] // SRA_FETCH receives: val(species_tag), val(taxonid) [only when --run_sra_fetch] // → emits: val(species_tag), path(norm_R1.fastq.gz), path(norm_R2.fastq.gz), path(se) // RNASEQ_PREPARE receives: tuple val(species_tag), val(meta), val(genome_fa), path(r1), path(r2), path(se) @@ -33,47 +31,6 @@ // FUNANNOTATE_PREDICT receives: tuple val(meta), val(genome_fa) // → emits metadata: tuple val(meta) -// Soft-mask each assembly using funannotate mask with tantan. -// storeDir caches the masked FASTA alongside the clean genome. -process MASKREPEAT_TANTAN_RUN { - label 'funannotate' - tag "${meta.id}" - - storeDir "${launchDir}/input_clean_genomes" - - cpus 8 - memory '16 GB' - time '2h' - - input: - tuple val(meta), val(genome_fa) - - output: - tuple val(meta), path("${meta.asmid}.masked.fasta.gz"), emit: masked - - script: - def asmid = meta.asmid - """ - source /etc/profile.d/modules.sh 2>/dev/null || true - # Inflate a gzipped clean genome to a local uncompressed copy; funannotate cannot - # read a gzipped FASTA via -i. Plain (uncompressed) genomes pass through unchanged. - GENOME_FA="${genome_fa}" - case "\$GENOME_FA" in - *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; - *) GENOME_IN="\$GENOME_FA" ;; - esac - funannotate mask -i "\$GENOME_IN" -o ${asmid}.masked.fasta -m tantan --cpus ${task.cpus} - # Deliver the soft-masked genome gzip-compressed to save space; consumers inflate it. - pigz -f ${asmid}.masked.fasta - """ - - stub: - def asmid = meta.asmid - """ - echo ">stub_${asmid}_masked" | pigz -c > ${asmid}.masked.fasta.gz - """ -} - // Query NCBI SRA for available paired-end RNA-seq accessions per species. // Lightweight: runs the esearch/efetch query only — no downloading. // Records up to 5 candidates (sorted by spot count desc) in a per-species CSV. @@ -1416,6 +1373,7 @@ include { ASM_STATS } from './modules/local/asm_stats' include { INPUT_CHECK } from './subworkflows/local/input_check' include { SETUP_DBS } from './subworkflows/local/setup_dbs' include { CLEAN_GENOMES } from './subworkflows/local/clean_genomes' +include { MASK_GENOME } from './subworkflows/local/mask_genome' include { ANTISMASH_RUN } from './modules/local/antismash_run' include { INTERPROSCAN_RUN } from './modules/local/interproscan_run' include { SIGNALP_RUN } from './modules/local/signalp_run' @@ -1479,28 +1437,10 @@ workflow { } // ── Repeat masking ──────────────────────────────────────────────────────── - // predict_genome_ch carries the genome path to use for prediction — either - // the tantan soft-masked genome (default) or the clean unmasked genome - // (--run_repeatmasker false). - def predict_genome_ch - if (params.run_repeatmasker.toBoolean()) { - MASKREPEAT_TANTAN_RUN(clean_genome_ch) - predict_genome_ch = MASKREPEAT_TANTAN_RUN.out.masked - .map { meta, masked_fa -> - tuple(meta, masked_fa.toAbsolutePath().toString()) - } - } else { - // --run_repeatmasker false: use masked genome if a prior run produced it, else unmasked. - predict_genome_ch = clean_genome_ch - .map { meta, genome_fa -> - def masked = genomeFile("${launchDir}/input_clean_genomes/${meta.asmid}.masked.fasta") - def use_fa = masked.exists() ? masked.toString() : genome_fa - if (params.debug.toBoolean()) { - log.info "[DEBUG] ${meta.asmid}: genome_fa=${use_fa} (masked=${masked.exists()})" - } - tuple(meta, use_fa) - } - } + // predict_genome_ch: tantan soft-masked (default) or clean/prior-masked genome. + // MASK_GENOME handles the run_repeatmasker if/else and storeDir-cached masking. + MASK_GENOME(clean_genome_ch) + def predict_genome_ch = MASK_GENOME.out.genomes // Gate the predict chain on funannotate DB + augustus config being ready. // SETUP_DBS was already called above; its storeDir-cached outputs are free diff --git a/modules/local/maskrepeat_tantan_run.nf b/modules/local/maskrepeat_tantan_run.nf new file mode 100644 index 0000000..541bfc3 --- /dev/null +++ b/modules/local/maskrepeat_tantan_run.nf @@ -0,0 +1,41 @@ +// Soft-mask each assembly using funannotate mask with tantan. +// storeDir caches the masked FASTA alongside the clean genome at +// input_clean_genomes/.masked.fasta.gz, so this runs at most once per assembly. +// Downstream tools inflate it on the fly. +// Resources overridden per-profile by withName: '.*:MASKREPEAT_TANTAN_RUN' in +// conf/profile_annotate.config (short queue, 8 cpus, 16 GB, 2h). +process MASKREPEAT_TANTAN_RUN { + label 'funannotate' + label 'process_medium' + tag "${meta.id}" + + storeDir "${launchDir}/input_clean_genomes" + + input: + tuple val(meta), val(genome_fa) + + output: + tuple val(meta), path("${meta.asmid}.masked.fasta.gz"), emit: masked + + script: + def asmid = meta.asmid + """ + source /etc/profile.d/modules.sh 2>/dev/null || true + # Inflate a gzipped clean genome to a local uncompressed copy; funannotate cannot + # read a gzipped FASTA via -i. Plain (uncompressed) genomes pass through unchanged. + GENOME_FA="${genome_fa}" + case "\$GENOME_FA" in + *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; + *) GENOME_IN="\$GENOME_FA" ;; + esac + funannotate mask -i "\$GENOME_IN" -o ${asmid}.masked.fasta -m tantan --cpus ${task.cpus} + # Deliver the soft-masked genome gzip-compressed to save space; consumers inflate it. + pigz -f ${asmid}.masked.fasta + """ + + stub: + def asmid = meta.asmid + """ + echo ">stub_${asmid}_masked" | pigz -c > ${asmid}.masked.fasta.gz + """ +} diff --git a/subworkflows/local/mask_genome.nf b/subworkflows/local/mask_genome.nf new file mode 100644 index 0000000..7dc096e --- /dev/null +++ b/subworkflows/local/mask_genome.nf @@ -0,0 +1,47 @@ +/* + * MASK_GENOME — tantan soft-masking with run_repeatmasker fallback + * + * When params.run_repeatmasker is true: runs MASKREPEAT_TANTAN_RUN (storeDir-cached). + * When false: reuses the masked genome from a prior run if it exists, otherwise + * passes the clean genome through unmasked. + * + * Emits: + * genomes — channel: tuple(val(meta), val(genome_fa)) + * genome_fa is an absolute-path String consumed by the predict chain. + * + * Note: genomeFile() is a top-level def in funannotate.nf, not accessible here. + * The fallback path inlines the equivalent logic (prefer .masked.fasta.gz if + * non-empty, else .masked.fasta). + */ + +include { MASKREPEAT_TANTAN_RUN } from './../../modules/local/maskrepeat_tantan_run' + +workflow MASK_GENOME { + + take: + ch_clean // channel: tuple(val(meta), val(genome_fa)) + + main: + def ch_predict + if (params.run_repeatmasker.toBoolean()) { + MASKREPEAT_TANTAN_RUN(ch_clean) + ch_predict = MASKREPEAT_TANTAN_RUN.out.masked + .map { meta, masked_fa -> tuple(meta, masked_fa.toAbsolutePath().toString()) } + } else { + // --run_repeatmasker false: use masked genome from a prior run if available. + ch_predict = ch_clean + .map { meta, genome_fa -> + def fgz = file("${launchDir}/input_clean_genomes/${meta.asmid}.masked.fasta.gz") + def fa = file("${launchDir}/input_clean_genomes/${meta.asmid}.masked.fasta") + def masked = (fgz.exists() && fgz.size() > 0) ? fgz : fa + def use_fa = masked.exists() ? masked.toString() : genome_fa + if (params.debug.toBoolean()) { + log.info "[DEBUG] ${meta.asmid}: genome_fa=${use_fa} (masked=${masked.exists()})" + } + tuple(meta, use_fa) + } + } + + emit: + genomes = ch_predict +} From ffbc841b1aca64100466dacd57b61b6e4d6c04ac Mon Sep 17 00:00:00 2001 From: Jason Stajich Date: Mon, 29 Jun 2026 10:31:40 -0700 Subject: [PATCH 3/5] =?UTF-8?q?feat(#8):=20extract=20SRA=20query/fetch=20p?= =?UTF-8?q?rocesses=20=E2=86=92=20modules/local/=20+=20FETCH=5FRNASEQ=20su?= =?UTF-8?q?bworkflow?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Six processes moved out of funannotate.nf monolith: modules/local/sra_query.nf — single-species SRA query (storeDir) modules/local/sra_query_batch.nf — batched query (publishDir, maxForks 4) modules/local/collect_sra_query.nf — merge per-species CSVs into manifest modules/local/write_empty_reads.nf — zero-byte read placeholders for no-data species modules/local/sra_fetch.nf — PE download + normalize (EBI FTP fallback) modules/local/sra_fetch_se.nf — SE download + normalize Composed into subworkflows/local/fetch_rnaseq.nf (FETCH_RNASEQ): take: ch_species (species_tag, taxonid) - Step 1: SRA query (skip_sra_query reuses cached CSVs; else SRA_QUERY_BATCH) - Step 2: COLLECT_SRA_QUERY — always runs - Step 3: three-way branch (SRA_FETCH / SRA_FETCH_SE / WRITE_EMPTY_READS) skipped when stop_after_sra_query=true emit: reads (species_tag, r1, r2, se), sra_manifest funannotate.nf: ~110 lines of query/routing orchestration replaced by: FETCH_RNASEQ(sra_input) reads_ch = FETCH_RNASEQ.out.reads withName patterns updated to '.*:SRA_QUERY_BATCH' etc. in profile_annotate.config. SRA_QUERY (standalone, unused in current workflow) retains bare name pattern. Stub-run gate: 15 completed, 0 failed (run_sra_fetch=false default) Dual-stub gate: 17 completed, 0 failed (run_sra_fetch=true via params-file) FETCH_RNASEQ:SRA_QUERY_BATCH, COLLECT_SRA_QUERY, SRA_FETCH all correctly namespaced Co-Authored-By: Claude Sonnet 4.6 Claude-Session: https://claude.ai/code/session_01AbWwmm9v6dJn13sgAJn3Sa --- conf/profile_annotate.config | 11 +- funannotate.nf | 791 +---------------------------- modules/local/collect_sra_query.nf | 32 ++ modules/local/sra_fetch.nf | 267 ++++++++++ modules/local/sra_fetch_se.nf | 141 +++++ modules/local/sra_query.nf | 64 +++ modules/local/sra_query_batch.nf | 108 ++++ modules/local/write_empty_reads.nf | 32 ++ subworkflows/local/fetch_rnaseq.nf | 133 +++++ 9 files changed, 793 insertions(+), 786 deletions(-) create mode 100644 modules/local/collect_sra_query.nf create mode 100644 modules/local/sra_fetch.nf create mode 100644 modules/local/sra_fetch_se.nf create mode 100644 modules/local/sra_query.nf create mode 100644 modules/local/sra_query_batch.nf create mode 100644 modules/local/write_empty_reads.nf create mode 100644 subworkflows/local/fetch_rnaseq.nf diff --git a/conf/profile_annotate.config b/conf/profile_annotate.config index e130bbd..4fbabe3 100644 --- a/conf/profile_annotate.config +++ b/conf/profile_annotate.config @@ -145,29 +145,30 @@ process { queue = 'short'; cpus = 8; memory = '16 GB'; time = '2h' } + // SRA_QUERY (standalone, not in subworkflow) — matches bare name withName: 'SRA_QUERY' { queue = 'short'; cpus = 2; memory = '4 GB'; time = '60m' maxForks = 10; errorStrategy = 'retry'; maxRetries = 2 } - withName: 'SRA_QUERY_BATCH' { + withName: '.*:SRA_QUERY_BATCH' { queue = 'short'; cpus = 1; memory = '4 GB'; time = '2h' maxForks = 10; errorStrategy = 'retry'; maxRetries = 2 } - withName: 'COLLECT_SRA_QUERY' { + withName: '.*:COLLECT_SRA_QUERY' { queue = 'short'; cpus = 1; memory = '2 GB'; time = '20m' } - withName: 'WRITE_EMPTY_READS' { + withName: '.*:WRITE_EMPTY_READS' { queue = 'short'; cpus = 1; memory = '1 GB'; time = '5m' } - withName: 'SRA_FETCH' { + withName: '.*:SRA_FETCH' { errorStrategy = 'retry'; maxRetries = 3 queue = { task.attempt == 1 ? 'short' : 'epyc' } memory = { task.attempt == 1 ? '48 GB' : '192 GB' } time = { task.attempt == 1 ? '2h' : task.attempt == 2 ? '8h' : task.attempt == 3 ? '24h' : '48h' } cpus = { task.attempt == 1 ? 24 : 32 } } - withName: 'SRA_FETCH_SE' { + withName: '.*:SRA_FETCH_SE' { errorStrategy = 'retry'; maxRetries = 3 queue = { task.attempt == 1 ? 'short' : 'epyc' } memory = { task.attempt == 1 ? '48 GB' : '192 GB' } diff --git a/funannotate.nf b/funannotate.nf index 040315d..d9d879a 100644 --- a/funannotate.nf +++ b/funannotate.nf @@ -31,673 +31,6 @@ // FUNANNOTATE_PREDICT receives: tuple val(meta), val(genome_fa) // → emits metadata: tuple val(meta) -// Query NCBI SRA for available paired-end RNA-seq accessions per species. -// Lightweight: runs the esearch/efetch query only — no downloading. -// Records up to 5 candidates (sorted by spot count desc) in a per-species CSV. -// storeDir caches results so re-runs skip the network query. -// To invalidate the cache for a species, delete rnaseq_reads/sra_query/.sra_query.csv -process SRA_QUERY { - label 'edirect' - tag "$species_tag" - - storeDir "${launchDir}/rnaseq_reads/sra_query" - - cpus 1 - memory '4 GB' - time '30m' - - input: - tuple val(species_tag), val(taxonid) - - output: - tuple val(species_tag), path("${species_tag}.sra_query.csv"), emit: query_result - - script: - """ - set -euo pipefail - - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${species_tag}.sra_query.csv - - esearch -db sra \\ - -query "txid${taxonid}[Organism:noexp] AND RNA-Seq[Strategy] AND PAIRED[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND (BGISEQ[Platform] OR Illumina[Platform])" | \\ - efetch -format runinfo > _runinfo.tmp - - # col 1=Run, col 4=spots, col 13=LibraryStrategy, col 16=LibraryLayout, col 19=Platform - # Prepend a platform rank (0=Illumina, 1=BGI/other) so the top 5 prefer Illumina, - # then by spot count desc; BGI/other only fill remaining slots when Illumina runs out. - awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="PAIRED" && \$1~/^[SDE]RR/ && \$4+0>=250000 {rank=(\$19~/[Ii]llumina/)?0:1; printf "%d,%s,%s,%s\\n", rank, \$1, \$4, \$19}' _runinfo.tmp | \\ - sort -t',' -k1,1n -k3,3rn | \\ - head -n 5 | \\ - while IFS=',' read -r rank acc spots platform; do - printf '%s,%s,%s,%s,%s,PAIRED\\n' "${species_tag}" "${taxonid}" "\$acc" "\$spots" "\$platform" - done >> ${species_tag}.sra_query.csv - - rm -f _runinfo.tmp - NHITS=\$(awk 'END{print NR-1}' ${species_tag}.sra_query.csv) - echo "[INFO] Found \$NHITS paired-end SRA accessions for ${species_tag} (taxonid=${taxonid})" - - # SE fallback: if no PE hits found and enable_single_end is true, query SINGLE layout - if [ "${params.enable_single_end}" = "true" ] && [ "\$NHITS" -eq 0 ]; then - esearch -db sra \\ - -query "txid${taxonid}[Organism:noexp] AND RNA-Seq[Strategy] AND SINGLE[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND Illumina[Platform]" | \\ - efetch -format runinfo > _runinfo_se.tmp - awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="SINGLE" && \$1~/^[SDE]RR/ && \$4+0>=250000 {printf "%s,%s,%s\\n", \$1, \$4, \$19}' _runinfo_se.tmp | \\ - sort -t',' -k2 -rn | \\ - head -n ${params.max_rnaseq_se_runs} | \\ - while IFS=',' read -r acc spots platform; do - printf '%s,%s,%s,%s,%s,SINGLE\\n' "${species_tag}" "${taxonid}" "\$acc" "\$spots" "\$platform" - done >> ${species_tag}.sra_query.csv - rm -f _runinfo_se.tmp - NHITS=\$(awk 'END{print NR-1}' ${species_tag}.sra_query.csv) - echo "[INFO] SE fallback: found \$NHITS single-end accessions for ${species_tag}" - fi - """ - - stub: - """ - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${species_tag}.sra_query.csv - printf '%s,%s,SRR000001,1000000,ILLUMINA,PAIRED\n' "${species_tag}" "${taxonid}" >> ${species_tag}.sra_query.csv - echo "[STUB] SRA_QUERY for ${species_tag}" - """ -} - -// Batched SRA query: handles params.sra_query_batch_size species per SLURM job. -// maxForks 4 caps concurrent jobs to avoid overwhelming NCBI. -// Per-species esearch/efetch is retried up to 3 times inline with exponential -// backoff before writing an empty CSV. Existing per-species CSVs already -// present in rnaseq_reads/sra_query/ are reused without re-querying. -process SRA_QUERY_BATCH { - label 'edirect' - tag "${species_tags[0]}_+${species_tags.size() - 1}_more" - - publishDir "${launchDir}/rnaseq_reads/sra_query", mode: 'copy', overwrite: false - - maxForks 4 - cpus 1 - memory '4 GB' - time '4h' - - input: - tuple val(species_tags), val(taxonids) - - output: - path("*.sra_query.csv"), emit: query_results - - script: - def cache_dir = "${launchDir}/rnaseq_reads/sra_query" - def batch_args = [species_tags, taxonids].transpose() - .collect { st, tid -> "${st}\\t${tid}" } - .join('\\n') - """ - set -uo pipefail - - printf '${batch_args}\\n' > batch_input.tsv - - query_species() { - local stag="\$1" tid="\$2" attempt - - for attempt in 1 2 3; do - rm -f "_runinfo_\${stag}.tmp" - if timeout 120 bash -c \\ - "esearch -db sra -query 'txid\${tid}[Organism:noexp] AND RNA-Seq[Strategy] AND PAIRED[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND (BGISEQ[Platform] OR Illumina[Platform])' | efetch -format runinfo" \\ - < /dev/null > "_runinfo_\${stag}.tmp"; then - return 0 - fi - echo "[WARN] Attempt \${attempt}/3 failed or timed out for \${stag}" - [ "\${attempt}" -lt 3 ] && sleep \$((attempt * 30)) - done - rm -f "_runinfo_\${stag}.tmp" - return 1 - } - - while IFS=\$(printf '\\t') read -r species_tag taxonid; do - cached="${cache_dir}/\${species_tag}.sra_query.csv" - if [ -s "\$cached" ]; then - cp "\$cached" "\${species_tag}.sra_query.csv" - echo "[INFO] Reusing cached result for \${species_tag}" - continue - fi - - if query_species "\${species_tag}" "\${taxonid}"; then - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\\n' > "\${species_tag}.sra_query.csv" - # col 1=Run, col 4=spots, col 13=LibraryStrategy, col 16=LibraryLayout, col 19=Platform - # Prepend a platform rank (0=Illumina, 1=BGI/other) so the top 5 prefer Illumina, - # then by spot count desc; BGI/other only fill remaining slots when Illumina runs out. - awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="PAIRED" && \$1~/^[SDE]RR/ && \$4+0>=250000 {rank=(\$19~/[Ii]llumina/)?0:1; printf "%d,%s,%s,%s\\n", rank, \$1, \$4, \$19}' "_runinfo_\${species_tag}.tmp" | \\ - sort -t',' -k1,1n -k3,3rn | \\ - head -n 5 | \\ - while IFS=',' read -r rank acc spots platform; do - printf '%s,%s,%s,%s,%s,PAIRED\\n' "\${species_tag}" "\${taxonid}" "\$acc" "\$spots" "\$platform" - done >> "\${species_tag}.sra_query.csv" - rm -f "_runinfo_\${species_tag}.tmp" - NHITS=\$(awk 'END{print NR-1}' "\${species_tag}.sra_query.csv") - echo "[INFO] Found \$NHITS paired-end accessions for \${species_tag} (taxonid=\${taxonid})" - # SE fallback: if no PE hits and enable_single_end, query SINGLE layout - if [ "${params.enable_single_end}" = "true" ] && [ "\$NHITS" -eq 0 ]; then - rm -f "_runinfo_se_\${species_tag}.tmp" - if timeout 120 bash -c \\ - "esearch -db sra -query 'txid\${taxonid}[Organism:noexp] AND RNA-Seq[Strategy] AND SINGLE[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND Illumina[Platform]' | efetch -format runinfo" \\ - < /dev/null > "_runinfo_se_\${species_tag}.tmp"; then - awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="SINGLE" && \$1~/^[SDE]RR/ && \$4+0>=250000 {printf "%s,%s,%s\\n", \$1, \$4, \$19}' "_runinfo_se_\${species_tag}.tmp" | \\ - sort -t',' -k2 -rn | \\ - head -n ${params.max_rnaseq_se_runs} | \\ - while IFS=',' read -r acc spots platform; do - printf '%s,%s,%s,%s,%s,SINGLE\\n' "\${species_tag}" "\${taxonid}" "\$acc" "\$spots" "\$platform" - done >> "\${species_tag}.sra_query.csv" - fi - rm -f "_runinfo_se_\${species_tag}.tmp" - NHITS=\$(awk 'END{print NR-1}' "\${species_tag}.sra_query.csv") - echo "[INFO] SE fallback: \$NHITS single-end accessions for \${species_tag}" - fi - else - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\\n' > "\${species_tag}.sra_query.csv" - echo "[WARN] All 3 attempts failed for \${species_tag}; writing empty CSV" - fi - done < batch_input.tsv - """ - - stub: - def stub_args = [species_tags, taxonids].transpose() - .collect { st, tid -> "${st}\\t${tid}" } - .join('\\n') - """ - printf '${stub_args}\\n' > batch_input.tsv - while IFS=\$(printf '\\t') read -r species_tag taxonid; do - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\\n' > "\${species_tag}.sra_query.csv" - printf '%s,%s,SRR000001,1000000,ILLUMINA,PAIRED\\n' "\${species_tag}" "\${taxonid}" >> "\${species_tag}.sra_query.csv" - done < batch_input.tsv - echo "[STUB] SRA_QUERY_BATCH (${species_tags.size()} species)" - """ -} - -// Merge all per-species SRA query CSVs into a single named manifest. -// Output: {stem}.rnaseq_sra.csv written alongside the input samples file. -// Columns: species_tag, taxonid, sra_accession, spots, platform -process COLLECT_SRA_QUERY { - label 'setup' - publishDir { file(params.samples).parent.toAbsolutePath().toString() }, mode: 'copy' - - cpus 1 - memory '1 GB' - time '10m' - - input: - path(query_csvs) - val(stem) - - output: - path("${stem}.rnaseq_sra.csv"), emit: manifest - - script: - """ - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${stem}.rnaseq_sra.csv - for f in ${query_csvs}; do - tail -n +2 "\$f" >> ${stem}.rnaseq_sra.csv - done - NSPECIES=\$(awk -F',' 'NR>1{print \$1}' ${stem}.rnaseq_sra.csv | sort -u | wc -l) - NACCESSIONS=\$(awk 'NR>1' ${stem}.rnaseq_sra.csv | wc -l) - echo "[INFO] ${stem}.rnaseq_sra.csv: \$NACCESSIONS accessions across \$NSPECIES species with RNA-seq data" - """ - - stub: - """ - printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${stem}.rnaseq_sra.csv - """ -} - -// Write zero-byte paired FASTQ placeholder files for species with no SRA data. -// Called only for species whose SRA_QUERY CSV has no data rows, avoiding a -// SLURM job allocation for what would be an immediate empty-file write. -process WRITE_EMPTY_READS { - label 'setup' - tag "$species_tag" - - storeDir "${launchDir}/rnaseq_reads" - - cpus 1 - memory '1 GB' - time '5m' - - input: - val(species_tag) - - output: - tuple val(species_tag), path("${species_tag}_norm_R1.fastq.gz"), path("${species_tag}_norm_R2.fastq.gz"), - path("${species_tag}_norm_SE.fastq.gz"), emit: reads - - script: - """ - : > ${species_tag}_norm_R1.fastq.gz - : > ${species_tag}_norm_R2.fastq.gz - : > ${species_tag}_norm_SE.fastq.gz - echo "[INFO] No SRA data for ${species_tag}; created empty read placeholders" - """ - - stub: - """ - : > ${species_tag}_norm_R1.fastq.gz - : > ${species_tag}_norm_R2.fastq.gz - : > ${species_tag}_norm_SE.fastq.gz - """ -} - -// Download and normalize up to params.max_rnaseq_runs SRA accessions for species -// that have RNA-seq data. Accessions are read from the pre-queried per-species CSV -// produced by SRA_QUERY, so no NCBI network call is made here. -// Only invoked for species with data rows in their SRA_QUERY CSV; WRITE_EMPTY_READS -// handles the no-data case at the channel level. -process SRA_FETCH { - label 'sra' - tag "$species_tag" - - storeDir "${launchDir}/rnaseq_reads" - - cpus 32 - memory '96 GB' - time '2h' - - input: - tuple val(species_tag), path(sra_query_csv) - - output: - tuple val(species_tag), path("${species_tag}_norm_R1.fastq.gz"), path("${species_tag}_norm_R2.fastq.gz"), - path("${species_tag}_norm_SE.fastq.gz"), emit: reads - path("${species_tag}.se_candidates.csv"), optional: true, emit: se_candidates - path("${species_tag}.blacklist_candidates.csv"), optional: true, emit: blacklist_candidates - - script: - """ - - # Output files must always exist (storeDir requirement). - : > ${species_tag}_norm_R1.fastq.gz - : > ${species_tag}_norm_R2.fastq.gz - : > ${species_tag}_norm_SE.fastq.gz - - # Read pre-queried accessions from SRA_QUERY CSV (up to max_rnaseq_runs). - # Consult the project override file for per-accession actions: - # skip — exclude entirely - # rename_headers — parallel-fastq-dump as normal, then seqkit renumbers headers - # with sequential integers so R1/R2 names match exactly. - # Fixes BGI/MGISEQ read-name mismatch from block-parallel splits. - # Accessions absent from the override file use the default parallel-fastq-dump path. - BLACKLIST="${launchDir}/rnaseq_blacklist.csv" - RAW_ACCESSIONS=\$(awk -F',' 'NR>1 {print \$3}' ${sra_query_csv} | head -n ${params.max_rnaseq_runs}) - TAXONID=\$(awk -F',' 'NR==2 {print \$2; exit}' ${sra_query_csv}) - - # Helper: record an accession that yielded a non-empty _1 but no _2 (single-end data - # mislabeled PAIRED in SRA). Writes a row in blacklist format so it can be pasted - # straight into rnaseq_blacklist.csv (acc,species_tag,taxonid,SE_trinity); a trailing - # spots column (SRA-reported, col 4 of the query CSV; empty if unknown) rides along as - # info and is ignored by the blacklist parser. Collected to rnaseq_se_candidates.csv. - flag_se_candidate() { - local acc="\$1" spots - spots=\$(awk -F',' -v a="\$acc" 'NR>1 && \$3==a {print \$4; exit}' ${sra_query_csv}) - echo "[SE_CANDIDATE] \$acc ${species_tag} spots=\${spots:-NA} — add to rnaseq_blacklist.csv as SE_trinity" - echo "\$acc,${species_tag},\$TAXONID,SE_trinity,\${spots}" >> ${species_tag}.se_candidates.csv - } - - # Helper: record an accession whose download failed outright (parallel-fastq-dump and the - # EBI FTP fallback both produced nothing usable). Writes a row in rnaseq_blacklist.csv - # column order (acc,species_tag,taxonid,skip) with action=skip so it can be pasted - # straight into the blacklist after review. Trailing spots column is info only and is - # ignored by the blacklist parser. Collected to rnaseq_blacklist_candidates.csv. - flag_blacklist_candidate() { - local acc="\$1" spots - spots=\$(awk -F',' -v a="\$acc" 'NR>1 && \$3==a {print \$4; exit}' ${sra_query_csv}) - echo "[BLACKLIST_CANDIDATE] \$acc ${species_tag} spots=\${spots:-NA} — download failed; add to rnaseq_blacklist.csv as skip" - echo "\$acc,${species_tag},\$TAXONID,skip,\${spots}" >> ${species_tag}.blacklist_candidates.csv - } - - # Helper: look up the explicit override action for an accession from the blacklist - # (col 4 = action: skip | rename_headers). Returns empty string if not listed. - acc_action() { - local acc="\$1" - [ -f "\$BLACKLIST" ] && awk -F',' -v a="\$acc" 'NR>1 && \$1==a {print \$4; exit}' "\$BLACKLIST" || true - } - - # Helper: look up sequencing platform from the per-species sra_query CSV (col 5). - # Returns the raw SRA platform string, e.g. BGISEQ, ILLUMINA. - acc_platform() { - local acc="\$1" - awk -F',' -v a="\$acc" 'NR>1 && \$3==a {print \$5; exit}' ${sra_query_csv} - } - - # Resolve the effective download strategy for an accession: - # explicit blacklist action takes precedence; - # BGISEQ platform auto-maps to rename_headers; - # everything else is empty (default parallel-fastq-dump path). - acc_strategy() { - local acc="\$1" action platform - action=\$(acc_action "\$acc") - if [ -n "\$action" ]; then - echo "\$action" - else - platform=\$(acc_platform "\$acc") - case "\$platform" in - BGISEQ|BGI*|MGI*) echo "rename_headers" ;; - *) echo "" ;; - esac - fi - } - - # Compute the EBI FTP directory URL for an SRA accession. - # Path formula: vol1/fastq/// - # subdir is omitted for ≤9-char accessions; for longer ones it zero-pads - # chars 10-onward of the accession string into a 3-char field. - ebi_ftp_dir() { - local acc="\$1" len prefix base - len="\${#acc}" - prefix="\${acc:0:6}" - base="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/\${prefix}" - if [ "\$len" -le 9 ]; then printf '%s/%s' "\$base" "\$acc" - elif [ "\$len" -eq 10 ]; then printf '%s/00%s/%s' "\$base" "\${acc:9:1}" "\$acc" - elif [ "\$len" -eq 11 ]; then printf '%s/0%s/%s' "\$base" "\${acc:9:2}" "\$acc" - else printf '%s/%s/%s' "\$base" "\${acc:9:3}" "\$acc" - fi - } - - # Partition accessions into skip / active. - # SE_trinity entries are treated as skip here: they route to SRA_FETCH_SE instead. - ACCESSIONS="" - for ACC in \$RAW_ACCESSIONS; do - STRATEGY=\$(acc_strategy "\$ACC") - ACTION=\$(acc_action "\$ACC") - if [ "\$STRATEGY" = "skip" ] || [ "\$ACTION" = "SE_trinity" ]; then - echo "[INFO] Skipping accession \$ACC for ${species_tag} (action: \${ACTION:-\$STRATEGY})" - else - ACCESSIONS="\$ACCESSIONS \$ACC" - fi - done - ACCESSIONS=\$(echo "\$ACCESSIONS" | xargs) # trim whitespace - - if [ -z "\$ACCESSIONS" ]; then - echo "[INFO] No paired-end RNA-seq runs found for ${species_tag} (no accessions in query CSV or all skipped)" - else - echo "[INFO] SRA accessions for ${species_tag}: \$ACCESSIONS" - TMPDIR=\${SCRATCH:-/tmp} - mkdir -p reads - - # Download and concatenate in accession order so R1/R2 stay matched. - for ACC in \$ACCESSIONS; do - STRATEGY=\$(acc_strategy "\$ACC") - echo "[INFO] Downloading \$ACC (strategy: \${STRATEGY:-default}) ..." - maxspot="" - if [ ${params.max_rnaseq_spot} -gt "0" ]; then - maxspot="-X ${params.max_rnaseq_spot}" - fi - - if [ "\$STRATEGY" = "rename_headers" ]; then - # BGI/MGISEQ: parallel-fastq-dump block-splitting desynchronises spot IDs - # between R1 and R2. Download identically to the default path, then pipe - # each file through seqkit to replace every header with a sequential integer - # before handing off — guarantees R1 read N and R2 read N have matching names. - parallel-fastq-dump --sra-id \$ACC --threads ${task.cpus} \ - --outdir reads/ --split-files --gzip --tmpdir \$TMPDIR || { - echo "[WARN] Download failed for \$ACC (\$maxspot), skipping" - flag_blacklist_candidate "\$ACC" - continue - } - if [ -f reads/\${ACC}_1.fastq.gz ] && [ -f reads/\${ACC}_2.fastq.gz ]; then - seqkit replace -j ${task.cpus} -p '.+' -r '{nr}' reads/\${ACC}_1.fastq.gz \ - | ${params.fastq_hdr_script} --read 1 /dev/stdin \ - --max-reads ${params.max_rnaseq_reads} \ - | pigz -c >> \$TMPDIR/${species_tag}_R1.fastq.gz - seqkit replace -j ${task.cpus} -p '.+' -r '{nr}' reads/\${ACC}_2.fastq.gz \ - | ${params.fastq_hdr_script} --read 2 /dev/stdin \ - --max-reads ${params.max_rnaseq_reads} \ - | pigz -c >> \$TMPDIR/${species_tag}_R2.fastq.gz - rm reads/\${ACC}_[12].fastq.gz - elif [ -s reads/\${ACC}_1.fastq.gz ]; then - flag_se_candidate "\$ACC" - rm -f reads/\${ACC}_1.fastq.gz - else - echo "[WARN] Missing pair for \$ACC after download, skipping" - flag_blacklist_candidate "\$ACC" - fi - else - parallel-fastq-dump --sra-id \$ACC --threads ${task.cpus} \ - --outdir reads/ --split-files \$maxspot --gzip --tmpdir \$TMPDIR || { - echo "[WARN] Download failed for \$ACC (\$maxspot), skipping" - flag_blacklist_candidate "\$ACC" - continue - } - if [ -f reads/\${ACC}_1.fastq.gz ] && [ -f reads/\${ACC}_2.fastq.gz ]; then - parallel -j 2 ${params.fastq_hdr_script} --read {} reads/\${ACC}_{}.fastq.gz \ - --max-reads ${params.max_rnaseq_reads} \ - \\| pigz -c \\>\\> \$TMPDIR/${species_tag}_R{}.fastq.gz ::: 1 2 - rm reads/\${ACC}_[12].fastq.gz - elif [ -s reads/\${ACC}_1.fastq.gz ]; then - flag_se_candidate "\$ACC" - rm -f reads/\${ACC}_1.fastq.gz - else - echo "[WARN] Missing pair for \$ACC after download, skipping" - flag_blacklist_candidate "\$ACC" - fi - fi - done - rm -rf reads - ENFORCE="${params.readlen_script}" - [[ -x "\$ENFORCE" ]] || { echo "[ERROR] enforce_seqpair_readlen not found or not executable at \$ENFORCE"; exit 1; } - if ! "\$ENFORCE" in=\$TMPDIR/${species_tag}_R1.fastq.gz \ - in2=\$TMPDIR/${species_tag}_R2.fastq.gz \ - out=\$TMPDIR/${species_tag}_trunc_R1.fastq.gz \ - out2=\$TMPDIR/${species_tag}_trunc_R2.fastq.gz minlen=75; then - echo "[WARN] enforce_seqpair_readlen failed for ${species_tag} (likely parallel-fastq-dump read-name mismatch); retrying via EBI FTP..." - rm -f \$TMPDIR/${species_tag}_R1.fastq.gz \$TMPDIR/${species_tag}_R2.fastq.gz \ - \$TMPDIR/${species_tag}_trunc_R1.fastq.gz \$TMPDIR/${species_tag}_trunc_R2.fastq.gz - mkdir -p reads_ebi - for ACC in \$ACCESSIONS; do - EBI_DIR=\$(ebi_ftp_dir "\$ACC") - aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \ - "\${EBI_DIR}/\${ACC}_1.fastq.gz" -d reads_ebi/ -o "\${ACC}_1.fastq.gz" || true - aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \ - "\${EBI_DIR}/\${ACC}_2.fastq.gz" -d reads_ebi/ -o "\${ACC}_2.fastq.gz" || true - if [ ! -s reads_ebi/\${ACC}_2.fastq.gz ]; then - if [ -s reads_ebi/\${ACC}_1.fastq.gz ]; then - echo "[WARN] \$ACC: no paired-end R2 at EBI but R1 present (single-end data)" - flag_se_candidate "\$ACC" - else - echo "[WARN] \$ACC: no paired-end R2 at EBI (single-end or accession absent); skipping" - flag_blacklist_candidate "\$ACC" - fi - rm -f reads_ebi/\${ACC}_1.fastq.gz reads_ebi/\${ACC}_2.fastq.gz - continue - fi - if [ ! -s reads_ebi/\${ACC}_1.fastq.gz ]; then - echo "[WARN] \$ACC: R1 missing at EBI; skipping" - flag_blacklist_candidate "\$ACC" - rm -f reads_ebi/\${ACC}_2.fastq.gz - continue - fi - parallel -j 2 ${params.fastq_hdr_script} --read {} reads_ebi/\${ACC}_{}.fastq.gz \ - --max-reads ${params.max_rnaseq_reads} \ - \\| pigz -c \\>\\> \$TMPDIR/${species_tag}_R{}.fastq.gz ::: 1 2 - rm -f reads_ebi/\${ACC}_[12].fastq.gz - done - rm -rf reads_ebi - "\$ENFORCE" in=\$TMPDIR/${species_tag}_R1.fastq.gz \ - in2=\$TMPDIR/${species_tag}_R2.fastq.gz \ - out=\$TMPDIR/${species_tag}_trunc_R1.fastq.gz \ - out2=\$TMPDIR/${species_tag}_trunc_R2.fastq.gz minlen=75 \ - || { echo "[ERROR] enforce_seqpair_readlen failed for ${species_tag} even after EBI FTP fallback"; exit 1; } - fi - bbnorm.sh in=\$TMPDIR/${species_tag}_trunc_R1.fastq.gz in2=\$TMPDIR/${species_tag}_trunc_R2.fastq.gz \ - out1=\$TMPDIR/${species_tag}_norm_R1.fastq.gz \ - out2=\$TMPDIR/${species_tag}_norm_R2.fastq.gz target=30 ecc=t - - fastp --in1 \$TMPDIR/${species_tag}_norm_R1.fastq.gz \ - --in2 \$TMPDIR/${species_tag}_norm_R2.fastq.gz \ - --out1 ${species_tag}_norm_R1.fastq.gz --out2 ${species_tag}_norm_R2.fastq.gz \ - --thread ${task.cpus} --detect_adapter_for_pe \ - --cut_front --cut_front_window_size 1 --cut_front_mean_quality 5 \ - --cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 5 \ - --cut_right --cut_right_window_size 4 --cut_right_mean_quality 5 \ - --length_required 25 - - rm \$TMPDIR/${species_tag}_* - - NPAIRS=\$(zcat ${species_tag}_norm_R1.fastq.gz 2>/dev/null | awk 'NR%4==1' | wc -l || echo 0) - echo "[INFO] Combined \$NPAIRS normalized read pairs for ${species_tag}" - - mkdir -p "${launchDir}/rnaseq_reads" - fi - """ - - stub: - """ - : > ${species_tag}_norm_R1.fastq.gz - : > ${species_tag}_norm_R2.fastq.gz - : > ${species_tag}_norm_SE.fastq.gz - echo "[STUB] SRA_FETCH for ${species_tag}" - """ -} - -// Download and normalize single-end RNA-seq for species with no paired-end data, -// or where rnaseq_blacklist.csv lists accessions with action=SE_trinity (SRA metadata -// says PAIRED but the data is actually single-end). -// Outputs zero-byte PE stubs so all read channels carry the same 4-tuple shape. -// storeDir note: if a species was previously fetched via SRA_FETCH (PE), re-routing -// it to SRA_FETCH_SE requires deleting rnaseq_reads/${species_tag}_norm_* files first. -process SRA_FETCH_SE { - label 'sra' - tag "$species_tag" - - storeDir "${launchDir}/rnaseq_reads" - - cpus 32 - memory '96 GB' - time '2h' - - input: - tuple val(species_tag), path(sra_query_csv) - - output: - tuple val(species_tag), - path("${species_tag}_norm_R1.fastq.gz"), - path("${species_tag}_norm_R2.fastq.gz"), - path("${species_tag}_norm_SE.fastq.gz"), emit: reads - - script: - """ - - # Zero-byte PE stubs (this process produces SE output only). - : > ${species_tag}_norm_R1.fastq.gz - : > ${species_tag}_norm_R2.fastq.gz - : > ${species_tag}_norm_SE.fastq.gz - - BLACKLIST="${launchDir}/rnaseq_blacklist.csv" - - acc_action() { - local acc="\$1" - [ -f "\$BLACKLIST" ] && awk -F',' -v a="\$acc" 'NR>1 && \$1==a {print \$4; exit}' "\$BLACKLIST" || true - } - - ebi_ftp_dir() { - local acc="\$1" len prefix base - len="\${#acc}" - prefix="\${acc:0:6}" - base="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/\${prefix}" - if [ "\$len" -le 9 ]; then printf '%s/%s' "\$base" "\$acc" - elif [ "\$len" -eq 10 ]; then printf '%s/00%s/%s' "\$base" "\${acc:9:1}" "\$acc" - elif [ "\$len" -eq 11 ]; then printf '%s/0%s/%s' "\$base" "\${acc:9:2}" "\$acc" - else printf '%s/%s/%s' "\$base" "\${acc:9:3}" "\$acc" - fi - } - - # Collect SE accessions from the query CSV: - # SE_trinity entries: col 6 (layout) may be PAIRED in SRA but blacklist says SE_trinity. - # Download with --split-files and take _1 only (real SE data). - # SINGLE layout entries: col 6 == SINGLE; pfd gives ACC.fastq.gz or ACC_1.fastq.gz. - SE_ACCESSIONS="" - while IFS=',' read -r stag tid acc spots platform layout rest; do - [ "\$stag" = "species_tag" ] && continue - ACTION=\$(acc_action "\$acc") - if [ "\$ACTION" = "SE_trinity" ] || [ "\${layout:-PAIRED}" = "SINGLE" ]; then - SE_ACCESSIONS="\$SE_ACCESSIONS \$acc" - fi - done < ${sra_query_csv} - SE_ACCESSIONS=\$(echo "\$SE_ACCESSIONS" | tr ' ' '\\n' | grep -v '^\$' | head -n ${params.max_rnaseq_se_runs} | tr '\\n' ' ' | xargs) - - if [ -z "\$SE_ACCESSIONS" ]; then - echo "[INFO] No SE accessions found for ${species_tag}" - exit 0 - fi - - echo "[INFO] SE accessions for ${species_tag}: \$SE_ACCESSIONS" - TMPDIR=\${SCRATCH:-/tmp} - mkdir -p reads - - for ACC in \$SE_ACCESSIONS; do - ACTION=\$(acc_action "\$ACC") - echo "[INFO] Downloading \$ACC (SE mode, action: \${ACTION:-default}) ..." - - # Download with --split-files. For SE_trinity (mislabeled PAIRED) this yields _1/_2; - # we take only _1 (the actual SE reads) and discard _2. For genuine SINGLE layout, - # pfd typically produces ACC_1.fastq.gz or ACC.fastq.gz. - parallel-fastq-dump --sra-id "\$ACC" --threads ${task.cpus} \\ - --outdir reads/ --split-files --gzip --tmpdir "\$TMPDIR" || { - echo "[WARN] pfd failed for \$ACC; trying EBI FTP..." - EBI_DIR=\$(ebi_ftp_dir "\$ACC") - # SE_trinity (mislabeled PAIRED): EBI has ACC_1.fastq.gz - # Genuine SINGLE: EBI has ACC.fastq.gz - aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \\ - "\${EBI_DIR}/\${ACC}_1.fastq.gz" -d reads/ -o "\${ACC}_1.fastq.gz" 2>/dev/null || true - if [ ! -s "reads/\${ACC}_1.fastq.gz" ]; then - aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \\ - "\${EBI_DIR}/\${ACC}.fastq.gz" -d reads/ -o "\${ACC}.fastq.gz" 2>/dev/null || true - fi - } - - SE_FILE="" - if [ -s "reads/\${ACC}_1.fastq.gz" ]; then SE_FILE="reads/\${ACC}_1.fastq.gz" - elif [ -s "reads/\${ACC}.fastq.gz" ]; then SE_FILE="reads/\${ACC}.fastq.gz" - fi - - if [ -n "\$SE_FILE" ]; then - ${params.fastq_hdr_script} --read 1 "\$SE_FILE" \\ - --max-reads ${params.max_rnaseq_reads} \\ - | pigz -c >> "\$TMPDIR/${species_tag}_SE.fastq.gz" - rm -f "reads/\${ACC}_1.fastq.gz" "reads/\${ACC}_2.fastq.gz" "reads/\${ACC}.fastq.gz" - else - echo "[WARN] No SE file found for \$ACC after download; skipping" - rm -f "reads/\${ACC}"*.fastq.gz - fi - done - rm -rf reads - - if [ ! -s "\$TMPDIR/${species_tag}_SE.fastq.gz" ]; then - echo "[WARN] No SE reads downloaded for ${species_tag}; leaving empty SE placeholder" - exit 0 - fi - - bbnorm.sh in="\$TMPDIR/${species_tag}_SE.fastq.gz" \\ - out="\$TMPDIR/${species_tag}_norm_SE_bbn.fastq.gz" target=30 ecc=f - - fastp --in1 "\$TMPDIR/${species_tag}_norm_SE_bbn.fastq.gz" \\ - --out1 ${species_tag}_norm_SE.fastq.gz \\ - --thread ${task.cpus} \\ - --cut_front --cut_front_window_size 1 --cut_front_mean_quality 5 \\ - --cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 5 \\ - --cut_right --cut_right_window_size 4 --cut_right_mean_quality 5 \\ - --length_required 25 - - rm -f "\$TMPDIR/${species_tag}_SE.fastq.gz" "\$TMPDIR/${species_tag}_norm_SE_bbn.fastq.gz" - - NREADS=\$(zcat ${species_tag}_norm_SE.fastq.gz 2>/dev/null | awk 'NR%4==1' | wc -l || echo 0) - echo "[INFO] \$NREADS normalized SE reads for ${species_tag}" - - mkdir -p "${launchDir}/rnaseq_reads" - """ - - stub: - """ - : > ${species_tag}_norm_R1.fastq.gz - : > ${species_tag}_norm_R2.fastq.gz - : > ${species_tag}_norm_SE.fastq.gz - echo "[STUB] SRA_FETCH_SE for ${species_tag}" - """ -} - // Run funannotate train on the representative (first) assembly of each species, then // archive the Trinity-GG transcripts (normalized reads are in rnaseq_reads) // reads into rnaseq_data/ so all other strains can skip those expensive steps. @@ -1374,6 +707,7 @@ include { INPUT_CHECK } from './subworkflows/local/input_check' include { SETUP_DBS } from './subworkflows/local/setup_dbs' include { CLEAN_GENOMES } from './subworkflows/local/clean_genomes' include { MASK_GENOME } from './subworkflows/local/mask_genome' +include { FETCH_RNASEQ } from './subworkflows/local/fetch_rnaseq' include { ANTISMASH_RUN } from './modules/local/antismash_run' include { INTERPROSCAN_RUN } from './modules/local/interproscan_run' include { SIGNALP_RUN } from './modules/local/signalp_run' @@ -1454,129 +788,25 @@ workflow { .map { row -> row[0..-3] } // FUNANNOTATE_PREDICT input tuple drops taxonid (not needed after masking/clean). - // When SRA is enabled: SRA_FETCH fetches reads once per species; RNASEQ_PREPARE runs - // funannotate train on the representative assembly and archives Trinity-GG, trimmed, and - // normalized reads to rnaseq_data/; all other strains run FUNANNOTATE_TRAIN --trinity. + // When SRA is enabled: FETCH_RNASEQ fetches reads once per species; RNASEQ_PREPARE runs + // funannotate train on the representative assembly and archives Trinity-GG transcripts + // to rnaseq_data/; all other strains run FUNANNOTATE_TRAIN --trinity. def predict_input_ch def reads_ch = Channel.empty() if (params.run_sra_fetch.toBoolean()) { // Build per-species input: group assemblies, keep first taxonid per species. def sra_input = predict_genome_ch - .map { meta, genome_fa -> + .map { meta, _genome_fa -> def species_tag = meta.species.replaceAll(/\s+/, '_') tuple(species_tag, meta.taxonid) } .groupTuple(by: 0) .map { species_tag, taxonids -> tuple(species_tag, taxonids[0]) } - // Step 1: query or reuse cached per-species SRA query results. - // skip_sra_query=true reads existing CSVs from rnaseq_reads/sra_query/ directly, - // bypassing SRA_QUERY_BATCH entirely (no SLURM jobs submitted). - def sra_query_results - if (params.skip_sra_query.toBoolean()) { - sra_query_results = sra_input - .map { species_tag, _taxonid -> - def csv = file("${launchDir}/rnaseq_reads/sra_query/${species_tag}.sra_query.csv") - if (!csv.exists()) { - log.warn "skip_sra_query: no cached CSV for ${species_tag} — skipping this species" - return null - } - tuple(species_tag, csv) - } - .filter { it != null } - } else { - def sra_batched = sra_input - .collate(params.sra_query_batch_size) - .map { batch -> tuple(batch.collect { it[0] }, batch.collect { it[1] }) } - SRA_QUERY_BATCH(sra_batched) - sra_query_results = SRA_QUERY_BATCH.out.query_results - .flatten() - .map { csv -> tuple(csv.baseName.replaceAll(/\.sra_query$/, ''), csv) } - } - - // Step 2: Collect all per-species results into {stem}.rnaseq_sra.csv - def stem = file(params.samples).baseName - COLLECT_SRA_QUERY( - sra_query_results.map { _stag, csv -> csv }.collect(), - stem - ) - - if (!params.stop_after_sra_query.toBoolean()) { - // Step 3: Classify each species CSV for routing. - // Read blacklist once so closures below can check SE_trinity accessions. - // Uses a Map for O(1) lookup. - def blPath = file("${launchDir}/rnaseq_blacklist.csv") - def blMap = blPath.exists() - ? blPath.readLines().drop(1) - .findAll { it.trim() && !it.startsWith('#') } - .collectEntries { line -> - def cols = line.split(',') - cols.size() >= 4 ? [(cols[0].trim()): cols[3].trim()] : [:] - } - : [:] - - // csvHasPE: CSV has at least one PAIRED accession not blocked or overridden to SE. - def csvHasPE = { csv -> - csv.readLines().drop(1).findAll { it.trim() }.any { line -> - def cols = line.split(',') - if (cols.size() < 3) return false - def layout = cols.size() > 5 ? cols[5].trim() : 'PAIRED' - def action = blMap.get(cols[2].trim(), '') - layout == 'PAIRED' && action != 'skip' && action != 'SE_trinity' - } - } - - // csvHasSEtrinity: CSV has PAIRED accessions overridden to SE via SE_trinity blacklist. - // These bypass the enable_single_end gate — they are a manual per-accession override. - def csvHasSEtrinity = { csv -> - csv.readLines().drop(1).findAll { it.trim() }.any { line -> - def cols = line.split(',') - cols.size() >= 3 && blMap.get(cols[2].trim(), '') == 'SE_trinity' - } - } - - // csvHasSingleLayout: CSV has at least one genuine SINGLE-layout accession. - // Only active when enable_single_end=true. - def csvHasSingleLayout = { csv -> - csv.readLines().drop(1).findAll { it.trim() }.any { line -> - def cols = line.split(',') - cols.size() > 5 && cols[5].trim() == 'SINGLE' && blMap.get(cols[2].trim(), '') != 'skip' - } - } - - // Three-way branch: - // has_pe → SRA_FETCH (PE wins; SE_trinity entries ignored here, handled by SRA_FETCH) - // has_se → SRA_FETCH_SE (SE_trinity always; SINGLE layout only if enable_single_end) - // no_data → WRITE_EMPTY_READS - def branched_sra = sra_query_results - .branch { - has_pe: csvHasPE.call(it[1]) - has_se: csvHasSEtrinity.call(it[1]) || - (params.enable_single_end.toBoolean() && csvHasSingleLayout.call(it[1])) - no_data: true - } + FETCH_RNASEQ(sra_input) + reads_ch = FETCH_RNASEQ.out.reads - SRA_FETCH(branched_sra.has_pe) - SRA_FETCH_SE(branched_sra.has_se) - WRITE_EMPTY_READS(branched_sra.no_data.map { stag, _csv -> stag }) - - // Accessions found to be single-end (non-empty _1, no _2) during the PE fetch are - // recorded as blacklist-ready rows; merge all per-task notes into one reviewable - // file at the project root. Add these to rnaseq_blacklist.csv as SE_trinity and - // rerun to route them through SRA_FETCH_SE. - SRA_FETCH.out.se_candidates - .collectFile(name: 'rnaseq_se_candidates.csv', storeDir: launchDir, newLine: false) - - // Accessions whose download failed outright (parallel-fastq-dump + EBI FTP both - // produced nothing) are recorded in rnaseq_blacklist.csv column order so they can be - // reviewed and pasted straight into the blacklist as skip entries; merged at root. - SRA_FETCH.out.blacklist_candidates - .collectFile(name: 'rnaseq_blacklist_candidates.csv', storeDir: launchDir, newLine: false) - reads_ch = SRA_FETCH.out.reads - .mix(SRA_FETCH_SE.out.reads) - .mix(WRITE_EMPTY_READS.out.reads) - - if (!params.stop_after_sra_fetch.toBoolean()) { + if (!params.stop_after_sra_fetch.toBoolean() && !params.stop_after_sra_query.toBoolean()) { // Build per-assembly channel keyed by species_tag with SRA reads joined. // reads_ch is now a 4-tuple: (species_tag, r1, r2, se) def assembly_with_reads = predict_genome_ch @@ -1658,13 +888,12 @@ workflow { } FUNANNOTATE_TRAIN(train_todo) predict_input_ch = FUNANNOTATE_TRAIN.out.mix(train_done).mix(predict_no_rnaseq) - } // end if (!params.stop_after_sra_fetch) - } // end if (!params.stop_after_sra_query) + } // end if (!stop_after_sra_fetch && !stop_after_sra_query) } else { predict_input_ch = predict_genome_ch } - if ((!params.stop_after_sra_fetch.toBoolean() && !params.stop_after_sra_query.toBoolean()) || !params.run_sra_fetch.toBoolean()) { + if (!params.run_sra_fetch.toBoolean() || (!params.stop_after_sra_fetch.toBoolean() && !params.stop_after_sra_query.toBoolean())) { def predict_ch = predict_input_ch .filter { meta, _gfa -> gbkResult("${params.target}/${meta.id}/predict_results", meta.id as String) == null || staleRnaseq(meta.id as String, meta.species as String) diff --git a/modules/local/collect_sra_query.nf b/modules/local/collect_sra_query.nf new file mode 100644 index 0000000..f8ef2b9 --- /dev/null +++ b/modules/local/collect_sra_query.nf @@ -0,0 +1,32 @@ +// Merge all per-species SRA query CSVs into a single named manifest. +// Output: {stem}.rnaseq_sra.csv written alongside the input samples file. +// Columns: species_tag, taxonid, sra_accession, spots, platform, layout +process COLLECT_SRA_QUERY { + label 'setup' + label 'process_single' + + publishDir { file(params.samples).parent.toAbsolutePath().toString() }, mode: 'copy' + + input: + path(query_csvs) + val(stem) + + output: + path("${stem}.rnaseq_sra.csv"), emit: manifest + + script: + """ + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${stem}.rnaseq_sra.csv + for f in ${query_csvs}; do + tail -n +2 "\$f" >> ${stem}.rnaseq_sra.csv + done + NSPECIES=\$(awk -F',' 'NR>1{print \$1}' ${stem}.rnaseq_sra.csv | sort -u | wc -l) + NACCESSIONS=\$(awk 'NR>1' ${stem}.rnaseq_sra.csv | wc -l) + echo "[INFO] ${stem}.rnaseq_sra.csv: \$NACCESSIONS accessions across \$NSPECIES species with RNA-seq data" + """ + + stub: + """ + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${stem}.rnaseq_sra.csv + """ +} diff --git a/modules/local/sra_fetch.nf b/modules/local/sra_fetch.nf new file mode 100644 index 0000000..1dcc514 --- /dev/null +++ b/modules/local/sra_fetch.nf @@ -0,0 +1,267 @@ +// Download and normalize up to params.max_rnaseq_runs paired-end SRA accessions. +// Reads accessions from the pre-queried per-species CSV (SRA_QUERY/SRA_QUERY_BATCH); +// no NCBI network call is made here. storeDir caches normalized reads. +// Resources overridden per-profile by withName: '.*:SRA_FETCH' in +// conf/profile_annotate.config (retry-escalating queue/memory/time). +process SRA_FETCH { + label 'sra' + label 'process_high' + tag "$species_tag" + + storeDir "${launchDir}/rnaseq_reads" + + input: + tuple val(species_tag), path(sra_query_csv) + + output: + tuple val(species_tag), path("${species_tag}_norm_R1.fastq.gz"), path("${species_tag}_norm_R2.fastq.gz"), + path("${species_tag}_norm_SE.fastq.gz"), emit: reads + path("${species_tag}.se_candidates.csv"), optional: true, emit: se_candidates + path("${species_tag}.blacklist_candidates.csv"), optional: true, emit: blacklist_candidates + + script: + """ + + # Output files must always exist (storeDir requirement). + : > ${species_tag}_norm_R1.fastq.gz + : > ${species_tag}_norm_R2.fastq.gz + : > ${species_tag}_norm_SE.fastq.gz + + # Read pre-queried accessions from SRA_QUERY CSV (up to max_rnaseq_runs). + # Consult the project override file for per-accession actions: + # skip — exclude entirely + # rename_headers — parallel-fastq-dump as normal, then seqkit renumbers headers + # with sequential integers so R1/R2 names match exactly. + # Fixes BGI/MGISEQ read-name mismatch from block-parallel splits. + # Accessions absent from the override file use the default parallel-fastq-dump path. + BLACKLIST="${launchDir}/rnaseq_blacklist.csv" + RAW_ACCESSIONS=\$(awk -F',' 'NR>1 {print \$3}' ${sra_query_csv} | head -n ${params.max_rnaseq_runs}) + TAXONID=\$(awk -F',' 'NR==2 {print \$2; exit}' ${sra_query_csv}) + + # Helper: record an accession that yielded a non-empty _1 but no _2 (single-end data + # mislabeled PAIRED in SRA). Writes a row in blacklist format so it can be pasted + # straight into rnaseq_blacklist.csv (acc,species_tag,taxonid,SE_trinity); a trailing + # spots column (SRA-reported, col 4 of the query CSV; empty if unknown) rides along as + # info and is ignored by the blacklist parser. Collected to rnaseq_se_candidates.csv. + flag_se_candidate() { + local acc="\$1" spots + spots=\$(awk -F',' -v a="\$acc" 'NR>1 && \$3==a {print \$4; exit}' ${sra_query_csv}) + echo "[SE_CANDIDATE] \$acc ${species_tag} spots=\${spots:-NA} — add to rnaseq_blacklist.csv as SE_trinity" + echo "\$acc,${species_tag},\$TAXONID,SE_trinity,\${spots}" >> ${species_tag}.se_candidates.csv + } + + # Helper: record an accession whose download failed outright (parallel-fastq-dump and the + # EBI FTP fallback both produced nothing usable). Writes a row in rnaseq_blacklist.csv + # column order (acc,species_tag,taxonid,skip) with action=skip so it can be pasted + # straight into the blacklist after review. Trailing spots column is info only and is + # ignored by the blacklist parser. Collected to rnaseq_blacklist_candidates.csv. + flag_blacklist_candidate() { + local acc="\$1" spots + spots=\$(awk -F',' -v a="\$acc" 'NR>1 && \$3==a {print \$4; exit}' ${sra_query_csv}) + echo "[BLACKLIST_CANDIDATE] \$acc ${species_tag} spots=\${spots:-NA} — download failed; add to rnaseq_blacklist.csv as skip" + echo "\$acc,${species_tag},\$TAXONID,skip,\${spots}" >> ${species_tag}.blacklist_candidates.csv + } + + # Helper: look up the explicit override action for an accession from the blacklist + # (col 4 = action: skip | rename_headers). Returns empty string if not listed. + acc_action() { + local acc="\$1" + [ -f "\$BLACKLIST" ] && awk -F',' -v a="\$acc" 'NR>1 && \$1==a {print \$4; exit}' "\$BLACKLIST" || true + } + + # Helper: look up sequencing platform from the per-species sra_query CSV (col 5). + # Returns the raw SRA platform string, e.g. BGISEQ, ILLUMINA. + acc_platform() { + local acc="\$1" + awk -F',' -v a="\$acc" 'NR>1 && \$3==a {print \$5; exit}' ${sra_query_csv} + } + + # Resolve the effective download strategy for an accession: + # explicit blacklist action takes precedence; + # BGISEQ platform auto-maps to rename_headers; + # everything else is empty (default parallel-fastq-dump path). + acc_strategy() { + local acc="\$1" action platform + action=\$(acc_action "\$acc") + if [ -n "\$action" ]; then + echo "\$action" + else + platform=\$(acc_platform "\$acc") + case "\$platform" in + BGISEQ|BGI*|MGI*) echo "rename_headers" ;; + *) echo "" ;; + esac + fi + } + + # Compute the EBI FTP directory URL for an SRA accession. + # Path formula: vol1/fastq/// + # subdir is omitted for ≤9-char accessions; for longer ones it zero-pads + # chars 10-onward of the accession string into a 3-char field. + ebi_ftp_dir() { + local acc="\$1" len prefix base + len="\${#acc}" + prefix="\${acc:0:6}" + base="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/\${prefix}" + if [ "\$len" -le 9 ]; then printf '%s/%s' "\$base" "\$acc" + elif [ "\$len" -eq 10 ]; then printf '%s/00%s/%s' "\$base" "\${acc:9:1}" "\$acc" + elif [ "\$len" -eq 11 ]; then printf '%s/0%s/%s' "\$base" "\${acc:9:2}" "\$acc" + else printf '%s/%s/%s' "\$base" "\${acc:9:3}" "\$acc" + fi + } + + # Partition accessions into skip / active. + # SE_trinity entries are treated as skip here: they route to SRA_FETCH_SE instead. + ACCESSIONS="" + for ACC in \$RAW_ACCESSIONS; do + STRATEGY=\$(acc_strategy "\$ACC") + ACTION=\$(acc_action "\$ACC") + if [ "\$STRATEGY" = "skip" ] || [ "\$ACTION" = "SE_trinity" ]; then + echo "[INFO] Skipping accession \$ACC for ${species_tag} (action: \${ACTION:-\$STRATEGY})" + else + ACCESSIONS="\$ACCESSIONS \$ACC" + fi + done + ACCESSIONS=\$(echo "\$ACCESSIONS" | xargs) # trim whitespace + + if [ -z "\$ACCESSIONS" ]; then + echo "[INFO] No paired-end RNA-seq runs found for ${species_tag} (no accessions in query CSV or all skipped)" + else + echo "[INFO] SRA accessions for ${species_tag}: \$ACCESSIONS" + TMPDIR=\${SCRATCH:-/tmp} + mkdir -p reads + + # Download and concatenate in accession order so R1/R2 stay matched. + for ACC in \$ACCESSIONS; do + STRATEGY=\$(acc_strategy "\$ACC") + echo "[INFO] Downloading \$ACC (strategy: \${STRATEGY:-default}) ..." + maxspot="" + if [ ${params.max_rnaseq_spot} -gt "0" ]; then + maxspot="-X ${params.max_rnaseq_spot}" + fi + + if [ "\$STRATEGY" = "rename_headers" ]; then + # BGI/MGISEQ: parallel-fastq-dump block-splitting desynchronises spot IDs + # between R1 and R2. Download identically to the default path, then pipe + # each file through seqkit to replace every header with a sequential integer + # before handing off — guarantees R1 read N and R2 read N have matching names. + parallel-fastq-dump --sra-id \$ACC --threads ${task.cpus} \ + --outdir reads/ --split-files --gzip --tmpdir \$TMPDIR || { + echo "[WARN] Download failed for \$ACC (\$maxspot), skipping" + flag_blacklist_candidate "\$ACC" + continue + } + if [ -f reads/\${ACC}_1.fastq.gz ] && [ -f reads/\${ACC}_2.fastq.gz ]; then + seqkit replace -j ${task.cpus} -p '.+' -r '{nr}' reads/\${ACC}_1.fastq.gz \ + | ${params.fastq_hdr_script} --read 1 /dev/stdin \ + --max-reads ${params.max_rnaseq_reads} \ + | pigz -c >> \$TMPDIR/${species_tag}_R1.fastq.gz + seqkit replace -j ${task.cpus} -p '.+' -r '{nr}' reads/\${ACC}_2.fastq.gz \ + | ${params.fastq_hdr_script} --read 2 /dev/stdin \ + --max-reads ${params.max_rnaseq_reads} \ + | pigz -c >> \$TMPDIR/${species_tag}_R2.fastq.gz + rm reads/\${ACC}_[12].fastq.gz + elif [ -s reads/\${ACC}_1.fastq.gz ]; then + flag_se_candidate "\$ACC" + rm -f reads/\${ACC}_1.fastq.gz + else + echo "[WARN] Missing pair for \$ACC after download, skipping" + flag_blacklist_candidate "\$ACC" + fi + else + parallel-fastq-dump --sra-id \$ACC --threads ${task.cpus} \ + --outdir reads/ --split-files \$maxspot --gzip --tmpdir \$TMPDIR || { + echo "[WARN] Download failed for \$ACC (\$maxspot), skipping" + flag_blacklist_candidate "\$ACC" + continue + } + if [ -f reads/\${ACC}_1.fastq.gz ] && [ -f reads/\${ACC}_2.fastq.gz ]; then + parallel -j 2 ${params.fastq_hdr_script} --read {} reads/\${ACC}_{}.fastq.gz \ + --max-reads ${params.max_rnaseq_reads} \ + \\| pigz -c \\>\\> \$TMPDIR/${species_tag}_R{}.fastq.gz ::: 1 2 + rm reads/\${ACC}_[12].fastq.gz + elif [ -s reads/\${ACC}_1.fastq.gz ]; then + flag_se_candidate "\$ACC" + rm -f reads/\${ACC}_1.fastq.gz + else + echo "[WARN] Missing pair for \$ACC after download, skipping" + flag_blacklist_candidate "\$ACC" + fi + fi + done + rm -rf reads + ENFORCE="${params.readlen_script}" + [[ -x "\$ENFORCE" ]] || { echo "[ERROR] enforce_seqpair_readlen not found or not executable at \$ENFORCE"; exit 1; } + if ! "\$ENFORCE" in=\$TMPDIR/${species_tag}_R1.fastq.gz \ + in2=\$TMPDIR/${species_tag}_R2.fastq.gz \ + out=\$TMPDIR/${species_tag}_trunc_R1.fastq.gz \ + out2=\$TMPDIR/${species_tag}_trunc_R2.fastq.gz minlen=75; then + echo "[WARN] enforce_seqpair_readlen failed for ${species_tag} (likely parallel-fastq-dump read-name mismatch); retrying via EBI FTP..." + rm -f \$TMPDIR/${species_tag}_R1.fastq.gz \$TMPDIR/${species_tag}_R2.fastq.gz \ + \$TMPDIR/${species_tag}_trunc_R1.fastq.gz \$TMPDIR/${species_tag}_trunc_R2.fastq.gz + mkdir -p reads_ebi + for ACC in \$ACCESSIONS; do + EBI_DIR=\$(ebi_ftp_dir "\$ACC") + aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \ + "\${EBI_DIR}/\${ACC}_1.fastq.gz" -d reads_ebi/ -o "\${ACC}_1.fastq.gz" || true + aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \ + "\${EBI_DIR}/\${ACC}_2.fastq.gz" -d reads_ebi/ -o "\${ACC}_2.fastq.gz" || true + if [ ! -s reads_ebi/\${ACC}_2.fastq.gz ]; then + if [ -s reads_ebi/\${ACC}_1.fastq.gz ]; then + echo "[WARN] \$ACC: no paired-end R2 at EBI but R1 present (single-end data)" + flag_se_candidate "\$ACC" + else + echo "[WARN] \$ACC: no paired-end R2 at EBI (single-end or accession absent); skipping" + flag_blacklist_candidate "\$ACC" + fi + rm -f reads_ebi/\${ACC}_1.fastq.gz reads_ebi/\${ACC}_2.fastq.gz + continue + fi + if [ ! -s reads_ebi/\${ACC}_1.fastq.gz ]; then + echo "[WARN] \$ACC: R1 missing at EBI; skipping" + flag_blacklist_candidate "\$ACC" + rm -f reads_ebi/\${ACC}_2.fastq.gz + continue + fi + parallel -j 2 ${params.fastq_hdr_script} --read {} reads_ebi/\${ACC}_{}.fastq.gz \ + --max-reads ${params.max_rnaseq_reads} \ + \\| pigz -c \\>\\> \$TMPDIR/${species_tag}_R{}.fastq.gz ::: 1 2 + rm -f reads_ebi/\${ACC}_[12].fastq.gz + done + rm -rf reads_ebi + "\$ENFORCE" in=\$TMPDIR/${species_tag}_R1.fastq.gz \ + in2=\$TMPDIR/${species_tag}_R2.fastq.gz \ + out=\$TMPDIR/${species_tag}_trunc_R1.fastq.gz \ + out2=\$TMPDIR/${species_tag}_trunc_R2.fastq.gz minlen=75 \ + || { echo "[ERROR] enforce_seqpair_readlen failed for ${species_tag} even after EBI FTP fallback"; exit 1; } + fi + bbnorm.sh in=\$TMPDIR/${species_tag}_trunc_R1.fastq.gz in2=\$TMPDIR/${species_tag}_trunc_R2.fastq.gz \ + out1=\$TMPDIR/${species_tag}_norm_R1.fastq.gz \ + out2=\$TMPDIR/${species_tag}_norm_R2.fastq.gz target=30 ecc=t + + fastp --in1 \$TMPDIR/${species_tag}_norm_R1.fastq.gz \ + --in2 \$TMPDIR/${species_tag}_norm_R2.fastq.gz \ + --out1 ${species_tag}_norm_R1.fastq.gz --out2 ${species_tag}_norm_R2.fastq.gz \ + --thread ${task.cpus} --detect_adapter_for_pe \ + --cut_front --cut_front_window_size 1 --cut_front_mean_quality 5 \ + --cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 5 \ + --cut_right --cut_right_window_size 4 --cut_right_mean_quality 5 \ + --length_required 25 + + rm \$TMPDIR/${species_tag}_* + + NPAIRS=\$(zcat ${species_tag}_norm_R1.fastq.gz 2>/dev/null | awk 'NR%4==1' | wc -l || echo 0) + echo "[INFO] Combined \$NPAIRS normalized read pairs for ${species_tag}" + + mkdir -p "${launchDir}/rnaseq_reads" + fi + """ + + stub: + """ + : > ${species_tag}_norm_R1.fastq.gz + : > ${species_tag}_norm_R2.fastq.gz + : > ${species_tag}_norm_SE.fastq.gz + echo "[STUB] SRA_FETCH for ${species_tag}" + """ +} diff --git a/modules/local/sra_fetch_se.nf b/modules/local/sra_fetch_se.nf new file mode 100644 index 0000000..568a867 --- /dev/null +++ b/modules/local/sra_fetch_se.nf @@ -0,0 +1,141 @@ +// Download and normalize single-end RNA-seq accessions (SE_trinity blacklist +// overrides or genuine SINGLE-layout SRA entries). Produces zero-byte PE stubs +// so all read channels carry the same 4-tuple shape. +// Resources overridden per-profile by withName: '.*:SRA_FETCH_SE' in +// conf/profile_annotate.config (retry-escalating queue/memory/time). +process SRA_FETCH_SE { + label 'sra' + label 'process_high' + tag "$species_tag" + + storeDir "${launchDir}/rnaseq_reads" + + input: + tuple val(species_tag), path(sra_query_csv) + + output: + tuple val(species_tag), + path("${species_tag}_norm_R1.fastq.gz"), + path("${species_tag}_norm_R2.fastq.gz"), + path("${species_tag}_norm_SE.fastq.gz"), emit: reads + + script: + """ + + # Zero-byte PE stubs (this process produces SE output only). + : > ${species_tag}_norm_R1.fastq.gz + : > ${species_tag}_norm_R2.fastq.gz + : > ${species_tag}_norm_SE.fastq.gz + + BLACKLIST="${launchDir}/rnaseq_blacklist.csv" + + acc_action() { + local acc="\$1" + [ -f "\$BLACKLIST" ] && awk -F',' -v a="\$acc" 'NR>1 && \$1==a {print \$4; exit}' "\$BLACKLIST" || true + } + + ebi_ftp_dir() { + local acc="\$1" len prefix base + len="\${#acc}" + prefix="\${acc:0:6}" + base="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/\${prefix}" + if [ "\$len" -le 9 ]; then printf '%s/%s' "\$base" "\$acc" + elif [ "\$len" -eq 10 ]; then printf '%s/00%s/%s' "\$base" "\${acc:9:1}" "\$acc" + elif [ "\$len" -eq 11 ]; then printf '%s/0%s/%s' "\$base" "\${acc:9:2}" "\$acc" + else printf '%s/%s/%s' "\$base" "\${acc:9:3}" "\$acc" + fi + } + + # Collect SE accessions from the query CSV: + # SE_trinity entries: col 6 (layout) may be PAIRED in SRA but blacklist says SE_trinity. + # Download with --split-files and take _1 only (real SE data). + # SINGLE layout entries: col 6 == SINGLE; pfd gives ACC.fastq.gz or ACC_1.fastq.gz. + SE_ACCESSIONS="" + while IFS=',' read -r stag tid acc spots platform layout rest; do + [ "\$stag" = "species_tag" ] && continue + ACTION=\$(acc_action "\$acc") + if [ "\$ACTION" = "SE_trinity" ] || [ "\${layout:-PAIRED}" = "SINGLE" ]; then + SE_ACCESSIONS="\$SE_ACCESSIONS \$acc" + fi + done < ${sra_query_csv} + SE_ACCESSIONS=\$(echo "\$SE_ACCESSIONS" | tr ' ' '\\n' | grep -v '^\$' | head -n ${params.max_rnaseq_se_runs} | tr '\\n' ' ' | xargs) + + if [ -z "\$SE_ACCESSIONS" ]; then + echo "[INFO] No SE accessions found for ${species_tag}" + exit 0 + fi + + echo "[INFO] SE accessions for ${species_tag}: \$SE_ACCESSIONS" + TMPDIR=\${SCRATCH:-/tmp} + mkdir -p reads + + for ACC in \$SE_ACCESSIONS; do + ACTION=\$(acc_action "\$ACC") + echo "[INFO] Downloading \$ACC (SE mode, action: \${ACTION:-default}) ..." + + # Download with --split-files. For SE_trinity (mislabeled PAIRED) this yields _1/_2; + # we take only _1 (the actual SE reads) and discard _2. For genuine SINGLE layout, + # pfd typically produces ACC_1.fastq.gz or ACC.fastq.gz. + parallel-fastq-dump --sra-id "\$ACC" --threads ${task.cpus} \\ + --outdir reads/ --split-files --gzip --tmpdir "\$TMPDIR" || { + echo "[WARN] pfd failed for \$ACC; trying EBI FTP..." + EBI_DIR=\$(ebi_ftp_dir "\$ACC") + # SE_trinity (mislabeled PAIRED): EBI has ACC_1.fastq.gz + # Genuine SINGLE: EBI has ACC.fastq.gz + aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \\ + "\${EBI_DIR}/\${ACC}_1.fastq.gz" -d reads/ -o "\${ACC}_1.fastq.gz" 2>/dev/null || true + if [ ! -s "reads/\${ACC}_1.fastq.gz" ]; then + aria2c --max-connection-per-server=4 --min-split-size=1M --max-tries=3 --retry-wait=5 \\ + "\${EBI_DIR}/\${ACC}.fastq.gz" -d reads/ -o "\${ACC}.fastq.gz" 2>/dev/null || true + fi + } + + SE_FILE="" + if [ -s "reads/\${ACC}_1.fastq.gz" ]; then SE_FILE="reads/\${ACC}_1.fastq.gz" + elif [ -s "reads/\${ACC}.fastq.gz" ]; then SE_FILE="reads/\${ACC}.fastq.gz" + fi + + if [ -n "\$SE_FILE" ]; then + ${params.fastq_hdr_script} --read 1 "\$SE_FILE" \\ + --max-reads ${params.max_rnaseq_reads} \\ + | pigz -c >> "\$TMPDIR/${species_tag}_SE.fastq.gz" + rm -f "reads/\${ACC}_1.fastq.gz" "reads/\${ACC}_2.fastq.gz" "reads/\${ACC}.fastq.gz" + else + echo "[WARN] No SE file found for \$ACC after download; skipping" + rm -f "reads/\${ACC}"*.fastq.gz + fi + done + rm -rf reads + + if [ ! -s "\$TMPDIR/${species_tag}_SE.fastq.gz" ]; then + echo "[WARN] No SE reads downloaded for ${species_tag}; leaving empty SE placeholder" + exit 0 + fi + + bbnorm.sh in="\$TMPDIR/${species_tag}_SE.fastq.gz" \\ + out="\$TMPDIR/${species_tag}_norm_SE_bbn.fastq.gz" target=30 ecc=f + + fastp --in1 "\$TMPDIR/${species_tag}_norm_SE_bbn.fastq.gz" \\ + --out1 ${species_tag}_norm_SE.fastq.gz \\ + --thread ${task.cpus} \\ + --cut_front --cut_front_window_size 1 --cut_front_mean_quality 5 \\ + --cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 5 \\ + --cut_right --cut_right_window_size 4 --cut_right_mean_quality 5 \\ + --length_required 25 + + rm -f "\$TMPDIR/${species_tag}_SE.fastq.gz" "\$TMPDIR/${species_tag}_norm_SE_bbn.fastq.gz" + + NREADS=\$(zcat ${species_tag}_norm_SE.fastq.gz 2>/dev/null | awk 'NR%4==1' | wc -l || echo 0) + echo "[INFO] \$NREADS normalized SE reads for ${species_tag}" + + mkdir -p "${launchDir}/rnaseq_reads" + """ + + stub: + """ + : > ${species_tag}_norm_R1.fastq.gz + : > ${species_tag}_norm_R2.fastq.gz + : > ${species_tag}_norm_SE.fastq.gz + echo "[STUB] SRA_FETCH_SE for ${species_tag}" + """ +} diff --git a/modules/local/sra_query.nf b/modules/local/sra_query.nf new file mode 100644 index 0000000..6b6f933 --- /dev/null +++ b/modules/local/sra_query.nf @@ -0,0 +1,64 @@ +// Single-species NCBI SRA query. storeDir caches the per-species CSV so re-runs +// skip the network query. To invalidate, delete rnaseq_reads/sra_query/.sra_query.csv. +// The batched variant (SRA_QUERY_BATCH) is preferred for many species at once. +process SRA_QUERY { + label 'edirect' + label 'process_single' + tag "$species_tag" + + storeDir "${launchDir}/rnaseq_reads/sra_query" + + input: + tuple val(species_tag), val(taxonid) + + output: + tuple val(species_tag), path("${species_tag}.sra_query.csv"), emit: query_result + + script: + """ + set -euo pipefail + + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${species_tag}.sra_query.csv + + esearch -db sra \\ + -query "txid${taxonid}[Organism:noexp] AND RNA-Seq[Strategy] AND PAIRED[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND (BGISEQ[Platform] OR Illumina[Platform])" | \\ + efetch -format runinfo > _runinfo.tmp + + # col 1=Run, col 4=spots, col 13=LibraryStrategy, col 16=LibraryLayout, col 19=Platform + # Prepend a platform rank (0=Illumina, 1=BGI/other) so the top 5 prefer Illumina, + # then by spot count desc; BGI/other only fill remaining slots when Illumina runs out. + awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="PAIRED" && \$1~/^[SDE]RR/ && \$4+0>=250000 {rank=(\$19~/[Ii]llumina/)?0:1; printf "%d,%s,%s,%s\\n", rank, \$1, \$4, \$19}' _runinfo.tmp | \\ + sort -t',' -k1,1n -k3,3rn | \\ + head -n 5 | \\ + while IFS=',' read -r rank acc spots platform; do + printf '%s,%s,%s,%s,%s,PAIRED\\n' "${species_tag}" "${taxonid}" "\$acc" "\$spots" "\$platform" + done >> ${species_tag}.sra_query.csv + + rm -f _runinfo.tmp + NHITS=\$(awk 'END{print NR-1}' ${species_tag}.sra_query.csv) + echo "[INFO] Found \$NHITS paired-end SRA accessions for ${species_tag} (taxonid=${taxonid})" + + # SE fallback: if no PE hits found and enable_single_end is true, query SINGLE layout + if [ "${params.enable_single_end}" = "true" ] && [ "\$NHITS" -eq 0 ]; then + esearch -db sra \\ + -query "txid${taxonid}[Organism:noexp] AND RNA-Seq[Strategy] AND SINGLE[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND Illumina[Platform]" | \\ + efetch -format runinfo > _runinfo_se.tmp + awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="SINGLE" && \$1~/^[SDE]RR/ && \$4+0>=250000 {printf "%s,%s,%s\\n", \$1, \$4, \$19}' _runinfo_se.tmp | \\ + sort -t',' -k2 -rn | \\ + head -n ${params.max_rnaseq_se_runs} | \\ + while IFS=',' read -r acc spots platform; do + printf '%s,%s,%s,%s,%s,SINGLE\\n' "${species_tag}" "${taxonid}" "\$acc" "\$spots" "\$platform" + done >> ${species_tag}.sra_query.csv + rm -f _runinfo_se.tmp + NHITS=\$(awk 'END{print NR-1}' ${species_tag}.sra_query.csv) + echo "[INFO] SE fallback: found \$NHITS single-end accessions for ${species_tag}" + fi + """ + + stub: + """ + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\n' > ${species_tag}.sra_query.csv + printf '%s,%s,SRR000001,1000000,ILLUMINA,PAIRED\n' "${species_tag}" "${taxonid}" >> ${species_tag}.sra_query.csv + echo "[STUB] SRA_QUERY for ${species_tag}" + """ +} diff --git a/modules/local/sra_query_batch.nf b/modules/local/sra_query_batch.nf new file mode 100644 index 0000000..218f6a8 --- /dev/null +++ b/modules/local/sra_query_batch.nf @@ -0,0 +1,108 @@ +// Batched SRA query: handles params.sra_query_batch_size species per SLURM job. +// maxForks 4 caps concurrent jobs to avoid overwhelming NCBI. +// Per-species esearch/efetch is retried up to 3 times inline with exponential +// backoff before writing an empty CSV. Existing per-species CSVs already +// present in rnaseq_reads/sra_query/ are reused without re-querying. +// Resources overridden per-profile by withName: '.*:SRA_QUERY_BATCH' in +// conf/profile_annotate.config (short queue, retry on failure). +process SRA_QUERY_BATCH { + label 'edirect' + label 'process_single' + tag "${species_tags[0]}_+${species_tags.size() - 1}_more" + + publishDir "${launchDir}/rnaseq_reads/sra_query", mode: 'copy', overwrite: false + + maxForks 4 + + input: + tuple val(species_tags), val(taxonids) + + output: + path("*.sra_query.csv"), emit: query_results + + script: + def cache_dir = "${launchDir}/rnaseq_reads/sra_query" + def batch_args = [species_tags, taxonids].transpose() + .collect { st, tid -> "${st}\\t${tid}" } + .join('\\n') + """ + set -uo pipefail + + printf '${batch_args}\\n' > batch_input.tsv + + query_species() { + local stag="\$1" tid="\$2" attempt + + for attempt in 1 2 3; do + rm -f "_runinfo_\${stag}.tmp" + if timeout 120 bash -c \\ + "esearch -db sra -query 'txid\${tid}[Organism:noexp] AND RNA-Seq[Strategy] AND PAIRED[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND (BGISEQ[Platform] OR Illumina[Platform])' | efetch -format runinfo" \\ + < /dev/null > "_runinfo_\${stag}.tmp"; then + return 0 + fi + echo "[WARN] Attempt \${attempt}/3 failed or timed out for \${stag}" + [ "\${attempt}" -lt 3 ] && sleep \$((attempt * 30)) + done + rm -f "_runinfo_\${stag}.tmp" + return 1 + } + + while IFS=\$(printf '\\t') read -r species_tag taxonid; do + cached="${cache_dir}/\${species_tag}.sra_query.csv" + if [ -s "\$cached" ]; then + cp "\$cached" "\${species_tag}.sra_query.csv" + echo "[INFO] Reusing cached result for \${species_tag}" + continue + fi + + if query_species "\${species_tag}" "\${taxonid}"; then + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\\n' > "\${species_tag}.sra_query.csv" + # col 1=Run, col 4=spots, col 13=LibraryStrategy, col 16=LibraryLayout, col 19=Platform + # Prepend a platform rank (0=Illumina, 1=BGI/other) so the top 5 prefer Illumina, + # then by spot count desc; BGI/other only fill remaining slots when Illumina runs out. + awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="PAIRED" && \$1~/^[SDE]RR/ && \$4+0>=250000 {rank=(\$19~/[Ii]llumina/)?0:1; printf "%d,%s,%s,%s\\n", rank, \$1, \$4, \$19}' "_runinfo_\${species_tag}.tmp" | \\ + sort -t',' -k1,1n -k3,3rn | \\ + head -n 5 | \\ + while IFS=',' read -r rank acc spots platform; do + printf '%s,%s,%s,%s,%s,PAIRED\\n' "\${species_tag}" "\${taxonid}" "\$acc" "\$spots" "\$platform" + done >> "\${species_tag}.sra_query.csv" + rm -f "_runinfo_\${species_tag}.tmp" + NHITS=\$(awk 'END{print NR-1}' "\${species_tag}.sra_query.csv") + echo "[INFO] Found \$NHITS paired-end accessions for \${species_tag} (taxonid=\${taxonid})" + # SE fallback: if no PE hits and enable_single_end, query SINGLE layout + if [ "${params.enable_single_end}" = "true" ] && [ "\$NHITS" -eq 0 ]; then + rm -f "_runinfo_se_\${species_tag}.tmp" + if timeout 120 bash -c \\ + "esearch -db sra -query 'txid\${taxonid}[Organism:noexp] AND RNA-Seq[Strategy] AND SINGLE[Layout] AND 00000000075[ReadLength] : 00000000300[ReadLength] AND Illumina[Platform]' | efetch -format runinfo" \\ + < /dev/null > "_runinfo_se_\${species_tag}.tmp"; then + awk -F',' 'NR>1 && \$13=="RNA-Seq" && \$16=="SINGLE" && \$1~/^[SDE]RR/ && \$4+0>=250000 {printf "%s,%s,%s\\n", \$1, \$4, \$19}' "_runinfo_se_\${species_tag}.tmp" | \\ + sort -t',' -k2 -rn | \\ + head -n ${params.max_rnaseq_se_runs} | \\ + while IFS=',' read -r acc spots platform; do + printf '%s,%s,%s,%s,%s,SINGLE\\n' "\${species_tag}" "\${taxonid}" "\$acc" "\$spots" "\$platform" + done >> "\${species_tag}.sra_query.csv" + fi + rm -f "_runinfo_se_\${species_tag}.tmp" + NHITS=\$(awk 'END{print NR-1}' "\${species_tag}.sra_query.csv") + echo "[INFO] SE fallback: \$NHITS single-end accessions for \${species_tag}" + fi + else + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\\n' > "\${species_tag}.sra_query.csv" + echo "[WARN] All 3 attempts failed for \${species_tag}; writing empty CSV" + fi + done < batch_input.tsv + """ + + stub: + def stub_args = [species_tags, taxonids].transpose() + .collect { st, tid -> "${st}\\t${tid}" } + .join('\\n') + """ + printf '${stub_args}\\n' > batch_input.tsv + while IFS=\$(printf '\\t') read -r species_tag taxonid; do + printf 'species_tag,taxonid,sra_accession,spots,platform,layout\\n' > "\${species_tag}.sra_query.csv" + printf '%s,%s,SRR000001,1000000,ILLUMINA,PAIRED\\n' "\${species_tag}" "\${taxonid}" >> "\${species_tag}.sra_query.csv" + done < batch_input.tsv + echo "[STUB] SRA_QUERY_BATCH (${species_tags.size()} species)" + """ +} diff --git a/modules/local/write_empty_reads.nf b/modules/local/write_empty_reads.nf new file mode 100644 index 0000000..8ab0685 --- /dev/null +++ b/modules/local/write_empty_reads.nf @@ -0,0 +1,32 @@ +// Write zero-byte paired FASTQ placeholder files for species with no SRA data. +// Runs in the driver process (no SLURM job) for the no-data branch from SRA branching. +// storeDir caches results so placeholders are not re-created on resume. +process WRITE_EMPTY_READS { + label 'setup' + label 'process_single' + tag "$species_tag" + + storeDir "${launchDir}/rnaseq_reads" + + input: + val(species_tag) + + output: + tuple val(species_tag), path("${species_tag}_norm_R1.fastq.gz"), path("${species_tag}_norm_R2.fastq.gz"), + path("${species_tag}_norm_SE.fastq.gz"), emit: reads + + script: + """ + : > ${species_tag}_norm_R1.fastq.gz + : > ${species_tag}_norm_R2.fastq.gz + : > ${species_tag}_norm_SE.fastq.gz + echo "[INFO] No SRA data for ${species_tag}; created empty read placeholders" + """ + + stub: + """ + : > ${species_tag}_norm_R1.fastq.gz + : > ${species_tag}_norm_R2.fastq.gz + : > ${species_tag}_norm_SE.fastq.gz + """ +} diff --git a/subworkflows/local/fetch_rnaseq.nf b/subworkflows/local/fetch_rnaseq.nf new file mode 100644 index 0000000..d0d73da --- /dev/null +++ b/subworkflows/local/fetch_rnaseq.nf @@ -0,0 +1,133 @@ +/* + * FETCH_RNASEQ — NCBI SRA query + download + normalization + * + * Steps: + * 1. Query SRA per species (SRA_QUERY_BATCH or reuse cached CSVs via skip_sra_query). + * 2. Merge all per-species CSVs into {stem}.rnaseq_sra.csv (COLLECT_SRA_QUERY). + * 3. Route each species to SRA_FETCH (PE), SRA_FETCH_SE (SE_trinity / SINGLE), + * or WRITE_EMPTY_READS (no data). Skipped when stop_after_sra_query=true. + * + * Emits: + * reads — channel: tuple(val(species_tag), path(r1), path(r2), path(se)) + * Empty channel when stop_after_sra_query=true. + * sra_manifest — path to {stem}.rnaseq_sra.csv (always emitted after step 2). + */ + +include { SRA_QUERY_BATCH } from './../../modules/local/sra_query_batch' +include { COLLECT_SRA_QUERY } from './../../modules/local/collect_sra_query' +include { WRITE_EMPTY_READS } from './../../modules/local/write_empty_reads' +include { SRA_FETCH } from './../../modules/local/sra_fetch' +include { SRA_FETCH_SE } from './../../modules/local/sra_fetch_se' + +workflow FETCH_RNASEQ { + + take: + ch_species // channel: tuple(val(species_tag), val(taxonid)) + + main: + // ── Step 1: SRA query ──────────────────────────────────────────────────── + def sra_query_results + if (params.skip_sra_query.toBoolean()) { + // Reuse pre-existing per-species CSVs; skip all NCBI network calls. + sra_query_results = ch_species + .map { species_tag, _taxonid -> + def csv = file("${launchDir}/rnaseq_reads/sra_query/${species_tag}.sra_query.csv") + if (!csv.exists()) { + log.warn "skip_sra_query: no cached CSV for ${species_tag} — skipping this species" + return null + } + tuple(species_tag, csv) + } + .filter { it != null } + } else { + def sra_batched = ch_species + .collate(params.sra_query_batch_size) + .map { batch -> tuple(batch.collect { it[0] }, batch.collect { it[1] }) } + SRA_QUERY_BATCH(sra_batched) + sra_query_results = SRA_QUERY_BATCH.out.query_results + .flatten() + .map { csv -> tuple(csv.baseName.replaceAll(/\.sra_query$/, ''), csv) } + } + + // ── Step 2: collect manifest ───────────────────────────────────────────── + def stem = file(params.samples).baseName + COLLECT_SRA_QUERY( + sra_query_results.map { _stag, csv -> csv }.collect(), + stem + ) + + // ── Step 3: fetch reads (skipped when stop_after_sra_query=true) ───────── + def ch_reads + if (params.stop_after_sra_query.toBoolean()) { + log.info "[FETCH_RNASEQ] stop_after_sra_query=true; skipping SRA download" + ch_reads = Channel.empty() + } else { + // Build the blacklist map once (O(1) per-accession lookup in the closures). + def blPath = file("${launchDir}/rnaseq_blacklist.csv") + def blMap = blPath.exists() + ? blPath.readLines().drop(1) + .findAll { it.trim() && !it.startsWith('#') } + .collectEntries { line -> + def cols = line.split(',') + cols.size() >= 4 ? [(cols[0].trim()): cols[3].trim()] : [:] + } + : [:] + + // csvHasPE: at least one PAIRED accession not blocked or overridden to SE. + def csvHasPE = { csv -> + csv.readLines().drop(1).findAll { it.trim() }.any { line -> + def cols = line.split(',') + if (cols.size() < 3) return false + def layout = cols.size() > 5 ? cols[5].trim() : 'PAIRED' + def action = blMap.get(cols[2].trim(), '') + layout == 'PAIRED' && action != 'skip' && action != 'SE_trinity' + } + } + + // csvHasSEtrinity: at least one PAIRED accession overridden to SE via blacklist. + def csvHasSEtrinity = { csv -> + csv.readLines().drop(1).findAll { it.trim() }.any { line -> + def cols = line.split(',') + cols.size() >= 3 && blMap.get(cols[2].trim(), '') == 'SE_trinity' + } + } + + // csvHasSingleLayout: at least one genuine SINGLE-layout accession (when enable_single_end=true). + def csvHasSingleLayout = { csv -> + csv.readLines().drop(1).findAll { it.trim() }.any { line -> + def cols = line.split(',') + cols.size() > 5 && cols[5].trim() == 'SINGLE' && blMap.get(cols[2].trim(), '') != 'skip' + } + } + + // Three-way branch: + // has_pe → SRA_FETCH (PE wins; SE_trinity entries are handled inside SRA_FETCH) + // has_se → SRA_FETCH_SE (SE_trinity always; SINGLE layout only if enable_single_end) + // no_data → WRITE_EMPTY_READS + def branched_sra = sra_query_results + .branch { + has_pe: csvHasPE.call(it[1]) + has_se: csvHasSEtrinity.call(it[1]) || + (params.enable_single_end.toBoolean() && csvHasSingleLayout.call(it[1])) + no_data: true + } + + SRA_FETCH(branched_sra.has_pe) + SRA_FETCH_SE(branched_sra.has_se) + WRITE_EMPTY_READS(branched_sra.no_data.map { stag, _csv -> stag }) + + // Merge candidate files for human review (not consumed by downstream processes). + SRA_FETCH.out.se_candidates + .collectFile(name: 'rnaseq_se_candidates.csv', storeDir: launchDir, newLine: false) + SRA_FETCH.out.blacklist_candidates + .collectFile(name: 'rnaseq_blacklist_candidates.csv', storeDir: launchDir, newLine: false) + + ch_reads = SRA_FETCH.out.reads + .mix(SRA_FETCH_SE.out.reads) + .mix(WRITE_EMPTY_READS.out.reads) + } + + emit: + reads = ch_reads + sra_manifest = COLLECT_SRA_QUERY.out.manifest +} From 9cfc749a351a16e439fc573ceb1a36afbb4a6c52 Mon Sep 17 00:00:00 2001 From: Jason Stajich Date: Mon, 29 Jun 2026 10:45:03 -0700 Subject: [PATCH 4/5] =?UTF-8?q?feat(#9):=20extract=20RNASEQ=5FPREPARE=20+?= =?UTF-8?q?=20FUNANNOTATE=5FTRAIN=20+=20FUNANNOTATE=5FPREDICT=20=E2=86=92?= =?UTF-8?q?=20TRAIN=5FPREDICT=20subworkflow?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three processes moved out of funannotate.nf monolith: modules/local/rnaseq_prepare.nf — Trinity-GG once per species (storeDir rnaseq_data/) modules/local/funannotate_train.nf — PASA alignment per assembly (training_target//) modules/local/funannotate_predict.nf — gene prediction (target//, Option B persistence) Composed into subworkflows/local/train_predict.nf (TRAIN_PREDICT): take: ch_genomes (meta, genome_fa), ch_reads (species_tag, r1, r2, se) - When run_sra_fetch+reads available: RNASEQ_PREPARE (repr only) → FUNANNOTATE_TRAIN (all) - Species with no reads bypass TRAIN entirely; predict still runs - FUNANNOTATE_PREDICT: skipped when gbk exists and not stale (staleRnaseq logic inlined) emit: metadata (val meta), done (*.predict.done marker) funannotate.nf: ~130 lines of orchestration replaced by TRAIN_PREDICT(predict_genome_ch, reads_ch) staleRnaseq/gbkResult inlined in train_predict.nf as top-level function declarations (same pattern as genomeFile in clean_genomes.nf and mask_genome.nf — these top-level `def` functions from funannotate.nf are not accessible from subworkflow files). withName patterns updated: RNASEQ_PREPARE/FUNANNOTATE_TRAIN/FUNANNOTATE_PREDICT → '.*:RNASEQ_PREPARE' etc. in conf/profile_annotate.config. Stub-run gate: 15 completed, 0 failed (run_sra_fetch=false) Dual-stub gate: 17 completed, 0 failed (run_sra_fetch=true) TRAIN_PREDICT:FUNANNOTATE_PREDICT correctly namespaced in both Co-Authored-By: Claude Sonnet 4.6 Claude-Session: https://claude.ai/code/session_01AbWwmm9v6dJn13sgAJn3Sa --- conf/profile_annotate.config | 6 +- funannotate.nf | 612 +-------------------------- modules/local/funannotate_predict.nf | 153 +++++++ modules/local/funannotate_train.nf | 170 ++++++++ modules/local/rnaseq_prepare.nf | 103 +++++ subworkflows/local/train_predict.nf | 143 +++++++ 6 files changed, 581 insertions(+), 606 deletions(-) create mode 100644 modules/local/funannotate_predict.nf create mode 100644 modules/local/funannotate_train.nf create mode 100644 modules/local/rnaseq_prepare.nf create mode 100644 subworkflows/local/train_predict.nf diff --git a/conf/profile_annotate.config b/conf/profile_annotate.config index 4fbabe3..6d1ebba 100644 --- a/conf/profile_annotate.config +++ b/conf/profile_annotate.config @@ -176,13 +176,13 @@ process { cpus = { task.attempt == 1 ? 24 : 32 } } - withName: 'RNASEQ_PREPARE' { + withName: '.*:RNASEQ_PREPARE' { queue = 'epyc'; cpus = 16 memory = { task.attempt == 1 ? '52 GB' : (48 + (task.attempt - 1) * 48) + ' GB' } time = '120h'; maxRetries = 3 } - withName: 'FUNANNOTATE_TRAIN' { + withName: '.*:FUNANNOTATE_TRAIN' { maxRetries = 3 queue = { task.attempt <= 2 ? 'epyc' : 'highmem' } memory = { task.attempt == 1 ? '48 GB' : task.attempt == 2 ? '192 GB' : '256 GB' } @@ -191,7 +191,7 @@ process { clusterOptions = '-N 1 -n 1 --job-name=nf-TRAIN' } - withName: 'FUNANNOTATE_PREDICT' { + withName: '.*:FUNANNOTATE_PREDICT' { queue = 'epyc' memory = { task.attempt == 1 ? '32 GB' : '64 GB' } time = '72h' diff --git a/funannotate.nf b/funannotate.nf index d9d879a..c4e2dfc 100644 --- a/funannotate.nf +++ b/funannotate.nf @@ -19,492 +19,7 @@ // meta.asmid, meta.species, meta.strain, meta.locustag, meta.busco, // meta.transl_table, meta.taxonid carry payload used inside process scripts. // header_length is NOT in meta — it comes from params.header_length (default 24). -// -// GENOME_CLEAN receives: tuple val(meta), path(genome_gz), val(taxondb) -// → emits: tuple val(meta), path(genome_fa) [storeDir writes input_clean_genomes/.fa.gz] -// SRA_FETCH receives: val(species_tag), val(taxonid) [only when --run_sra_fetch] -// → emits: val(species_tag), path(norm_R1.fastq.gz), path(norm_R2.fastq.gz), path(se) -// RNASEQ_PREPARE receives: tuple val(species_tag), val(meta), val(genome_fa), path(r1), path(r2), path(se) -// → emits: val(species_tag), path(trinity-GG.fasta) [storeDir caches in rnaseq_data/] -// FUNANNOTATE_TRAIN receives: tuple val(meta), val(genome_fa), path(r1), path(r2), path(se), path(trinity_fa) -// → emits: tuple val(meta), val(genome_fa) -// FUNANNOTATE_PREDICT receives: tuple val(meta), val(genome_fa) -// → emits metadata: tuple val(meta) -// Run funannotate train on the representative (first) assembly of each species, then -// archive the Trinity-GG transcripts (normalized reads are in rnaseq_reads) -// reads into rnaseq_data/ so all other strains can skip those expensive steps. -// storeDir skips this process entirely if all five output files already exist. -process RNASEQ_PREPARE { - label 'funannotate' - tag "$species_tag" - - storeDir "${launchDir}/rnaseq_data" - - cpus 16 - memory '96 GB' - time '120h' - - input: - tuple val(species_tag), val(meta), val(genome_fa), path(r1), path(r2), path(se) - - output: - tuple val(species_tag), - path("${species_tag}.trinity-GG.fasta"), emit: shared - - script: - def out = meta.id - def species = meta.species - def strain = meta.strain - def header_length = params.header_length - """ - # ── Empty-reads sentinel: no RNA-seq found by SRA_FETCH / SRA_FETCH_SE ── - if [ ! -s "${r1}" ] && [ ! -s "${se}" ]; then - echo "[INFO] No RNAseq reads for ${species_tag}; writing empty shared markers" - touch ${species_tag}.trinity-GG.fasta - exit 0 - fi - - # ── If representative was already trained, just extract shared files ────── - TRAIN_GFF3="${params.training_target}/${out}/training/funannotate_train.pasa.gff3" - if [ -f "\$TRAIN_GFF3" ]; then - echo "[INFO] Training already complete for ${out}; extracting shared files to rnaseq_data" - TRAINDIR="${params.training_target}/${out}/training" - TRINITY_FA=\$(find \$TRAINDIR -maxdepth 1 -name "trinity.fasta" | head -1) - if [ -n "\$TRINITY_FA" ]; then - cp "\$TRINITY_FA" ${species_tag}.trinity-GG.fasta - else - touch ${species_tag}.trinity-GG.fasta - fi - exit 0 - fi - - source /etc/profile.d/modules.sh 2>/dev/null || true - - export AUGUSTUS_CONFIG_PATH=${params.augustus_config} - export FUNANNOTATE_DB=${params.funannotate_db} - TMPDIR=\${SCRATCH:-/tmp} - - # ── Run full funannotate train on the representative genome ─────────────── - # Use SCRATCH for the funannotate output dir so Trinity/HISAT2/normalize - # intermediates land on fast local storage and don't consume project quota. - echo "[INFO] RNASEQ_PREPARE: running funannotate train for representative ${out} (species: ${species_tag})" - - # Inflate a gzipped clean genome to a local uncompressed copy; funannotate cannot - # read a gzipped FASTA via -i. Plain (uncompressed) genomes pass through unchanged. - GENOME_FA="${genome_fa}" - case "\$GENOME_FA" in - *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; - *) GENOME_IN="\$GENOME_FA" ;; - esac - - if [ -s "${r1}" ]; then - funannotate train -i "\$GENOME_IN" -o \$SCRATCH/${out} \\ - --left_norm ${r1} --right_norm ${r2} --aligners minimap2 \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --jaccard_clip --no-progress --min_coverage 4 \\ - --max_intronlen ${params.max_intronlen} \\ - --stop_after_trinity --no_trimmomatic - else - echo "[INFO] RNASEQ_PREPARE: using single-end reads for ${out}" - funannotate train -i "\$GENOME_IN" -o \$SCRATCH/${out} \\ - --single_norm ${se} --aligners minimap2 \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --no-progress --min_coverage 4 \\ - --max_intronlen ${params.max_intronlen} \\ - --stop_after_trinity --no_trimmomatic - fi - - # ── Copy shared outputs to rnaseq_data/ ────────────────────────────────── - TRAINDIR="\$SCRATCH/${out}/training" - TRINITY_FA=\$(find \$TRAINDIR -maxdepth 1 -name "trinity.fasta" | head -1) - if [ -n "\$TRINITY_FA" ]; then - cp "\$TRINITY_FA" ${species_tag}.trinity-GG.fasta - else - echo "[WARN] No trinity.fasta found under \$TRAINDIR for ${out}" - touch ${species_tag}.trinity-GG.fasta - fi - - # ── Clean up scratch output dir (all intermediates were temporary) ──────── - rm -rf "\$SCRATCH/${out}" - echo "[INFO] RNASEQ_PREPARE complete for ${species_tag}" - """ - - stub: - def out = meta.id - """ - echo ">stub_trinity_${species_tag}" > ${species_tag}.trinity-GG.fasta - mkdir -p ${params.training_target}/${out}/training - touch ${params.training_target}/${out}/training/funannotate_train.pasa.gff3 - """ -} - -// For non-representative strains: funannotate train --trinity runs only -// PASA (skips Trimmomatic, normalization, HISAT2, and Trinity-GG assembly). -// Falls back to a full train when no shared Trinity is available (e.g. species with -// a single strain or when run_sra_fetch is false). -process FUNANNOTATE_TRAIN { - label 'funannotate' - tag "${meta.id}" - - cpus 16 - memory '96 GB' - time '120h' - - input: - tuple val(meta), val(genome_fa), path(r1), path(r2), path(se), path(trinity_fa) - - output: - tuple val(meta), val(genome_fa) - - script: - def out = meta.id - def asmid = meta.asmid - def species = meta.species - def strain = meta.strain - def locustag = meta.locustag - def header_length = params.header_length - def pasa_db_arg = "--pasa_db sqlite" - """ - # ── Skip if no RNA-seq data at all ──────────────────────────────────────── - if [ ! -s "${r1}" ] && [ ! -s "${se}" ] && [ ! -s "${trinity_fa}" ]; then - echo "[INFO] No RNAseq data for ${out}, skipping funannotate train" - exit 0 - fi - - # ── Skip if training output already present and rnaseq is not newer than GBK ── - # Accept a compressed prediction (.gbk.gz) as "done" so folders can be space-saved. - TRAIN_GFF3="${params.training_target}/${out}/training/funannotate_train.pasa.gff3" - PREDICT_GBK="${params.target}/${out}/predict_results/${out}.gbk" - [ -f "\$PREDICT_GBK" ] || PREDICT_GBK="${params.target}/${out}/predict_results/${out}.gbk.gz" - if [ -f "\$TRAIN_GFF3" ]; then - RETRAIN=0 - if [ -f "\$PREDICT_GBK" ]; then - # Re-train if the rnaseq reads are newer than the existing prediction GBK. - if [ -s "${r1}" ] && [ "${r1}" -nt "\$PREDICT_GBK" ]; then - echo "[INFO] RNAseq R1 reads newer than predict GBK for ${out}; retraining" -# rm -rf "${params.training_target}/${out}/training" - RETRAIN=1 - elif [ -s "${se}" ] && [ "${se}" -nt "\$PREDICT_GBK" ]; then - echo "[INFO] RNAseq SE reads newer than predict GBK for ${out}; retraining" -# rm -rf "${params.training_target}/${out}/training" - RETRAIN=1 - fi - fi - if [ \$RETRAIN -eq 0 ]; then - echo "[INFO] Training already complete for ${out}; skipping" - exit 0 - fi - fi - - source /etc/profile.d/modules.sh 2>/dev/null || true - module load funannotate/dev-1.8.18 - export AUGUSTUS_CONFIG_PATH=${params.augustus_config} - export FUNANNOTATE_DB=${params.funannotate_db} - TMPDIR=\${SCRATCH:-/tmp} - export PASACONF="" - pasa_db_arg="--pasa_db sqlite" - # ── Optional per-task MariaDB for PASA ──────────────────────────────────── - if [ "${params.pasa_mysql}" = "true" ]; then - MYSQL_SCRATCH=${params.training_target}/${out}/training/mysql_db - if [ ! -f \$MYSQL_SCRATCH/mysql/conf/my.cnf ]; then - echo "[INFO] Setting up temporary MariaDB for PASA at \$MYSQL_SCRATCH" - mkdir -p \$MYSQL_SCRATCH/db \$MYSQL_SCRATCH/conf - rsync -a ${params.mysql_datadir}/mysql \$MYSQL_SCRATCH/db/ || \ - { echo "ERROR: Failed to copy mysql data from ${params.mysql_datadir}" >&2; exit 1; } - cp ${params.pasa_conf_dir}/my.cnf \$MYSQL_SCRATCH/conf/my.cnf || \ - { echo "ERROR: Failed to copy my.cnf" >&2; exit 1; } - fi - MYHOSTNAME=\$(hostname -s) - PORT=\$(shuf -i3000-4999 -n1) - export PASACONF=\$MYSQL_SCRATCH/conf/pasa-local-\${MYHOSTNAME}.config.txt - cp ${params.pasa_conf_dir}/conf.txt \$PASACONF - sed -i "s/^MYSQLSERVER.*\$/MYSQLSERVER=\${MYHOSTNAME}:\${PORT}/" \$PASACONF - perl -i -p -e "s/port = \\d+/port = \${PORT}/" \$MYSQL_SCRATCH/conf/my.cnf - # ── may be unnecessary if overridden by -B option later? ── - export SINGULARITY_BINDPATH=\$TMPDIR,\$MYSQL_SCRATCH/mysql_db - stop_mysqldb() { singularity instance stop mysqldb_${asmid} 2>/dev/null || true; } - trap "stop_mysqldb; exit 130" SIGHUP SIGINT SIGTERM - trap "stop_mysqldb" EXIT - module load singularity - singularity instance start --writable-tmpfs \\ - -B \$MYSQL_SCRATCH/conf/my.cnf:/etc/mysql/my.cnf,\$MYSQL_SCRATCH/db/:/var/lib/mysql,\$MYSQL_SCRATCH/conf:/usr/conf \\ - ${params.mariadb_sif} mysqldb_${asmid} /usr/bin/mysqld_safe - pasa_db_arg="--pasa_db mysql" - sleep 5 - fi - - # Inflate a gzipped clean genome to a local uncompressed copy; funannotate cannot - # read a gzipped FASTA via -i. Plain (uncompressed) genomes pass through unchanged. - GENOME_FA="${genome_fa}" - case "\$GENOME_FA" in - *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; - *) GENOME_IN="\$GENOME_FA" ;; - esac - - # ── Use shared Trinity transcripts (PASA only) or run full train ────────── - if [ -s "${trinity_fa}" ]; then - if [ -s "${r1}" ]; then - echo "[INFO] Running funannotate train (PASA+PE) for ${out} using shared Trinity" - funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ - --trinity ${trinity_fa} --left_norm ${r1} --right_norm ${r2} \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --jaccard_clip --no-progress \\ - --max_intronlen ${params.max_intronlen} \\ - \$pasa_db_arg - elif [ -s "${se}" ]; then - echo "[INFO] Running funannotate train (PASA+SE) for ${out} using shared Trinity" - funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ - --trinity ${trinity_fa} --single_norm ${se} \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --no-progress \\ - --max_intronlen ${params.max_intronlen} \\ - \$pasa_db_arg - else - echo "[INFO] Running funannotate train (PASA only, no reads) for ${out} using shared Trinity" - funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ - --trinity ${trinity_fa} --left_norm ${r1} --right_norm ${r2} \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --jaccard_clip --no-progress \\ - --max_intronlen ${params.max_intronlen} \\ - \$pasa_db_arg - fi - elif [ -s "${r1}" ]; then - echo "[INFO] Running funannotate train (full PE, no shared Trinity) for ${out}" - funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ - --left_norm ${r1} --right_norm ${r2} --aligners minimap2 \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --jaccard_clip --no-progress --min_coverage 4 \\ - --max_intronlen ${params.max_intronlen} \\ - \$pasa_db_arg - else - echo "[INFO] Running funannotate train (full SE, no shared Trinity) for ${out}" - funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ - --single_norm ${se} --aligners minimap2 \\ - --species "${species}" --strain "${strain}" \\ - --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ - --header_length ${header_length} \\ - --no-progress --min_coverage 4 \\ - --max_intronlen ${params.max_intronlen} \\ - \$pasa_db_arg - fi - - # ── Remove large intermediates not needed for predict or update ───────────── - # Keeps: *.bam, *.bai, *.pasa.gff3, *.stringtie.gtf, *.transcripts.gff3 - TRAINDIR="${params.training_target}/${out}/training" - echo "[INFO] Removing large training intermediates in \$TRAINDIR" - rm -rf "\$TRAINDIR/hisat2" - rm -rf "\$TRAINDIR/trinity_gg" - echo "[INFO] Training cleanup complete for ${out}" - echo "mysql is ${params.pasa_mysql}" - if [ "${params.pasa_mysql}" = "true" ]; then stop_mysqldb; fi - echo "[INFO] stopped mysql" - """ - - stub: - def out = meta.id - """ - echo "[STUB] FUNANNOTATE_TRAIN stub for ${out}" - mkdir -p ${params.training_target}/${out}/training - touch ${params.training_target}/${out}/training/funannotate_train.pasa.gff3 - """ -} - -// Option B persistence model: funannotate predict computes DIRECTLY into the persistent -// per-genome dir (${params.target}/${out}), symmetric with FUNANNOTATE_TRAIN writing to -// training_target. funannotate checkpoints into predict_misc/, so a restart after an -// OOM/timeout/orchestrator death resumes completed steps in place rather than starting -// over. There is no publishDir copy and no work-dir<->target rsync: the durable output is -// written where downstream steps already read it. Large intermediates still go to the -// node-local --tmpdir. The Nextflow output is a small marker file (nothing consumes the -// predict dir as a channel; downstream rebuilds metadata from the CSV and gates on the -// on-disk GBK), so emitting a marker keeps the DAG edge without copying the result tree. -process FUNANNOTATE_PREDICT { - label 'funannotate' - tag "${meta.id}" - - cpus 16 - memory '32 GB' - time '32h' - - input: - tuple val(meta), val(genome_fa) - - output: - val meta, emit: metadata - path("${meta.id}.predict.done"), emit: done - - script: - def out = meta.id - def asmid = meta.asmid - def species = meta.species - def strain = meta.strain - def locustag = meta.locustag - def busco_lineage = meta.busco - def header_length = params.header_length - def transl_table = meta.transl_table - """ - source /etc/profile.d/modules.sh 2>/dev/null || true - module load funannotate/dev-1.8.18 - export AUGUSTUS_CONFIG_PATH=${params.augustus_config} - export FUNANNOTATE_DB=${params.funannotate_db} - TMPDIR=\${SCRATCH:-/tmp} - - PREDICTDIR="${params.target}/${out}" - PREDICT_GBK="\$PREDICTDIR/predict_results/${out}.gbk" - - if [ "${params.debug.toBoolean()}" = "true" ]; then - echo "[DEBUG] out=${out} asmid=${asmid} species=${species} strain=${strain}" - echo "[DEBUG] locustag=${locustag} busco=${busco_lineage} transl_table=${transl_table}" - echo "[DEBUG] proteins=${params.proteins} genome_fa=${genome_fa}" - echo "[DEBUG] PREDICTDIR=\$PREDICTDIR TMPDIR=\$TMPDIR pwd=\$(pwd)" - fi - - # ── Skip vs. refresh decision ───────────────────────────────────────────── - # The workflow schedules this process when the GBK is missing OR stale (rnaseq/trinity - # newer than the GBK, per staleRnaseq()). Re-derive staleness here from the same on-disk - # timestamps so a current GBK short-circuits, but a stale one forces a clean re-predict. - if [ -s "\$PREDICT_GBK" ]; then - SPECIES_TAG=\$(printf '%s' "${species}" | sed -E 's/[[:space:]]+/_/g') - STALE=0 - for f in "${launchDir}/rnaseq_reads/\${SPECIES_TAG}_norm_R1.fastq.gz" \\ - "${launchDir}/rnaseq_reads/\${SPECIES_TAG}_norm_SE.fastq.gz" \\ - "${launchDir}/rnaseq_data/\${SPECIES_TAG}.trinity-GG.fasta"; do - if [ -s "\$f" ] && [ "\$f" -nt "\$PREDICT_GBK" ]; then STALE=1; fi - done - if [ "\$STALE" -eq 0 ]; then - echo "[INFO] Prediction already complete and current for ${out}; nothing to do" - touch ${out}.predict.done - exit 0 - fi - echo "[INFO] Stale prediction for ${out}: rnaseq/trinity newer than GBK — clearing predict outputs for a fresh run" - rm -rf "\$PREDICTDIR/predict_results" "\$PREDICTDIR/predict_misc" - fi - - mkdir -p "\$PREDICTDIR" - - # ── Guard against a corrupt partial from a previous attempt ─────────────── - # funannotate resumes from predict_misc/. If predict_results/ exists without a - # predict_misc/ (a half-written tree with no checkpoints and no GBK), clear it so - # predict starts the consensus/output step from a clean state instead of choking on it. - if [ ! -d "\$PREDICTDIR/predict_misc" ] && [ -d "\$PREDICTDIR/predict_results" ]; then - echo "[WARN] predict_results/ present without predict_misc/ for ${out}; clearing stale partial" - rm -rf "\$PREDICTDIR/predict_results" - fi - - # funannotate predict expects training data at /training; point it at the - # persistent training dir. The symlink lives in the persistent project tree (no - # publishDir to recursively copy the target), so it is left in place. - if [ -d "${params.training_target}/${out}/training" ]; then - ln -sfn "${params.training_target}/${out}/training" "\$PREDICTDIR/training" - fi - - TBL2ASN_PARAMS="-l paired-ends" - - # Inflate a gzipped clean/masked genome to a local uncompressed copy; funannotate - # cannot read a gzipped FASTA via -i, and the pre-flight awk below also needs plain - # text. Plain (uncompressed) genomes pass through unchanged. - GENOME_FA="${genome_fa}" - case "\$GENOME_FA" in - *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; - *) GENOME_IN="\$GENOME_FA" ;; - esac - - # ── Too-small-genome pre-flight guard ──────────────────────────────────── - # Assemblies that are both small AND fragmented cannot yield funannotate's - # required 30 training models; predict would run for hours then abort with - # "Not enough gene models N to train Augustus (30 required), exiting". Detect - # that up front from cheap contig stats and skip cleanly (flag, no crash). - # Requires BOTH gates so complete small genomes (e.g. Malassezia) are unaffected. - # Disabled when predict_min_asm_bp=0. - SKIP_REPORT="${params.target}/predict_skipped_too_small.tsv" - if [ "${params.predict_min_asm_bp}" -gt 0 ]; then - # Per-contig lengths -> sort descending -> N50 (portable; no gawk asort). - read ASM_BP ASM_CTG ASM_N50 < <( - awk '/^>/{if(len)print len;len=0;next}{len+=length(\$0)}END{if(len)print len}' "\$GENOME_IN" \\ - | sort -rn \\ - | awk '{L[NR]=\$1;tot+=\$1}END{half=tot/2;run=0;n50=0;for(i=1;i<=NR;i++){run+=L[i];if(run>=half){n50=L[i];break}}print tot, NR, n50}') - echo "[INFO] Pre-flight assembly stats for ${out}: \${ASM_BP} bp, \${ASM_CTG} contigs, N50 \${ASM_N50}" - SMALL=0; FRAG=0 - [ "\$ASM_BP" -lt "${params.predict_min_asm_bp}" ] && SMALL=1 - { [ "\$ASM_N50" -lt "${params.predict_frag_max_n50}" ] || [ "\$ASM_CTG" -gt "${params.predict_frag_max_contigs}" ]; } && FRAG=1 - if [ "\$SMALL" -eq 1 ] && [ "\$FRAG" -eq 1 ]; then - echo "[WARN] ${out} is too small/fragmented for funannotate training (\${ASM_BP} bp, \${ASM_CTG} contigs, N50 \${ASM_N50}); skipping predict" >&2 - mkdir -p "${params.target}" - [ -s "\$SKIP_REPORT" ] || printf 'out\tasmid\tlocustag\treason\ttotal_bp\tcontigs\tN50\n' > "\$SKIP_REPORT" - printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "${out}" "${asmid}" "${locustag}" "preflight_small_fragmented" "\$ASM_BP" "\$ASM_CTG" "\$ASM_N50" >> "\$SKIP_REPORT" - touch "\$PREDICTDIR/${out}.predict.skipped_too_small" - touch ${out}.predict.done - exit 0 - fi - fi - - funannotate predict --name ${locustag} -i "\$GENOME_IN" --strain "${strain}" \\ - -o "\$PREDICTDIR" -s "${species}" --cpu ${task.cpus} --busco_db ${busco_lineage} \\ - --AUGUSTUS_CONFIG_PATH \$AUGUSTUS_CONFIG_PATH -w codingquarry:0 glimmerhmm:0 \\ - --min_training_models 30 --tmpdir \$TMPDIR --SeqCenter ${params.seqcenter} \\ - --keep_no_stops --header_length ${header_length} --protein_evidence ${params.proteins} \\ - --max_intronlen ${params.max_intronlen} --min_intronlen ${params.min_intronlen} \\ - --tbl2asn "\$TBL2ASN_PARAMS" --table ${transl_table} --auto-skip-genemark || true - - # ── Post-predict catch ──────────────────────────────────────────────────── - # If predict produced no GBK, distinguish the known "too few training models" - # outcome (an unfixable property of the assembly) from a genuine error. The - # former is flagged and skipped so it does not abort the batch; anything else - # still hard-fails so real problems surface. - if [ ! -s "\$PREDICT_GBK" ]; then - PLOG="\$PREDICTDIR/logfiles/funannotate-predict.log" - if [ -f "\$PLOG" ] && grep -q "Not enough gene models .* to train Augustus" "\$PLOG"; then - NMODELS=\$(grep -oE "Not enough gene models [0-9]+" "\$PLOG" | grep -oE "[0-9]+" | tail -1) - echo "[WARN] ${out}: funannotate found only \${NMODELS:-&2 - mkdir -p "${params.target}" - [ -s "\$SKIP_REPORT" ] || printf 'out\tasmid\tlocustag\treason\ttotal_bp\tcontigs\tN50\n' > "\$SKIP_REPORT" - printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "${out}" "${asmid}" "${locustag}" "funannotate_too_few_models:\${NMODELS:-NA}" "" "" "" >> "\$SKIP_REPORT" - touch "\$PREDICTDIR/${out}.predict.skipped_too_small" - touch ${out}.predict.done - exit 0 - fi - echo "ERROR: funannotate predict did not produce expected GBK: \$PREDICT_GBK" >&2 - exit 1 - fi - if [ -d "\$PREDICTDIR/predict_misc/ab_initio_parameters" ]; then - mv "\$PREDICTDIR/predict_misc/ab_initio_parameters" "\$PREDICTDIR" - mv "\$PREDICTDIR/predict_misc/trnascan.no-overlaps.gff3" "\$PREDICTDIR" - rm -rf "\$PREDICTDIR/predict_misc" - mkdir -p "\$PREDICTDIR/predict_misc" - mv "\$PREDICTDIR/ab_initio_parameters" "\$PREDICTDIR/trnascan.no-overlaps.gff3" "\$PREDICTDIR/predict_misc" - fi - find "\$PREDICTDIR/predict_results/" -maxdepth 1 \\( -name "*.txt" -o -name "*.mrna-transcripts.fa" \\) -print0 \ - | xargs -0 --no-run-if-empty pigz - sync - touch ${out}.predict.done - echo "[INFO] Prediction complete for ${out} at \$PREDICTDIR" - """ - - stub: - def out = meta.id - """ - echo "[STUB] Would run funannotate predict for ${out} using ${genome_fa}" - [ -f "${genome_fa}" ] || [ -f "${genome_fa}.gz" ] || { echo "ERROR: genome not found at ${genome_fa}[.gz]" >&2; exit 1; } - mkdir -p ${params.target}/${out}/predict_results ${params.target}/${out}/predict_misc - # non-empty so downstream size>0 gating (predict_ch / postpredict) is exercised - echo "LOCUS stub_${out}" > ${params.target}/${out}/predict_results/${out}.gbk - echo ">stub_${out}_p1" > ${params.target}/${out}/predict_results/${out}.proteins.fa - touch ${out}.predict.done - """ -} process FUNANNOTATE_ANNOTATE { label 'funannotate' @@ -708,6 +223,7 @@ include { SETUP_DBS } from './subworkflows/local/setup_dbs' include { CLEAN_GENOMES } from './subworkflows/local/clean_genomes' include { MASK_GENOME } from './subworkflows/local/mask_genome' include { FETCH_RNASEQ } from './subworkflows/local/fetch_rnaseq' +include { TRAIN_PREDICT } from './subworkflows/local/train_predict' include { ANTISMASH_RUN } from './modules/local/antismash_run' include { INTERPROSCAN_RUN } from './modules/local/interproscan_run' include { SIGNALP_RUN } from './modules/local/signalp_run' @@ -787,14 +303,11 @@ workflow { .combine(SETUP_DBS.out.config) .map { row -> row[0..-3] } - // FUNANNOTATE_PREDICT input tuple drops taxonid (not needed after masking/clean). - // When SRA is enabled: FETCH_RNASEQ fetches reads once per species; RNASEQ_PREPARE runs - // funannotate train on the representative assembly and archives Trinity-GG transcripts - // to rnaseq_data/; all other strains run FUNANNOTATE_TRAIN --trinity. - def predict_input_ch + // SRA read fetching (FETCH_RNASEQ) + training (RNASEQ_PREPARE + FUNANNOTATE_TRAIN) + // + prediction (FUNANNOTATE_PREDICT). All three are composed in TRAIN_PREDICT. + // When run_sra_fetch=false, reads channel is empty and genomes go straight to predict. def reads_ch = Channel.empty() if (params.run_sra_fetch.toBoolean()) { - // Build per-species input: group assemblies, keep first taxonid per species. def sra_input = predict_genome_ch .map { meta, _genome_fa -> def species_tag = meta.species.replaceAll(/\s+/, '_') @@ -805,131 +318,24 @@ workflow { FETCH_RNASEQ(sra_input) reads_ch = FETCH_RNASEQ.out.reads - - if (!params.stop_after_sra_fetch.toBoolean() && !params.stop_after_sra_query.toBoolean()) { - // Build per-assembly channel keyed by species_tag with SRA reads joined. - // reads_ch is now a 4-tuple: (species_tag, r1, r2, se) - def assembly_with_reads = predict_genome_ch - .map { meta, genome_fa -> - def species_tag = meta.species.replaceAll(/\s+/, '_') - tuple(species_tag, meta, genome_fa) - } - .combine(reads_ch, by: 0) - // assembly_with_reads: (species_tag, meta, genome_fa, r1, r2, se) - - // RNASEQ_PREPARE: run funannotate train --stop_after_trinity once per species on - // the representative (first) assembly, then cache the Trinity-GG FASTA in rnaseq_data/ - // so all other strains share it. Normalized reads stay in rnaseq_reads/ (SRA_FETCH storeDir). - // pasa.gff3 is NOT produced here (--stop_after_trinity stops before PASA); - // it is produced by FUNANNOTATE_TRAIN for every strain including the representative. - // Species whose representative r1 and se are both zero-length skip RNASEQ_PREPARE - // entirely; an empty trinity FASTA is written locally without submitting a SLURM job. - def repr_ch = assembly_with_reads - .groupTuple(by: 0) - .map { species_tag, metas, genomes, r1s, r2s, ses -> - tuple(species_tag, metas[0], genomes[0], r1s[0], r2s[0], ses[0]) - } - - def repr_branched = repr_ch.branch { - has_reads: it[3].size() > 0 || it[5].size() > 0 // r1=[3] or se=[5] - no_reads: true - } - - RNASEQ_PREPARE(repr_branched.has_reads) - - // For species with no RNA-seq reads, write an empty trinity FASTA to rnaseq_data/ - // in the driver process (no SLURM job) and emit it directly as a shared channel item. - def empty_shared_ch = repr_branched.no_reads - .map { species_tag, _meta, _gfa, _r1, _r2, _se -> - def empty_fa = file("${launchDir}/rnaseq_data/${species_tag}.trinity-GG.fasta") - if (!empty_fa.exists()) { - empty_fa.parent.mkdirs() - empty_fa.text = '' - } - tuple(species_tag, empty_fa) - } - - def shared_ch = RNASEQ_PREPARE.out.shared.mix(empty_shared_ch) - - // Join shared Trinity from rnaseq_data back to every assembly for FUNANNOTATE_TRAIN. - // Normalized reads (r1/r2/se) come from SRA_FETCH/SRA_FETCH_SE via assembly_with_reads. - def train_input = assembly_with_reads - .combine(shared_ch, by: 0) - .map { species_tag, meta, genome_fa, r1, r2, se, trinity_fa -> - tuple(meta, genome_fa, r1, r2, se, trinity_fa) - } - // train_input: meta=0, genome_fa=1, r1=2, r2=3, se=4, trinity_fa=5 - - // Branch on r1 (idx 2), se (idx 4), or trinity_fa (idx 5) sizes. - // Assemblies with no RNA-seq bypass FUNANNOTATE_TRAIN entirely. - def branched = train_input.branch { - has_rnaseq: it[2].size() > 0 || it[4].size() > 0 || it[5].size() > 0 - no_rnaseq: true - } - def predict_no_rnaseq = branched.no_rnaseq - .map { meta, genome_fa, _r1, _r2, _se, _tf -> - tuple(meta, genome_fa) - } - - // Skip TRAIN at the channel level when pasa.gff3 already exists and is non-empty, - // UNLESS the rnaseq reads or trinity FASTA is newer than the existing prediction GBK - // (staleRnaseq), in which case we re-run training so predict can be refreshed too. - def train_todo = branched.has_rnaseq.filter { meta, _gfa, _r1, _r2, _se, _tf -> - def gff3 = file("${params.training_target}/${meta.id}/training/funannotate_train.pasa.gff3") - !gff3.exists() || gff3.size() == 0 || staleRnaseq(meta.id as String, meta.species as String) - } - def train_done = branched.has_rnaseq - .filter { meta, _gfa, _r1, _r2, _se, _tf -> - def gff3 = file("${params.training_target}/${meta.id}/training/funannotate_train.pasa.gff3") - gff3.exists() && gff3.size() > 0 && !staleRnaseq(meta.id as String, meta.species as String) - } - .map { meta, genome_fa, _r1, _r2, _se, _tf -> - tuple(meta, genome_fa) - } - FUNANNOTATE_TRAIN(train_todo) - predict_input_ch = FUNANNOTATE_TRAIN.out.mix(train_done).mix(predict_no_rnaseq) - } // end if (!stop_after_sra_fetch && !stop_after_sra_query) - } else { - predict_input_ch = predict_genome_ch } if (!params.run_sra_fetch.toBoolean() || (!params.stop_after_sra_fetch.toBoolean() && !params.stop_after_sra_query.toBoolean())) { - def predict_ch = predict_input_ch - .filter { meta, _gfa -> - gbkResult("${params.target}/${meta.id}/predict_results", meta.id as String) == null || staleRnaseq(meta.id as String, meta.species as String) - } - FUNANNOTATE_PREDICT(predict_ch) + TRAIN_PREDICT(predict_genome_ch, reads_ch) // ── Post-predict steps and annotation ──────────────────────────────────── // postpredict: all samples with a completed predict_results/*.gbk, whether // produced in this run or a prior one. This is the source for all optional // pre-annotate steps and for FUNANNOTATE_ANNOTATE itself. - // INPUT_CHECK.out.samples already has taxon/asmid/suppress/n_test filters applied; - // we just add the predict-results existence check on top. def postpredict = INPUT_CHECK.out.samples - // Only genomes whose prediction was already complete AND current in a PRIOR run. - // This is the exact logical complement of the predict_ch filter, so this set is - // disjoint from the genomes (re)predicted in THIS run (which arrive via - // FUNANNOTATE_PREDICT.out.metadata below). Keeping them disjoint means no genome - // is fed downstream twice and stale genomes correctly wait for the fresh predict. + // Only genomes predicted in a PRIOR run — complement of what TRAIN_PREDICT runs now. .filter { meta -> gbkResult("${params.target}/${meta.id}/predict_results", meta.id as String) != null && !staleRnaseq(meta.id as String, meta.species as String) } - // annotate_ready_ch threads through optional pre-annotate steps. Each optional - // step splits the channel into "needs to run" vs "already done", processes the - // former, then mixes the freshly-completed items back. FUNANNOTATE_ANNOTATE only - // fires once all requested optional steps are complete for a given sample. - // Joining ANTISMASH/INTERPRO/SIGNALP output back through predict_meta reconstructs - // the metadata tuple while encoding the dependency edge in the channel DAG. - // - // Same-run completion gate: genomes predicted in THIS run flow in via - // FUNANNOTATE_PREDICT.out.metadata (a real channel edge, so downstream waits for - // predict to finish), while prior-run genomes flow in via postpredict (available - // immediately). The two sets are disjoint by the filters above, so a plain mix - // needs no dedup. (The optional steps below are still each gated behind their - // params — run_antismash/interpro/signalp/update/annotate, all default false.) - def predict_meta = postpredict.mix(FUNANNOTATE_PREDICT.out.metadata) + // Same-run completion gate: new predictions arrive via TRAIN_PREDICT.out.metadata; + // prior-run completions arrive via postpredict (available immediately). + def predict_meta = postpredict.mix(TRAIN_PREDICT.out.metadata) def annotate_ready_ch = predict_meta if (params.run_antismash.toBoolean()) { diff --git a/modules/local/funannotate_predict.nf b/modules/local/funannotate_predict.nf new file mode 100644 index 0000000..d6cb639 --- /dev/null +++ b/modules/local/funannotate_predict.nf @@ -0,0 +1,153 @@ +// Run funannotate predict for one assembly. Writes directly into params.target// +// (Option B persistence: no publishDir copy). Emits a small marker file to carry the +// DAG edge without transferring the full predict tree through Nextflow's work/ directory. +// Resources overridden by withName: '.*:FUNANNOTATE_PREDICT' in conf/profile_annotate.config. +process FUNANNOTATE_PREDICT { + label 'funannotate' + label 'process_high' + tag "${meta.id}" + + input: + tuple val(meta), val(genome_fa) + + output: + val meta, emit: metadata + path("${meta.id}.predict.done"), emit: done + + script: + def out = meta.id + def asmid = meta.asmid + def species = meta.species + def strain = meta.strain + def locustag = meta.locustag + def busco_lineage = meta.busco + def header_length = params.header_length + def transl_table = meta.transl_table + """ + source /etc/profile.d/modules.sh 2>/dev/null || true + module load funannotate/dev-1.8.18 + export AUGUSTUS_CONFIG_PATH=${params.augustus_config} + export FUNANNOTATE_DB=${params.funannotate_db} + TMPDIR=\${SCRATCH:-/tmp} + + PREDICTDIR="${params.target}/${out}" + PREDICT_GBK="\$PREDICTDIR/predict_results/${out}.gbk" + + if [ "${params.debug.toBoolean()}" = "true" ]; then + echo "[DEBUG] out=${out} asmid=${asmid} species=${species} strain=${strain}" + echo "[DEBUG] locustag=${locustag} busco=${busco_lineage} transl_table=${transl_table}" + echo "[DEBUG] proteins=${params.proteins} genome_fa=${genome_fa}" + echo "[DEBUG] PREDICTDIR=\$PREDICTDIR TMPDIR=\$TMPDIR pwd=\$(pwd)" + fi + + # ── Skip vs. refresh decision ───────────────────────────────────────────── + if [ -s "\$PREDICT_GBK" ]; then + SPECIES_TAG=\$(printf '%s' "${species}" | sed -E 's/[[:space:]]+/_/g') + STALE=0 + for f in "${launchDir}/rnaseq_reads/\${SPECIES_TAG}_norm_R1.fastq.gz" \\ + "${launchDir}/rnaseq_reads/\${SPECIES_TAG}_norm_SE.fastq.gz" \\ + "${launchDir}/rnaseq_data/\${SPECIES_TAG}.trinity-GG.fasta"; do + if [ -s "\$f" ] && [ "\$f" -nt "\$PREDICT_GBK" ]; then STALE=1; fi + done + if [ "\$STALE" -eq 0 ]; then + echo "[INFO] Prediction already complete and current for ${out}; nothing to do" + touch ${out}.predict.done + exit 0 + fi + echo "[INFO] Stale prediction for ${out}: rnaseq/trinity newer than GBK — clearing predict outputs for a fresh run" + rm -rf "\$PREDICTDIR/predict_results" "\$PREDICTDIR/predict_misc" + fi + + mkdir -p "\$PREDICTDIR" + + # ── Guard against a corrupt partial from a previous attempt ─────────────── + if [ ! -d "\$PREDICTDIR/predict_misc" ] && [ -d "\$PREDICTDIR/predict_results" ]; then + echo "[WARN] predict_results/ present without predict_misc/ for ${out}; clearing stale partial" + rm -rf "\$PREDICTDIR/predict_results" + fi + + # Point funannotate at the persistent training dir via symlink. + if [ -d "${params.training_target}/${out}/training" ]; then + ln -sfn "${params.training_target}/${out}/training" "\$PREDICTDIR/training" + fi + + TBL2ASN_PARAMS="-l paired-ends" + + # Inflate a gzipped clean/masked genome to a local uncompressed copy. + GENOME_FA="${genome_fa}" + case "\$GENOME_FA" in + *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; + *) GENOME_IN="\$GENOME_FA" ;; + esac + + # ── Too-small-genome pre-flight guard ──────────────────────────────────── + SKIP_REPORT="${params.target}/predict_skipped_too_small.tsv" + if [ "${params.predict_min_asm_bp}" -gt 0 ]; then + read ASM_BP ASM_CTG ASM_N50 < <( + awk '/^>/{if(len)print len;len=0;next}{len+=length(\$0)}END{if(len)print len}' "\$GENOME_IN" \\ + | sort -rn \\ + | awk '{L[NR]=\$1;tot+=\$1}END{half=tot/2;run=0;n50=0;for(i=1;i<=NR;i++){run+=L[i];if(run>=half){n50=L[i];break}}print tot, NR, n50}') + echo "[INFO] Pre-flight assembly stats for ${out}: \${ASM_BP} bp, \${ASM_CTG} contigs, N50 \${ASM_N50}" + SMALL=0; FRAG=0 + [ "\$ASM_BP" -lt "${params.predict_min_asm_bp}" ] && SMALL=1 + { [ "\$ASM_N50" -lt "${params.predict_frag_max_n50}" ] || [ "\$ASM_CTG" -gt "${params.predict_frag_max_contigs}" ]; } && FRAG=1 + if [ "\$SMALL" -eq 1 ] && [ "\$FRAG" -eq 1 ]; then + echo "[WARN] ${out} is too small/fragmented for funannotate training (\${ASM_BP} bp, \${ASM_CTG} contigs, N50 \${ASM_N50}); skipping predict" >&2 + mkdir -p "${params.target}" + [ -s "\$SKIP_REPORT" ] || printf 'out\tasmid\tlocustag\treason\ttotal_bp\tcontigs\tN50\n' > "\$SKIP_REPORT" + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "${out}" "${asmid}" "${locustag}" "preflight_small_fragmented" "\$ASM_BP" "\$ASM_CTG" "\$ASM_N50" >> "\$SKIP_REPORT" + touch "\$PREDICTDIR/${out}.predict.skipped_too_small" + touch ${out}.predict.done + exit 0 + fi + fi + + funannotate predict --name ${locustag} -i "\$GENOME_IN" --strain "${strain}" \\ + -o "\$PREDICTDIR" -s "${species}" --cpu ${task.cpus} --busco_db ${busco_lineage} \\ + --AUGUSTUS_CONFIG_PATH \$AUGUSTUS_CONFIG_PATH -w codingquarry:0 glimmerhmm:0 \\ + --min_training_models 30 --tmpdir \$TMPDIR --SeqCenter ${params.seqcenter} \\ + --keep_no_stops --header_length ${header_length} --protein_evidence ${params.proteins} \\ + --max_intronlen ${params.max_intronlen} --min_intronlen ${params.min_intronlen} \\ + --tbl2asn "\$TBL2ASN_PARAMS" --table ${transl_table} --auto-skip-genemark || true + + # ── Post-predict catch ──────────────────────────────────────────────────── + if [ ! -s "\$PREDICT_GBK" ]; then + PLOG="\$PREDICTDIR/logfiles/funannotate-predict.log" + if [ -f "\$PLOG" ] && grep -q "Not enough gene models .* to train Augustus" "\$PLOG"; then + NMODELS=\$(grep -oE "Not enough gene models [0-9]+" "\$PLOG" | grep -oE "[0-9]+" | tail -1) + echo "[WARN] ${out}: funannotate found only \${NMODELS:-&2 + mkdir -p "${params.target}" + [ -s "\$SKIP_REPORT" ] || printf 'out\tasmid\tlocustag\treason\ttotal_bp\tcontigs\tN50\n' > "\$SKIP_REPORT" + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "${out}" "${asmid}" "${locustag}" "funannotate_too_few_models:\${NMODELS:-NA}" "" "" "" >> "\$SKIP_REPORT" + touch "\$PREDICTDIR/${out}.predict.skipped_too_small" + touch ${out}.predict.done + exit 0 + fi + echo "ERROR: funannotate predict did not produce expected GBK: \$PREDICT_GBK" >&2 + exit 1 + fi + if [ -d "\$PREDICTDIR/predict_misc/ab_initio_parameters" ]; then + mv "\$PREDICTDIR/predict_misc/ab_initio_parameters" "\$PREDICTDIR" + mv "\$PREDICTDIR/predict_misc/trnascan.no-overlaps.gff3" "\$PREDICTDIR" + rm -rf "\$PREDICTDIR/predict_misc" + mkdir -p "\$PREDICTDIR/predict_misc" + mv "\$PREDICTDIR/ab_initio_parameters" "\$PREDICTDIR/trnascan.no-overlaps.gff3" "\$PREDICTDIR/predict_misc" + fi + find "\$PREDICTDIR/predict_results/" -maxdepth 1 \\( -name "*.txt" -o -name "*.mrna-transcripts.fa" \\) -print0 \ + | xargs -0 --no-run-if-empty pigz + sync + touch ${out}.predict.done + echo "[INFO] Prediction complete for ${out} at \$PREDICTDIR" + """ + + stub: + def out = meta.id + """ + echo "[STUB] Would run funannotate predict for ${out} using ${genome_fa}" + [ -f "${genome_fa}" ] || [ -f "${genome_fa}.gz" ] || { echo "ERROR: genome not found at ${genome_fa}[.gz]" >&2; exit 1; } + mkdir -p ${params.target}/${out}/predict_results ${params.target}/${out}/predict_misc + echo "LOCUS stub_${out}" > ${params.target}/${out}/predict_results/${out}.gbk + echo ">stub_${out}_p1" > ${params.target}/${out}/predict_results/${out}.proteins.fa + touch ${out}.predict.done + """ +} diff --git a/modules/local/funannotate_train.nf b/modules/local/funannotate_train.nf new file mode 100644 index 0000000..41986a6 --- /dev/null +++ b/modules/local/funannotate_train.nf @@ -0,0 +1,170 @@ +// Run funannotate train (PASA alignment) for a single assembly, using shared Trinity-GG +// transcripts produced by RNASEQ_PREPARE (representative strain) or running a full +// Trinity+PASA train when no shared Trinity is available (fallback). +// Writes output directly to params.training_target//. No publishDir — the persistent +// training directory is the primary output, accessed by FUNANNOTATE_PREDICT via symlink. +// Resources overridden by withName: '.*:FUNANNOTATE_TRAIN' in conf/profile_annotate.config. +process FUNANNOTATE_TRAIN { + label 'funannotate' + label 'process_high' + tag "${meta.id}" + + input: + tuple val(meta), val(genome_fa), path(r1), path(r2), path(se), path(trinity_fa) + + output: + tuple val(meta), val(genome_fa) + + script: + def out = meta.id + def asmid = meta.asmid + def species = meta.species + def strain = meta.strain + def locustag = meta.locustag + def header_length = params.header_length + def pasa_db_arg = "--pasa_db sqlite" + """ + # ── Skip if no RNA-seq data at all ──────────────────────────────────────── + if [ ! -s "${r1}" ] && [ ! -s "${se}" ] && [ ! -s "${trinity_fa}" ]; then + echo "[INFO] No RNAseq data for ${out}, skipping funannotate train" + exit 0 + fi + + # ── Skip if training output already present and rnaseq is not newer than GBK ── + TRAIN_GFF3="${params.training_target}/${out}/training/funannotate_train.pasa.gff3" + PREDICT_GBK="${params.target}/${out}/predict_results/${out}.gbk" + [ -f "\$PREDICT_GBK" ] || PREDICT_GBK="${params.target}/${out}/predict_results/${out}.gbk.gz" + if [ -f "\$TRAIN_GFF3" ]; then + RETRAIN=0 + if [ -f "\$PREDICT_GBK" ]; then + if [ -s "${r1}" ] && [ "${r1}" -nt "\$PREDICT_GBK" ]; then + echo "[INFO] RNAseq R1 reads newer than predict GBK for ${out}; retraining" + RETRAIN=1 + elif [ -s "${se}" ] && [ "${se}" -nt "\$PREDICT_GBK" ]; then + echo "[INFO] RNAseq SE reads newer than predict GBK for ${out}; retraining" + RETRAIN=1 + fi + fi + if [ \$RETRAIN -eq 0 ]; then + echo "[INFO] Training already complete for ${out}; skipping" + exit 0 + fi + fi + + source /etc/profile.d/modules.sh 2>/dev/null || true + module load funannotate/dev-1.8.18 + export AUGUSTUS_CONFIG_PATH=${params.augustus_config} + export FUNANNOTATE_DB=${params.funannotate_db} + TMPDIR=\${SCRATCH:-/tmp} + export PASACONF="" + pasa_db_arg="--pasa_db sqlite" + # ── Optional per-task MariaDB for PASA ──────────────────────────────────── + if [ "${params.pasa_mysql}" = "true" ]; then + MYSQL_SCRATCH=${params.training_target}/${out}/training/mysql_db + if [ ! -f \$MYSQL_SCRATCH/mysql/conf/my.cnf ]; then + echo "[INFO] Setting up temporary MariaDB for PASA at \$MYSQL_SCRATCH" + mkdir -p \$MYSQL_SCRATCH/db \$MYSQL_SCRATCH/conf + rsync -a ${params.mysql_datadir}/mysql \$MYSQL_SCRATCH/db/ || \ + { echo "ERROR: Failed to copy mysql data from ${params.mysql_datadir}" >&2; exit 1; } + cp ${params.pasa_conf_dir}/my.cnf \$MYSQL_SCRATCH/conf/my.cnf || \ + { echo "ERROR: Failed to copy my.cnf" >&2; exit 1; } + fi + MYHOSTNAME=\$(hostname -s) + PORT=\$(shuf -i3000-4999 -n1) + export PASACONF=\$MYSQL_SCRATCH/conf/pasa-local-\${MYHOSTNAME}.config.txt + cp ${params.pasa_conf_dir}/conf.txt \$PASACONF + sed -i "s/^MYSQLSERVER.*\$/MYSQLSERVER=\${MYHOSTNAME}:\${PORT}/" \$PASACONF + perl -i -p -e "s/port = \\d+/port = \${PORT}/" \$MYSQL_SCRATCH/conf/my.cnf + export SINGULARITY_BINDPATH=\$TMPDIR,\$MYSQL_SCRATCH/mysql_db + stop_mysqldb() { singularity instance stop mysqldb_${asmid} 2>/dev/null || true; } + trap "stop_mysqldb; exit 130" SIGHUP SIGINT SIGTERM + trap "stop_mysqldb" EXIT + module load singularity + singularity instance start --writable-tmpfs \\ + -B \$MYSQL_SCRATCH/conf/my.cnf:/etc/mysql/my.cnf,\$MYSQL_SCRATCH/db/:/var/lib/mysql,\$MYSQL_SCRATCH/conf:/usr/conf \\ + ${params.mariadb_sif} mysqldb_${asmid} /usr/bin/mysqld_safe + pasa_db_arg="--pasa_db mysql" + sleep 5 + fi + + # Inflate a gzipped clean genome to a local uncompressed copy. + GENOME_FA="${genome_fa}" + case "\$GENOME_FA" in + *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; + *) GENOME_IN="\$GENOME_FA" ;; + esac + + # ── Use shared Trinity transcripts (PASA only) or run full train ────────── + if [ -s "${trinity_fa}" ]; then + if [ -s "${r1}" ]; then + echo "[INFO] Running funannotate train (PASA+PE) for ${out} using shared Trinity" + funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ + --trinity ${trinity_fa} --left_norm ${r1} --right_norm ${r2} \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --jaccard_clip --no-progress \\ + --max_intronlen ${params.max_intronlen} \\ + \$pasa_db_arg + elif [ -s "${se}" ]; then + echo "[INFO] Running funannotate train (PASA+SE) for ${out} using shared Trinity" + funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ + --trinity ${trinity_fa} --single_norm ${se} \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --no-progress \\ + --max_intronlen ${params.max_intronlen} \\ + \$pasa_db_arg + else + echo "[INFO] Running funannotate train (PASA only, no reads) for ${out} using shared Trinity" + funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ + --trinity ${trinity_fa} --left_norm ${r1} --right_norm ${r2} \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --jaccard_clip --no-progress \\ + --max_intronlen ${params.max_intronlen} \\ + \$pasa_db_arg + fi + elif [ -s "${r1}" ]; then + echo "[INFO] Running funannotate train (full PE, no shared Trinity) for ${out}" + funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ + --left_norm ${r1} --right_norm ${r2} --aligners minimap2 \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --jaccard_clip --no-progress --min_coverage 4 \\ + --max_intronlen ${params.max_intronlen} \\ + \$pasa_db_arg + else + echo "[INFO] Running funannotate train (full SE, no shared Trinity) for ${out}" + funannotate train -i "\$GENOME_IN" -o ${params.training_target}/${out} \\ + --single_norm ${se} --aligners minimap2 \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --no-progress --min_coverage 4 \\ + --max_intronlen ${params.max_intronlen} \\ + \$pasa_db_arg + fi + + # ── Remove large intermediates not needed for predict or update ───────────── + TRAINDIR="${params.training_target}/${out}/training" + echo "[INFO] Removing large training intermediates in \$TRAINDIR" + rm -rf "\$TRAINDIR/hisat2" + rm -rf "\$TRAINDIR/trinity_gg" + echo "[INFO] Training cleanup complete for ${out}" + echo "mysql is ${params.pasa_mysql}" + if [ "${params.pasa_mysql}" = "true" ]; then stop_mysqldb; fi + echo "[INFO] stopped mysql" + """ + + stub: + def out = meta.id + """ + echo "[STUB] FUNANNOTATE_TRAIN stub for ${out}" + mkdir -p ${params.training_target}/${out}/training + touch ${params.training_target}/${out}/training/funannotate_train.pasa.gff3 + """ +} diff --git a/modules/local/rnaseq_prepare.nf b/modules/local/rnaseq_prepare.nf new file mode 100644 index 0000000..3a9ccd5 --- /dev/null +++ b/modules/local/rnaseq_prepare.nf @@ -0,0 +1,103 @@ +// Run funannotate train --stop_after_trinity for the representative assembly of a species. +// Archives shared Trinity-GG FASTA to rnaseq_data/ so all non-representative strains can +// skip the Trinity assembly step. storeDir-cached. +// Resources overridden by withName: '.*:RNASEQ_PREPARE' in conf/profile_annotate.config. +process RNASEQ_PREPARE { + label 'funannotate' + label 'process_high' + tag "$species_tag" + + storeDir "${launchDir}/rnaseq_data" + + input: + tuple val(species_tag), val(meta), val(genome_fa), path(r1), path(r2), path(se) + + output: + tuple val(species_tag), + path("${species_tag}.trinity-GG.fasta"), emit: shared + + script: + def out = meta.id + def species = meta.species + def strain = meta.strain + def header_length = params.header_length + """ + # ── Empty-reads sentinel: no RNA-seq found by SRA_FETCH / SRA_FETCH_SE ── + if [ ! -s "${r1}" ] && [ ! -s "${se}" ]; then + echo "[INFO] No RNAseq reads for ${species_tag}; writing empty shared markers" + touch ${species_tag}.trinity-GG.fasta + exit 0 + fi + + # ── If representative was already trained, just extract shared files ────── + TRAIN_GFF3="${params.training_target}/${out}/training/funannotate_train.pasa.gff3" + if [ -f "\$TRAIN_GFF3" ]; then + echo "[INFO] Training already complete for ${out}; extracting shared files to rnaseq_data" + TRAINDIR="${params.training_target}/${out}/training" + TRINITY_FA=\$(find \$TRAINDIR -maxdepth 1 -name "trinity.fasta" | head -1) + if [ -n "\$TRINITY_FA" ]; then + cp "\$TRINITY_FA" ${species_tag}.trinity-GG.fasta + else + touch ${species_tag}.trinity-GG.fasta + fi + exit 0 + fi + + source /etc/profile.d/modules.sh 2>/dev/null || true + + export AUGUSTUS_CONFIG_PATH=${params.augustus_config} + export FUNANNOTATE_DB=${params.funannotate_db} + TMPDIR=\${SCRATCH:-/tmp} + + echo "[INFO] RNASEQ_PREPARE: running funannotate train for representative ${out} (species: ${species_tag})" + + # Inflate a gzipped clean genome to a local uncompressed copy. + GENOME_FA="${genome_fa}" + case "\$GENOME_FA" in + *.gz) echo "[INFO] Inflating compressed genome \$GENOME_FA"; pigz -dc "\$GENOME_FA" > genome_input.fa; GENOME_IN="\$(pwd)/genome_input.fa" ;; + *) GENOME_IN="\$GENOME_FA" ;; + esac + + if [ -s "${r1}" ]; then + funannotate train -i "\$GENOME_IN" -o \$SCRATCH/${out} \\ + --left_norm ${r1} --right_norm ${r2} --aligners minimap2 \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --jaccard_clip --no-progress --min_coverage 4 \\ + --max_intronlen ${params.max_intronlen} \\ + --stop_after_trinity --no_trimmomatic + else + echo "[INFO] RNASEQ_PREPARE: using single-end reads for ${out}" + funannotate train -i "\$GENOME_IN" -o \$SCRATCH/${out} \\ + --single_norm ${se} --aligners minimap2 \\ + --species "${species}" --strain "${strain}" \\ + --cpus ${task.cpus} --memory ${task.memory.toGiga()}G \\ + --header_length ${header_length} \\ + --no-progress --min_coverage 4 \\ + --max_intronlen ${params.max_intronlen} \\ + --stop_after_trinity --no_trimmomatic + fi + + # ── Copy shared outputs to rnaseq_data/ ────────────────────────────────── + TRAINDIR="\$SCRATCH/${out}/training" + TRINITY_FA=\$(find \$TRAINDIR -maxdepth 1 -name "trinity.fasta" | head -1) + if [ -n "\$TRINITY_FA" ]; then + cp "\$TRINITY_FA" ${species_tag}.trinity-GG.fasta + else + echo "[WARN] No trinity.fasta found under \$TRAINDIR for ${out}" + touch ${species_tag}.trinity-GG.fasta + fi + + rm -rf "\$SCRATCH/${out}" + echo "[INFO] RNASEQ_PREPARE complete for ${species_tag}" + """ + + stub: + def out = meta.id + """ + echo ">stub_trinity_${species_tag}" > ${species_tag}.trinity-GG.fasta + mkdir -p ${params.training_target}/${out}/training + touch ${params.training_target}/${out}/training/funannotate_train.pasa.gff3 + """ +} diff --git a/subworkflows/local/train_predict.nf b/subworkflows/local/train_predict.nf new file mode 100644 index 0000000..6a71051 --- /dev/null +++ b/subworkflows/local/train_predict.nf @@ -0,0 +1,143 @@ +/* + * TRAIN_PREDICT — RNA-seq training + gene prediction for one batch of assemblies. + * + * Takes: + * ch_genomes — channel: tuple(val(meta), val(genome_fa)) (from MASK_GENOME) + * ch_reads — channel: tuple(val(species_tag), path(r1), path(r2), path(se)) + * May be empty when run_sra_fetch=false or stop_after_sra_fetch=true. + * + * Steps: + * 1. When reads are available: + * a. RNASEQ_PREPARE: run once per species on the representative assembly (storeDir). + * b. FUNANNOTATE_TRAIN: run PASA alignment for every assembly using the shared Trinity. + * Species with no reads bypass training entirely. + * 2. FUNANNOTATE_PREDICT: run for all assemblies not yet predicted (or stale). + * + * Emits: + * metadata — val(meta) for each assembly that completed predict + * done — path to *.predict.done marker file + */ + +include { RNASEQ_PREPARE } from './../../modules/local/rnaseq_prepare' +include { FUNANNOTATE_TRAIN } from './../../modules/local/funannotate_train' +include { FUNANNOTATE_PREDICT } from './../../modules/local/funannotate_predict' + +// GBK may be stored compressed to save space; return the first existing non-empty file. +def _gbkResult(String dir, String id) { + def plain = file("${dir}/${id}.gbk") + if (plain.exists() && plain.size() > 0) return plain + def gz = file("${dir}/${id}.gbk.gz") + if (gz.exists() && gz.size() > 0) return gz + return null +} + +// True when rnaseq reads or shared Trinity are newer than the existing prediction GBK. +def _staleRnaseq(String id, String species) { + def species_tag = species.replaceAll(/\s+/, '_') + def gbk = _gbkResult("${params.target}/${id}/predict_results", id) + if (gbk == null) return false + def gbkMod = gbk.lastModified() + def r1 = file("${launchDir}/rnaseq_reads/${species_tag}_norm_R1.fastq.gz") + def se = file("${launchDir}/rnaseq_reads/${species_tag}_norm_SE.fastq.gz") + def trinity = file("${launchDir}/rnaseq_data/${species_tag}.trinity-GG.fasta") + return (r1.exists() && r1.size() > 0 && r1.lastModified() > gbkMod) || + (se.exists() && se.size() > 0 && se.lastModified() > gbkMod) || + (trinity.exists() && trinity.size() > 0 && trinity.lastModified() > gbkMod) +} + +workflow TRAIN_PREDICT { + + take: + ch_genomes // tuple(val(meta), val(genome_fa)) + ch_reads // tuple(val(species_tag), path(r1), path(r2), path(se)) — may be empty + + main: + def predict_input_ch + + if (params.run_sra_fetch.toBoolean() && !params.stop_after_sra_fetch.toBoolean()) { + // ── Build per-assembly channel keyed by species_tag with SRA reads joined ── + // reads_ch is a 4-tuple: (species_tag, r1, r2, se) + def assembly_with_reads = ch_genomes + .map { meta, genome_fa -> + def species_tag = meta.species.replaceAll(/\s+/, '_') + tuple(species_tag, meta, genome_fa) + } + .combine(ch_reads, by: 0) + // assembly_with_reads: (species_tag, meta, genome_fa, r1, r2, se) + + // ── RNASEQ_PREPARE: once per species on the representative (first) assembly ── + // storeDir-cached; subsequent runs reuse the Trinity-GG FASTA. + def repr_ch = assembly_with_reads + .groupTuple(by: 0) + .map { species_tag, metas, genomes, r1s, r2s, ses -> + tuple(species_tag, metas[0], genomes[0], r1s[0], r2s[0], ses[0]) + } + + def repr_branched = repr_ch.branch { + has_reads: it[3].size() > 0 || it[5].size() > 0 // r1=[3] or se=[5] + no_reads: true + } + + RNASEQ_PREPARE(repr_branched.has_reads) + + // For species with no RNA-seq, write an empty Trinity FASTA in the driver process + // (no SLURM job) and emit it directly as a shared channel item. + def empty_shared_ch = repr_branched.no_reads + .map { species_tag, _meta, _gfa, _r1, _r2, _se -> + def empty_fa = file("${launchDir}/rnaseq_data/${species_tag}.trinity-GG.fasta") + if (!empty_fa.exists()) { + empty_fa.parent.mkdirs() + empty_fa.text = '' + } + tuple(species_tag, empty_fa) + } + + def shared_ch = RNASEQ_PREPARE.out.shared.mix(empty_shared_ch) + + // ── Join shared Trinity back to every assembly for FUNANNOTATE_TRAIN ── + def train_input = assembly_with_reads + .combine(shared_ch, by: 0) + .map { species_tag, meta, genome_fa, r1, r2, se, trinity_fa -> + tuple(meta, genome_fa, r1, r2, se, trinity_fa) + } + // train_input: (meta, genome_fa, r1, r2, se, trinity_fa) + + // Assemblies with no RNA-seq bypass FUNANNOTATE_TRAIN entirely. + def branched = train_input.branch { + has_rnaseq: it[2].size() > 0 || it[4].size() > 0 || it[5].size() > 0 + no_rnaseq: true + } + def predict_no_rnaseq = branched.no_rnaseq + .map { meta, genome_fa, _r1, _r2, _se, _tf -> tuple(meta, genome_fa) } + + // Skip TRAIN when pasa.gff3 already exists and is not stale relative to reads. + def train_todo = branched.has_rnaseq.filter { meta, _gfa, _r1, _r2, _se, _tf -> + def gff3 = file("${params.training_target}/${meta.id}/training/funannotate_train.pasa.gff3") + !gff3.exists() || gff3.size() == 0 || _staleRnaseq(meta.id as String, meta.species as String) + } + def train_done = branched.has_rnaseq + .filter { meta, _gfa, _r1, _r2, _se, _tf -> + def gff3 = file("${params.training_target}/${meta.id}/training/funannotate_train.pasa.gff3") + gff3.exists() && gff3.size() > 0 && !_staleRnaseq(meta.id as String, meta.species as String) + } + .map { meta, genome_fa, _r1, _r2, _se, _tf -> tuple(meta, genome_fa) } + + FUNANNOTATE_TRAIN(train_todo) + predict_input_ch = FUNANNOTATE_TRAIN.out.mix(train_done).mix(predict_no_rnaseq) + } else { + // No RNA-seq: pass genomes straight to predict. + predict_input_ch = ch_genomes + } + + // ── FUNANNOTATE_PREDICT ───────────────────────────────────────────────── + def predict_ch = predict_input_ch + .filter { meta, _gfa -> + _gbkResult("${params.target}/${meta.id}/predict_results", meta.id as String) == null || + _staleRnaseq(meta.id as String, meta.species as String) + } + FUNANNOTATE_PREDICT(predict_ch) + + emit: + metadata = FUNANNOTATE_PREDICT.out.metadata + done = FUNANNOTATE_PREDICT.out.done +} From f8c206305d6b05be6d0c68952ce986ac657aba6f Mon Sep 17 00:00:00 2001 From: Jason Stajich Date: Mon, 29 Jun 2026 10:58:55 -0700 Subject: [PATCH 5/5] =?UTF-8?q?feat(#10):=20extract=20FUNANNOTATE=5FANNOTA?= =?UTF-8?q?TE=20+=20FUNANNOTATE=5FUPDATE=20=E2=86=92=20ANNOTATE=5FGENOME?= =?UTF-8?q?=20subworkflow?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two processes moved out of funannotate.nf monolith: modules/local/funannotate_annotate.nf — funannotate annotate (functional annotation) modules/local/funannotate_update.nf — funannotate update (PASA re-alignment post-predict) Composed with existing tool modules into subworkflows/local/annotate_genome.nf (ANNOTATE_GENOME): take: ch_predict_meta (val meta), ch_reads (species_tag, r1, r2, se) - Optional ANTISMASH_RUN → INTERPROSCAN_RUN → SIGNALP_RUN gating chain (unchanged logic) - Optional FUNANNOTATE_UPDATE (joins reads_ch; gated by gbk existence; update→annotate) - Optional FUNANNOTATE_ANNOTATE (gated behind all requested pre-annotate steps) emit: versions (ch_versions from ANTISMASH/INTERPRO/SIGNALP) funannotate.nf: ~100 lines of post-predict orchestration replaced by: ANNOTATE_GENOME(predict_meta, reads_ch) ch_versions = ch_versions.mix(ANNOTATE_GENOME.out.versions) withName patterns updated: ANTISMASH_RUN, INTERPROSCAN_RUN, SIGNALP_RUN, FUNANNOTATE_ANNOTATE, FUNANNOTATE_UPDATE → '.*:NAME' in conf/profile_annotate.config. _gbkResult inlined as top-level function in annotate_genome.nf (same pattern as other subworkflows). funannotate.nf now contains: only includes, 3 utility defs (gbkResult/genomeFile/staleRnaseq), and the workflow{} block. No inline process definitions remain. Stub-run gate: 15 completed, 0 failed Dual-stub gate: 17 completed, 0 failed ANNOTATE_GENOME:ANTISMASH_RUN, INTERPROSCAN_RUN, SIGNALP_RUN, FUNANNOTATE_ANNOTATE all correctly namespaced Co-Authored-By: Claude Sonnet 4.6 Claude-Session: https://claude.ai/code/session_01AbWwmm9v6dJn13sgAJn3Sa --- conf/profile_annotate.config | 10 +- funannotate.nf | 258 +------------------------- modules/local/funannotate_annotate.nf | 54 ++++++ modules/local/funannotate_update.nf | 95 ++++++++++ subworkflows/local/annotate_genome.nf | 119 ++++++++++++ 5 files changed, 280 insertions(+), 256 deletions(-) create mode 100644 modules/local/funannotate_annotate.nf create mode 100644 modules/local/funannotate_update.nf create mode 100644 subworkflows/local/annotate_genome.nf diff --git a/conf/profile_annotate.config b/conf/profile_annotate.config index 6d1ebba..984106c 100644 --- a/conf/profile_annotate.config +++ b/conf/profile_annotate.config @@ -199,25 +199,25 @@ process { clusterOptions = '-N 1 -n 1 --job-name=nf-PREDICT' } - withName: 'ANTISMASH_RUN' { + withName: '.*:ANTISMASH_RUN' { queue = 'epyc'; cpus = 8; memory = '24 GB'; time = '8h' } - withName: 'INTERPROSCAN_RUN' { + withName: '.*:INTERPROSCAN_RUN' { queue = 'epyc'; cpus = 8; memory = '32 GB'; time = '60h' } - withName: 'SIGNALP_RUN' { + withName: '.*:SIGNALP_RUN' { // signalp6-gpu: GPU partition (gpu allows >2h; short_gpu if always ≤2h). queue = 'gpu' cpus = 16; memory = '64 GB'; time = '8h' clusterOptions = '-N 1 -n 1 --gres=gpu:1' } - withName: 'FUNANNOTATE_ANNOTATE' { + withName: '.*:FUNANNOTATE_ANNOTATE' { queue = 'epyc'; cpus = 16; memory = '64 GB'; time = '24h' } - withName: 'FUNANNOTATE_UPDATE' { + withName: '.*:FUNANNOTATE_UPDATE' { queue = 'epyc'; cpus = 16; memory = '64 GB'; time = '48h' } } diff --git a/funannotate.nf b/funannotate.nf index c4e2dfc..3c435c0 100644 --- a/funannotate.nf +++ b/funannotate.nf @@ -21,160 +21,6 @@ // header_length is NOT in meta — it comes from params.header_length (default 24). -process FUNANNOTATE_ANNOTATE { - label 'funannotate' - tag "${meta.id}" - - cpus 16 - memory '32 GB' - time '48h' - - input: - val(meta) - - output: - tuple val(meta), path("${meta.id}.annotate.done"), emit: marker - - script: - def out = meta.id - def species = meta.species - def strain = meta.strain - def locustag = meta.locustag - def busco_lineage = meta.busco - def header_length = params.header_length - def antiSm = file("${params.target}/${meta.id}/antismash_local/${meta.id}.gbk") - def antiSmArg = antiSm.exists() ? "--antismash ${antiSm}" : "" - """ - source /etc/profile.d/modules.sh 2>/dev/null || true - - export AUGUSTUS_CONFIG_PATH=${params.augustus_config} - export FUNANNOTATE_DB=${params.funannotate_db} - TMPDIR=\${SCRATCH:-/tmp} - - funannotate annotate -i ${params.target}/${out} -o ${params.target}/${out} \\ - --species "${species}" --strain "${strain}" \\ - --busco_db ${busco_lineage} --rename ${locustag} \\ - --sbt ${params.sbt_template} \\ - --header_length ${header_length} \\ - ${antiSmArg} \\ - --cpu ${task.cpus} --tmpdir \$TMPDIR - - EXPECTED_GBK="${params.target}/${out}/annotate_results/${out}.gbk" - if [ ! -f "\$EXPECTED_GBK" ]; then - echo "ERROR: funannotate annotate did not produce expected GBK: \$EXPECTED_GBK" >&2 - exit 1 - fi - touch ${out}.annotate.done - """ - - stub: - def out = meta.id - """ - echo "[STUB] Would run funannotate annotate for ${out}" - mkdir -p ${params.target}/${out}/annotate_results ${params.target}/${out}/annotate_misc - touch ${params.target}/${out}/annotate_results/${out}.gbk - touch ${out}.annotate.done - """ -} - -process FUNANNOTATE_UPDATE { - label 'funannotate' - tag "${meta.id}" - - cpus 16 - memory '96 GB' - time '48h' - - input: - tuple val(meta), path(r1), path(r2) - - output: - val meta - - script: - def out = meta.id - def asmid = meta.asmid - def species = meta.species - def strain = meta.strain - def locustag = meta.locustag - def busco_lineage = meta.busco - def header_length = params.header_length - def pasa_db_arg = "--pasa_db sqlite" - """ - # ── Skip if no reads (empty marker file from SRA_FETCH) ────────────────── - if [ ! -s "${r1}" ]; then - echo "[INFO] No RNAseq reads for ${out}, skipping funannotate update" - exit 0 - fi - - source /etc/profile.d/modules.sh 2>/dev/null || true - - export AUGUSTUS_CONFIG_PATH=${params.augustus_config} - export FUNANNOTATE_DB=${params.funannotate_db} - TMPDIR=\${SCRATCH:-/tmp} - export PASACONF="" - pasa_db_arg="--pasa_db sqlite" - # ── Optional per-task MariaDB for PASA ──────────────────────────────────── - if [ "${params.pasa_mysql}" = "true" ]; then - MYSQL_SCRATCH=${params.training_target}/${out}/training/mysql_db - if [ ! -f \$MYSQL_SCRATCH/mysql/conf/my.cnf ]; then - echo "[INFO] Setting up temporary MariaDB for PASA at \$MYSQL_SCRATCH" - mkdir -p \$MYSQL_SCRATCH/db \$MYSQL_SCRATCH/conf - rsync -a ${params.mysql_datadir}/mysql \$MYSQL_SCRATCH/db/ || \ - { echo "ERROR: Failed to copy mysql data from ${params.mysql_datadir}" >&2; exit 1; } - cp ${params.pasa_conf_dir}/my.cnf \$MYSQL_SCRATCH/conf/my.cnf || \ - { echo "ERROR: Failed to copy my.cnf" >&2; exit 1; } - fi - MYHOSTNAME=\$(hostname -s) - PORT=\$(shuf -i3000-4999 -n1) - export PASACONF=\$MYSQL_SCRATCH/conf/pasa-local-\${MYHOSTNAME}.config.txt - cp ${params.pasa_conf_dir}/conf.txt \$PASACONF - sed -i "s/^MYSQLSERVER.*\$/MYSQLSERVER=\${MYHOSTNAME}:\${PORT}/" \$PASACONF - perl -i -p -e "s/port = \\d+/port = \${PORT}/" \$MYSQL_SCRATCH/conf/my.cnf - export SINGULARITY_BINDPATH=\$TMPDIR,\$MYSQL_SCRATCH/mysql_db - stop_mysqldb() { singularity instance stop mysqldb_${asmid} 2>/dev/null || true; } - trap "stop_mysqldb; exit 130" SIGHUP SIGINT SIGTERM - trap "stop_mysqldb" EXIT - module load singularity - singularity instance start --writable-tmpfs \\ - -B \$MYSQL_SCRATCH/conf/my.cnf:/etc/mysql/my.cnf,\$MYSQL_SCRATCH/db/:/var/lib/mysql,\$MYSQL_SCRATCH/conf:/usr/conf \\ - ${params.mariadb_sif} mysqldb_${asmid} /usr/bin/mysqld_safe - pasa_db_arg="--pasa_db mysql" - sleep 5 - fi - - # Link training data into work dir so funannotate update finds it at the relative path it expects. - mkdir -p ${out} - if [ -d "${params.training_target}/${out}/training" ]; then - ln -sfn "${params.training_target}/${out}/training" "${out}/training" - fi - - # r1/r2 are pre-normalized reads from SRA_FETCH (fastp-trimmed + bbnorm-normalized). - # funannotate update will still run its internal alignment step against these. - echo "[INFO] Running funannotate update for ${out}" - funannotate update -i ${params.target}/${out} \\ - --left ${r1} --right ${r2} \\ - --cpus ${task.cpus} \\ - \$pasa_db_arg - if [ "${params.pasa_mysql}" = "true" ]; then stop_mysqldb; fi - echo "[INFO] stopped mysql" - EXPECTED="${params.target}/${out}/update_results/${out}.gbk" - if [ ! -f "\$EXPECTED" ]; then - echo "ERROR: funannotate update did not produce expected GBK: \$EXPECTED" >&2 - exit 1 - fi - """ - - stub: - def out = meta.id - """ - echo "[STUB] FUNANNOTATE_UPDATE stub for ${out} (r1=${r1}, r2=${r2})" - mkdir -p ${params.target}/${out}/update_results - touch ${params.target}/${out}/update_results/${out}.tbl - touch ${params.target}/${out}/update_results/${out}.gbk - touch ${params.target}/${out}/update_results/${out}.gff3 - """ -} // A funannotate step's GenBank output may be stored uncompressed (.gbk) or // gzip-compressed (.gbk.gz) so completed folders can be compressed to save space. @@ -224,9 +70,7 @@ include { CLEAN_GENOMES } from './subworkflows/local/clean_genomes' include { MASK_GENOME } from './subworkflows/local/mask_genome' include { FETCH_RNASEQ } from './subworkflows/local/fetch_rnaseq' include { TRAIN_PREDICT } from './subworkflows/local/train_predict' -include { ANTISMASH_RUN } from './modules/local/antismash_run' -include { INTERPROSCAN_RUN } from './modules/local/interproscan_run' -include { SIGNALP_RUN } from './modules/local/signalp_run' +include { ANNOTATE_GENOME } from './subworkflows/local/annotate_genome' workflow { // `--help` prints schema-driven parameter help (grouped, with types/defaults) and exits. @@ -323,105 +167,17 @@ workflow { if (!params.run_sra_fetch.toBoolean() || (!params.stop_after_sra_fetch.toBoolean() && !params.stop_after_sra_query.toBoolean())) { TRAIN_PREDICT(predict_genome_ch, reads_ch) - // ── Post-predict steps and annotation ──────────────────────────────────── - // postpredict: all samples with a completed predict_results/*.gbk, whether - // produced in this run or a prior one. This is the source for all optional - // pre-annotate steps and for FUNANNOTATE_ANNOTATE itself. + // Genomes predicted in a PRIOR run (complement of what TRAIN_PREDICT runs this run). def postpredict = INPUT_CHECK.out.samples - // Only genomes predicted in a PRIOR run — complement of what TRAIN_PREDICT runs now. .filter { meta -> - gbkResult("${params.target}/${meta.id}/predict_results", meta.id as String) != null && !staleRnaseq(meta.id as String, meta.species as String) + gbkResult("${params.target}/${meta.id}/predict_results", meta.id as String) != null && + !staleRnaseq(meta.id as String, meta.species as String) } - - // Same-run completion gate: new predictions arrive via TRAIN_PREDICT.out.metadata; - // prior-run completions arrive via postpredict (available immediately). + // Merge prior-run and same-run predictions; feed into annotation chain. def predict_meta = postpredict.mix(TRAIN_PREDICT.out.metadata) - def annotate_ready_ch = predict_meta - - if (params.run_antismash.toBoolean()) { - def as_todo = annotate_ready_ch.filter { meta -> - def asDir = file("${params.target}/${meta.id}/antismash_local") - !(asDir.isDirectory() && asDir.list()?.any { it.endsWith('.json') || it.endsWith('.json.gz') }) - } - def as_done = annotate_ready_ch.filter { meta -> - def asDir = file("${params.target}/${meta.id}/antismash_local") - asDir.isDirectory() && asDir.list()?.any { it.endsWith('.json') || it.endsWith('.json.gz') } - } - ANTISMASH_RUN(as_todo) - ch_versions = ch_versions.mix(ANTISMASH_RUN.out.versions) - def as_completed = ANTISMASH_RUN.out.results - .map { meta, _files -> meta } - annotate_ready_ch = as_completed.mix(as_done) - } - - if (params.run_interpro.toBoolean()) { - def ipr_todo = annotate_ready_ch.filter { meta -> - !file("${params.target}/${meta.id}/annotate_misc/iprscan.xml").exists() - } - def ipr_done = annotate_ready_ch.filter { meta -> - file("${params.target}/${meta.id}/annotate_misc/iprscan.xml").exists() - } - INTERPROSCAN_RUN(ipr_todo) - ch_versions = ch_versions.mix(INTERPROSCAN_RUN.out.versions) - def ipr_completed = INTERPROSCAN_RUN.out.results - .map { meta, _xml -> meta } - annotate_ready_ch = ipr_completed.mix(ipr_done) - } - if (params.run_signalp.toBoolean()) { - def sp_todo = annotate_ready_ch.filter { meta -> - !file("${params.target}/${meta.id}/annotate_misc/signalp.results.txt").exists() - } - def sp_done = annotate_ready_ch.filter { meta -> - file("${params.target}/${meta.id}/annotate_misc/signalp.results.txt").exists() - } - SIGNALP_RUN(sp_todo) - ch_versions = ch_versions.mix(SIGNALP_RUN.out.versions) - def sp_completed = SIGNALP_RUN.out.results - .map { meta, _txt -> meta } - annotate_ready_ch = sp_completed.mix(sp_done) - } - - if (params.run_update.toBoolean()) { - if (params.run_sra_fetch.toBoolean()) { - // UPDATE runs from predict results in parallel with antismash/interpro/signalp. - // Reads are joined from SRA_FETCH (storeDir-cached, so prior-run reads are reused). - // The join on upd_signal gates annotate_ready_ch so ANNOTATE waits for UPDATE. - def upd_input = predict_meta - .map { meta -> - def species_tag = meta.species.replaceAll(/\s+/, '_') - tuple(species_tag, meta) - } - .combine(reads_ch, by: 0) - .map { _st, meta, r1, r2 -> - tuple(meta, r1, r2) - } - def upd_todo = upd_input.filter { meta, _r1, _r2 -> - gbkResult("${params.target}/${meta.id}/update_results", meta.id as String) == null - } - def upd_done_signal = upd_input - .filter { meta, _r1, _r2 -> - gbkResult("${params.target}/${meta.id}/update_results", meta.id as String) != null - } - .map { meta, _r1, _r2 -> tuple(meta.id, 'upd') } - FUNANNOTATE_UPDATE(upd_todo) - def upd_signal = FUNANNOTATE_UPDATE.out - .map { meta -> tuple(meta.id, 'upd') } - .mix(upd_done_signal) - annotate_ready_ch = annotate_ready_ch - .map { meta -> tuple(meta.id, meta) } - .join(upd_signal) - .map { _id, meta, _flag -> meta } - } else { - log.warn "run_update=true but run_sra_fetch=false; funannotate update skipped (no reads available)" - } - } - - if (params.run_annotate.toBoolean()) { - FUNANNOTATE_ANNOTATE(annotate_ready_ch.filter { meta -> - gbkResult("${params.target}/${meta.id}/annotate_results", meta.id as String) == null - }) - } + ANNOTATE_GENOME(predict_meta, reads_ch) + ch_versions = ch_versions.mix(ANNOTATE_GENOME.out.versions) } // end if (!params.stop_after_sra_fetch || !params.run_sra_fetch) } diff --git a/modules/local/funannotate_annotate.nf b/modules/local/funannotate_annotate.nf new file mode 100644 index 0000000..d23b85c --- /dev/null +++ b/modules/local/funannotate_annotate.nf @@ -0,0 +1,54 @@ +// Run funannotate annotate for one assembly. Writes output to params.target//annotate_results/. +// Resources overridden by withName: '.*:FUNANNOTATE_ANNOTATE' in conf/profile_annotate.config. +process FUNANNOTATE_ANNOTATE { + label 'funannotate' + label 'process_high' + tag "${meta.id}" + + input: + val(meta) + + output: + tuple val(meta), path("${meta.id}.annotate.done"), emit: marker + + script: + def out = meta.id + def species = meta.species + def strain = meta.strain + def locustag = meta.locustag + def busco_lineage = meta.busco + def header_length = params.header_length + def antiSm = file("${params.target}/${meta.id}/antismash_local/${meta.id}.gbk") + def antiSmArg = antiSm.exists() ? "--antismash ${antiSm}" : "" + """ + source /etc/profile.d/modules.sh 2>/dev/null || true + + export AUGUSTUS_CONFIG_PATH=${params.augustus_config} + export FUNANNOTATE_DB=${params.funannotate_db} + TMPDIR=\${SCRATCH:-/tmp} + + funannotate annotate -i ${params.target}/${out} -o ${params.target}/${out} \\ + --species "${species}" --strain "${strain}" \\ + --busco_db ${busco_lineage} --rename ${locustag} \\ + --sbt ${params.sbt_template} \\ + --header_length ${header_length} \\ + ${antiSmArg} \\ + --cpu ${task.cpus} --tmpdir \$TMPDIR + + EXPECTED_GBK="${params.target}/${out}/annotate_results/${out}.gbk" + if [ ! -f "\$EXPECTED_GBK" ]; then + echo "ERROR: funannotate annotate did not produce expected GBK: \$EXPECTED_GBK" >&2 + exit 1 + fi + touch ${out}.annotate.done + """ + + stub: + def out = meta.id + """ + echo "[STUB] Would run funannotate annotate for ${out}" + mkdir -p ${params.target}/${out}/annotate_results ${params.target}/${out}/annotate_misc + touch ${params.target}/${out}/annotate_results/${out}.gbk + touch ${out}.annotate.done + """ +} diff --git a/modules/local/funannotate_update.nf b/modules/local/funannotate_update.nf new file mode 100644 index 0000000..831ea30 --- /dev/null +++ b/modules/local/funannotate_update.nf @@ -0,0 +1,95 @@ +// Run funannotate update (PASA re-alignment with RNA-seq reads post-predict). +// Resources overridden by withName: '.*:FUNANNOTATE_UPDATE' in conf/profile_annotate.config. +process FUNANNOTATE_UPDATE { + label 'funannotate' + label 'process_high' + tag "${meta.id}" + + input: + tuple val(meta), path(r1), path(r2) + + output: + val meta + + script: + def out = meta.id + def asmid = meta.asmid + def species = meta.species + def strain = meta.strain + def locustag = meta.locustag + def busco_lineage = meta.busco + def header_length = params.header_length + def pasa_db_arg = "--pasa_db sqlite" + """ + # ── Skip if no reads (empty marker file from SRA_FETCH) ────────────────── + if [ ! -s "${r1}" ]; then + echo "[INFO] No RNAseq reads for ${out}, skipping funannotate update" + exit 0 + fi + + source /etc/profile.d/modules.sh 2>/dev/null || true + + export AUGUSTUS_CONFIG_PATH=${params.augustus_config} + export FUNANNOTATE_DB=${params.funannotate_db} + TMPDIR=\${SCRATCH:-/tmp} + export PASACONF="" + pasa_db_arg="--pasa_db sqlite" + # ── Optional per-task MariaDB for PASA ──────────────────────────────────── + if [ "${params.pasa_mysql}" = "true" ]; then + MYSQL_SCRATCH=${params.training_target}/${out}/training/mysql_db + if [ ! -f \$MYSQL_SCRATCH/mysql/conf/my.cnf ]; then + echo "[INFO] Setting up temporary MariaDB for PASA at \$MYSQL_SCRATCH" + mkdir -p \$MYSQL_SCRATCH/db \$MYSQL_SCRATCH/conf + rsync -a ${params.mysql_datadir}/mysql \$MYSQL_SCRATCH/db/ || \ + { echo "ERROR: Failed to copy mysql data from ${params.mysql_datadir}" >&2; exit 1; } + cp ${params.pasa_conf_dir}/my.cnf \$MYSQL_SCRATCH/conf/my.cnf || \ + { echo "ERROR: Failed to copy my.cnf" >&2; exit 1; } + fi + MYHOSTNAME=\$(hostname -s) + PORT=\$(shuf -i3000-4999 -n1) + export PASACONF=\$MYSQL_SCRATCH/conf/pasa-local-\${MYHOSTNAME}.config.txt + cp ${params.pasa_conf_dir}/conf.txt \$PASACONF + sed -i "s/^MYSQLSERVER.*\$/MYSQLSERVER=\${MYHOSTNAME}:\${PORT}/" \$PASACONF + perl -i -p -e "s/port = \\d+/port = \${PORT}/" \$MYSQL_SCRATCH/conf/my.cnf + export SINGULARITY_BINDPATH=\$TMPDIR,\$MYSQL_SCRATCH/mysql_db + stop_mysqldb() { singularity instance stop mysqldb_${asmid} 2>/dev/null || true; } + trap "stop_mysqldb; exit 130" SIGHUP SIGINT SIGTERM + trap "stop_mysqldb" EXIT + module load singularity + singularity instance start --writable-tmpfs \\ + -B \$MYSQL_SCRATCH/conf/my.cnf:/etc/mysql/my.cnf,\$MYSQL_SCRATCH/db/:/var/lib/mysql,\$MYSQL_SCRATCH/conf:/usr/conf \\ + ${params.mariadb_sif} mysqldb_${asmid} /usr/bin/mysqld_safe + pasa_db_arg="--pasa_db mysql" + sleep 5 + fi + + # Link training data into work dir so funannotate update finds it at the relative path it expects. + mkdir -p ${out} + if [ -d "${params.training_target}/${out}/training" ]; then + ln -sfn "${params.training_target}/${out}/training" "${out}/training" + fi + + echo "[INFO] Running funannotate update for ${out}" + funannotate update -i ${params.target}/${out} \\ + --left ${r1} --right ${r2} \\ + --cpus ${task.cpus} \\ + \$pasa_db_arg + if [ "${params.pasa_mysql}" = "true" ]; then stop_mysqldb; fi + echo "[INFO] stopped mysql" + EXPECTED="${params.target}/${out}/update_results/${out}.gbk" + if [ ! -f "\$EXPECTED" ]; then + echo "ERROR: funannotate update did not produce expected GBK: \$EXPECTED" >&2 + exit 1 + fi + """ + + stub: + def out = meta.id + """ + echo "[STUB] FUNANNOTATE_UPDATE stub for ${out} (r1=${r1}, r2=${r2})" + mkdir -p ${params.target}/${out}/update_results + touch ${params.target}/${out}/update_results/${out}.tbl + touch ${params.target}/${out}/update_results/${out}.gbk + touch ${params.target}/${out}/update_results/${out}.gff3 + """ +} diff --git a/subworkflows/local/annotate_genome.nf b/subworkflows/local/annotate_genome.nf new file mode 100644 index 0000000..6c44a1e --- /dev/null +++ b/subworkflows/local/annotate_genome.nf @@ -0,0 +1,119 @@ +/* + * ANNOTATE_GENOME — optional pre-annotate steps + funannotate update + funannotate annotate. + * + * Accepts all assemblies with a completed predict GBK (both prior-run and same-run). + * Chains optional tools (ANTISMASH, INTERPRO, SIGNALP) before annotation; each tool + * splits the channel into todo/done, processes todo, then re-mixes so FUNANNOTATE_ANNOTATE + * only fires once all requested tools are complete. + * + * Takes: + * ch_predict_meta — val(meta) for every assembly ready for annotation + * ch_reads — tuple(val(species_tag), path(r1), path(r2), path(se)) + * May be empty when run_sra_fetch=false; only used by UPDATE. + * + * Emits: + * versions — mixed versions.yml paths from ANTISMASH / INTERPRO / SIGNALP + */ + +include { ANTISMASH_RUN } from './../../modules/local/antismash_run' +include { INTERPROSCAN_RUN } from './../../modules/local/interproscan_run' +include { SIGNALP_RUN } from './../../modules/local/signalp_run' +include { FUNANNOTATE_UPDATE } from './../../modules/local/funannotate_update' +include { FUNANNOTATE_ANNOTATE } from './../../modules/local/funannotate_annotate' + +// GBK may be stored compressed to save space; return the first non-empty file or null. +def _gbkResult(String dir, String id) { + def plain = file("${dir}/${id}.gbk") + if (plain.exists() && plain.size() > 0) return plain + def gz = file("${dir}/${id}.gbk.gz") + if (gz.exists() && gz.size() > 0) return gz + return null +} + +workflow ANNOTATE_GENOME { + + take: + ch_predict_meta // val(meta) + ch_reads // tuple(val(species_tag), path(r1), path(r2), path(se)) — may be empty + + main: + def ch_versions = Channel.empty() + def annotate_ready_ch = ch_predict_meta + + if (params.run_antismash.toBoolean()) { + def as_todo = annotate_ready_ch.filter { meta -> + def asDir = file("${params.target}/${meta.id}/antismash_local") + !(asDir.isDirectory() && asDir.list()?.any { it.endsWith('.json') || it.endsWith('.json.gz') }) + } + def as_done = annotate_ready_ch.filter { meta -> + def asDir = file("${params.target}/${meta.id}/antismash_local") + asDir.isDirectory() && asDir.list()?.any { it.endsWith('.json') || it.endsWith('.json.gz') } + } + ANTISMASH_RUN(as_todo) + ch_versions = ch_versions.mix(ANTISMASH_RUN.out.versions) + annotate_ready_ch = ANTISMASH_RUN.out.results.map { meta, _files -> meta }.mix(as_done) + } + + if (params.run_interpro.toBoolean()) { + def ipr_todo = annotate_ready_ch.filter { meta -> + !file("${params.target}/${meta.id}/annotate_misc/iprscan.xml").exists() + } + def ipr_done = annotate_ready_ch.filter { meta -> + file("${params.target}/${meta.id}/annotate_misc/iprscan.xml").exists() + } + INTERPROSCAN_RUN(ipr_todo) + ch_versions = ch_versions.mix(INTERPROSCAN_RUN.out.versions) + annotate_ready_ch = INTERPROSCAN_RUN.out.results.map { meta, _xml -> meta }.mix(ipr_done) + } + + if (params.run_signalp.toBoolean()) { + def sp_todo = annotate_ready_ch.filter { meta -> + !file("${params.target}/${meta.id}/annotate_misc/signalp.results.txt").exists() + } + def sp_done = annotate_ready_ch.filter { meta -> + file("${params.target}/${meta.id}/annotate_misc/signalp.results.txt").exists() + } + SIGNALP_RUN(sp_todo) + ch_versions = ch_versions.mix(SIGNALP_RUN.out.versions) + annotate_ready_ch = SIGNALP_RUN.out.results.map { meta, _txt -> meta }.mix(sp_done) + } + + if (params.run_update.toBoolean()) { + if (params.run_sra_fetch.toBoolean()) { + def upd_input = ch_predict_meta + .map { meta -> + def species_tag = meta.species.replaceAll(/\s+/, '_') + tuple(species_tag, meta) + } + .combine(ch_reads, by: 0) + .map { _st, meta, r1, r2, _se -> tuple(meta, r1, r2) } + def upd_todo = upd_input.filter { meta, _r1, _r2 -> + _gbkResult("${params.target}/${meta.id}/update_results", meta.id as String) == null + } + def upd_done_signal = upd_input + .filter { meta, _r1, _r2 -> + _gbkResult("${params.target}/${meta.id}/update_results", meta.id as String) != null + } + .map { meta, _r1, _r2 -> tuple(meta.id, 'upd') } + FUNANNOTATE_UPDATE(upd_todo) + def upd_signal = FUNANNOTATE_UPDATE.out + .map { meta -> tuple(meta.id, 'upd') } + .mix(upd_done_signal) + annotate_ready_ch = annotate_ready_ch + .map { meta -> tuple(meta.id, meta) } + .join(upd_signal) + .map { _id, meta, _flag -> meta } + } else { + log.warn "run_update=true but run_sra_fetch=false; funannotate update skipped (no reads available)" + } + } + + if (params.run_annotate.toBoolean()) { + FUNANNOTATE_ANNOTATE(annotate_ready_ch.filter { meta -> + _gbkResult("${params.target}/${meta.id}/annotate_results", meta.id as String) == null + }) + } + + emit: + versions = ch_versions +}