@Donovan Parks: I ran my pipeline with this command:

nextflow run biodiversitydata-se/sbdi-phylomarkercheck -r v1.0.2 -profile docker -resume -params-file bac-ssu-r232.params.yaml

on a machine with a total of 256 virtual cores and 1 TiB RAM. The param file looked like this:

markername: 'bac-ssu-r232'
input: 'input/bac-ssu-r232.fna'
hmm: 'bac.hmm'
hmmkey: '16S_rRNA'
gtdb_metadata: '<https://data.gtdb.ecogenomic.org/releases/release232/232.0/bac120_metadata_r232.tsv.gz>'
n_per_species: '1,5,20'
outdir: 'r232-1'
max_cpus: 12
max_memory: 256.GB
non_gap_prop: 0.8
phylogeny: '<https://data.gtdb.ecogenomic.org/releases/release232/232.0/bac120_r232.tree>'
model: 'GTR+F+I+G4'

(As you can see, I didn't use all the resources available.) Where the input file was made with

make

rules like this:

ssu_all_r232.fna.gz:
        wget <https://data.gtdb.ecogenomic.org/releases/release232/232.0/genomic_files_all/$@>

arc-ssu-r232.fna: ssu_all_r232.fna.gz
        zgrep -A 1 'd__Arc' $< | grep -v '^--' > $@

bac-ssu-r232.fna: ssu_all_r232.fna.gz
        zgrep -A 1 'd__Bac' $< | grep -v '^--' > $@

I.e. I subset the official file by domain. The corresponding param file for Archaea looked like this:

markername: 'arc-ssu-r232'
input: 'input/arc-ssu-r232.fna'
hmm: 'arc.hmm'
hmmkey: '16S_rRNA'
gtdb_metadata: '<https://data.gtdb.ecogenomic.org/releases/release232/232.0/ar53_metadata_r232.tsv.gz>'
n_per_species: '1,5,20'
outdir: 'r232-1'
max_cpus: 12
non_gap_prop: 0.8
phylogeny: '<https://data.gtdb.ecogenomic.org/releases/release232/232.0/ar53_r232.tree>'
model: 'GTR+F+I+G4'

I also had a

nextflow.config

like this:

$ cat nextflow.config
process {
  memory = 120.GB
  time = 30.d
  cpus = 8

  withName: HMMER_HMMALIGN {
    memory = 60.GB
    time = 30.d
    cpus = 4
  }

  withName: HMMER_ESLALIMASK {
    memory = 60.GB
    time = 30.d
    cpus = 4
  }

  withName: SATIVA {
    cpus = 24
    memory = 240.GB
    time = 30.d
  }

  withName: EXTRACTSEQNAMES {
    cpus = 2
    time = 30.d
  }

  withName: IQTREE_OPTIMIZE {
    cpus = 12
    memory = 120.GB
    time = 30.d
  }

  withName: ALIGNINPUT {
    cpus = 12
    memory = 240.GB
    time = 30.d
  }

  withName: MASKINPUT {
    cpus = 12
    memory = 240.GB
    time = 30.d
  }
}

This increases the resources for various steps. I hope you can replicate with this info! 🙂

Hi Daniel. Thank you for the details. Can you please email me the bac.hmm and arc.hmm file you are using? I think a first goal would be to verify I can reproduce your results on my system then we can bat around ideas on how to improve computational efficiency.

The

bac.hmm

and

arc.hmm

are the Barrnap files (also attached here). (They weren't available for a while so I used local copies. Now I believe they're back in the Barrnap repo.)

Great. Thanks! I'm re-running the archaeal data now. Hopefully the observed difference was simply due to running different versions of the workflow (head of master vs 1.0.2).

That head of master is not the same as 1.0.2 is an oversight of mine I believe... I'm not always so careful when nobody's watching.

FWIW, I have a method that tries to clean marker gene data sets that can run the GTDB R232 data in a day. What is still unclear is how it compares to something like SATIVA. So one goal here is to see if I can understand the difference between this pipeline and the SBDI SATIVA pipeline. Unfortunately, my gut tells me SATIVA will be better, but we will see. Either way, it would be great to make the SBDI SATIVA pipeline more efficient since a potential good, conservative marker set would be to run both (multiple) cleaning pipelines.

Sometimes I get these errors too for various images. Usually temporary. However, I didn't run the pipeline with singularity but Docker.

I can work around this by adding the following to the nextflow.config file:

withName: '.*SATIVA' {
        container = '<docker://quay.io/biocontainers/sativa:0.9.3--py312h031d066_0>'
    }

Yes, but I'm not sure if that will work. As I said, I've had problems with Singularity. I actually had a lot of problems to get the bioconda recipe for Sativa to work in the first place and this is yet another reason to reverse engineer and rely on the individual tools instead.

Sounds good. I don't have a direct way to run Docker at the moment so hopefully this will work.

Yes, this file contains up to 20 correct sequences per species. Sounds strange. The idea is of course that it should be reproducible.

I think it is reproducible as this sequence was missing in all my previous runs as well (which used the master branch, but presumably this isn't that different). Can't see why this sequence would be filtered out since it is from s__Methanocaldococcus lauensis which only has 4 16S sequences to select from and this is a near full length 16S sequence on a long contig.

I was more thinking of the difference between my Docker run and your Singularity run that worries me plus of course that a correct sequence is filtered before entering Sativa. Prefiltering does have its problems and I'm not sure if I do the best I could, particularly when it comes to length

I'll see if I can get things running with Docker. This should remove one big variable. I think we may have "fixed" the fwd/rev issue on our end so that the sequences are all oriented correctly.

That has been my guess but I didn't trust that so I didn't remove the code