Segmentation fault error

[blasarre]

Hello, and thanks for developing this tool. I installed laca on my university HPC following the "Installation from GitHub repository" instructions. Unfortunately, I keep hitting a Segmentation fault error, both when analyzing my own data and when using the test data you provided in raw.fastq.gz. I've included the log text from the run with your test data below for reference. Do you have any insights/suggestions as to the cause of this error and how to remedy it?

Thanks for your help,
Breah

�[38;20m2025-02-25 12:59:49,102 - root - INFO - LACA version: 0.3+2.g5fcefb2.dirty (laca.py:61)�[0m
�[38;20m2025-02-25 12:59:49,390 - root - DEBUG - Executing: snakemake all --directory '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2' --snakefile '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/Snakefile' --configfile '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/config.yaml' --use-conda --conda-prefix '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Database/conda_envs' --use-singularity --singularity-prefix '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Database/singularity_envs' --singularity-args '--bind /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/resources/guppy_barcoding/:/opt/ont/guppy/data/barcoding/,/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Reads'  --rerun-triggers mtime --rerun-incomplete --scheduler greedy --jobs 128 --nolock   --resources mem=717 mem_mb=734877 java_mem=610       (laca.py:105)�[0m
Config file config.yaml is extended by additional config specified via the command line.
Building DAG of jobs...
Pulling singularity image docker://genomicpariscentre/guppy:3.3.3.
Creating conda environment ../envs/laca/laca/workflow/envs/cutadapt.yaml...
Downloading and installing remote packages.
Environment for /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/rules/../envs/cutadapt.yaml created (location: Database/conda_envs/e4b3be83aa2be3b1302cbdf4cab39043_)
Creating conda environment ../envs/laca/laca/workflow/envs/isONcorCon.yaml...
Downloading and installing remote packages.
Environment for /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/rules/../envs/isONcorCon.yaml created (location: Database/conda_envs/aefa100ed5426368adafb85ee09974d6_)
Creating conda environment ../envs/laca/laca/workflow/envs/mmseqs2.yaml...
Downloading and installing remote packages.
Environment for /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/rules/../envs/mmseqs2.yaml created (location: Database/conda_envs/d0b0384d02d03d94aee67487f393ebd1_)
Creating conda environment ../envs/laca/laca/workflow/envs/q2plugs.yaml...
Downloading and installing remote packages.
Environment for /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/rules/../envs/q2plugs.yaml created (location: Database/conda_envs/6eae29c2b8398a61e8ca33b90dd729dc_)
Using shell: /usr/bin/bash
Provided cores: 128
Rules claiming more threads will be scaled down.
Provided resources: mem=717, mem_mb=734877, java_mem=610
Job stats:
job                      count    min threads    max threads
---------------------  -------  -------------  -------------
all                          1              1              1
check_primers_repseqs        1              2              2
cls_isONclust                1              1              1
cls_kmerCon                  1              1              1
cls_meshclust                1              1              1
col_q2blast_batch            1              1              1
collect_consensus            1              1              1
combine_cls                  1              1              1
combine_fastq                1              1              1
count_matrix                 1              1              1
demux_check                  1              1              1
drep_consensus               1              6              6
exclude_empty_fqs            1              1              1
get_taxonomy                 1              1              1
get_tree                     1              1              1
guppy                        1              6              6
isONclust                    1              6              6
matrix_seqid                 1              1              1
q2_fasttree                  1              6              6
q2_repseqs                   1              1              1
q2export_tree                1              1              1
rename_drep_seqs             1              1              1
repseqs_split                1              1              1
total                       23              1              6

Select jobs to execute...

[Tue Feb 25 13:04:24 2025]
localrule guppy:
    output: demux_guppy
    log: logs/demultiplex/guppy.log
    jobid: 9
    benchmark: benchmarks/demultiplex/guppy.txt
    reason: Missing output files: demux_guppy
    threads: 6
    resources: tmpdir=/scratch/blasarre/6062851, mem=50

Activating singularity image /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Database/singularity_envs/be79a9f6f5e87678ce46ad686c92cb19.simg
INFO:    Converting SIF file to temporary sandbox...
ONT Guppy barcoding software version 3.3.3+fa743a6
input path:         /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Reads
save path:          demux_guppy
arrangement files:  barcode_arrs_16S-GXO192.cfg 
min. score front:   60
min. score rear:    60


Found 1 fastq files.

0%   10   20   30   40   50   60   70   80   90   100%
|----|----|----|----|----|----|----|----|----|----|
***************************************************
Done in 182432 ms.
INFO:    Cleaning up image...
[Tue Feb 25 13:07:28 2025]
Finished job 9.
1 of 23 steps (4%) done
Select jobs to execute...

[Tue Feb 25 13:07:28 2025]
localcheckpoint demux_check:
    input: demux_guppy
    output: demultiplexed
    log: logs/demultiplex/check.log
    jobid: 8
    benchmark: benchmarks/demultiplex/check.txt
    reason: Missing output files: demultiplexed; Input files updated by another job: demux_guppy
    resources: tmpdir=/scratch/blasarre/6062851
DAG of jobs will be updated after completion.

[Tue Feb 25 13:07:29 2025]
Finished job 8.
2 of 23 steps (9%) done
Creating conda environment ../envs/laca/laca/workflow/envs/yacrd.yaml...
Downloading and installing remote packages.
Environment for /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/rules/../envs/yacrd.yaml created (location: Database/conda_envs/0f356a0f79553ccf51ae2f747bef105c_)
Select jobs to execute...

[Tue Feb 25 13:07:31 2025]
localrule collect_fastq:
    input: demultiplexed/BRK13
    output: qc/BRK13.fastq
    log: logs/demultiplex/collect_fastq/BRK13.log
    jobid: 40
    benchmark: benchmarks/demultiplex/collect_fastq/BRK13.txt
    reason: Missing output files: qc/BRK13.fastq; Updated input files: demultiplexed/BRK13
    wildcards: barcode=BRK13
    resources: tmpdir=/scratch/blasarre/6062851

[Tue Feb 25 13:07:31 2025]
Finished job 40.
3 of 30 steps (10%) done
Select jobs to execute...

[Tue Feb 25 13:07:31 2025]
rule check_primers:
    input: qc/BRK13.fastq
    output: qc/primers_passed/BRK13F.fastq, qc/primers_unpassed/BRK13F.fastq
    log: logs/qc/check_primersF/BRK13.log
    jobid: 39
    benchmark: benchmarks/qc/check_primersF/BRK13.txt
    reason: Missing output files: qc/primers_passed/BRK13F.fastq, qc/primers_unpassed/BRK13F.fastq; Input files updated by another job: qc/BRK13.fastq
    wildcards: barcode=BRK13
    threads: 6
    resources: tmpdir=/scratch/blasarre/6062851, mem=10, time=1

Activating conda environment: Database/conda_envs/e4b3be83aa2be3b1302cbdf4cab39043_
[Tue Feb 25 13:07:34 2025]
Finished job 39.
4 of 30 steps (13%) done
Removing temporary output qc/BRK13.fastq.
Select jobs to execute...

[Tue Feb 25 13:07:34 2025]
rule check_primersR:
    input: qc/primers_unpassed/BRK13F.fastq
    output: qc/primers_passed/BRK13R.fastq, qc/primers_unpassed/BRK13.fastq
    log: logs/qc/check_primersR/BRK13.log
    jobid: 41
    benchmark: benchmarks/qc/check_primersR/BRK13.txt
    reason: Missing output files: qc/primers_passed/BRK13R.fastq; Input files updated by another job: qc/primers_unpassed/BRK13F.fastq
    wildcards: barcode=BRK13
    threads: 6
    resources: tmpdir=/scratch/blasarre/6062851, mem=10, time=1

Activating conda environment: Database/conda_envs/e4b3be83aa2be3b1302cbdf4cab39043_
[Tue Feb 25 13:07:35 2025]
Finished job 41.
5 of 30 steps (17%) done
Removing temporary output qc/primers_unpassed/BRK13F.fastq.
Removing temporary output qc/primers_unpassed/BRK13.fastq.
Select jobs to execute...

[Tue Feb 25 13:07:35 2025]
rule revcomp_fq_combine:
    input: qc/primers_passed/BRK13F.fastq, qc/primers_passed/BRK13R.fastq
    output: qc/primers_passed/BRK13R_revcomp.fastq, qc/primers_passed/BRK13.fastq
    jobid: 38
    reason: Missing output files: qc/primers_passed/BRK13.fastq; Input files updated by another job: qc/primers_passed/BRK13F.fastq, qc/primers_passed/BRK13R.fastq
    wildcards: barcode=BRK13
    threads: 2
    resources: tmpdir=/scratch/blasarre/6062851, mem=10, time=1

[Tue Feb 25 13:07:35 2025]
Finished job 38.
6 of 30 steps (20%) done
Removing temporary output qc/primers_passed/BRK13F.fastq.
Removing temporary output qc/primers_passed/BRK13R.fastq.
Removing temporary output qc/primers_passed/BRK13R_revcomp.fastq.
Select jobs to execute...

[Tue Feb 25 13:07:35 2025]
rule minimap2ava_yacrd:
    input: qc/primers_passed/BRK13.fastq
    output: qc/yacrd/BRK13.paf
    log: logs/qc/yacrd/BRK13_ava.log
    jobid: 42
    benchmark: benchmarks/qc/yacrd/BRK13_ava.txt
    reason: Missing output files: qc/yacrd/BRK13.paf; Input files updated by another job: qc/primers_passed/BRK13.fastq
    wildcards: barcode=BRK13
    threads: 6
    resources: tmpdir=/scratch/blasarre/6062851, mem=50, time=1

Activating conda environment: Database/conda_envs/0f356a0f79553ccf51ae2f747bef105c_
/usr/bin/bash: line 1: 580155 Segmentation fault      (core dumped) minimap2 -x ava-ont -g 500 -f 10000 -t 6 qc/primers_passed/BRK13.fastq qc/primers_passed/BRK13.fastq > qc/yacrd/BRK13.paf 2> logs/qc/yacrd/BRK13_ava.log
[Tue Feb 25 13:07:35 2025]
Error in rule minimap2ava_yacrd:
    jobid: 42
    input: qc/primers_passed/BRK13.fastq
    output: qc/yacrd/BRK13.paf
    log: logs/qc/yacrd/BRK13_ava.log (check log file(s) for error details)
    conda-env: /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Database/conda_envs/0f356a0f79553ccf51ae2f747bef105c_
    shell:
        minimap2 -x ava-ont -g 500 -f 10000 -t 6 qc/primers_passed/BRK13.fastq qc/primers_passed/BRK13.fastq > qc/yacrd/BRK13.paf 2> logs/qc/yacrd/BRK13_ava.log
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)

Removing output files of failed job minimap2ava_yacrd since they might be corrupted:
qc/yacrd/BRK13.paf
Shutting down, this might take some time.
Exiting because a job execution failed. Look above for error message
Complete log: .snakemake/log/2025-02-25T125949.996297.snakemake.log
�[31;1m2025-02-25 13:07:35,704 - root - CRITICAL - Command 'snakemake all --directory '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2' --snakefile '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/Snakefile' --configfile '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/config.yaml' --use-conda --conda-prefix '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Database/conda_envs' --use-singularity --singularity-prefix '/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Database/singularity_envs' --singularity-args '--bind /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/envs/laca/laca/workflow/resources/guppy_barcoding/:/opt/ont/guppy/data/barcoding/,/lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_2/Reads'  --rerun-triggers mtime --rerun-incomplete --scheduler greedy --jobs 128 --nolock   --resources mem=717 mem_mb=734877 java_mem=610      ' returned non-zero exit status 1. (laca.py:113)�[0m

[blasarre]

I was able to get laca working with a few modifications of the yaml files.

First, I updated the yacrd.yaml file to specify minimap2 v2.18:

name: yacrd
channels:
  - conda-forge
  - bioconda
  - defaults
dependencies:
  - yacrd
  - minimap2=2.18

This fixed the issue I reported yesterday, but I hit a new issue with medaka (despite using medaka v1.12.1):

TF_CPP_MIN_LOG_LEVEL is set to '3'
Checking program versions
This is medaka 1.12.1
samtools: error while loading shared libraries: libncurses.so.5: cannot open shared object file: No such file or directory

From what I could find online, this appears to be a common issue having to do with the version of samtools and/or channel priorities for installation. Based on this information, I modified the medaka.yaml file to prioritize conda-forge and to specify samtools v1.11:

name: medaka
channels:
  - conda-forge
  - bioconda
  - defaults
dependencies:
  - python=3.8
  - minimap2
  - samtools=1.11
  - bcftools
  - pyabpoa
  - cudatoolkit=11.2
  - cudnn=8.1
  - pip:
    - medaka>=1.8, <2.0

I ran laca successfully using these modified yaml files on the raw.fastq.gz test data and got the expected 4 OTUs, so I think it is running correctly, but please let me know if you anticipate any issues with the validity of the laca output using these modifications.

Thanks,
Breah

[yanhui09]

Hi Breah,

Glad that you fixed the problem. It shall be the annoying problem of software inconsistency. If you get the expected OTUs, it shall be good with the installation, etc.

[blasarre]

Although the analysis did run to completion, and did give me four OTUs, the OTU order and counts aren't an exact match for what was reported in an earlier issue.

Is it okay and/or expected for the OTU order to vary? And can you tell me what the count # for each OTU is supposed to be?

Here are the contents of my output files for reference:

taxonomy.tsv

OTU_1	k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Pseudomonadales;f__Pseudomonadaceae;g__Pseudomonas;s__Pseudomonas_aeruginosa
OTU_2	k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Salmonella;s__Salmonella_enterica
OTU_3	k__Bacteria;p__Firmicutes;c__Bacilli;o__Lactobacillales;f__Listeriaceae;g__Listeria;s__Listeria_monocytogenes
OTU_4	k__Bacteria;p__Firmicutes;c__Bacilli;o__Lactobacillales;f__Lactobacillaceae;g__Limosilactobacillus;s__Lactobacillus_fermentum

count_matrix.tsv

#OTU ID	BRK13
OTU_1	30
OTU_2	298
OTU_3	785
OTU_4	72

rep_seqs.fasta

>OTU_1
ATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTAAGTTGGGAGGAAGGGCAGTAAGTTAATACCTTGCTGTTTTGACGTTACCAACAGAATAAGCACCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCAGCAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTACTGAGCTAGAGTACGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCGTAAACGATGTCGACTAGCCGTTGGGATCCTTGAGATCTTAGTGGCGCAGCTAACGCGATAAGTCGACCGCCTGGGGAGTACGGCCGCAAGGTTAAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGAAGCAACGCGAAGAACCTTACCTGGCCTTGACATGCTGAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTCAGACACAGGTGCTGCATGGCTGTCGTCAGCTCGTGTCGTGAGATGTTGGGTTAAGTCCCGTAACGAGCGCAACCCTTGTCCTTAGTTACCAGCACCTCGGGTGGGCACTCTAAGGAGACTGCCGGTGACAAACCGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACGGCCAGGGCTACACACGTGCTACAATGGTCGGTACAAAGGGTTGCCAAGCCGCGAGGTGGAGCTAATCCCATAAAACCGATCGTAGTCCGGATCGCAGTCTGCAACTCGACTGCGTGAAGTCGGAATCGCTAGTAATCGTGAATCAGAATGTCACGGTGAATACGTTCCCGGGCCTTGTACACACCGCCCGTC
>OTU_2
CACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTATGAAGAAGGCCTTCGGGTTGTAAAGTACTTTCAGCGGGGAGGAAGGTGTTGTGGTTAATACCTTCGCGCAATTGACGTTACCCGCAGAAGAAGCACCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAGGCTTGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCGTAAACGATGTCTACTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTAACGCGTTAAGTCGACCGCCTGGGGAGTACGGCCGCAAGGTTAAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGCAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAAGTTTCCAGAGATGAGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGTGTTGTGAAATGTTGGGTTAAGTCCCGCAACGAGCGCAACCCTTATCCTTTGTTGCCAGCGGTTCGGCCGGGAACTCAAAGGAGACTGCCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACGACCAGGGCTACACACGTGCTACAATGGCGCATACAAAGAGAAGCGACCTCGCGAGAGCAAGCGGACCTCATAAAGTGCGTCGTAGTCCGGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTGGATCAGAATGCCACGGTGAATACGTTCCCGGGCCTTGTACACACCGCCCGTC
>OTU_3
GAGGCAGCAGTAGGGAATCTTCCGCAATGGACGAAAGTCTGACGGAGCAACGCCGCGTGTAGTGAAGAAGGTTTTCGGATCGTAAAGCTCTGTTGTTAGAGAAGAACAAGGATACGAGTAACTGCTTGTCCCTTGACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGCGCAGGCGGTTTTTAAGTCTGATGTGAAAGCCCCCGGCTCAACCGGGGAGGGTCATTGGAAACTGGAAGACTTGAGTGCAGAAGAGGAGAGTGGAATTCCACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAGGCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCGTAAACGATGAGTGCTAAGTGTTAGGGGGTTTCCGCCCCTTAGTGCTGCAGCTAACGCATTAAGCACTCCGCCTGGGGAGTACGACCGCAAGGTTGAAACTCAAAGGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGAAGCAACGCGAAGAACCTTACCAGGTCTTGACATCCTTTGACCACTCTAGAGATAGAGCTTTCCCTTCGGGGACAAAGTGACAGGTGGTGCATGGTTGTCGTCAGCTCGTGTCGTGAGATGTTGGGTTAAGTCCCGCAACGAGCGCAACCCTTGATTTTAGTTGCCAGCATTTAGTTGGGCACTCTAAGGTGACTGCCGGTGACAAACCGGAGGAAGGTGGGGATGACGTCAAATCATCATGCCCCTTATGACCTGGGCTACACACGTGCTACAATGGATAGTACAAAGGGTCGCGAAACCGCGAGGTGAAGCTAATCCCATAAAACTGTTCTCAGTTCGGATTGTAGGCTGCAACTCGCCTACATGAAGCCGGAATCGCTAGTAATCGTGGATCAGCATGCCACGGTGAATACGTTCCCGGGCCTTGTACACACCGCCCGTC
>OTU_4
AATCTTCCACAATGGGCGCAAGCCTGATGGAGCAACACCGCGTGAGTGAAGAAGGGTTTCGGCTCGTAAAGCTCTGTTGTTAAAGAAGAACACGTATGAGAGTAACTGTTCATACGTTGACGGTATTTAACCAGAAAGTCACGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGTGGCAAGCGTTATCCGGATTTATTGGGCGTAAAGAGAGTGCAGGCGGTTTTCTAAGTCTGATGTGAAAGCCTTCGGCTTAACCGGAGAAGTGCATCGGAAACTGGATAACTTGAGTGCAGAAGAGGGTAGTGGAACTCCATGTGTAGCGGTGGAATGCGTAGATATATGGAAGAACACCAGTGGCGAAGGCGGCTACCTGGTCTGCAACTGACGCTGAGACTCGAAAGCATGGGTAGCGAACAGGATTAGATACCCTGGTAGTCCATGCCGTAAACGATGAGTGCTAGGTGTTGGAGGGTTTCCGCCCTTCAGTGCCGGAGCTAACGCATTAAGCACTCCGCCTGGGGAGTACGACCGCAAGGTTGAAACTCAAAGGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGAAGCTACGCGAAGAACCTTACCAGGTCTTGACATCTTGCGCCAACCCTAGAGATAGGGCGTTTCCTTCGGGAACGCAATGACAGGTGGTGCATGGTCGTCGTCAGCTCGTGTCGTGAGATGTTGGGTTAAGTCCCGCAACGAGCGCAACCCTTGTTACTAGTTGCCAGCATTAAGTTGGGCACTCTAGTGAGACTGCCGGTGACAAACCGGAGGAAGGTGGGGACGACGTCAGATCATCATGCCCCTTATGACCTGGGCTACACACGTGCTACAATGGACGGTACAACGAGTCGCGAACTCGCGAGGGCAAGCAAATCTCTTAAAACCGTTCTCAGTTCGGACTGCAGGCTGCAACTCGCCTGCACGAAGTCGGAATCGCTAGTAATCGCGGATCAGCATGCCGCGGTGAATACGTTCCCGGGCCTTGTACACACCGCCCGTC

Thanks,
Breah

[yanhui09]

can I see your config file? there're several modes for LACA. The order of the OTUs has no meaning as long as the OTU sequence match the taxonomy

[blasarre]

Here's the config.yaml.

# variables to configure the workflow on a global scale
basecalled_dir: /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_3/Reads
database_dir: /lustre/hdd/LAS/gbeattie-lab/blasarre/Microbiome_ONT/laca_test_3/Database
demultiplexed_dir:
merge_runs:

threads:
  normal: 2 # common jobs
  large: 6 # threads-dependent jobs
mem:
  normal: 10
  large: 50
# runtime (h) only for cluster excution
runtime:
  default: 5
  long: 30
  simple: 1

# Pooling samples for denoise [true|false]
pool: true
##############################################################
# Demultiplex [nanopore]
demuxer: guppy # guppy | minibar 
guppy:
  barcode_kits: SQK16S-GXO192
minibar:
  args: -p 0.75 -l 160
nreads_m: 1000 # rm barcoded samples in shallow sequencing
##############################################################
# Quality control
# subsample reads if with inadequate memory
subsample: false
# filter reads by quality and length range
seqkit:
  n: 5000 # subsampling number [int]
  min_qual: 7
  min_len: 800
  max_len: 1600

# primer-pattern check
primer_check: true
# order matters
fprimer:
  27Fa: AGAGTTTGATYMTGGCTYAG
  27Fb: AGGGTTCGATTCTGGCTCAG
  338Fa: ACWCCTACGGGWGGCAGCAG
  338Fb: GACTCCTACGGGAGGCWGCAG
rprimer:
  1540R: TACGGYTACCTTGTTACGACT
  1391R: GACGGGCGGTGTGTRCA

cutadapt:
  max_errors: 0.2
  min_overlap: 6
  # https://cutadapt.readthedocs.io/en/stable/guide.html#action
  # trim|retain|none
  action: retain

# filter chimeric reads with yacrd
chimera_filt: true
yacrd:
  minimap2:
    # ont: -g 500; pacbio: -g 5000
    g: 500
  # ont|pacbio P6-C4: -c 4; pacbio Sequal: -c 3
  c: 4
  n: 0.4
##############################################################
#  global read clustering
min_cluster_size: 20
cluster:
- isONclust
- meshclust
#  - isONclust
#  - umapclust
#  - meshclust

isONclust:
  # --ont: --k13 --w20; --isoseq: --k15 --w50
  k: 13
  w: 20

umapclust:
  # To mitagate the memory bottleneck of umap, batch mode is induced. 
  # default: -1, no batch use
  # or use "auto" estimation;
  max_batch_size: -1
  umap:
    n_neighbors: 15 # Increase this for global structure
    min_dist: 0.1 # Lower this to increase density for clustering (0.0, 0.99)
  hdbscan:
    min_bin_size: 50 # Define the smallest size to bin
    min_samples: 50 # Increase this for more conservative clustering - more points as noise
    epsilon: 0.5 # Lower this for more micro-clusters

meshclust:
  max_batch_size: -1 # default: -1, no batch use
  t: #keep it empty for auto-estimation
##############################################################
# consensus mode [kmerCon|miniCon|isoCon|umiCon]
min_support_reads: 20
# kmerCon: consensus from direct kmer bins
# miniCon: refined kmer consensus based on pairwise alignment with minimap2 [orignal method in https://genome.cshlp.org/content/30/3/437]
# isoCon: IsoCon (close similarity)
# umiCon: suspend primer-check with cutadapt to preseve UMIs in QC
consensus:
- kmerCon
  #- kmerCon # the first is the final kOTUs 
  #- isoCon
  #- miniCon 
  #- umiCon

miniCon:
  max_batch_size: 5000 # -1, no batch use, may have low yeild for large clusters
  min_score_frac: 0.85

isoCon:
  run_isONcorrect: true
  # speed up with small values, e.g. 500, 1000, 
  # but may affect quality if nearest neighbors not found, default: 2**32
  # https://github.com/ksahlin/IsoCon/issues/3
  neighbor_search_depth: default
  p_value_threshold: 0.01
  # evenly paritioned to run in batch if the number of reads bigger than the set value
  # default: -1, no batch use
  max_batch_size: -1

minimap2:
  # choose by reads type; https://github.com/lh3/minimap2
  # ont: ava-ont; pb: ava-pb
  x_ava: ava-ont # [yacrd & miniCon]
  # ont: map-ont; pb-clr: map-pb, pb-css: asm20 (< v2.18), map-hifi (> v2.19)
  x_map: map-ont # [polish & quant]
  # a high threshold when mapping against redunadnt 16S sequences.
  # https://github.com/lh3/minimap2/issues/897
  f: 10000 # [yacrd & miniCon & polish & quant]
  # maximally allowed gap-compressed per-base sequence divergence as valid mapping, 
  # filter mismappings due to partial sequence similarity when creating abundance matrix with minimap2
  max_de: 0.075 # [quant]

racon:
  iter: 2
  m: 8
  x: -6
  g: -8
  w: 500

medaka:
  iter: 1
  m: r941_min_hac_g507
##############################################################
# create abundance matrix by sequence id or minimap2 mapping, the first is the final
# seqid|minimap2: seqid, recommended in pool mode; minimap2, recommended in single sample mode
quant:
- seqid
uchime: false
##############################################################
# Taxonomy assignment
classifier:
- q2blast   # the first to be the final
  #- kraken2
  #- mmseqs2

mmseqs2:
  min_id: 0.99
  c: 0.9 # the fraction of aligned residues
  taxdb: SILVA # predefined taxonomy database name in `mmseqs databases`; The provided NT database do not have taxonomy info yet.
  # make a custom seqTaxDB using local blast databases;
  blastdb_alias: refseq_rna   # alias for blast database; leave this as empty string [""] if you want to use the default database by mmseqs2
  lca: 4 #LCA Mode 1: single search LCA, 2/3: approximate 2bLCA (not supported in NT-NT match ), 4: top hit [4]

kraken2:
  # Prebuilt database: https://benlangmead.github.io/aws-indexes/k2
  prebuilt: https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20220607.tar.gz
  classify_cmd: --confidence 0.05   # use --confidence 0.05 for standard database 

q2blast:
  split_by_nseq: 20 # split by N sequence for parallel
##############################################################
# Build a de novo phylogenetic tree
# shortest amplified region in multi-primer PCR design
fprimer_max:
  338Fa: ACWCCTACGGGWGGCAGCAG
  338Fb: GACTCCTACGGGAGGCWGCAG
rprimer_min:
  1391R: GACGGGCGGTGTGTRCA

phylogeny:
- FastTree   # the first to be the final
  #- IQ-TREE
  #- RAxML

q2phylo:
  # supplementary parameters for q2_phylogen pipeline `align-to-tree-mafft-{fasttree, iqtree, raxml}`
  # https://docs.qiime2.org/2021.11/plugins/available/phylogeny/
  fasttree: ''
  iqtree: ''
  raxml: ''
##############################################################
# UMI extraction with seqkit and cutadapt
umi:
  # keep it loose for putative UMIs in tolerated error
  loc: 200 # robust location containing the adaptor-UMI-primer pattern, leaving only the start/end regions to find UMIs
  len: 15 # base length of UMI on 5' and 3' side. Their length shall be the same.
  # concatenated UMI regions on both sides, whitespaces will be trimmed automatically
  # use IUPAC code https://www.bioinformatics.org/sms/iupac.html
  pattern: NNNNN NNNNN NNNNN NNNNN NNNNN NNNNN
  # extract potential UMIs with cutadapt
  cutadapt:
    max_errors: 0.2 # -e, high error tolerance 
    min_overlap: 8 # -O
  cl_identity: 0.9 # umi cluster similarity
  # keep it conservative for trusted UMI bins
  s: 90 # Check start of read up to s bp for validated UMIs.
  e: 90 # Check end of read up to f bp for validated UMIs.
  # params to filter umi bins, inherited from 
  # https://github.com/SorenKarst/longread_umi/blob/master/scripts/nanopore_pipeline.sh
  # https://github.com/SorenKarst/longread_umi/blob/master/scripts/umi_binning.sh
  u: 3.5 # UMI match error filter
  U: 30  # UMI match error SD filter
  O: 0.20 # Min read orientation fraction
  N: 10000 # Maximum number of reads +/-
  S: 10 # UMI bin size/UMI cluster size cutoff.
# linker and primer for umi extraction
flinker: GTCTCGTGGGCTCGG
rlinker: GTCTCGTGGGCTCGG
#########################################################################################
# in silico test with Badreads
simulate:
  min_id: 90,95,97,98,99
  max_id: 95,97,98,99,100
  min_len: -1 # -1 infinite 
  max_len: -1
  n_refs: 10 # n sequences as reference
  # n depth values for the reference sequences, following its order of ref_1,ref_2,...ref_n
  # appended in fasta header, used by Badread
  depth: 1,1,1,1,1,1,1,1,1,1
  badread:
    # an absolute value (e.g. 250M) or a relative depth (e.g. 25x) 
    quantity: 10x,20x,40x,60x,80x,100x,200x,400x,800x
    identity: 87.5,97.5,5
    error_model: nanopore2020
    qscore_model: nanopore2020
    start_adapter: 100,100
    end_adapter: 100,100
    start_adapter_seq:
    end_adapter_seq:
    junk_reads: 0
    random_reads: 0
    chimeras: 1
    glitches: 0,0,0