From 4ccd5479ccd93e3404ad25d5b59a0f322758cc5f Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Mon, 22 Jun 2026 15:46:13 +0200 Subject: [PATCH 01/37] Merge changes from main into develop (#703) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Víctor López <98259577+victor5lm@users.noreply.github.com> From b15805e2944bfd3d4aca1491b0b0eadc91b274d7 Mon Sep 17 00:00:00 2001 From: ESapenaVentura Date: Tue, 23 Jun 2026 13:48:35 +0200 Subject: [PATCH 02/37] Enabled explicit local access for wkhtmltopdf --- buisciii/bioinfo_doc.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/buisciii/bioinfo_doc.py b/buisciii/bioinfo_doc.py index f3387344..ae9d2f95 100644 --- a/buisciii/bioinfo_doc.py +++ b/buisciii/bioinfo_doc.py @@ -552,7 +552,10 @@ def convert_to_pdf(self, html_file): pdf_file = html_file.replace(".html", ".pdf") try: pdfkit.from_file( - html_file, output_path=pdf_file, configuration=self.config_pdfkit + html_file, + output_path=pdf_file, + configuration=self.config_pdfkit, + options={"enable-local-file-access": ""}, ) log.info(f"PDF file created successfully: {pdf_file}") except OSError as e: From abf31f16ef08d0091ae5f9fd38004ff0fa2cf618 Mon Sep 17 00:00:00 2001 From: ESapenaVentura Date: Tue, 23 Jun 2026 13:52:53 +0200 Subject: [PATCH 03/37] Updated CHANGELOG --- CHANGELOG.md | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index cf26f603..ffe292ca 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,30 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [2.3.4dev] - 2026-XX-XX : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.4dev + +### Credits + +- [Enrique Sapena](https://github.com/ESapenaVentura) + +### Template fixes and updates + +### Modules + +#### Implementation + +#### Added enhancements + +#### Fixes + +- Fixed local file access for PDF generation via wkhtmltopdf [#705](https://github.com/BU-ISCIII/buisciii-tools/pull/705) + +#### Changed + +#### Removed + +### Requirements + ## [2.3.3] - 2026-06-22 : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.3 ### Credits From 8a82362fce5fbd2fa9016ff042fb02f512dbd971 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Wed, 1 Jul 2026 10:11:12 +0200 Subject: [PATCH 04/37] Updated wkhtmltopdf required version --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index cdad4196..f1903464 100644 --- a/environment.yml +++ b/environment.yml @@ -1,7 +1,7 @@ channels: - conda-forge dependencies: -- wkhtmltopdf>=0.12.4 +- wkhtmltopdf>=0.12.6.1 - tree>=2.0.2 - pip>22.0.2 - python>3.9 From a2d7a4b352fde98b8e2dcbe9af2c7cdd6e0bb86d Mon Sep 17 00:00:00 2001 From: victor5lm Date: Wed, 1 Jul 2026 10:14:52 +0200 Subject: [PATCH 05/37] Fixed linting for github actions --- .github/workflows/python_lint.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python_lint.yml b/.github/workflows/python_lint.yml index 6a2cde04..d961f8ee 100644 --- a/.github/workflows/python_lint.yml +++ b/.github/workflows/python_lint.yml @@ -4,7 +4,7 @@ on: push: branches: "**" pull_request: - types: [opened, reopened, synchronize, closed] + types: [opened, reopened, synchronize] branches: "**" jobs: From fd1d76ee6eabd08426a33d9ff47c0fb81638e01d Mon Sep 17 00:00:00 2001 From: victor5lm Date: Wed, 1 Jul 2026 10:15:02 +0200 Subject: [PATCH 06/37] Updated assembly config for multiqc --- buisciii/templates/assembly/DOC/hpc_slurm_assembly.config | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config b/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config index 89faddc4..5561630c 100644 --- a/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config +++ b/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config @@ -221,7 +221,8 @@ process { saveAs: { filename -> filename.equals('versions.yml') ? null : filename } ] } - withName: 'MULTIQC_CUSTOM' { + withName: 'CUSTOM_MULTIQC' { + ext.args = '--interactive' publishDir = [ [ path: { "${params.outdir}/99-stats/multiqc" }, From 142ad81ffc6c5c7a08a0bdd7217b895a4a212ea3 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Wed, 1 Jul 2026 10:15:19 +0200 Subject: [PATCH 07/37] Updated assembly results lablog --- buisciii/templates/assembly/RESULTS/lablog_assembly_results | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/buisciii/templates/assembly/RESULTS/lablog_assembly_results b/buisciii/templates/assembly/RESULTS/lablog_assembly_results index 18cea185..e9f71a8a 100644 --- a/buisciii/templates/assembly/RESULTS/lablog_assembly_results +++ b/buisciii/templates/assembly/RESULTS/lablog_assembly_results @@ -6,8 +6,8 @@ mkdir -p $DELIVERY_FOLDER/assembly cd $DELIVERY_FOLDER/assembly # Links to reports -ln -s ../../../ANALYSIS/*ASSEMBLY/multiqc/multiqc_report.html . -ln -s ../../../ANALYSIS/*ASSEMBLY/multiqc/summary_assembly_metrics_mqc.csv . +ln -s ../../../ANALYSIS/*ASSEMBLY/99-stats/multiqc/multiqc_report.html . +ln -s ../../../ANALYSIS/*ASSEMBLY/99-stats/summary_assembly_metrics_mqc.csv . ln -s ../../../ANALYSIS/*ASSEMBLY/Kmerfinder/kmerfinder_summary.csv . ln -s ../../../ANALYSIS/*ASSEMBLY/03-assembly/quast/global_report/report.html quast_global_report.html From 969e2e84a531546e3776d202edbc5cf815195cb9 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Wed, 1 Jul 2026 10:20:14 +0200 Subject: [PATCH 08/37] Updated CHANGELOG.md --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index ffe292ca..5ab694fd 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,9 +9,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ### Credits - [Enrique Sapena](https://github.com/ESapenaVentura) +- [Victor Lopez](https://github.com/victor5lm) ### Template fixes and updates +- Updated assembly config for MultiQC and RESULTS lablog [#708](https://github.com/BU-ISCIII/buisciii-tools/pull/708) + ### Modules #### Implementation @@ -28,6 +31,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ### Requirements +- Updated wkhtmltopdf version requirement from 0.12.4 to 0.12.6.1 [#708](https://github.com/BU-ISCIII/buisciii-tools/pull/708) + ## [2.3.3] - 2026-06-22 : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.3 ### Credits From 1656275ac07c463edf9c9d0f5b2c6388d1801880 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 12:18:25 +0200 Subject: [PATCH 09/37] Updated exomeeb template --- .../exomeeb/ANALYSIS/ANALYSIS01_EXOME/lablog | 2 +- .../exomeeb/DOC/hpc_slurm_sarek.config | 30 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/buisciii/templates/exomeeb/ANALYSIS/ANALYSIS01_EXOME/lablog b/buisciii/templates/exomeeb/ANALYSIS/ANALYSIS01_EXOME/lablog index ccfadeab..c3e71ab2 100755 --- a/buisciii/templates/exomeeb/ANALYSIS/ANALYSIS01_EXOME/lablog +++ b/buisciii/templates/exomeeb/ANALYSIS/ANALYSIS01_EXOME/lablog @@ -38,7 +38,7 @@ nextflow run /data/ucct/bi/pipelines/nf-core-sarek/nf-core-sarek_3.4.4/3_4_4/mai --dbsnp '/data/ucct/bi/references/eukaria/homo_sapiens/hg19/1000genomes_b37/variants/20131208/dbsnp_138_mod.b37.vcf.gz' \\ --dbsnp_tbi '/data/ucct/bi/references/eukaria/homo_sapiens/hg19/1000genomes_b37/variants/20131208/dbsnp_138_mod.b37.vcf.gz.tbi' \\ --igenomes_ignore \\ - --igenomes_base '/data/ucct/bi/references/igenomes/' \\ + --igenomes_base '/data/ucct/bi/references/refgenie/alias/' \\ --tools 'haplotypecaller' \\ --joint_germline \\ -resume diff --git a/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config b/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config index a0ce9ec9..2c80c9ad 100644 --- a/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config +++ b/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config @@ -26,23 +26,23 @@ params { max_time = '48.h' genomes { 'GATK.GRCh37' { - ascat_alleles = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/G1000_alleles_hg19.zip" + ascat_alleles = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/G1000_alleles_hg19.zip" ascat_genome = 'hg19' - ascat_loci = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/G1000_loci_hg19.zip" - ascat_loci_gc = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/GC_G1000_hg19.zip" - ascat_loci_rt = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/RT_G1000_hg19.zip" - bwa = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/" - chr_dir = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/Chromosomes" - dict = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.dict" - fasta = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta" - fasta_fai = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta.fai" - germline_resource = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/af-only-gnomad.raw.sites.vcf.gz" - germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/af-only-gnomad.raw.sites.vcf.gz.tbi" + ascat_loci = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/G1000_loci_hg19.zip" + ascat_loci_gc = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/GC_G1000_hg19.zip" + ascat_loci_rt = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/RT_G1000_hg19.zip" + bwa = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/bwa_index/v20230419/" + chr_dir = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/chromosomes/v20230419/" + dict = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.dict" + fasta = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta" + fasta_fai = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta.fai" + germline_resource = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz" + germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz.tbi" intervals = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/intervals/wgs_calling_regions_Sarek.list" - known_snps = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/1000G_phase1.snps.high_confidence.b37.vcf.gz" - known_snps_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" - known_indels = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" - known_indels_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi" + known_snps = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz" + known_snps_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" + known_indels = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" + known_indels_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi" snpeff_db = 'GRCh37.87' snpeff_genome = 'GRCh37' snpeff_version = '5.1' From 308540206dd3a0918c7c391c3057d5c6ba142c98 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 12:18:35 +0200 Subject: [PATCH 10/37] Updated DROP template --- .../DATE_ANALYSIS01_DROP/AberrantExpression/config.yaml | 2 +- .../ANALYSIS/DATE_ANALYSIS01_DROP/AberrantSplicing/config.yaml | 2 +- .../DROP/ANALYSIS/DATE_ANALYSIS01_DROP/MAE/config.yaml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantExpression/config.yaml b/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantExpression/config.yaml index 36cea8e6..a9acbf59 100644 --- a/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantExpression/config.yaml +++ b/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantExpression/config.yaml @@ -9,7 +9,7 @@ sampleAnnotation: sample_annotation_AE.tsv geneAnnotation: gencode29: /data/ucct/bi/references/refgenie/alias/hg19/gtf/gencode29/gencode.v29lift37.annotation.gtf genomeAssembly: hg19 -genome: /data/ucct/bi/references/igenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/WholeGenomeFasta/genome.fa +genome: /data/ucct/bi/references/refgenie/alias/Homo_sapiens_Ensembl_GRCh37/fasta/v20230419/genome.fa # You can define multiple reference genomes in yaml format, ncbi: path/to/ncbi, ucsc: path/to/ucsc # the keywords that define the path should be in the GENOME column of the sample annotation table diff --git a/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantSplicing/config.yaml b/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantSplicing/config.yaml index 31574896..7f3cfa8f 100644 --- a/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantSplicing/config.yaml +++ b/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/AberrantSplicing/config.yaml @@ -9,7 +9,7 @@ sampleAnnotation: sample_annotation_AS.tsv geneAnnotation: gencode29: /data/ucct/bi/references/refgenie/alias/hg19/gtf/gencode29/gencode.v29lift37.annotation.gtf genomeAssembly: hg19 -genome: /data/ucct/bi/references/igenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/WholeGenomeFasta/genome.fa +genome: /data/ucct/bi/references/refgenie/alias/Homo_sapiens_Ensembl_GRCh37/fasta/v20230419/genome.fa # You can define multiple reference genomes in yaml format, ncbi: path/to/ncbi, ucsc: path/to/ucsc # the keywords that define the path should be in the GENOME column of the sample annotation table diff --git a/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/MAE/config.yaml b/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/MAE/config.yaml index 5167d9c2..9b4aa13a 100644 --- a/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/MAE/config.yaml +++ b/buisciii/templates/DROP/ANALYSIS/DATE_ANALYSIS01_DROP/MAE/config.yaml @@ -9,7 +9,7 @@ sampleAnnotation: sample_annotation_MAE.tsv geneAnnotation: gencode29: /data/ucct/bi/references/refgenie/alias/hg19/gtf/gencode29/gencode.v29lift37.annotation.gtf genomeAssembly: hg19 -genome: /data/ucct/bi/references/igenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/WholeGenomeFasta/genome.fa +genome: /data/ucct/bi/references/refgenie/alias/Homo_sapiens_Ensembl_GRCh37/fasta/v20230419/genome.fa # You can define multiple reference genomes in yaml format, ncbi: path/to/ncbi, ucsc: path/to/ucsc # the keywords that define the path should be in the GENOME column of the sample annotation table From 68647ddad3bc638374767a2b2c69cf51dd4abe6c Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 12:18:47 +0200 Subject: [PATCH 11/37] Updated exometrio template --- .../ANALYSIS/ANALYSIS01_EXOME/lablog | 2 +- .../exometrio/DOC/hpc_slurm_sarek.config | 30 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/buisciii/templates/exometrio/ANALYSIS/ANALYSIS01_EXOME/lablog b/buisciii/templates/exometrio/ANALYSIS/ANALYSIS01_EXOME/lablog index fa13243c..a05ce981 100755 --- a/buisciii/templates/exometrio/ANALYSIS/ANALYSIS01_EXOME/lablog +++ b/buisciii/templates/exometrio/ANALYSIS/ANALYSIS01_EXOME/lablog @@ -58,7 +58,7 @@ nextflow run /data/ucct/bi/pipelines/nf-core-sarek/nf-core-sarek_3.4.4/3_4_4/mai --save_reference \\ --joint_germline \\ --igenomes_ignore \\ - --igenomes_base '/data/ucct/bi/references/igenomes/' \\ + --igenomes_base '/data/ucct/bi/references/refgenie/alias/' \\ -resume EOF diff --git a/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config b/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config index a0ce9ec9..2c80c9ad 100644 --- a/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config +++ b/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config @@ -26,23 +26,23 @@ params { max_time = '48.h' genomes { 'GATK.GRCh37' { - ascat_alleles = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/G1000_alleles_hg19.zip" + ascat_alleles = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/G1000_alleles_hg19.zip" ascat_genome = 'hg19' - ascat_loci = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/G1000_loci_hg19.zip" - ascat_loci_gc = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/GC_G1000_hg19.zip" - ascat_loci_rt = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/RT_G1000_hg19.zip" - bwa = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/" - chr_dir = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/Chromosomes" - dict = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.dict" - fasta = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta" - fasta_fai = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta.fai" - germline_resource = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/af-only-gnomad.raw.sites.vcf.gz" - germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/af-only-gnomad.raw.sites.vcf.gz.tbi" + ascat_loci = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/G1000_loci_hg19.zip" + ascat_loci_gc = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/GC_G1000_hg19.zip" + ascat_loci_rt = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/RT_G1000_hg19.zip" + bwa = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/bwa_index/v20230419/" + chr_dir = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/chromosomes/v20230419/" + dict = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.dict" + fasta = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta" + fasta_fai = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta.fai" + germline_resource = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz" + germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz.tbi" intervals = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/intervals/wgs_calling_regions_Sarek.list" - known_snps = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/1000G_phase1.snps.high_confidence.b37.vcf.gz" - known_snps_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" - known_indels = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" - known_indels_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi" + known_snps = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz" + known_snps_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" + known_indels = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" + known_indels_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi" snpeff_db = 'GRCh37.87' snpeff_genome = 'GRCh37' snpeff_version = '5.1' From 863e25cc8eb3dd9eb61c94639e7c2d0fa0607d2d Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 12:18:55 +0200 Subject: [PATCH 12/37] Updated wgstrio template --- .../wgstrio/ANALYSIS/ANALYSIS01_GENOME/lablog | 2 +- .../wgstrio/DOC/hpc_slurm_sarek.config | 30 +++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/buisciii/templates/wgstrio/ANALYSIS/ANALYSIS01_GENOME/lablog b/buisciii/templates/wgstrio/ANALYSIS/ANALYSIS01_GENOME/lablog index 83af4fb5..3234ee44 100755 --- a/buisciii/templates/wgstrio/ANALYSIS/ANALYSIS01_GENOME/lablog +++ b/buisciii/templates/wgstrio/ANALYSIS/ANALYSIS01_GENOME/lablog @@ -52,7 +52,7 @@ nextflow run /data/ucct/bi/pipelines/nf-core-sarek/nf-core-sarek_3.4.4/3_4_4/mai --dbsnp '/data/ucct/bi/references/eukaria/homo_sapiens/hg19/1000genomes_b37/variants/20131208/dbsnp_138_mod.b37.vcf.gz' \\ --dbsnp_tbi '/data/ucct/bi/references/eukaria/homo_sapiens/hg19/1000genomes_b37/variants/20131208/dbsnp_138_mod.b37.vcf.gz.tbi' \\ --igenomes_ignore \\ - --igenomes_base '/data/ucct/bi/references/igenomes/' \\ + --igenomes_base '/data/ucct/bi/references/refgenie/alias/' \\ --tools 'haplotypecaller' \\ --joint_germline \\ -resume diff --git a/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config b/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config index a0ce9ec9..2c80c9ad 100644 --- a/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config +++ b/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config @@ -26,23 +26,23 @@ params { max_time = '48.h' genomes { 'GATK.GRCh37' { - ascat_alleles = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/G1000_alleles_hg19.zip" + ascat_alleles = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/G1000_alleles_hg19.zip" ascat_genome = 'hg19' - ascat_loci = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/G1000_loci_hg19.zip" - ascat_loci_gc = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/GC_G1000_hg19.zip" - ascat_loci_rt = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/ASCAT/RT_G1000_hg19.zip" - bwa = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/BWAIndex/" - chr_dir = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/Chromosomes" - dict = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.dict" - fasta = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta" - fasta_fai = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Sequence/WholeGenomeFasta/human_g1k_v37_decoy.fasta.fai" - germline_resource = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/af-only-gnomad.raw.sites.vcf.gz" - germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/af-only-gnomad.raw.sites.vcf.gz.tbi" + ascat_loci = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/G1000_loci_hg19.zip" + ascat_loci_gc = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/GC_G1000_hg19.zip" + ascat_loci_rt = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-ASCAT/v20230419/RT_G1000_hg19.zip" + bwa = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/bwa_index/v20230419/" + chr_dir = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/chromosomes/v20230419/" + dict = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.dict" + fasta = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta" + fasta_fai = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta.fai" + germline_resource = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz" + germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz.tbi" intervals = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/intervals/wgs_calling_regions_Sarek.list" - known_snps = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/1000G_phase1.snps.high_confidence.b37.vcf.gz" - known_snps_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" - known_indels = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" - known_indels_tbi = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/GATKBundle/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi" + known_snps = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz" + known_snps_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" + known_indels = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" + known_indels_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz.tbi" snpeff_db = 'GRCh37.87' snpeff_genome = 'GRCh37' snpeff_version = '5.1' From 676f7d0591025548803710432c7023579c357501 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 12:23:25 +0200 Subject: [PATCH 13/37] Updated missing paths --- buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config | 2 +- buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config | 2 +- buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config b/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config index 2c80c9ad..728c1f2b 100644 --- a/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config +++ b/buisciii/templates/exomeeb/DOC/hpc_slurm_sarek.config @@ -38,7 +38,7 @@ params { fasta_fai = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta.fai" germline_resource = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz" germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz.tbi" - intervals = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/intervals/wgs_calling_regions_Sarek.list" + intervals = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-intervals/v20230419/wgs_calling_regions_Sarek.list" known_snps = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz" known_snps_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" known_indels = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" diff --git a/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config b/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config index 2c80c9ad..728c1f2b 100644 --- a/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config +++ b/buisciii/templates/exometrio/DOC/hpc_slurm_sarek.config @@ -38,7 +38,7 @@ params { fasta_fai = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta.fai" germline_resource = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz" germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz.tbi" - intervals = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/intervals/wgs_calling_regions_Sarek.list" + intervals = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-intervals/v20230419/wgs_calling_regions_Sarek.list" known_snps = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz" known_snps_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" known_indels = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" diff --git a/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config b/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config index 2c80c9ad..728c1f2b 100644 --- a/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config +++ b/buisciii/templates/wgstrio/DOC/hpc_slurm_sarek.config @@ -38,7 +38,7 @@ params { fasta_fai = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/fasta/v20230419/human_g1k_v37_decoy.fasta.fai" germline_resource = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz" germline_resource_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/af-only-gnomad.raw.sites.vcf.gz.tbi" - intervals = "${params.igenomes_base}/Homo_sapiens/GATK/GRCh37/Annotation/intervals/wgs_calling_regions_Sarek.list" + intervals = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-intervals/v20230419/wgs_calling_regions_Sarek.list" known_snps = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz" known_snps_tbi = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/1000G_phase1.snps.high_confidence.b37.vcf.gz.tbi" known_indels = "${params.igenomes_base}/Homo_sapiens_GATK_GRCh37/annotation-GATKBundle/v20230419/{1000G_phase1,Mills_and_1000G_gold_standard}.indels.b37.vcf.gz" From 8831064b332d35dbcef5a95475c41e47cfecc011 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 12:25:06 +0200 Subject: [PATCH 14/37] Updated CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5ab694fd..e1a02936 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ### Template fixes and updates - Updated assembly config for MultiQC and RESULTS lablog [#708](https://github.com/BU-ISCIII/buisciii-tools/pull/708) +- Updated iGenomes paths [#710](https://github.com/BU-ISCIII/buisciii-tools/pull/710) ### Modules From 91dcf0c7bcadea23e14179a3a110ea84a9992bee Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 13:17:05 +0200 Subject: [PATCH 15/37] Updated assembly template --- buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog | 1 + 1 file changed, 1 insertion(+) diff --git a/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog b/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog index d8f80540..740d4b88 100644 --- a/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog +++ b/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog @@ -123,6 +123,7 @@ nextflow run /data/ucct/bi/pipelines/nf-core-bacass/nf-core-bacass-2.5.0/main.nf --save_trimmed ${SAVETRIMMED} \\ --skip_kraken2 true \\ --skip_kmerfinder false \\ + --skip_busco true \\ --kmerfinderdb /data/ucct/bi/references/kmerfinder/latest/bacteria \\ --ncbi_assembly_metadata /data/ucct/bi/references/bacteria/20240626/assembly_summary_refseq.txt \\ ${PROKKA_ARGS} \\ From 84e6c68379acbc24a692772056cf23d4dfd42837 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 13:17:14 +0200 Subject: [PATCH 16/37] Updated characterization template --- .../ANALYSIS/ANALYSIS01_CHARACTERIZATION/02-ariba/lablog | 2 +- .../ANALYSIS/ANALYSIS01_CHARACTERIZATION/05-mlva/lablog | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/02-ariba/lablog b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/02-ariba/lablog index 9414e334..a8447f25 100644 --- a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/02-ariba/lablog +++ b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/02-ariba/lablog @@ -1,3 +1,3 @@ # module load singularity -cp /data/ucct/bi/references/ariba/databases.txt . +cp /data/ucct/bi/references/refgenie/alias/ariba/databases.txt . diff --git a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/05-mlva/lablog b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/05-mlva/lablog index 742bf830..3d79e7fa 100644 --- a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/05-mlva/lablog +++ b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/05-mlva/lablog @@ -6,7 +6,7 @@ mkdir MLVA_output scratch_dir=$(pwd | sed 's|/data/ucct/bi/scratch_tmp|/scratch|g') cp ../../*_ASSEMBLY/03-assembly/unicycler/*.fasta* assemblies/ gzip -d assemblies/*.fasta.gz -available_primers=$(ls /data/ucct/bi/references/MLVA/*primer* | rev | cut -d "/" -f1 | rev | cut -d "_" -f1) +available_primers=$(ls /data/ucct/bi/references/refgenie/alias/MLVA/primers/v20231004/*primer* | rev | cut -d "/" -f1 | rev | cut -d "_" -f1) echo "Available primers:" select primer in $available_primers; do @@ -18,6 +18,6 @@ select primer in $available_primers; do fi done -primer_file=$(ls /data/ucct/bi/references/MLVA/${primer}*) +primer_file=$(ls /data/ucct/bi/references/refgenie/alias/MLVA/primers/v20231004/${primer}*) echo "srun --partition short_idx --chdir ${scratch_dir} --output logs/MLVA.log --job-name MLVA python /data/ucct/bi/pipelines/mlva/MLVA_finder.py -c -i assemblies -o MLVA_output -p ${primer_file} --full-locus-name --predicted-PCR-size-table --flanking-seq 20 &" > _01_mlva.sh From a6955558b41b99bc5eac4a14a875238d746b9444 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 13:17:27 +0200 Subject: [PATCH 17/37] Updated pikavirus template --- buisciii/templates/pikavirus/DOC/hpc_slurm_pikavirus.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/buisciii/templates/pikavirus/DOC/hpc_slurm_pikavirus.config b/buisciii/templates/pikavirus/DOC/hpc_slurm_pikavirus.config index a4ac7caa..c44c3330 100644 --- a/buisciii/templates/pikavirus/DOC/hpc_slurm_pikavirus.config +++ b/buisciii/templates/pikavirus/DOC/hpc_slurm_pikavirus.config @@ -10,7 +10,7 @@ process { params { config_profile_name = 'ISCIII HPC profile' config_profile_description = 'Profile designed for the High Performance Computer in the ISCIII' - kraken2_db = "/data/ucct/bi/references/kraken/minikraken_8GB_20200312" + kraken2_db = "/data/ucct/bi/references/refgenie/alias/kraken/minikraken_8GB/v20200312" vir_ref_dir = "/data/ucct/bi/references/PikaVirus/viral_assemblies_for_pikavirus" vir_dir_repo = "/data/ucct/bi/references/PikaVirus/viral_assemblies.tsv" bact_ref_dir = "/data/ucct/bi/references/PikaVirus/bacteria_assemblies_for_pikavirus" From 46c717e5dffba40a77d1614fcd9f9469018aa3ec Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 13:17:34 +0200 Subject: [PATCH 18/37] Updated viralrecon template --- .../templates/viralrecon/DOC/viralrecon_sars_nanopore.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/buisciii/templates/viralrecon/DOC/viralrecon_sars_nanopore.config b/buisciii/templates/viralrecon/DOC/viralrecon_sars_nanopore.config index 06875ddb..648e93e3 100644 --- a/buisciii/templates/viralrecon/DOC/viralrecon_sars_nanopore.config +++ b/buisciii/templates/viralrecon/DOC/viralrecon_sars_nanopore.config @@ -10,7 +10,7 @@ process { withName: 'ARTIC_MINION' { ext.args = [ '--normalise 500', - '--scheme-directory /data/ucct/bi/references/virus/2019-nCoV/primer_schemes/', + '--scheme-directory /data/ucct/bi/references/refgenie/alias/coronaviridae/primer_schemes/', '--medaka' ].join(' ').trim() } From 9e5ba0ce653f6c0234e0f73201f3a55b99a4a4b1 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 13:21:51 +0200 Subject: [PATCH 19/37] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e1a02936..f69308dc 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Updated assembly config for MultiQC and RESULTS lablog [#708](https://github.com/BU-ISCIII/buisciii-tools/pull/708) - Updated iGenomes paths [#710](https://github.com/BU-ISCIII/buisciii-tools/pull/710) +- Updated the assembly template to skip busco and some paths in correspondence to refgenie system [#711](https://github.com/BU-ISCIII/buisciii-tools/pull/711) ### Modules From a746a3d48619adc2ea7e76686bae5aad4494e378 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 15:19:29 +0200 Subject: [PATCH 20/37] Updated CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f69308dc..10961129 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [2.3.4dev] - 2026-XX-XX : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.4dev +## [2.3.4] - 2026-07-02 : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.4 ### Credits From 2d7074c2841dbfbd68fe736536e9051fe8bb8424 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 15:19:35 +0200 Subject: [PATCH 21/37] Updated pyproject.toml --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 4e567274..f8bd066e 100755 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "buisciii-tools" -version = "2.3.3" +version = "2.3.4" dynamic = ["dependencies"] authors = [ From 9c6bcacefb43c19c06c3a3dc2d2538ac35be19a7 Mon Sep 17 00:00:00 2001 From: victor5lm Date: Thu, 2 Jul 2026 15:19:43 +0200 Subject: [PATCH 22/37] Updated __main__.py --- buisciii/__main__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/buisciii/__main__.py b/buisciii/__main__.py index 490de2b9..8052a15b 100755 --- a/buisciii/__main__.py +++ b/buisciii/__main__.py @@ -66,7 +66,7 @@ def run_buisciii(): highlight=False, ) - __version__ = "2.3.3" + __version__ = "2.3.4" stderr.print( "[grey39] BUISCIII-tools version {}".format(__version__), highlight=False ) From 72774ba0b27f4b865e59f84d3fbf0b887972e23a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Thu, 9 Jul 2026 14:37:26 +0200 Subject: [PATCH 23/37] updated bacass to version 2.6.1, no need to skip busco --- .../ANALYSIS/ANALYSIS01_ASSEMBLY/lablog | 12 ++-- .../assembly/DOC/hpc_slurm_assembly.config | 68 +++++++++++-------- 2 files changed, 45 insertions(+), 35 deletions(-) diff --git a/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog b/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog index 740d4b88..730ce832 100644 --- a/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog +++ b/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog @@ -26,7 +26,7 @@ assembly_options=("short" "long" "hybrid") print_color "Indicate the preferred assembly mode:" 'blue' select ASSEMBLY_MODE in "${assembly_options[@]}"; do if [ -n "$ASSEMBLY_MODE" ]; then - if [ $ASSEMBLY_MODE == "short" ]; then + if [ $ASSEMBLY_MODE == "short" ]; then ASSEMBLER="unicycler" elif [ "$ASSEMBLY_MODE" == "long" ] || [ "$ASSEMBLY_MODE" == "hybrid" ]; then ASSEMBLER="dragonflye" @@ -46,7 +46,7 @@ select TRIMMED in "${trim_options[@]}"; do # rename trimmed if [ "$TRIMMED" == "Yes" ] || [ "$TRIMMED" == "y" ]; then SAVETRIMMED="true" - else + else SAVETRIMMED="false" fi @@ -112,18 +112,16 @@ cat < assembly.sbatch # module load Nextflow/25.04.6 singularity export NXF_OPTS="-Xms500M -Xmx8G" -nextflow run /data/ucct/bi/pipelines/nf-core-bacass/nf-core-bacass-2.5.0/main.nf \\ +nextflow run /data/ucct/bi/pipelines/nf-core-bacass/nf-core-bacass-2.6.1/2_6_1/main.nf \\ -c ../../DOC/hpc_slurm_assembly.config \\ -profile singularity \\ --input samplesheet.csv \\ --outdir ./ \\ --assembly_type ${ASSEMBLY_MODE} \\ --assembler ${ASSEMBLER} \\ - --skip_polish true \\ + --skip_polish \\ --save_trimmed ${SAVETRIMMED} \\ - --skip_kraken2 true \\ - --skip_kmerfinder false \\ - --skip_busco true \\ + --skip_kraken2 \\ --kmerfinderdb /data/ucct/bi/references/kmerfinder/latest/bacteria \\ --ncbi_assembly_metadata /data/ucct/bi/references/bacteria/20240626/assembly_summary_refseq.txt \\ ${PROKKA_ARGS} \\ diff --git a/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config b/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config index 5561630c..fced18cc 100644 --- a/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config +++ b/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config @@ -13,7 +13,7 @@ process { queue = 'middle_idx' jobName = { "$task.name - $task.hash" } conda = null - + errorStrategy = { task.exitStatus in [140,143,137,138,104,134,139] ? 'retry' : 'finish'; task.exitStatus in [1,4,255] ? 'ignore' : 'finish' } maxRetries = 1 maxErrors = '-1' @@ -37,6 +37,7 @@ params { */ params { publish_dir_mode = 'copy' } process { + withName: '.*:.*:FASTQ_TRIM_FASTP_FASTQC:FASTQC_RAW' { maxRetries = 2 memory = {12.GB * task.attempt} @@ -118,34 +119,33 @@ process { mode: params.publish_dir_mode, ] ] - } - withName: '.*:.*:KMERFINDER_SUBWORKFLOW:FIND_DOWNLOAD_REFERENCE' { - publishDir = [ - path: { "${params.outdir}/../../REFERENCES" }, - pattern: "*.{fna,gff}.gz", - mode: params.publish_dir_mode, - saveAs: { filename -> - if (filename.equals('versions.yml')){ - null - } else { - "${refmeta.toString().replace(' ', '_')}/${filename}" - } - } - ] } - withName: '.*:.*:KMERFINDER_SUBWORKFLOW:KMERFINDER' { + withName: '.*:.*:KMERFINDER_SUMMARY_DOWNLOAD:NCBI_DATASETS_DOWNLOAD' { + publishDir = [ + path: { "${params.outdir}/../../REFERENCES" }, + mode: params.publish_dir_mode, + pattern: "*.{fna,gff}", + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + + withName: '.*:.*:KMERFINDER_SUMMARY_DOWNLOAD:KMERFINDER_KMERFINDER' { maxRetries = 2 - memory = {12.GB * task.attempt} - errorStrategy = { task.exitStatus in [1] ? 'retry' : 'finish'} + memory = {12.GB * task.attempt} + errorStrategy = { task.exitStatus in [1] ? 'retry' : 'finish'} publishDir = [ path: { "${params.outdir}/02-taxonomy_contamination/kmerfinder/${meta.id}" }, - mode: params.publish_dir_mode + mode: params.publish_dir_mode, + pattern: "*.{txt,json}", + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } ] } - withName: '.*:.*:KMERFINDER_SUBWORKFLOW:KMERFINDER_SUMMARY' { + withName: '.*:.*:KMERFINDER_SUMMARY_DOWNLOAD:KMERFINDER_SUMMARY' { publishDir = [ path: { "${params.outdir}/99-stats" }, - mode: params.publish_dir_mode + mode: params.publish_dir_mode, + pattern: "*.csv", + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } ] } withName: 'KRAKEN2|KRAKEN2_LONG' { @@ -161,9 +161,9 @@ process { path: { "${params.outdir}/03-assembly/${params.assembler}" }, mode: params.publish_dir_mode, saveAs: { filename -> - if (filename.endsWith('.scaffolds.fa.gz') || - filename.endsWith('.contigs.fasta.gz') || - filename.endsWith('.contigs.fa') || + if (filename.endsWith('.scaffolds.fa.gz') || + filename.endsWith('.contigs.fasta.gz') || + filename.endsWith('.contigs.fa') || filename.endsWith('.fasta.gz')) { "${meta.id}.fasta.gz" } else { @@ -221,13 +221,25 @@ process { saveAs: { filename -> filename.equals('versions.yml') ? null : filename } ] } - withName: 'CUSTOM_MULTIQC' { - ext.args = '--interactive' + withName: 'BUSCO_BUSCO' { publishDir = [ - [ + path: { "${params.outdir}/05-annotation/busco" }, + mode: params.publish_dir_mode, + saveAs: { filename -> + if (filename.equals('versions.yml') || filename.startsWith('busco_downloads/')) { + null + } else { + filename + } + } + ] + } + withName: 'MULTIQC_CUSTOM' { + publishDir = [ + [ path: { "${params.outdir}/99-stats/multiqc" }, mode: params.publish_dir_mode, - saveAs: { filename -> + saveAs: { filename -> if (filename.equals('versions.yml') || filename.endsWith('.csv')) { null } else { From 0559972635986c7d58efcdd6692241c80f01be81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Thu, 9 Jul 2026 14:40:27 +0200 Subject: [PATCH 24/37] updated changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 10961129..c989456f 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,12 +10,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - [Enrique Sapena](https://github.com/ESapenaVentura) - [Victor Lopez](https://github.com/victor5lm) +- [Sara Monzon](https://github.com/saramonzon) ### Template fixes and updates - Updated assembly config for MultiQC and RESULTS lablog [#708](https://github.com/BU-ISCIII/buisciii-tools/pull/708) - Updated iGenomes paths [#710](https://github.com/BU-ISCIII/buisciii-tools/pull/710) - Updated the assembly template to skip busco and some paths in correspondence to refgenie system [#711](https://github.com/BU-ISCIII/buisciii-tools/pull/711) +- Updated assembly template to nf-core/bacass 2.6.1 and restored BUSCO execution with matching config changes [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) ### Modules From 462b6265042351a002afda26d34124df8d71aaa5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 12:23:35 +0200 Subject: [PATCH 25/37] include R1 and R2 in assembly sample names for cleaning --- buisciii/templates/services.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/buisciii/templates/services.json b/buisciii/templates/services.json index b4a7a472..21f84cc6 100755 --- a/buisciii/templates/services.json +++ b/buisciii/templates/services.json @@ -9,7 +9,7 @@ "description": "nf-core/bacass: Simple bacterial assembly and annotation pipeline", "clean": { "folders":[], - "files":["01-processing/fastp/sample_name_1.fastp.fastq.gz", "01-processing/fastp/sample_name_2.fastp.fastq.gz"] + "files":["01-processing/fastp/sample_name_R1.fastp.fastq.gz", "01-processing/fastp/sample_name_R2.fastp.fastq.gz"] }, "clean_script": "", "no_copy": ["RAW", "TMP", "latest"], From fd3ebd83caa8406613c96ccd83ad00c3ca84e4a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 13:52:22 +0200 Subject: [PATCH 26/37] fix assembly Bacass output handling --- .../templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog | 7 ++++++- buisciii/templates/assembly/DOC/hpc_slurm_assembly.config | 2 +- .../templates/assembly/RESULTS/lablog_assembly_results | 4 ++-- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog b/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog index 730ce832..0048aca3 100644 --- a/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog +++ b/buisciii/templates/assembly/ANALYSIS/ANALYSIS01_ASSEMBLY/lablog @@ -57,6 +57,11 @@ select TRIMMED in "${trim_options[@]}"; do done print_color "Selected trimmed file option: $TRIMMED save trimmed" 'green' +SAVE_TRIMMED_ARG="" +if [ "$SAVETRIMMED" = "true" ]; then + SAVE_TRIMMED_ARG="--save_trimmed" +fi + # Select Prokka gram type gram_options=("+" "-" "skip") @@ -120,7 +125,7 @@ nextflow run /data/ucct/bi/pipelines/nf-core-bacass/nf-core-bacass-2.6.1/2_6_1/m --assembly_type ${ASSEMBLY_MODE} \\ --assembler ${ASSEMBLER} \\ --skip_polish \\ - --save_trimmed ${SAVETRIMMED} \\ + ${SAVE_TRIMMED_ARG} \\ --skip_kraken2 \\ --kmerfinderdb /data/ucct/bi/references/kmerfinder/latest/bacteria \\ --ncbi_assembly_metadata /data/ucct/bi/references/bacteria/20240626/assembly_summary_refseq.txt \\ diff --git a/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config b/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config index fced18cc..881b11f2 100644 --- a/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config +++ b/buisciii/templates/assembly/DOC/hpc_slurm_assembly.config @@ -234,7 +234,7 @@ process { } ] } - withName: 'MULTIQC_CUSTOM' { + withName: 'CUSTOM_MULTIQC' { publishDir = [ [ path: { "${params.outdir}/99-stats/multiqc" }, diff --git a/buisciii/templates/assembly/RESULTS/lablog_assembly_results b/buisciii/templates/assembly/RESULTS/lablog_assembly_results index e9f71a8a..6623df96 100644 --- a/buisciii/templates/assembly/RESULTS/lablog_assembly_results +++ b/buisciii/templates/assembly/RESULTS/lablog_assembly_results @@ -7,8 +7,8 @@ cd $DELIVERY_FOLDER/assembly # Links to reports ln -s ../../../ANALYSIS/*ASSEMBLY/99-stats/multiqc/multiqc_report.html . -ln -s ../../../ANALYSIS/*ASSEMBLY/99-stats/summary_assembly_metrics_mqc.csv . -ln -s ../../../ANALYSIS/*ASSEMBLY/Kmerfinder/kmerfinder_summary.csv . +ln -s ../../../ANALYSIS/*ASSEMBLY/99-stats/multiqc/summary_assembly_metrics_mqc.csv . +ln -s ../../../ANALYSIS/*ASSEMBLY/99-stats/kmerfinder_summary.csv . ln -s ../../../ANALYSIS/*ASSEMBLY/03-assembly/quast/global_report/report.html quast_global_report.html # Links to per reference reports From de91d01949789c4acf4ae2d45b98a8acbf536b8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 13:52:29 +0200 Subject: [PATCH 27/37] support Bacass numeric sample IDs downstream --- .../01-preprocessing/lablog | 38 +++++++++++++------ .../03-amrfinderplus/lablog | 36 ++++++++++++++---- .../ANALYSIS/ANALYSIS02_PLASMIDID/lablog | 18 +++++++-- .../ANALYSIS01_SNIPPY/02-preprocessing/lablog | 28 ++++++++++---- 4 files changed, 91 insertions(+), 29 deletions(-) diff --git a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/01-preprocessing/lablog b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/01-preprocessing/lablog index d176e321..9932320a 100644 --- a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/01-preprocessing/lablog +++ b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/01-preprocessing/lablog @@ -1,23 +1,41 @@ # module load singularity +# Bacass rejects all-numeric IDs and appends _S in its samplesheet. Keep the +# service-wide samples_id.txt canonical and translate only when reading Bacass +# output files. +bacass_id() { + if [[ "$1" =~ ^[0-9]+$ ]]; then + printf '%s_S\n' "$1" + else + printf '%s\n' "$1" + fi + +} + # if assembly pipeline was performed first and the trimmed sequences were saved, this should work: read -p $'\e[1;37mDid you save the trimmed reads from previous assembly pipeline? [y/N]: \e[1;38;5;220m' -n 1 answer; tput sgr0; echo if [ "$answer" == "y" ]; then echo "Creating links to trimmed reads..." - while read in; do - mkdir ${in} - cd ${in} - ln -s ../../../*/01-processing/fastp/${in}_1.fastp.fastq.gz ${in}_R1_filtered.fastq.gz - ln -s ../../../*/01-processing/fastp/${in}_2.fastp.fastq.gz ${in}_R2_filtered.fastq.gz - cd - + shopt -s nullglob + while IFS= read -r in; do + assembly_id=$(bacass_id "$in") + r1=(../../*_ASSEMBLY/01-processing/fastp/${assembly_id}_R1.fastp.fastq.gz) + r2=(../../*_ASSEMBLY/01-processing/fastp/${assembly_id}_R2.fastp.fastq.gz) + if (( ${#r1[@]} != 1 || ${#r2[@]} != 1 )); then + echo "ERROR: expected one Bacass read pair for ${in} (${assembly_id}); found R1=${#r1[@]}, R2=${#r2[@]}" >&2 + exit 1 + fi + mkdir -p "${in}" + ln -s -f "$(readlink -f "${r1[0]}")" "${in}/${in}_R1_filtered.fastq.gz" + ln -s -f "$(readlink -f "${r2[0]}")" "${in}/${in}_R2_filtered.fastq.gz" done < ../samples_id.txt echo -e "\e[32mLinks for $(cat ../samples_id.txt | wc -l) samples succesfully created.\e[0m" - - else + + else echo "Preparing _01_fastp.sh file for trimming..." - + mkdir logs scratch_dir=$(echo $(pwd) | sed 's@/data/ucct/bi/scratch_tmp/@/scratch/@g') cat ../samples_id.txt | xargs -I @@ echo "mkdir @@; srun --chdir ${scratch_dir} --mem 10G --time 1:00:00 --job-name FP.@@ --output logs/FP.@@.%j.log --partition short_idx --cpus-per-task 5 singularity exec -B ${scratch_dir}/../../../ -B /srv/fastq_repo/ /data/ucct/bi/pipelines/singularity-images/fastp:1.2.0--h43da1c4_0 fastp --in1 ${scratch_dir}/../00-reads/@@_R1.fastq.gz --in2 ${scratch_dir}/../00-reads/@@_R2.fastq.gz --thread 5 --cut_front --cut_tail --cut_mean_quality 15 --qualified_quality_phred 15 --trim_poly_x --length_required 50 --detect_adapter_for_pe --json ${scratch_dir}/@@/@@_fastp.json --html ${scratch_dir}/@@/@@_fastp.html --out1 ${scratch_dir}/@@/@@_R1_filtered.fastq.gz --out2 ${scratch_dir}/@@/@@_R2_filtered.fastq.gz --unpaired1 ${scratch_dir}/@@/@@_R1_unpaired.fastq.gz --unpaired2 ${scratch_dir}/@@/@@_R2_unpaired.fastq.gz &" > _01_fastp.sh @@ -25,5 +43,3 @@ read -p $'\e[1;37mDid you save the trimmed reads from previous assembly pipeline echo -e "\e[32mFile _01_fastp.sh ready.\e[0m" fi - - diff --git a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/03-amrfinderplus/lablog b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/03-amrfinderplus/lablog index 3670f82f..0c9d9f5a 100644 --- a/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/03-amrfinderplus/lablog +++ b/buisciii/templates/characterization/ANALYSIS/ANALYSIS01_CHARACTERIZATION/03-amrfinderplus/lablog @@ -2,9 +2,18 @@ # Only works with assemblies, prior to this you must have executed assembly pipeline +bacass_id() { + if [[ "$1" =~ ^[0-9]+$ ]]; then + printf '%s_S\n' "$1" + else + printf '%s\n' "$1" + fi +} + scratch_dir=$(echo $(pwd) | sed 's@/data/ucct/bi/scratch_tmp/@/scratch/@g') -mkdir logs +mkdir logs +: > _01_run_amrfinder.sh possible_organisms=("Acinetobacter_baumannii" "Burkholderia_cepacia" "Burkholderia_pseudomallei" "Campylobacter" "Citrobacter_freundii" "Clostridioides_difficile" "Enterobacter_asburiae" "Enterobacter_cloacae" "Enterococcus_faecalis" "Enterococcus_faecium" "Escherichia" "Klebsiella_oxytoca" "Klebsiella_pneumoniae" "Neisseria_gonorrhoeae" "Neisseria_meningitidis" "Pseudomonas_aeruginosa" "Salmonella" "Serratia_marcescens" "Staphylococcus_aureus" "Staphylococcus_pseudintermedius" "Streptococcus_agalactiae" "Streptococcus_pneumoniae" "Streptococcus_pyogenes" "Vibrio_cholerae" "Vibrio_parahaemolyticus" "Vibrio_vulnificus" "OTHER") echo @@ -13,11 +22,23 @@ PS3=$(echo -e "\n\033[1;37mSelect number:\033[0m ") select bacteria in "${possible_organisms[@]}"; do if [[ -n "$bacteria" ]]; then echo -e "\033[0;32mOrganism selected: ${bacteria}\033[0m" - if [ $bacteria = "OTHER" ]; then - cat ../samples_id.txt | while read in; do echo "srun --chdir $scratch_dir --partition middle_idx --output logs/AMRFINDER_${in}.%j.log --job-name AMRFINDER_${in} amrfinder -n $(ls ../../*ANALYSIS*ASSEMBLY/03-assembly/unicycler/${in}.fasta.gz) --name ${in} --plus -o ${in}_out.tsv &" >> _01_run_amrfinder.sh; done - else - cat ../samples_id.txt | while read in; do echo "srun --chdir $scratch_dir --partition middle_idx --output logs/AMRFINDER_${in}.%j.log --job-name AMRFINDER_${in} amrfinder -n $(ls ../../*ANALYSIS*ASSEMBLY/03-assembly/unicycler/${in}.fasta.gz) --organism ${bacteria} --name ${in} --plus -o ${in}_out.tsv &" >> _01_run_amrfinder.sh; done - fi + shopt -s nullglob + while IFS= read -r in; do + assembly_id=$(bacass_id "$in") + assemblies=(../../*ANALYSIS*ASSEMBLY/03-assembly/unicycler/${assembly_id}-unicycler.fasta.gz) + if (( ${#assemblies[@]} == 0 )); then + assemblies=(../../*ANALYSIS*ASSEMBLY/03-assembly/unicycler/${assembly_id}-unicycler.fasta) + fi + if (( ${#assemblies[@]} != 1 )); then + echo "ERROR: expected one Bacass assembly for ${in} (${assembly_id}); found ${#assemblies[@]}" >&2 + exit 1 + fi + organism_arg="" + if [ "$bacteria" != "OTHER" ]; then + organism_arg="--organism ${bacteria}" + fi + echo "srun --chdir $scratch_dir --partition middle_idx --output logs/AMRFINDER_${in}.%j.log --job-name AMRFINDER_${in} amrfinder -n ${assemblies[0]} ${organism_arg} --name ${in} --plus -o ${in}_out.tsv &" >> _01_run_amrfinder.sh + done < ../samples_id.txt break else echo -e "\n\033[0;31mInvalid input.\033[0m" @@ -27,11 +48,10 @@ done cat << "EOF" > _02_amrfinder_summary.sh #!/bin/bash output_file="amrfinder_summary.tsv" -{ +{ for file in *_out.tsv; do [ "$header_added" ] || { head -n 1 "$file"; header_added=true; } tail -n +2 "$file" done } > "$output_file" EOF - diff --git a/buisciii/templates/plasmidid/ANALYSIS/ANALYSIS02_PLASMIDID/lablog b/buisciii/templates/plasmidid/ANALYSIS/ANALYSIS02_PLASMIDID/lablog index 33c588cf..e5b6cd64 100644 --- a/buisciii/templates/plasmidid/ANALYSIS/ANALYSIS02_PLASMIDID/lablog +++ b/buisciii/templates/plasmidid/ANALYSIS/ANALYSIS02_PLASMIDID/lablog @@ -1,9 +1,21 @@ # module load singularity -mkdir logs +mkdir -p logs scratch_dir=$(echo $PWD | sed 's/\/data\/ucct\/bi\/scratch_tmp/\/scratch/g') -cat ../samples_id.txt | xargs -I @@ echo "srun --chdir $scratch_dir --output logs/PLASMIDID.@@.%j.log --job-name PLASMIDID_@@ --partition short_idx --time 05:00:00 env - PATH="$PATH" singularity exec --bind /data/ucct/bi/references/ --bind ${scratch_dir} --bind ${scratch_dir}/../../ /data/ucct/bi/pipelines/singularity-images/depot.galaxyproject.org-singularity-plasmidid-1.6.5--hdfd78af_0.img plasmidID -1 ${scratch_dir}/../*ASSEMBLY/01-processing/fastp/@@_1.fastp.fastq.gz -2 ${scratch_dir}/../*ASSEMBLY/01-processing/fastp/@@_2.fastp.fastq.gz -d /data/ucct/bi/references/plasmidID/plasmid_ddbb/20200203/20200203_plasmids.fasta -s @@ -c ${scratch_dir}/../*ASSEMBLY/03-assembly/unicycler/@@.fasta -a ${scratch_dir}/plasmidID_annotation_config_file.txt --no-trim -o ${scratch_dir} &" > _01_plasmidID.sh +bacass_id() { + if [[ "$1" =~ ^[0-9]+$ ]]; then + printf '%s_S\n' "$1" + else + printf '%s\n' "$1" + fi +} + +: > _01_plasmidID.sh +while IFS= read -r sample; do + assembly_id=$(bacass_id "$sample") + echo "srun --chdir $scratch_dir --output logs/PLASMIDID.${sample}.%j.log --job-name PLASMIDID_${sample} --partition short_idx --time 05:00:00 env - PATH=\"\$PATH\" singularity exec --bind /data/ucct/bi/references/ --bind ${scratch_dir} --bind ${scratch_dir}/../../ /data/ucct/bi/pipelines/singularity-images/depot.galaxyproject.org-singularity-plasmidid-1.6.5--hdfd78af_0.img plasmidID -1 ${scratch_dir}/../*ASSEMBLY/01-processing/fastp/${assembly_id}_R1.fastp.fastq.gz -2 ${scratch_dir}/../*ASSEMBLY/01-processing/fastp/${assembly_id}_R2.fastp.fastq.gz -d /data/ucct/bi/references/plasmidID/plasmid_ddbb/20200203/20200203_plasmids.fasta -s ${sample} -c ${scratch_dir}/../*ASSEMBLY/03-assembly/unicycler/${assembly_id}-unicycler.fasta -a ${scratch_dir}/plasmidID_annotation_config_file.txt --no-trim -o ${scratch_dir} &" >> _01_plasmidID.sh +done < ../samples_id.txt # micromamba activate outbreakinfo -echo "python3 summary_report_pid.py -i NO_GROUP" > _02_summary_table.sh \ No newline at end of file +echo "python3 summary_report_pid.py -i NO_GROUP" > _02_summary_table.sh diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/02-preprocessing/lablog b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/02-preprocessing/lablog index ff91dff4..8e38fb60 100644 --- a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/02-preprocessing/lablog +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/02-preprocessing/lablog @@ -1,16 +1,31 @@ # module load singularity +# Translate canonical numeric IDs only when consuming Bacass outputs. +bacass_id() { + if [[ "$1" =~ ^[0-9]+$ ]]; then + printf '%s_S\n' "$1" + else + printf '%s\n' "$1" + fi +} + # if assembly pipeline was performed first and the trimmed sequences were saved, this should work: read -p $'\e[1;37mDid you save the trimmed reads from previous assembly pipeline? [y/N]: \e[1;38;5;220m' -n 1 answer; tput sgr0; echo if [ "$answer" == "y" ]; then echo "Creating links to trimmed reads..." - while read in; do - mkdir ${in} - cd ${in} - ln -s ../../../*/01-processing/fastp/${in}_1.fastp.fastq.gz ${in}_R1_filtered.fastq.gz - ln -s ../../../*/01-processing/fastp/${in}_2.fastp.fastq.gz ${in}_R2_filtered.fastq.gz - cd - + shopt -s nullglob + while IFS= read -r in; do + assembly_id=$(bacass_id "$in") + r1=(../../*_ASSEMBLY/01-processing/fastp/${assembly_id}_R1.fastp.fastq.gz) + r2=(../../*_ASSEMBLY/01-processing/fastp/${assembly_id}_R2.fastp.fastq.gz) + if (( ${#r1[@]} != 1 || ${#r2[@]} != 1 )); then + echo "ERROR: expected one Bacass read pair for ${in} (${assembly_id}); found R1=${#r1[@]}, R2=${#r2[@]}" >&2 + exit 1 + fi + mkdir -p "${in}" + ln -s -f "$(readlink -f "${r1[0]}")" "${in}/${in}_R1_filtered.fastq.gz" + ln -s -f "$(readlink -f "${r2[0]}")" "${in}/${in}_R2_filtered.fastq.gz" done < ../samples_id.txt echo -e "\e[32mLinks for $(cat ../samples_id.txt | wc -l) samples succesfully created.\e[0m" @@ -24,4 +39,3 @@ read -p $'\e[1;37mDid you save the trimmed reads from previous assembly pipeline echo -e "\e[32mFile _01_fastp.sh ready.\e[0m" fi - From 7e65929c11a27edbdffc518ae467da873290e24e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 13:52:35 +0200 Subject: [PATCH 28/37] fix Snippy IQ-TREE model reuse workflow --- .../snippy/ANALYSIS/ANALYSIS01_SNIPPY/06-iqtree/lablog | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/06-iqtree/lablog b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/06-iqtree/lablog index 0f146ca3..d780b275 100644 --- a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/06-iqtree/lablog +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/06-iqtree/lablog @@ -3,9 +3,9 @@ scratch_dir=$(echo $PWD | sed 's/\/data\/ucct\/bi\/scratch_tmp/\/scratch/g') #1 -mkdir logs -echo "srun --chdir ${scratch_dir} --output logs/IQTREEMFP.%j.log --job-name IQTREEMFP --cpus-per-task 20 --mem 5G --partition short_idx --time 00:30:00 singularity exec -B ${scratch_dir}/../../../ /data/ucct/bi/pipelines/singularity-images/iqtree:2.1.4_beta--hdcc8f71_0 iqtree -s ${scratch_dir}/../05-snippy/phylo.aln -m MFP &" > _00_iqtreemfp.sh +mkdir -p logs +echo "srun --chdir ${scratch_dir} --output logs/IQTREEMFP.%j.log --job-name IQTREEMFP --cpus-per-task 20 --mem 5G --partition short_idx --time 02:30:00 singularity exec -B ${scratch_dir}/../../../ /data/ucct/bi/pipelines/singularity-images/iqtree:2.1.4_beta--hdcc8f71_0 iqtree -s ${scratch_dir}/../05-snippy/phylo.aln -m MFP &" > _00_iqtreemfp.sh #2 #bestfitmodel=$(grep Best logs/IQTREEMFP* | awk -F": " '{print $2}' | awk '{print $1}') -#echo "srun --chdir ${scratch_dir} --output logs/IQTREEFULLALIGN.%j.log --job-name IQTREEFULLALIGN --cpus-per-task 20 --mem 15G --partition short_idx --time 08:00:00 singularity exec -B ${scratch_dir}/../../../ /data/ucct/bi/pipelines/singularity-images/iqtree:2.1.4_beta--hdcc8f71_0 iqtree -s ${scratch_dir}/../05-snippy/phylo.aln -m $bestfitmodel -T 20 -B 1000 -pre phylo.iqtree.bootstrap &" > _01_iqtreeall.sh \ No newline at end of file +#echo "srun --chdir ${scratch_dir} --output logs/IQTREEFULLALIGN.%j.log --job-name IQTREEFULLALIGN --cpus-per-task 20 --mem 15G --partition short_idx --time 08:00:00 singularity exec -B ${scratch_dir}/../../../ /data/ucct/bi/pipelines/singularity-images/iqtree:2.1.4_beta--hdcc8f71_0 iqtree -s ${scratch_dir}/../05-snippy/phylo.aln -m $bestfitmodel -T 20 -B 1000 -pre phylo.iqtree.bootstrap &" > _01_iqtreeall.sh From 02ab65110ecf58b2ca1de57493b637054aa313bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 13:52:41 +0200 Subject: [PATCH 29/37] enhance outbreak summary workbook reporting --- .../RESULTS/generate_summary_outbreak.py | 546 +++++++++++++++++- .../RESULTS/summary_outbreak_template.xlsx | Bin 46139 -> 37706 bytes 2 files changed, 523 insertions(+), 23 deletions(-) diff --git a/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py b/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py index 73711d7b..5ffe6de2 100644 --- a/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py +++ b/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py @@ -1,7 +1,12 @@ import os import glob +import re +import shlex +import numpy as np import pandas as pd from openpyxl import load_workbook +from openpyxl.styles import Alignment, Border, Font, PatternFill +from openpyxl.utils import get_column_letter # ------------------------------------------------------ # File Paths @@ -11,7 +16,7 @@ xlsx_template = "summary_outbreak_template.xlsx" samples_file = "../ANALYSIS/samples_id.txt" tsv_file = glob.glob("../ANALYSIS/*_CHARACTERIZATION/99-stats/ariba_mlst_full.tsv") -csv_file = glob.glob("../ANALYSIS/*_ASSEMBLY/Kmerfinder/kmerfinder_summary.csv") +csv_file = glob.glob("../ANALYSIS/*_ASSEMBLY/99-stats/kmerfinder_summary.csv") mapping_file = glob.glob("../ANALYSIS/*_SNIPPY/99-stats/mapping_stats_summary.txt") wgs_metrics_file = glob.glob( "../ANALYSIS/*_SNIPPY/99-stats/wgs_metrics_all_filtered.txt" @@ -32,12 +37,35 @@ mlva_file = glob.glob( "../ANALYSIS/*_CHARACTERIZATION/05-mlva/MLVA_output/MLVA_analysis_assemblies.csv" ) +phylo_alignment_file = glob.glob("../ANALYSIS/*_SNIPPY/05-snippy/phylo.aln") +nonrecombinant_alignment_file = glob.glob( + "../ANALYSIS/*_SNIPPY/05-snippy/clean.core.aln" +) +snippy_commands_file = glob.glob("../ANALYSIS/*_SNIPPY/05-snippy/commands.out") def get_first_match(file_list): return file_list[0] if file_list else None +def normalize_sample_id(sample_id): + """Return the sample identifier used in ``samples_id.txt``. + + Analysis programs add pipeline-specific suffixes to the sample name. Keeping + the normalization in one place prevents otherwise valid results from being + skipped when they are matched to workbook rows. + """ + sample_id = os.path.basename(str(sample_id).strip()) + sample_id = re.sub( + r"(?:-unicycler)?\.scaffolds(?:\.(?:fa|fasta|fna)(?:\.gz)?)?$", + "", + sample_id, + flags=re.IGNORECASE, + ) + sample_id = re.sub(r"_S$", "", sample_id, flags=re.IGNORECASE) + return sample_id + + tsv_file = get_first_match(tsv_file) csv_file = get_first_match(csv_file) mapping_file = get_first_match(mapping_file) @@ -50,6 +78,9 @@ def get_first_match(file_list): amrfinder_dir = get_first_match(amrfinder_dir) plasmid_file = get_first_match(plasmid_file) mlva_file = get_first_match(mlva_file) +phylo_alignment_file = get_first_match(phylo_alignment_file) +nonrecombinant_alignment_file = get_first_match(nonrecombinant_alignment_file) +snippy_commands_file = get_first_match(snippy_commands_file) # ------------------------------------------------------ # File Reading Functions @@ -69,7 +100,7 @@ def read_samples(file_path): samples = {} with open(file_path, "r") as f: for index, line in enumerate(f, start=3): # Start at row 3 in Excel - sample_id = line.strip() + sample_id = normalize_sample_id(line) samples[sample_id] = index return samples @@ -124,7 +155,7 @@ def read_kmerfinder(csv_file): raise KeyError(f"Missing expected column: {col}") return { - str(row["sample_name"]): { + normalize_sample_id(row["sample_name"]): { "colE": row.iloc[1], "colF": row.iloc[2], "colG": row.iloc[4], @@ -292,17 +323,30 @@ def read_amrfinder_results(directory): """ resistance_dict = {} - for filename in os.listdir(directory): + for filename in sorted(os.listdir(directory)): if filename.endswith("_out.tsv"): - sample_id = filename.replace("_out.tsv", "") + sample_id = normalize_sample_id(filename.replace("_out.tsv", "")) file_path = os.path.join(directory, filename) - df = pd.read_csv(file_path, sep="\t") + df = pd.read_csv(file_path, sep="\t", keep_default_na=False) - # Check if required columns exist - if "Name" in df.columns and "Gene symbol" in df.columns: - # Extract unique resistance genes - filtered_genes = df["Gene symbol"].dropna().astype(str).unique() + # AMRFinderPlus renamed these fields in recent releases. + gene_column = next( + (column for column in ("Element symbol", "Gene symbol") if column in df.columns), + None, + ) + if gene_column: + # The output can also contain virulence/stress hits; only AMR + # elements belong in the resistance summary. + type_column = next( + (column for column in ("Type", "Element type") if column in df.columns), + None, + ) + if type_column: + df = df[df[type_column].astype(str).str.upper() == "AMR"] + filtered_genes = ( + df[gene_column].replace("", pd.NA).dropna().astype(str).unique() + ) gene_list = ", ".join(filtered_genes) resistance_dict[sample_id] = gene_list @@ -321,23 +365,26 @@ def read_amrfinderplus_resistance(directory): """ amr_data = {} - for filename in os.listdir(directory): + for filename in sorted(os.listdir(directory)): if filename.endswith("_out.tsv"): - sample_id = filename.replace("_out.tsv", "") + sample_id = normalize_sample_id(filename.replace("_out.tsv", "")) file_path = os.path.join(directory, filename) - df = pd.read_csv(file_path, sep="\t") + df = pd.read_csv(file_path, sep="\t", keep_default_na=False) # Remove unnecessary columns columns_to_remove = [ "Name", + "Protein id", "Protein identifier", + "HMM accession", "HMM id", "HMM description", ] df_filtered = df.drop( columns=[col for col in columns_to_remove if col in df.columns] ) + df_filtered = df_filtered.where(pd.notna(df_filtered), None) # Store the filtered DataFrame in the dictionary amr_data[sample_id] = df_filtered @@ -400,7 +447,7 @@ def read_quast_report(quast_file): raise KeyError(f"Missing required column: {col}") df_quast = df_quast.copy() - df_quast.loc[:, "Sample"] = df_quast["Assembly"].str.replace(".scaffolds", "") + df_quast.loc[:, "Sample"] = df_quast["Assembly"].map(normalize_sample_id) return df_quast.set_index("Sample")[required_columns].to_dict(orient="index") except Exception as e: @@ -446,9 +493,7 @@ def read_quast_per_reference(directory): ) # Limpiar el nombre de la muestra eliminando `.scaffolds` - df["Sample"] = df["Assembly"].str.replace( - ".scaffolds", "", regex=False - ) + df["Sample"] = df["Assembly"].map(normalize_sample_id) # Crear diccionario de resultados for _, row in df.iterrows(): @@ -527,6 +572,157 @@ def read_mlva_results(mlva_file): return [], {} +def read_fasta_alignment(file_path): + """Read an aligned FASTA file while preserving its sequence order.""" + names = [] + sequences = [] + current_name = None + current_sequence = [] + + with open(file_path, "r") as alignment: + for line_number, raw_line in enumerate(alignment, start=1): + line = raw_line.strip() + if not line: + continue + if line.startswith(">"): + if current_name is not None: + names.append(current_name) + sequences.append("".join(current_sequence).upper()) + current_name = line[1:].split()[0] + if not current_name: + raise ValueError(f"Empty FASTA identifier at {file_path}:{line_number}") + current_sequence = [] + else: + if current_name is None: + raise ValueError( + f"Sequence data before the first FASTA header at " + f"{file_path}:{line_number}" + ) + current_sequence.append(line) + + if current_name is not None: + names.append(current_name) + sequences.append("".join(current_sequence).upper()) + + if not names: + raise ValueError(f"No sequences found in alignment: {file_path}") + if len(names) != len(set(names)): + raise ValueError(f"Duplicate sequence identifiers in alignment: {file_path}") + + lengths = {len(sequence) for sequence in sequences} + if len(lengths) != 1: + raise ValueError( + f"Sequences have different lengths in alignment {file_path}: " + f"{sorted(lengths)}" + ) + + return names, sequences + + +def read_snippy_mapping_reference(commands_file): + """Read the actual ``--ref`` used by the generated Snippy commands. + + The lablog supplies a wildcard placeholder to ``snippy-multi``; commands.out + is therefore the authoritative record of the reference used for every run. + """ + result = { + "commands_file": commands_file or "Not found", + "references": [], + "display": "Not detected", + "notes": "Snippy commands.out was not found.", + } + if not commands_file: + return result + + references = [] + with open(commands_file, "r") as commands: + for line in commands: + try: + tokens = shlex.split(line) + except ValueError: + continue + if not tokens or os.path.basename(tokens[0]) != "snippy": + continue + for index, token in enumerate(tokens): + if token == "--ref" and index + 1 < len(tokens): + references.append(tokens[index + 1]) + elif token.startswith("--ref="): + references.append(token.split("=", 1)[1]) + + references = list(dict.fromkeys(references)) + result["references"] = references + if not references: + result["notes"] = f"No --ref argument found in {commands_file}." + return result + + descriptions = [] + resolved_files = [] + for reference in references: + if os.path.isabs(reference): + candidates = [ + reference, + re.sub( + r"^/scratch/", + "/data/ucct/bi/scratch_tmp/", + reference, + ), + ] + else: + candidates = [ + os.path.normpath(os.path.join(os.path.dirname(commands_file), reference)), + os.path.normpath(os.path.join("../REFERENCES", os.path.basename(reference))), + ] + resolved = next((path for path in candidates if os.path.isfile(path)), None) + resolved_files.append(resolved or reference) + + fasta_header = "" + if resolved: + with open(resolved, "r") as fasta: + for line in fasta: + if line.startswith(">"): + fasta_header = line[1:].strip() + break + descriptions.append(fasta_header or os.path.basename(reference)) + + result["display"] = "; ".join(descriptions) + result["notes"] = ( + f"Reference file(s): {', '.join(resolved_files)}. " + f"Detected from --ref in {commands_file}." + ) + if len(references) > 1: + result["notes"] += " Warning: multiple distinct references were detected." + return result + + +def calculate_pairwise_nucleotide_differences(sequences, chunk_size=100000): + """Calculate MEGA-style absolute differences with pairwise deletion. + + A position is compared for a pair only when both sequences contain an + unambiguous A, C, G, or T. Gaps, Ns, and other ambiguity codes are excluded + independently for each pair. + """ + sequence_array = np.vstack( + [np.frombuffer(sequence.encode("ascii"), dtype=np.uint8) for sequence in sequences] + ) + sequence_count, alignment_length = sequence_array.shape + differences = np.zeros((sequence_count, sequence_count), dtype=np.int64) + valid_bases = np.frombuffer(b"ACGT", dtype=np.uint8) + + # Work in position chunks so whole-genome alignments do not require large + # temporary arrays for every pair. + for start in range(0, alignment_length, chunk_size): + block = sequence_array[:, start : start + chunk_size] + valid = np.isin(block, valid_bases) + for first in range(sequence_count - 1): + pair_valid = valid[first + 1 :] & valid[first] + pair_differences = (block[first + 1 :] != block[first]) & pair_valid + counts = np.count_nonzero(pair_differences, axis=1) + differences[first, first + 1 :] += counts + + differences += differences.T + return differences + + # ------------------------------------------------------ # Excel Writing Functions # ------------------------------------------------------ @@ -639,9 +835,9 @@ def add_quast_stats(ws, samples, quast_dict): """ try: for sample, row_num in samples.items(): - sample_without_scaffolds = sample.replace(".scaffolds", "") - if sample_without_scaffolds in quast_dict: - quast_stats = quast_dict[sample_without_scaffolds] + normalized_sample = normalize_sample_id(sample) + if normalized_sample in quast_dict: + quast_stats = quast_dict[normalized_sample] ws[f"L{row_num}"] = quast_stats.get("# contigs (>= 1000 bp)", "NA") ws[f"N{row_num}"] = quast_stats.get("GC (%)", "NA") ws[f"P{row_num}"] = quast_stats.get("L50", "NA") @@ -780,15 +976,194 @@ def add_mlva_results(ws, mlva_headers, mlva_dict): mlva_dict (dict): Dictionary containing MLVA results per sample. """ try: - ws.append(["Sample ID"] + mlva_headers) + for column, header in enumerate(["Sample ID"] + mlva_headers, start=1): + ws.cell(1, column, header) - for sample_id, values in mlva_dict.items(): - ws.append([sample_id] + values) + for row_number, (sample_id, values) in enumerate(mlva_dict.items(), start=2): + for column, value in enumerate([sample_id] + values, start=1): + ws.cell(row_number, column, value) except Exception as e: print(f"Error adding MLVA results to xlsx: {e}") +def reset_worksheet(wb, sheet_name): + """Clear a template worksheet while retaining its name and position.""" + if sheet_name not in wb.sheetnames: + raise KeyError(f"Missing worksheet in template: {sheet_name}") + index = wb.sheetnames.index(sheet_name) + wb.remove(wb[sheet_name]) + return wb.create_sheet(sheet_name, index) + + +def add_readme_sheet(wb, entries): + """Populate the first worksheet with descriptions and data provenance.""" + ws = reset_worksheet(wb, "README") + headers = [ + "Sheet", + "Description", + "Source", + "Included samples/sequences", + "Alignment sites", + "Notes", + ] + header_fill = PatternFill("solid", fgColor="D9EAF7") + for column, header in enumerate(headers, start=1): + cell = ws.cell(1, column, header) + cell.font = Font(bold=True) + cell.fill = header_fill + cell.alignment = Alignment(wrap_text=True, vertical="top") + + for row_number, entry in enumerate(entries, start=2): + values = [ + entry.get("sheet"), + entry.get("description"), + entry.get("source"), + entry.get("count"), + entry.get("sites"), + entry.get("notes"), + ] + for column, value in enumerate(values, start=1): + ws.cell(row_number, column, value) + ws.cell(row_number, column).alignment = Alignment( + wrap_text=True, vertical="top" + ) + + widths = [30, 55, 42, 26, 18, 55] + for column, width in enumerate(widths, start=1): + ws.column_dimensions[get_column_letter(column)].width = width + ws.freeze_panes = "A2" + ws.auto_filter.ref = f"A1:F{len(entries) + 1}" + + +def apply_homogeneous_formatting(wb): + """Apply a consistent borderless, filterable layout to every worksheet.""" + header_fill = PatternFill("solid", fgColor="4472C4") + header_font = Font(bold=True, color="FFFFFF") + no_border = Border() + + for ws in wb.worksheets: + header_row = 2 if ws.title == "summary" else 1 + + # Determine the real populated range rather than using template-only + # formatting, which can extend hundreds of blank rows. + populated_cells = [ + cell + for row in ws.iter_rows() + for cell in row + if cell.value is not None + ] + if populated_cells: + last_row = max(cell.row for cell in populated_cells) + last_column = max(cell.column for cell in populated_cells) + else: + last_row = header_row + last_column = 1 + + # Remove borders from every instantiated template cell as well as all + # newly populated cells. + for row in ws.iter_rows(): + for cell in row: + cell.border = no_border + + style_rows = (1, 2) if ws.title == "summary" else (header_row,) + for row_number in style_rows: + for column in range(1, last_column + 1): + cell = ws.cell(row_number, column) + cell.fill = header_fill + cell.font = header_font + cell.alignment = Alignment( + horizontal="center", vertical="center", wrap_text=True + ) + + ws.freeze_panes = f"B{header_row + 1}" + header_has_values = any( + ws.cell(header_row, column).value is not None + for column in range(1, last_column + 1) + ) + if header_has_values: + ws.auto_filter.ref = ( + f"A{header_row}:{get_column_letter(last_column)}{last_row}" + ) + else: + ws.auto_filter.ref = None + + +def add_snp_distance_sheets( + wb, + samples, + alignment_file, + matrix_sheet, + pairs_sheet, + description, +): + """Create a symmetric distance matrix and a three-column pair table.""" + ws_matrix = reset_worksheet(wb, matrix_sheet) + ws_pairs = reset_worksheet(wb, pairs_sheet) + + if not alignment_file: + print(f"Warning: alignment not found for {description}") + return { + "source": "Not found", + "count": 0, + "sites": 0, + "missing": sorted(samples), + } + + names, sequences = read_fasta_alignment(alignment_file) + distances = calculate_pairwise_nucleotide_differences(sequences) + present_samples = {normalize_sample_id(name) for name in names} + missing_samples = sorted(set(samples) - present_samples) + site_count = len(sequences[0]) + + header_fill = PatternFill("solid", fgColor="D9EAF7") + table_row = 1 + ws_matrix.cell(table_row, 1, "Sample") + for column, name in enumerate(names, start=2): + ws_matrix.cell(table_row, column, name) + for row, name in enumerate(names, start=table_row + 1): + ws_matrix.cell(row, 1, name) + for column, value in enumerate(distances[row - table_row - 1], start=2): + ws_matrix.cell(row, column, int(value)) + for cell in ws_matrix[table_row]: + cell.font = Font(bold=True) + cell.fill = header_fill + cell.alignment = Alignment(text_rotation=90 if cell.column > 1 else 0) + for row in range(table_row + 1, table_row + 1 + len(names)): + ws_matrix.cell(row, 1).font = Font(bold=True) + ws_matrix.freeze_panes = "B2" + ws_matrix.column_dimensions["A"].width = 28 + for column in range(2, len(names) + 2): + # Leave enough room for five-digit distances without requiring users + # to resize every matrix column manually in Excel. + ws_matrix.column_dimensions[get_column_letter(column)].width = 12 + + pair_headers = ["Sample 1", "Sample 2", "Nucleotide differences"] + for column, header in enumerate(pair_headers, start=1): + cell = ws_pairs.cell(table_row, column, header) + cell.font = Font(bold=True) + cell.fill = header_fill + pair_row = table_row + 1 + for first in range(len(names) - 1): + for second in range(first + 1, len(names)): + ws_pairs.cell(pair_row, 1, names[first]) + ws_pairs.cell(pair_row, 2, names[second]) + ws_pairs.cell(pair_row, 3, int(distances[first, second])) + pair_row += 1 + ws_pairs.freeze_panes = "A2" + ws_pairs.auto_filter.ref = f"A{table_row}:C{pair_row - 1}" + ws_pairs.column_dimensions["A"].width = 28 + ws_pairs.column_dimensions["B"].width = 28 + ws_pairs.column_dimensions["C"].width = 23 + + return { + "source": alignment_file, + "count": len(names), + "sites": site_count, + "missing": missing_samples, + } + + # ------------------------------------------------------ # Main Functions # ------------------------------------------------------ @@ -814,6 +1189,7 @@ def main(): ) plasmid_data = read_plasmid_data(plasmid_file) if plasmid_file else pd.DataFrame() mlva_headers, mlva_dict = read_mlva_results(mlva_file) if mlva_file else ([], {}) + mapping_reference = read_snippy_mapping_reference(snippy_commands_file) wb = load_workbook(xlsx_template) @@ -824,6 +1200,11 @@ def main(): ws_mlva = wb["MLVA"] ws_amrfinder_resistance = wb["AMRFinderPlus Resistance result"] + if mapping_reference["display"] == "Not detected": + ws_summary["H1"] = "MAPPING: reference not detected" + else: + ws_summary["H1"] = f"MAPPING: {mapping_reference['display']}" + add_samples(ws_summary, samples) add_ariba_mlst_stats(ws_summary, samples, mlst_dict, mlst_header) add_kmerfinder_stats(ws_summary, samples, kmerfinder_dict) @@ -836,11 +1217,130 @@ def main(): add_amrfinder_results(ws_resistance, samples, amrfinder_dict) add_amrfinderplus_resistance(ws_amrfinder_resistance, samples, amr_resistance_data) + core_all_sites = add_snp_distance_sheets( + wb, + samples, + phylo_alignment_file, + "SNP core all sites", + "SNP core all sites pairs", + "Core-genome invariant and variable positions called in all included sequences", + ) + variable_no_recomb = add_snp_distance_sheets( + wb, + samples, + nonrecombinant_alignment_file, + "SNP variable no recomb", + "SNP variable no recomb pairs", + "Variable positions after removal of recombinant regions", + ) + if not plasmid_data.empty: add_plasmid_data(ws_plasmids, plasmid_data) if mlva_headers and mlva_dict: add_mlva_results(ws_mlva, mlva_headers, mlva_dict) + else: + ws_mlva["A1"] = "Sample ID" + + distance_method = ( + "Absolute nucleotide differences with pairwise deletion; only A/C/G/T " + "positions are compared, and gaps or ambiguous bases are excluded per pair." + ) + + def distance_notes(metadata): + missing = metadata.get("missing", []) + missing_note = ( + f" Absent outbreak samples: {', '.join(missing)}." + if missing + else " All outbreak samples are present." + ) + return distance_method + missing_note + + plasmid_count = ( + plasmid_data["Sample"].astype(str).nunique() + if not plasmid_data.empty and "Sample" in plasmid_data.columns + else 0 + ) + readme_entries = [ + { + "sheet": "README", + "description": "Workbook contents, sources, counts, and calculation notes.", + }, + { + "sheet": "summary", + "description": "Per-sample typing, taxonomy, mapping, variants, and assembly statistics.", + "source": f"Analysis summary files; {mapping_reference['commands_file']}", + "count": len(samples), + "notes": ( + f"Snippy mapping reference: {mapping_reference['display']}. " + f"{mapping_reference['notes']}" + ), + }, + { + "sheet": "SNP core all sites", + "description": "Symmetric SNP-distance matrix across invariant and variable core-genome positions called in all included sequences.", + "source": core_all_sites["source"], + "count": core_all_sites["count"], + "sites": core_all_sites["sites"], + "notes": distance_notes(core_all_sites), + }, + { + "sheet": "SNP core all sites pairs", + "description": "Long-format pair list corresponding to the SNP core all sites matrix.", + "source": core_all_sites["source"], + "count": core_all_sites["count"], + "sites": core_all_sites["sites"], + "notes": distance_notes(core_all_sites), + }, + { + "sheet": "SNP variable no recomb", + "description": "Symmetric SNP-distance matrix across variable positions after recombinant regions were removed.", + "source": variable_no_recomb["source"], + "count": variable_no_recomb["count"], + "sites": variable_no_recomb["sites"], + "notes": distance_notes(variable_no_recomb), + }, + { + "sheet": "SNP variable no recomb pairs", + "description": "Long-format pair list corresponding to the variable non-recombinant matrix.", + "source": variable_no_recomb["source"], + "count": variable_no_recomb["count"], + "sites": variable_no_recomb["sites"], + "notes": distance_notes(variable_no_recomb), + }, + { + "sheet": "plasmids", + "description": "Per-sample PlasmidID matches.", + "source": plasmid_file or "Not found", + "count": plasmid_count, + }, + { + "sheet": "virulence", + "description": "Per-sample virulence genes detected with ARIBA/VFDB.", + "source": virulence_file or "Not found", + "count": len(virulence_dict), + }, + { + "sheet": "Resistance result", + "description": "Per-sample antimicrobial-resistance gene summary from CARD and AMRFinderPlus.", + "source": "ARIBA/CARD and AMRFinderPlus", + "count": len(set(resistance_dic) | set(amrfinder_dict)), + }, + { + "sheet": "AMRFinderPlus Resistance result", + "description": "Detailed AMRFinderPlus element calls.", + "source": amrfinder_dir or "Not found", + "count": len(amr_resistance_data), + }, + { + "sheet": "MLVA", + "description": "Multiple-locus variable-number tandem-repeat analysis results.", + "source": mlva_file or "Not found", + "count": len(mlva_dict), + }, + ] + add_readme_sheet(wb, readme_entries) + apply_homogeneous_formatting(wb) output_xlsx = "summary_outbreak_filled.xlsx" wb.save(output_xlsx) diff --git a/buisciii/templates/characterization/RESULTS/summary_outbreak_template.xlsx b/buisciii/templates/characterization/RESULTS/summary_outbreak_template.xlsx index 9682810e09e8b9d1f3d8da5e82cb94565c95d653..f44205749ad472311a0bb476707f5f88d0c4d1ae 100644 GIT binary patch literal 37706 zcmeEv2UL?;*Dg||C`CX~siHKMF1^HtNRcX{B25HDN|4?IC?Z8r2p~1m1q75THPQtH zl-@%}dJPallKWD>aps#*=O6F-?_IN2)a1O+KD(S}@3T+xny4rc5Ypn|;gR5Z1kK*L zQgiKU2=J-__(KK!xd%2-u>xCK^Xgk#@;I2AU2D2_pq%gE{sFj}U4FPPQv$_hDc_jl zCLUX_JTCLc>vJ=NAPYywkMwaoKCJ9{FJ%^ZM$82TYRJ#GpIn!&uZ#7K9dXX!@AL3~ zcbqbTnTenx&b-6ji~peLr=k3}6RkW~jBT4jFWP5`tFl)N-8?7$j)OU^mw}9G+=)%v zk48Z(>5o;%U0UlmD z1s)y^5GDsRUYq-d=7!sU^KS;KK~o#t!$ON8LQ_6nZDW2HEc)QK&_tY+SwxOmw_xtO z2H$r?_q1^sKD)bS!+E#WQyjVOz??-_bq1R2r{F|H5I(wud!}Et-&5%w2_x0-I*#3V zhkpO~{MAb8%L-3A<_{fka77MMr{4=Wa1bg%b0+6W0ko~|;pI1fyuQD345~3N6j!Og^l`iZ2YC%!LJ9*!$K#?`$lP)=3X63KFVYX{}k97-| z2CC-gBWP$3TXQgmMb&C`SixEC1-0H9@i5 zK(ghc5ecHtRT@3P`%pIrX{*TWD@fynq;Biq%daUnvfp1xqz}%Kmz{p;G5OZdtrfBS zqE$8dadU}F!kXXbpv!V{S?X7$z9=Q1(kwWhSC6zjChK(2feUgde9roMilJo5MfTg1 ziteejsW&Wxv^YjpT)2a)Cmn07lNp{*Pfwm@m^tm0^{$bw?}XHssax!d`(oK$3rGy! zl|BuZsl2-)a45!^=kN=es%B?-FN5ncjI{2P8He^iV`aH}`Jq}z$^z)d>U2nqyb;I@ zL9Ig@mEJ6MsM>DIk1tO@fLih#te{=>9&au<BCg_*NL0Qz%3D( z??1Xi!mXdkRGau3i0MxxO=?}bi&>Yql;Wb4UzS#8`oJhPKX;Yw^c@FVd>Ymfu1M0s z>ap1@(f1kM@%J4ECx{QGTeaHJE(w4KoLqtoDq}+^;;qcMitFvhzZN!0kn=S7?9X~x zT-BdZOP^|FrxA zW|x7?W$W8O-zFCJoUu{|P%0gqKEGw=I{k( zmBf8{6kg^mVMJ%%l3#d~-sF3&=c4UHHYbZ!Wc)+ljVH{RPSZz04szpV z8ONJ{boA@&BRf=g#Q&gX27y9&x8undcMwlMJ3iOIKgJ#ESLLjqkX8Iq&1}x^6Y_P4 zSBY_~$d&alq5C2u8xb&q)l?IHqRTIAvpfg}@&oWjd=Drm9_)M-Wymxcm6-vyKgp`$ zqRIDIMA(nn+Gdtt++L%RPSQ$CugO~}?Ok7enokD$ppeD+si@4bs9Nea%1#>Aaar%L zipPVWT=AhzD?yXmG2oq35xUl%Xf-6>rgTa+qjQRPg2p9EKg#4VuA^n5#W);Cgu)3Q z4$CdgZJO$5lM1>MAn>7jf>=X8a)~lZL#1A4xP)H1OwM3Vhb=pV)9+J}2a3~!^-TPW z!f_YJOA#LnvQ&$Rm!UUWsr)WIHo`XvI!ZN4njceUegdjvyFje?O2k@Qdrsjco`tcy zeO{|8(cuHlHpe|l=35_epHO`68g)X;r0r(m#d{Y9PgPhv8B9YHaNH*t&kZ>9LaH#j znFPD49vE0bH$iVm(2^h~K9xp)JB2x4Lfpn}B0<)JEjv+`%;oGWmzn!dvuWCF9anvI z|3d~(U_W)FAMnUKgNJ9%jsIt45ZEfUecIO7jXg&eGUN2_vK#TX$H~tLi{cYdvpXJ& zvlm`Tz4QeC%-K`v*K=l2%UVulE5VJBkQU7zy1rMs(eaPf69~0$jh8izhk>*)qLB?9 zZYUckJLG(S;1a5;V`@rs&MnGSymWnX%6F1@R$2Qx=OVqh^*4HUc`8Pbq-8yD1+eRj(djJZFu7rp2a9~=m3 z1WUig&OGuP4$RKaM&Sx?GQ3@zSO};*v^6bRW^Yi`jwH>m~&dwq=y^>>a zB_oh7X0XPz10KcP*aRuW)~?Pm_ZLCt)~AiLb9MW5qGEd8%5~@5L`(98L2w6tG#^A# zLdw-)uIQPwY0cai2QGq9vC(qvv&=MV!DZncW!N}QN}`CIPa>L|spO;G z!=?GEf>|c|`IS2zFjG;&p+Y_hA-EO%mpAeBw#L6vjgcMMoG*Adwb+U zk9`doy`UATYo0FcGCn&W4M&+GFjjrDI4;~sS~Q$WZ?Pu?*OG|yJH)QbgB?+(vqbYU zU`8jeVQo;@$#I+&DF>wR2|Z+nP&Z=`xcm^Ed8eW_%uf_mr1@}hR6q#~zx@V*eNk!> zz)`a_3}(I)gDbXkw4QE&p&K%jQVwQVG!|U7#B|EqGv9oC`|FIR2mP9bdvNwtrYm?n z%x_rf2+RgM6CXbZcg|T^tm=VC-h}a%-+|a*`ufTOC*COEh(a=t!=uhVKrU`n)J~j@ zTAT14MlMc#!LDY;NJ^0r *pcp9IP?yHye_C~#W!td`b%Bc7mMg{ z90?9`<1RAu5d+T$$jemSkyXL#`Nk{Dv*A*hOQwT@&V?HygEI@serpV9i#CsiWALF@ z!>gcl<|baIbU1F*3Nt*B8Q9KiJjax-8x4o#NW(=*ja{;a{Fm(KnA6iUxv}F#JvhhW zJnSrHS}Aqhi3xGrf+ww(Pe(MapRsHa&~cKXq`UWvf22I()t`MfR{(9Q920}eulzT zpy(y;y&YD6fn!5AyS{1_Gsw`S88uuVYl>TVh;8T|!SW6=+Fobz4}4flY%mu^vLMvYdrV=P-E z56hPL$}*;lPBM=#4Y%aBbamBKEj^e3<30z_Uz4!tbs1PM;zt)M&%YcXEmQY&^#qnZ zpRdn$EhfN0>6OHLTykcFO{IdTCSfKi1!+GO1JTipbtCIGk%?3qNp^>e^H4q5fcrof zySAb;-TYS2K$r78VZgO*L(#M5!YZeUop_>{wc z_%N6X7SVK#%n&a#@Z=y=CS4(=Ac-+=(MUcSLvtS(W?TpwI9z7xdC60RP5jn04#rXbIvQv)Cz_L(vuXC!Mr@5iO#D~B z4OEj9(-QUQ;W|_3C8!9etXwrmd4dzG;=IrK#iOI7b?Q)8C;`;iJ}8F7=94w8>x=oz zuMBkJiZ41ovwFI?f0VO!?ec4EinBl!2Y#&ql!H@Qbw!mzC8kf#iS;volj^ceIgLkv zqt7DQs8HRBI(lgQsbL!Zcqdk~0UD2y-`J#5D(e`)3Z>)ZRn>_tz9M$n@sSnuo(t8s|y%XOztwCj9fKIgD%A? z$FhCGHV&fR1)J6NO%D+LL|4R6M6tC z^f+`Cde2`acbK;4(B`xn5=!qFReyA71M(=%uJ?o8s`lrd^bPVql=_FsN*}~MsJXT$>a!=h zXC?cdXyvXQk$ZN~?AdX5*N)#!Hk{@MXf59r)!Gw1uqUd&E6TlRhu5wh`Fk!L*|Xz! zla-QXKBzI@6MeQPy61xMo@n*19X)$?(C^t{uxrQfCL2!E09K)FS5#t8lz3Nkc5f9J z_w0DIYlr=w3lw{H{BE++)S?G9+IymYd!l)h%5# z2$*elbO-S z>y|+Ft&(gmZ()T{VTB~&w0z;wPNUjxquO30WPcJv@7Y_uXTy8n^!GyJBiYg-*|HC8a+(0t;*4-V$)<(Rh{8oh}aH9V;2&> z3mM;qhy#$}*$f%S#NA7iJC`~(F9FVsR_#jX?m}jFArd7H^ zhy}}CwFLevBCHYfzPm%G?rw)C(`P3_4>k5|&Dcj@u(UnLs?lt*_Lka`Scy%Ti9i3uZJTm#-|h9^QE(f&Gyjp%?dgGZ z*q+@~^#>u%N8^hheh!i#<}R@PJn2jvB0_oFYV{RcKIQs>Rm#&~mYflaiwYJeajLgnx z<@E4e`x_qqh}2?IsKV@nf|0z?pHpMtZ&7saZ-;-g{Qk)z{@sa-zpMvIr1KWEtpi#1 zr*Q>p#Z2>`Xx^;FKS-ZH+W-IF#N_!SKZ4ZD7zFkg!7-Hi8#XL;d;U}xfAL+2! zOMH{DK*8zZG2e;o&)c=N5&u{3fc+_7yL;g87Q(+UeKQktc3AO02^7xGeify6BKyx%I`~&ndWRMNJf*Mw3QBL*<{hg3 zg4+B|Ro%;@t%J|sB}_AwK^nm1h23pWYcCy-?4@JXUON6QSOR}_g6S91e9q0I>F>6N z|Ctk8ReHy|f1c8Qzlzd3k^Sc>&Hbw=y~B!sp3-^0tlik5>ObFZ4E^0y`wOmL@}B|v zjfTBE+c^`q#r5o^qwroj{ykXr(cs^_#q$rIUbY&K9qayiN}K;GO7BGW7f||FROzh_ zdWRMNJf+8e0j0NU^A1)2`Pv-$3u-etza4NWq}a=|irsXa-8&RY?4{%1gZ-C zc7w_wvxjG9@3iN!ZJjTpS{m$nw>-U|OI+e5Q z1G8(gTUyd`?Z;tFN5#viRa(_8FOIXW}f5e92);ynWI&YT5bvv)%`4asjnN98~JKgv^BZLtIvol$5YJkW?)%rLJ_Tvy0PGJ)++QvoVV8YijbAc5^T&cAY9l zt<7SRJ3yHEjhpDvWM(ETc5M)oj9#fv%Y@+87W%?)>q{8m{!2F(*VT1q%H=uS%GB@@ z0oZMQaH$C1N4XlsGk>z)YOYY7M4w&20=7fx=+G zy0{vcu4_>=TvBoaoiCCF18iPF<{RT|aFZbKCgy$!#2bgq1W8M~@XSFR56*&I))tVk zXC0aS7~}#5(wnc$;pTz~^RB7!4s7CGo-hP$`ZXQKONn#EO`s3F$64ZRkS(U@4Hs04 zHO6g*)D4mWo}0pz!Ng|FV#jTBbY5=NywUxIly#8+x*(7 z4+12$I;>n#6bJ+L;^xN|CNQqfBMrI`(AxMM@!V_-r5mivn1tk_=TjE;t8U);`OmQw z)!g@NSUEUJ*F}8q*zGvn*rNVrL*4*xQ;5BeKTM^s>JZkfb*!*owe~g7%9Hza~9U8 z$^(`h{UH(yUzJF<2S|eeO>vOXBd&!q+nVlrK4y{fs3?I1AJBCyNi32jRk;WyVY?!y zq6VE-3}B8i)qUBlNdmd@oAHBepx|NR+2#&O`~7YMSMi}#Y!`i0siCxNmugjg1;yO1 z^13eA8V`9xPn{=KWrZGPyCA130A)x#&($=tR#v)<9E|p2HVeXIyX>QK7fQ$`Q>$VH z#ZSDP(m1ka@_xB$Ci)p{|9`oxG6_uX$>ZAKHa@`d^beD(o`xPz6bl|so^2N4vR~6C zlXs_%yUfWi?->+#k5f$ELoMz)r;z-k=9t=rU?<@AW{Hg_O=Og=#E6_oR^{Ilj)Ih>60fda4m0bq=&#Mx;v&VBMv zsbecSN#%WmVn;ZMeO;gqU#7l>2j zTnc9(o>l7*d^;8Umu7OFta~gFe-<#JO2NG4L!QrftmUMs^Ur$q2EaEfxORLf_4_eA zXHPgo&8ph4@Y{#+KVlE(nYt%gu-b=jd2QVYBtYNySFp}AE6$7)4Y-YNzEbhw+5VrA z!A)-TC_wD7+~mSPBL=9_sjV!b{z;j}`pRg!pOfi^S9@!TK>wISDi5Ih6J;h;oS^#* zFNdobL5U1yx>QyK#VA%1XJUMMq<&IHJJzS@;!nz8;(X9DKPl4??-MEezflJJIxRoU zt-Vd_xY$x2oTEKx9zr#EU7a$Dv0OFHNZKhaNJ#7_NE*df4z%SJQ}SXP9sx zg$2{UGXMX?fUk45=gcF_%;Zwpw{5a4mx?@;Eb)qo%55k`qUwjeqvi-jsg9<%oaOT04{In)w`aeBch(nnS&mU7|haNM$a7$Gf zict)B{Zc;3H^u0|%iBo-DaOZNUPHo7k`!WM z(5*aSG{!q40OL&#uQ6v=|47SUrAVq?2jw4CBvP-v&flg;u3q2Dk5nXR`FM`MLGeIK zT@L?DvW?(USbJ1cSj(s3xL~e(zAx`4J?1j3duf&AnQ|}n<%6UrDTb3Tosv9E?nS&b zO7b!>?0&hTUKoy2p#MBxZ!-3fxQ1D6l*5;+?&GI!ltuJj*Yu7_jrFJ9p|JF!oa zY9D?k8o5%_jC1-g@PB5FGiUVY94P&0`Fe`bMAojWD^D_DtsJPrVD zdn)=Dn2vDi{mkliGpz)yzP-BdXE@uj^47lCTju97G0(eXSnBMrULw(@>K+9_+>F&r z$Da-O%||~Q2j1R_{wSwG_NlSf1NQqwl67vXgOy6&%Ju<+_RGNQHK{KRDRgccj*^C@ z^SChJ;ylbuBsrmrFY&)T-X?VY)g!U%?Vmp^6T1HDk=SMV=a2b>uI-1z#9Sa;jIE(s z0kX8vz-1z=vo0Y4R+7tJXhsOzOCpZQrJC#^F7A746(-|CTV=RY{BEn!{d?iBL!vbqm)9y|p5oUH#0nssf_7_Y;cISN)qgMiWk_rj zcqW0e^rG6(5DP4Mx%!>ZD1>gF=3jzvk(+;Kg=w4EjfAb z4jj$7(v2D5uA5BhUfkmHWiW0)nm4P|1=0@z?(LaJO+;@DrklUVZ8Wts;6?`*CS;;& zGgs%_N+C{t%Dn4g-rh9`+#0%wVHUVU2Gv#(wz>*q0`5tHtSt{uHBss^!Hb=7YwK;} z_OLnl!8sfnJ8=o=GK5yIYHC_+p`--iqNG9dEjXvug(Wyq6AqKi+`lok01n%5L#`@x zC)=YYQGqi{wm64WH10ia8S^RMerYCsZgD1Ew;zI9U*mOgThd)0TZ`!5SeYM1z(4~a z4D+}5R3-VH!x91CQrixEOYNV0ndZ5zd#b{|SO8yD`$V^Dk)YKhM^3BiMYo_vWX?5> zc#Do!X(bi&Vg0Ne=XHGvFP4YxCYnVA(rYfyQNwS=IWfI+RM!@pTg_p~5(b4+XQ%ua zh2uc#jCrBgD9W_cfvpXjkEIbJi1bV^)F3gzyfNB(139P1hq)p`6q_lw{-juj9*fHX1`A3)oSaW0|Sm26~ zL3G}Olm3ipIqP0C*YVo#8}~F9J}k6l9KSL!ay$k}>~crCjjzqdz;7DDu!_w<>E6J? zDI0oCrOe?avqcN%yO$PEIT30=!VQT~?3>v!C5DP_X@e0wPE(Sz~f4YEr%S;-H$@D5SE zA8(XO?xl4nA?Y1&bX6+Y8sm*Zds?aEv9eh5Mmb;D9t&$MN!Z4PHM+`ek7Y)oeYbJh z;gr+T5T)aT`_w@GTq%^Ud`G3Xr>q zdfrvk=sh<%e17V6S2wY6-JtXSkAW{qQW^_8MRhV3F&-9-Vk51*`t`41!5|6L+}V{Y zeYM@MR`Ny0W}S5(_t)hqKXCSYOP8CS=$gFlz{Mc3a#09-0+w;mz}k2&`n5M{sHXg} zj$>q&nq5fnx2aQ;(QKfIkB3)5@~5dIw3RxcJptW(R8Mr-775Hm&n9b7o5+!{m zUFrUC{*CQWE5BjBGQuZVh{do_Ez`-kNyHn`SZQaex}IBy=f4_SA4hh}ND(!HZj*i# zrAe9k7{69^B8on0GFUg{_*&Wi(PrjlCj-!UoSv$lIXxrMq*(%|UVO=9)c75E(}h%p zk51DS^Q=nY@AiYpiuR_ByMz!;veCtqSm5Gg`~zeOi~~wzt+i2V0$fz; zEU_Ov@wtdXfic{t37G2oZ^meE+PWL*V;kgZo312YJpOq$)q-*nxRtXGDMuT^vC#Oz z6s7M8d1|4HXXn|GEaet?pqxU`<(_*6<8$T^0SmPQmUE$BDZgdl4-{|Ve~b6O#rxml z{fFXh*eWx@Oo)fKNcX43TVyMBI<*5Muh0U&$Kd(g$7tBxA?G#EE}KHKX3oZqz9b#z({}W({Blt#$>nEwO^zQ1 z2cG`Ykq$zF55R_7x>^sESp>^+Am4mqF5rIaeW7EVD*5yogJ?%DaSIIPptG2%{G>Y8 zzb&^IJ*{D)lY_m@=puQ2CT%QyQ3tkkDzD>!>s&hSz74iW{CHf@y(v(m3vYMiTBK9f zVJ(~teSC`(%&vo7z;5Q|SIG5tv-50qG3NQ99mzd9GX)f8P8uES`|?K{m?ETH-$tIr zS!4&=VMFS0b25)+^ZFsBeHJ>P=PQnmKJqd)i<&Ncu_YM1?OR-y5J1=98F0OJ#(;PF2ULQ>@J!XKU zG{;o%>OlJ1Gb3eF3gVHEE}f*A5%wx|lF=}U)MTElX=?A8>Z)OBOx4XaR=*iw9@`}) zYZCMNnx0wq1Gdh5ufUUDj?^J0E^AjAbMGYxt9I2!F?!jUO5e>!``Qj&L7k-7JOmI5OE1h7?=_ns(GND*8#TN6BEH zHq6%E^@d1VDOKLBrK0{`4MM_wG9Rcc3GqE_KCpBa9}t!;qoN|D-FNB3ArivoeUJ3~ zcslWkE}f@RAw0NGmWApn@j*}SB0UA#gZnO=Kg37K;BlUXk&h&A(lcF!@7u{P(M|9B zon(5Ep6Tx-S27N)Q2Y+&)#Bz79(P}bgqFn-G$DT#V)8Wbk!o)fQ0fRAxKDpz!Kjm;-e=>UXdvt`iK`e+T+SM&H2T4d`&A{l>>R3CX<5Us-bE``X9_^SyV*(Q5mT zJ-mDU3%tuh=kSfWa*#|FSED?ys%(#exZ_*6jtsiISej29aU-)Nq9X-alkH29`daC& zTWhxsaohgFj*TM7zHj&AS=qHGxj;wgNA{;hMRdC(R0RVap*zkug!wZSrwaS|(;2DP zx|34x_how!yb*tB6xlr6F*oG?hz;?4C5P@U22~nWM(6C*I@UJ3kHMK&!_c^^zhGTG zrU%I+uMf5%KOB@UoM_o8RMgwPzp|P*!x-ibr;S%>NsJ8#5Zp}5JP>Caa9+Hiy)IuAam3Z@yvS9 zxs31~VDLRHrHSeY&*;VHX+ zlp>kkH@fxAmA*G08i_Btk9C#rz4y*D&iArA&qUwX>^#%9M5tPj%@aQD4IYS5*?uOI zd8U`ikb_kDMJHt8F6UK{jDv3_Hskz(tc!Bx>ZexhWo|#$kf)AeeFs0&>iOn+&^)aV zM*H27d$~c6ldn$Ck>#Vn4 z|2h2VQ3q(VyKmt@ zjn+sB_{Wf$E<=qs7P;JIi)92 z0VyYAB<+$nxToJznLsXhv_52m$S@d>%=!Oj0J3LxRCxVA05RnEKEgw#wFy_#oeaClPQI@APbUg6$ zn-a^5o_%r#7kTa)_shl32#V&ncClC)fLonZ?+~@qK9f_Kr>0*sqC+QXd zTcL=1`hKp4PZHk+jh}B*y4OLDuKf(E@QR|+X{Wa^psizg#@_cfj+#pS(OdD0_UR8% zl`py5Go|vTeaux~M};}7gosXL!0ygwxPLrMbxTh)>w}S&5b65kH?m*f$_7X+zH`0U ze$uAXpoR(EDScw^YKzai#TLvUHyEN2kXW&SRZ|%sm|Uq%L|(*-sY) z&d_8r#p<_slUFHSDSfDaz7HJsxldO2rTDxWF@B&}vD|*ajmFfH95Z5Nm!iN(!3N<# z2VT!q;#V;WCU=sRK4|Ao{OUc$pW_9s6B!?#7jzFX=Hnh`(?caih1 z1^3}1aWmPRj|o{AaQ9iRBEMwN~Jt_-b9A+<)54>1~6uaUUXX|HX zCp)}jl{qmZZbfkHnqx;VQ!i#^qNQzy8En@f0>>>6 zwfQA$+At#^?;#FV^)SeW71FI^ZnlOO?SWy|b!$Q)hm-rpivmSplos}LZcYoU4T!OZ zG_W8gSkSc>p=(Na!prgKqYIVV_581O9gQu?(JOW>}sa7z4H>Hs0ktM-iNfn zwW9-P>MJ7WvY8OV?ue|?(oCEea}C_}AVN9`J_q_(gmZD3Z(_&kvg2^2Enx`UZQRCU zi|*=Z#g$pQv9PqFzN!Y~#zF_05>z4zc{o2c+t&u}SgoHD_69{Uxk6m+mz=ZEd9IPc zbi`mtQ4i!iLU(CC2GM^9A+@nuHV6DB67H>-6EYcR7Bd#pU(rE?8g|01Q3UETxU8?C zlIhT2l{IlLYlyJ4H^B`RvP-=iUBeq)y5q{a7)v+o#>QY;jeA8rI2ko_K`Khy#ukEH z3d|RQGxth=a&sN=P6jbN2^gL|VU0WF;@SfTmEOf&m_ZJYXX2i(uI4*q-#3-w3gGL@ zqfV>C^=J$ly&e$;W+p~R3OYPoLScu|j{@QQjd)|6lUFtKS3L^8b`?bOdc0ffde@WI z6dSm-7+yAGfph#i8>kHLy+fRo-&E4m2f^l{H|*@^QIqOl9|w<*A!b1##AD6IV?$GY zSnGLA8HTxRVH&{zN6*lm8fi+eMz-~7>t#61?!&FrYx##Z`G*#vDMe!8LzCfUgOe*5 ztil9xb_tE9^gykLqbyb}?9uc2F00EO4R6N6-{3^wzei>E&0rO#${LC?bwRiwhNot* zXNr(MoS|DW&Hxv(gz%e=@!3ElYgS>*!jCC)AUs1eJ$-!@@Ke*w<2_T6MSXo3+{y{^ zrV$M^-zv!OqcX<@A3Gq5` zVVRr)M>JI9)1W)3Mx%9s^(s*~3mvVeT3|7CtM3EC^f>Wr%bpkb$k@vm+}orXT@lGA zXIfsY$WbwxzyflVBRa{bFP#2JSn!Gp>NFDa>c;hx_Yb$3#L*q}ROKTb!+W!0g0e3& zm}z;pDPjsE4D(E9Ih5e3`1}Cbg%nvG(g!j^o;OsZ*`3N*juHjyQG9p+Q!1Tn@mZO; zj6U^LDM#5;q3p(U615BLABntPwVuo*%#AyKA?>viO^F_rzo#4@v0konGGcY&I@%`H zmHKjWce1)`NlWE-WsCtCMYK&4AYVstLQ0)m*y&tD^(uvPPU#2W)bg@jy8pzdw2Mh=k@2LRu@iiB0rwFkpBA0=hLYz zmwAs8#_Lgk+IMrNWoiz+e0%e8>Zy(?40!y2%vhUM7fzlY$>|XicR7&y=IJ^%Lj^fb zF1Gj)sHfWE!PD{$IZv7e(uWn5k8z$fjQe!r)?AlSS2p>oS%p+EAki&Y`vVeIv8990 z)}LrjLhn~7Aj!dPV%dqOpwJ*yEBOmtj{6TM8r{BR65jqWDfuAug@EeV^Eazz1@{s2 zNv;)XcJeuch7zMeS>KbGRz~=j(XWwb4sd;W5;hN2s0$5JQIo$?%SDtpa_l-Z!h}50 z@ODjjjDX6Oy3660vih?}@48h)qR2DHHYJvRNTS@Wf*wttIlU>d@g(c5Kagnkm-h)P|smRM$3BKkcZ@DdBlM^FIV3?o?6+1O0rNx64DMlgYCt=j8 zL(h0Zc)7IxTFW|GHB5S7ir#B=TJ<6#=(G>#LNHW#bY3Bi7 zPN4CAv;05c$fcEiP;R!^1gL+II+?s&io<@PME%K#{~BGGvfczi zsoTpzDQD=17*$sJjZ+@RhaD27;{RH1)og+wpp%#ZPynpTQ`ZlNzQ6dP(gQod%RYgf z#435v$3bcW@=9EsY>7k1TA@!&xDyR-H-&o$sIk{E1R-NE6Fv?bijr=pnLu=uXm<`$ z7wM&zg*BAJ#$!phD@-7^Kv-~%a$SYAyRnM;@uXa6fQgjx{<^E-oYq&{o%5FW$xd$? zsj}sjlmYYky@K}V#Cb<4dUuOB!GtUR!1{il2$h>|Z4zNe8@0TAEvH!GD=5_WduF-p z$o`CyLXFY{v+q`cb@!7io9=ItZK&z{v0f_Uas7amneK`Q+y>HF7D(qGO0;ZBQKPMR z@xyqZ{-!klSGH}j2FT^?E!MpJ!Cd~23W&7(cU4R)D|gNrPp_)4pvBNzxx|pui3=uCm6Zl;wxYjS;iur>(ka zu(A(O_yfgb2Y86>hw}Q{NXy|nS4+G9b-BU1>L@V3g&cBT;q;J-2a;M2Yum)c_SX1& zS<3DDffVgd1C-6D)F3&%CLrOyp8_1y1WG@bnNe>kMfJCsM*4##U7y*qd@IH@(z;aN zCo|0LT_z+u6Ox+=$>#+e!rAHKoOE$+x;Rf=+!I}#uP!b?7Z;?vvhU9}!~ZD!!9Ap0 zdf@R^>DUR!PsGio97g%-N9q8uFu?4abzt!Ob)YRK`tQ|&C}6`i+*%#3A4yVR_fYC0 z281?u52))0t{>Vxc7T_Cl!5&G?(PBbk_+vwzQs!zpy@|~1)!6+f(0a5$2CCoJ3G{G zb^I{*0S$ojn+2v$dI0)E1!f1h*=%$Q{+wAq;Ch;rgx0z{#;-Gbv+QrSG>e(AAI<$^ zX!lKvo>IerUf@S!44{p-32DG=`Q5VApAB@G+dyDj+SK$T%~`edz~<&1SGlR_I}I?< zW^bZDY__kLeVW#mepbZ>>c?8l|9H`NT378y`1{+~0zvqJ{v`_NUjQJ`_yX>HZzt-` zUP`3wgg@Qu{WdHAk1V0}%D(ZPG=X3W9UH&5w`?t^ie(VGTa3#w1gV;7#ey|J`SU4k8cz_a|bmyC8j{|z%W;Jnnm^m7h_`b!4tQ@)mGzf8T4%;&;tNnxx&852Tq~H z*c@m;v-*?gZ-j@DRtlw%Ly;E>fGYdaO0G-oH~Wj4`f+)M2~H_Pl3FLhQ0JMB!d#b2 z-v>s1?I<+<5n$19N8zy_0rvEF6xKWcfNFGCN1^MF0O{L03VDA7DDk}^EP>B@_ zzUirDQn<(y6K$X|0*~*N9II#iYTECi_WcEde^b|p4(uu2x_R-Nk@4i{-!b0GC(rVHWvtS;O4r0z? z^S$3vo<_uzhKe4y0w?&HObY~L!^PgR#5upo0a?+(%wCLH!&a<&EH)1&~uv%U!29{Oi* zuy6m+x_ZxpUOpGZrRlp5`CDKA%6J5SrI?+qwA2i*p(00)ky!WfwC<9tWL= zt?RuM6STIuB7OSGLCSXoEDIxV0u!g7t{+Ie+vm~p&_l{0Bj;!oyHbTtm9)FX;`8hW zC3j-DlAncXym)u*H5tRl zY{0N~|Fa*SITx=cU09mr3kkR5M~BmHCK%d1rQ|9wrFI@#VTs8r5prRm1i>zrW9f5j z$lF5~A8>nR)Jm2=;gB!5#APK+SuOlE+Q)_EQ;B4kD9tp<5<1@T1Dw=h?YUe&rmNk=VrPCS1 zj5e5RGWi_A$XmWnPtV@JqSRmbB_WMd)pPV0t5 zpD~4mTQHRxL5J4mh+v!WCBaPW@78yLL+R{iU~Mx_;Ng*Pe}H0b<7j4Rz4@U@-|fE5 zkEn5)agx&yNeR`mysn<5Y0*iNo2k$t7p1E85v)6>5>LVxvgrURB7Oz zjF4BNLIfafj+V#c)lnJh(!%B~Z}CBK6Z)?^1d&f)`cs=pvI+=uoV{|nAzJ-ncv)!n zJF)2QN!4Lgc1+iqLh3tZQI%YV4_ce80u#NTRtg?T>i@2p12Wd|cDH(w9&{d==dXw}-tTs%)& zysi@!M#*akTR!|);L&t1VKVD=tw)o|jj&6KAzww4<+P7TX$j@1hSA%zy%HH|i6k9$ zUR!SHFOX-S4mm=udS>k9b7!t0WgAgA$%V*|pNXSk2W+cfgg`oIHInW>;Fpn2deMr-0JBz@A4nG2V9L@YR!%Zj%Oek@?<*(8{?fqw0Njj`n2>fwUk zX`3vQ%1DyxRtdeyFn**oOOl)E%Nx+Dri+o6h3E)wGMGgSPZKQIOg;({&!CI2=;M2T zSR^W9qf%wXUBm@~&CEjVn+Cg&Rbq8&78Qnw->#);;E}Wlof(}cxNAlEGRlH4MB#CF z5U(!fN_y*~$*lMECf85S7HU^KS1$b2aAhS>oL9h@Iv3?M$z)`$*iT$Zh?|fJ(^ve&!0B7rAXwiqV49@qDHu4 z9vJdpO2kRkz0v43gU%}J^u|ofk~>}l2}5diL$B6ivTYR|qL@4%z89FtAuj@N(KQV7Ty=bH%)*Kh<*VR2}F%%Ttv__rip{q?pXB?!GyWY$6{x zHIvBiTgmu(vTHi%1mQ8uGc^%Vu9cygH80QR-_bz|7Wn)J&u3B_>MEvG;`6`he0%BXf%iV+ z@!^#(3@Vpbo^^g;xqnCj7ZbysNqO#-9fNgGqFFUx<(;k@Lg{!}rd)b4594m07E(CR z9eeKiXJjYU*|8ZBUxxZ>C+OD<$0W}jOE$o7YT}fcF;`en7-14r9^|fISz=gxTVjpX{Au` zXud?B<_;NBcJNyi(YKDWq7&Qk0VnZE@$e|OYn(mU%JeQ6Y`R(HazjnQk$kjIaK7Hp z*V8FR0)no-Hp!8{Mwh4#zd^xpB(0RC=UyKiMpI6xZyWF9itpmc#9TXmI2)!`GLztP zX)cPCF{L11ggub2y_QMzOTGZH{7r?-W+Ym}`Kyi~9Ztg8tLYFGUM7w^5f+u>T7bZHP+{%4L``V zGQ`gf)B8qF7+vteXo^ZJy|LswXeM>5tTUPs&CbNeWepU~Hbm8s$hUNn;;;3c1!$5B z+|_b?n=ad_@mtbtrAb5EwE)laR73-5HtTptUd1U@hS2BJ@O~~E;=u2lNUi-kxikK-Q^idEr%XpO3?#bmL^65TW|>^N`GnJJ=cOq&*Z)XP-F(@k z;h0I~9W}GK%OO?kL-jfzC*^j_>;2Vl@6o<3vu4V1IrSx-78w`Y>Pl~j?pvqG66#l({ax-q zyXSlVEmc??;jk`=ac3Xv6E(33g*pn`H#Xck%2XuJeWE~s_vHS0Yu}&K{r`(E|JoUk z)b9+7{)W%0WeV_SWD;Q(0fsLJL*t}(F>c3wJ*EPAXMtD%h+&`s7!wQ(Df!6(MfnBA z`iTVvpa?+Mz3(Wy!w;ak(?HdNFs(qkf$=c1?&SQURBYNGJYDUf0cY6|#NabzEXOhMm# zfiPu41XfeP`{>Y3LEl7#Fy(YKR#U)RiqK8L?H1oStfqjsg`k^)z9t^wmX1WMrhr$- zqnm=h2o_<=-BhfmfS1Ign}WW+5Mj!R46LSrR~w?6g1$x!VTxBaR#ObYE5^``LGKA8 zj9J};-59K?7rny;H-@ph1#AqW`-ScT^adM3Yg-3O^9@}idIJcdF~1Y;bXao;T{C)h qh|oN*8>|^e<%n($dWC^7=SUCM>Lb9Ll?|jw7zqCWcWGGkfp`EVu+HBA literal 46139 zcmeFa2UrwY_C72qC|Lz0Dj7tg&cT2nk|k%M2?7F=V}oP~ zP0rBW|7x7sb#~QZ-TC&}-?z^5%y8${drzJ7o^Wqf->SNyaPkzg?Mb#3Z7aZj%jF5%xEGuI*|lpwnMQJE9}6N3#! zcxWsB%Tb>$c|LPbv`f})puLrL8>f<8Nkc*}x`nsrF+7Ni)^zFcS8NNTFAF>v(0|x0 ztrOO0b5KOao6}j9K@0P|oALIXcMwENUR_k^6Bq91o+7+y>huVqOsoz$7xj*;1*vuM zGYOszWvPsl3w0+GPEKtglm;G^@M7;=OQzfxpm}t`a35pG^iy``>0DLO8;593;0N3R z=>B{$SIQQdhh(qAMX&@AurRX-oBPD4afX+R*$V2HJ@er%#?ZapPa5lIs_h ze-)M7zo`7HsO0%Y^nde9YpUaGD!ezEh7x1vKs6p48>Re6IC%UZ!+oQvvW>*n$7^kIFR{?p7BM$enTR}u zB6nwZHUevF4rZ2?)|(Fa-8Pm-5^Htt*1-4sb6h#Vf}Ze6C?L3VtP05jO4h6LvEM_) zl7PlNcD5&%8bdWWJm8BtUI#9wNnqVy}#5x)68Ek zaIn=la}bJLj@7Uohp*Y>A#19kiy9mp9plMJr@g+mje*p$c{r1s`{L3HM`y0sN0B}8 zJUyY*&s%|kwTSKQ)Un!`k)6~@&|Ucc!a68VY+$Z*Z={qTTDr7Tgo7*g?qTq9S$H5 z7Mk1i#!3%wbEF<@Y>b1bMr$fOc4FshYby7rGWZXsocV35e70t%z_q&vePASFE$>5K z>*3~gA34Nk#Ni&p7R0P-~HdTnoNO<^G;wzrc*Zk{SR-P=QHL!}2)orc_t zIP9xJpa-esk)`VOENH-nV?Th{H(Pf&6x25paW*F{2SaEv=|3k5Bned5VZ& zUL4=#I7sa8;rDj(3T6?=$bl5 z5PEIuXt0pqSMCsd&?cPICj8GPasFn~N0EOvsq@o1`{H`;Kbu4eA!uG*C3V6D3;M}4 zukMnjP(pf)++1)wHSF9$hnzM#;{^W6cR%Hf4z)saylSI8waaUxPp7t$j5m#3B`AHN zG84Em4U3b}cCU(N(sr+iHq&-Tm$5-bVjH#grGRCnC4#K2>k7$s{=)1mYAjB&YTC+yP-j6|6z+F-{T6Dn@WwJB}?l zk9N15kZAB$UeiwZpvrGJWPlq-``-`A?{p7;IfQ)Pwa@S8xwq|-AUaQs?D&orr*h_K@ zI=+)R;zJtxDpw}3iEJ@iQ6jN4p1HBy*ak$(m!UUm`@RK4-l`|QJcb>G5_vvXd`3-Wp@5ceCa}LLvR&J^1Q6L@X@y(Y;Y3{aHR%6MsgD6xH!v&JiC6L@rEVJK19FqeK#DJo9Av zh%(%x3>~WNTZ7UPLi~9QdlRMQrCjlEiR7#zCPj&))A60ic{d1%B%i>>wZ&)xM0TIi z(060uL}@8qh7Qs8EdWIF)Dx4AVF#f^p38l=s;-Z6Pp}<^Un#l~pr!dEmT!tww8fZ3 ziNw{=_h8{fX=z%9-lOeXjuHtWjvK@NfD%cT`)*TRALX8QJB;*FbVrm(KNikv;y9F+ z+&aGVIqyaREfXfNzbUd!+qW1HSyWGaXACjmMD>rS@Nrh z@1R7o==jd$yc-5YYE59X+G6yhw8YlXcW22*8LnQ2{w6|nYqefIu@*Luu_lUfrW zWKEUagc;cl<7ceTep)Dl-8!I7DuHBH$(0G&2=|F$XTa>PM4p3U$`tV@1O1BfP86(N zV@?2`OR&p&RnqfF4!u-q9T80F2uuvc)c6goI5PlDHspG4h@Oc82DCH;U=^W&ER2M#jAsJ@U}2nMQJgp% z0E8NHO&Fq=qkw@e4S`sjC?E$TAt&S6AOHx7Q*4M6zX|{^8gbDWp|_xbK`jkISUhw9 zkeiW^hw*GM0PKlV?2QxWMd@h7WnhFpgaQV)Gz4Qgp@961gaV9bLjd4ZoZ@txxFAYL zBd+I0=<_IGNJ~QqRuKv)%t$E0cs3LOZp0~W#)*rfbTs0cFhbu$0Yh6FLa{bcKoBFL zIOEwc0C+N9@l?F{HI$CVTr|e$SXclstfe6gi-#TnN-z>iGM)_wfQ0diMDgNMC>@Qt z42;pqP{8n(hHxw=6i|kdP?qs*1OTLqSEP>@mqY1j%=O$Dof!p;XlaPRDnbDj7zq^_ z&qe}3?s!F>cyT3^j>cRQ#^{15U}Q@}B-SPhsKQ8igYj$>0F;PVl#CZwMd@h5MPq_4 zhXO{mG(=(XFaW@tjD)us&qf2loAHXb;>B;HbTr{IFhSQu0i#&Paq4fIu;aMhWh+XI97Vp}9*v3gO1br}ioGKR$h zgFWLFz2e3703h6i3(XYW7X^%Kk&MG4Vg!H&jD&`aVQ~O3FkUezUfc)(ikfm=Ger+a z0pnXFtr-bz7{d|)U}d~wRlK+@N=H*JG&A%%6fmhp zG6{=_2>{wN5;`!3B>}*;c%`=ZI*02wL)gst*vv57u#^N#l7@cVA?hj{y!ic302SL!BdRgn3_YxKakkU+j-uVicwjGFr)dY-iC<%j_>5Aq`gPlE^jB1E7cWc72=HM z0KFX@`<=SUO=EkHBDbmw%aF$T`N>VUHOmx^WikDv6i($((%M<(6dvVJ^4eL=sbxshy!q5-8E}AqSA(|X)3m9z5~kq(_g^9f}W!luMW=t9=GZ}4(ceYjoLJo+B8uMhhR1{ zUN$ofAoO@klKP60h=I`S8+7>hfJ(Dk^mhWQ=Zh{Xo8)2u)#pvYOm_NAHmdqG#S|vQ z6g~i|54m=Sg;obPP^K)b1hdQqS$+Y3(es}PiodT5{>MVkRI39MNH{Gk!9;UG3LxPu z3TCX*XRJ_}*%(t8A5-WHWab6!j(-pM)wqRtw9Hu9tRgxW*W5({aZlMukII=UqBU?=9t37m_pPQ^d;>M8?6p}pbB|d309g5 z3IkPWpUO|6 z?G6L24%FQeX)D2h4@mFGc0RlD%dZRr=su){c4>Ot{fqVmYZtywKpaRKhIn`?Do0ma4CN-)D*kO`He<}2ESm$zFYm<_7$Gq+LCM#|4rs^Jg zxy&sM=#4|&+z}&*J9Dr3YhUxTclhjY&(@~rt)vm;fi_oq26i$*vU-90dffw^dHXpa zl>;chEyvyJauG^E%owd7LOBrnhV4+;ZiXsUW$0${8l?WJTWqUpw>Eg0CGD2ai zS0!4xTdIL{-bMu4Z#3;SW)CpeZd=4U@tr3 znuAo1yR~pV*IdLd2e6sC5G*j=(?vz#zK;mv&qZv`_t}H1`f_?)5o?o;D|uswUQ?5c z8XzCfjY)qh(A{_Aog4^{EtA7#BpT|`*oL&A_NyExPC!b~e(-24-IqrgzS?%L=U`&> z_2zqa)c*%}Dqz=hp2hq`=;qJq8QK5%2@g#Akc`OxVF z(c>lHE~I!YAjtB38iVs~>wF9T`3|Uk^~no_TET^d^ajt*C{@2v<}qsJw@>miCM!^_ zkBb<`IoBtPN#~kkrR+h1n=aya&FvNmrL0weAxnJ*+E#3f78^}b#i@R_6hwjUxJm%N zvM%1x=G5blXe;G+7_Z_+(c2-0t056)7bJC`<^^mKqBzKTy{J?OSa=0^MEFl15B`dS zqel6J-G4tS!XNPPF10`h7-&^h2M`mMrkCVND`xw)oiP;R8qL_ z(HPY%N%Cps`L6mVS_A{38yng<`~Wl#-(IL)z9d#t4MxHj6K!j};S0;~&>c{i>xce&=;oY@%Coam6*`SH8~?xwL0_3F7X_MdmB#f` z)_)++r{Z95o2pin5qSJy46(Vp$mCdqHOe72jzpLgAlrNT0MlVpJs@P&V6Dy0eqUwj z3UU`YqT=7P=Opff+=#Ui(_7zYLX^RJV!B6~N2qE%%aHIsNKbIea{V6Trk+t`>N0HO z1rl!F*W9V$IoxQrM<9zXM?*&MPXul))k{0PnZ%)HpoC{m&6f6XGYpOEI^%m>O==e7 zhf856r_Tsp>*$*(Qc3E%a4z)GwyC!NZ`p3p$B4s57H-_)*Zmpyv@6(47MO<5- zb0{c~9lrA@_}Y0{>=tUlv?saNyfG!@Mp<7Z;4@R3cv)|x%I+g}*h*C7Wto}jX3`$q zig<%laaiAAyR@skMwRv`DuPQ*`HCzW;QrBBYYuMecnK47o)f(g$A zw{cKKg{i4Yx5{B;(%()?=nWgjsqm64vhCu(M0M`Qb$K~5W+vjah+E+nY7!Q5Pz+r$ z+2NkY6sZK~FVnnA3yO-wR+D=p%bQ74oEFln_i<(_T{N*AYq!kzyqbKgtY9W}bXvfz zNOCoW1zAz^vS)4nYMuyZ2$TLF6ixH%jSx~(A(x$M@0$zbAC7!aAr?r_kQNmc?x1$< ziX3q!BX3$vZ@5RdIJsO;d!Jc7|M1NT3NaGK?_dbcYs&W@jdi7b6!yCQ6;)T%lq_V~ zGU>L`e0d@myA`#PvN&WYVKYx4v2n8X{g`m+VY$=Y;ZC&4&@7sczUxVCXJ!7L82{vLSr&6T zSlY>XS)uTb<-W~1i}_BIeL9A3$JyT;)$Qvqg{8hUIu zRA>Ds+)&KEPG& zuKUMk^iXEp@9FaBAJ+|YpDal5irC5Xj~m=al{+U8^VqRUsRMiQ zU&v~2#l82(W_1E!)k2rY_}Hwj0<0FZlfOMKtB!8Pfqq&fc*iaof>py(5;=i>$9DrT zfj;k-k~ns7*F0-0eD2&8cI;S%0I~Yt$EvgzN$|0!H^GWwDY~3M>EpAi3wsgXFGY83 zRzrZ2A-nkE+_6iB3{W!CoV%VJ+YK)PH*{CrYdCgrlQ_Ba0Ttlbs|-h{;=unxR_`pn zxO!|>n*mm{ox9?W&FU1uYJbJO565P84Ja9wx;&=GE*ZNPB*8$*IQA;T(W^MH2?*|E zmyEp!?45r7Qbxy*)jXH3|2|edHgYaVr}le03*jo@D&V!@w&AhhvfyRmX5m@lTH-C@ zF5*$(Qs5=xCgN%1YU7RIj^GKCy6`L+5AJo>2=Qrgdu#a{K`bopSlltaV{^ygj`&iNU`zhSD{p73 zqB93Y_|1M~{rSE5AM!`@>+`$wTl0tWoAUefJMzc!>+-wuTk?nU8}s|}+w-CM4-I4> z%6`cluq>mOTSFPn0(`7Rey{vY{AT=!{i6MD`}O(p`W5>*`)&E*`33ta`L+A8`Q`Xo z`@#G!;tG?ya4mV6>s_Q^gNo%q3Jj`vS$SExS^2X$vw5<)viY*FX7grqXY*TeTJczM zS@Bt2wc@qnw&GvXAMDDHG?-I3SQxx_d&m4@2wr8bL8d{xLAk+WgI5O82E_)!1~~@F z2Gs@+4PF{V8Wb7?8e|zH8dMsbhJ^Sf;|l8yYcxya)Tk`lS;cwb{e)|Ad{C|>xB3@d z>*8)@(;)9_|k#`C9jq^*qR@jc73RK}z1?Vl1*#1F^+YGslA5mqP&E z@T%p#)>?!*!s^w*0CL->W<3uK5z}Mtd>?H#-OBASLsq^Xq~k7~D?+D{63*25KGJNu zg?oLP%zZsb#9cZ?gihu6?l$iA8M3|gATM|68WB3ZlyHvD_pxTvN8K3fLFn$%z9KKh zQ&g!tFNc{Cz329rcsP=p7iG3t?=J+>-m8~6M4k8Z<4~kmjao_<4i(<6*DH~iBw~C^ zd-Cqm9A0T3H=@^r?(UX_&#}yqRj;a^>eOfj#8|g}5u;5CIucOI?|me|*87Wqu?w&x z0aYg5M*`f0qrV7H$LZ8~!|iZ)S9|ccMrF)z4T{^pH9QG_YxH!X&Q=i)=u744d2|AB z+02pPL#CxK0;%4lpbWYm;#0Ez7<~^=pRSx+z`@LJy5rOP6^N3SIno%;Za_(J#D2{yOkOO3qAooPD266DBM{?dm7K zjO#l}hF;En%E9c{fge+HV;T4PDB1OL?neMuGaYAJlWC%aB`&-A{W|bNN}d1#*OE*V zAuK`f>Lg-gWt83^yYzjs&?(ODwfUw_P-Q7mOn^kHzC$?iNm$xTkOIJWO`<(~PDy4@` zz($&_Ie7ZZE^9U~}!>$bASLG6jK+G;9sr~RpHCy()!yX81#daLI zzA#~H3*JX;No@5{g^sW65P-oRxZvGd)Yc>=2IQ^V^e0Z-0sr*&tfLJ|-`bv4f^*J% zIQR@yUi4j5K;#6?^Z2{GJy2a6owGroH^I9tl0-@W) z(xo0Io{!a28Tqe;w|LW3`Z6+|%$$sPdduc+X(CI;Xc==#R>j*DsLtJu64}!9cyFQH zLsEZB(`7^<|7|eaSkLO)@(Baw!ugg>PjYzb(|xio&>Mr0&dK$H_4>qK6LO@q7@=l| z@XSlij;;Rw+jLY)-8RvlyDOz{`zxgu>fxtU`r3vsiI*gqAu5(XxiYxIy6hd^?`kd% zKw2c+9lF+r50)QSdR`dmx7d>K8GX>57M10tt^IrvblZF&#NJk9^-yl*a8fm|(6zpD zbhfR$tFnQTr_J47V0Zb^a-QTvNJYpf`nhI$N6b8fJcExX&QbQ@R+*k4r#3v>gIjHS z256A>1XY=y0vbd;LDih84CsVF{{XgY8bZjcV9o zKPBV~0U8e2@q{j@h7)!sVF{|?f?v$tA0lh_w}#h>*l<Gf`;shq{k2Wd(aHHl?DtxOF zs9%#zaG|?isz$&`-*dV|PImHTDL+Ny1DsfH&;}k@RWu_U+TveFjwRT6KU`IvAU)E~ zTzOG6|7_Pi${ZRjZj9Hhx9!u<8C77xKcOo;!8*TB3cdE^nMVSr-QxSSv|^w7IO%aL z|Al~yPnBX|kC1%NbmFr6&wJ%K>bEa_rWz@MJ^OsG&zhDHTy=JMijG);em~^JsgChV zm9iP25?nA;E$}+31TTy~M*Vap`0~k&=#h{&oJ4-9`n4w~?y{sXyfemm;m9Aq1aZ!| zWu6vQB(({1y=LvCF8iT|h04-%{(pD`d}~QuAiKs*$H>%ja9W-EqaabdCF?y>!|Q!#sb< z)x;Kh>L)2!Fm$$5bo+IA+K`!oTuf=rcp|-Q{r&Elw+S~t_rKk=4e<@`+TJy-0Z%qo z?=C6nwc+tOK2}$m?;Kk~Dvq6TALH*Mf>gVSUSRBR(CZmwauB9f{m2Y5CVHYplAT7? zt;)-v2-fVm&hy+7e;Lo~twgT&{TC;OYHZg0yQD?4lebC1SoT&NV}%*%&BH~WpG+Au zqc3>LpZ8e$Y{JO?6c(lf0S9x#eenx~qY z@^U{zqsOF=-+DTyE^^vX+)(_yINwH%Dvq!Pnl#$^lMyH0qS5Pa;^_@fC~#7Sk%pZL z!=q-V9x>5NiViMBw?FX|P42|di3s$!7{;aOXkC!DM7hzp+B70wFBTKVdT$!~G+v`JmweluNQV_#KfYgdb)9f z6sGduw|Qv<*Na#)^3Fc>e`x%NL4Q#fpfsP{dtWQC3ex-BEC?m z5-YQ+&!ht>Cf@KFW+kxz4y%0~8-Wjxi$A}r%^yN0HZ)x@GVbMhSLi0Des)1_-qC3D zqB)2*@izCmXMyHS^OdKCpFIttnv28^ck)ytQU^@ob=DIyhZT9b`qt_(%a$pJ6f$$` z%jNcPQe?+Q>0+Sm$lV_A94i|)@AlO@>h5n|j4;HRTfNk^Q@oPGR6A*1oY#IykPrr1 zFJw;%PM9}HiE1Br8I(yiPqTON#AR=v+S?8hPAIdUnz~SFUEqw9c$jH6yscPy$&BZq z$fUpIgHu~CSwTyNsPf3vsEi{hl+V_4w!&;yr*pILQYp-SO{gFCF7e>wU>20cTF(&QQNa*H=MBy$<91P`aH&!ppVvAuopRU#?yl8d(q3S21 zQy8HS(mvh15$UUc^P-jf*Ihln`rH?-q`&Se^wkf+pZ%y$B!F?P?j<$zE>RTsS*zTK z?dF`=Ne`ZVQVk>W!nP1Q`=PH+hdc;9gZi=?5j{rg1BcTO*qrh)@T7fjo+T2&j(DJQ zjScfbfjw}guK3c|?kw7*uk2YWNvHy70xEA97KVNL9QCCqr`E3ZR&NrzwB4`cJ263x z-f%LUTI`9Tgo3ZLwr`!**IgsN&WpZvK3{iPK5+JYQ1{~Nu7d~8I1lTZzU~Tq=$!bl zZsF@Libu{YkLvKzzB*&pBj?ITbzEO}2{$gg1mFo{JbHNMEP>>MUFuf9WJ9R;)X;@G z!H&DZhGYYJ4uQ7X{k`*4_IzKS6mxhOJ_Z3kpm^r|Pp^xP%E7m)VZYRW>{f(}ruvzm zzZEf@y5aTH)$q$Fy<(s?FFDUXPSqbr+pD>XPXngCrLEDoo%^Pwa)4bR_;wK8$fHfz zfIu&nun92?59A+ji$8n+y-n!b0DK*a=L>kWk)%>p?6;T~O7T4Q_UsBJ!wa(6HD6iO zT%#I&Hgc%le+6#ik~H!j&7$y{u;q?rkW#@^MEuj?jbYuPD{S%38T)vN)w9`~N^k3X z1X7GzB!iV771#JZFVR+A^}Nj!F|d8AtIDu{h+nvXKb6P*lmIh!PSEhM+;d|U2Yl=g zp%^=UF()P7mNkVN?^h3v3rvrdyAem{Q;li7qvU(twN+s@dWu4~bGTiP%6LoDi@tEt z>tUnLMBh7jcp6fM_A28vM{iJFaoN}f7IeA0{#Yt$p|UCN(~%wAeS`A7hjKFM;Wlfa zWw%$ucn)>2$7X9Kg~|6faQ(%=s|u_AI+QZHL3+_1E6SJJuJxF}9)dsFn0Aq7O%wQ# zU%PQ@oz3!!!tE{XD3`gSdE=B*aLwZ90?BS3&#z3CiY^rrNU4=+*ez@?ecrJKoAjMq z19?Q$Fky-IG#ZERdkO|`xY~EJwsGeS%eq$M4_@?gCE7tmkOd?uOifFqnC^bu+GcRo zi=to4_w%_eaUIXa2(n}zCo0P8o|U00{yetzMC2(&bWZW-Eu^LwIJxhmVJ8P0bl;al z!@SIM5%S@rDL!y)=2wSWlW_0bOlZKU~pQM|F`y!geW#N$)?CJl%cP zIUJ0xD*oc6`4ge$_;(@_Pb0b{BOTMQKJY|cCR5@SY)*KmgX%PnbS%R9z=!JO7i>;^ z=OyvsF=N(ee=w4=DQe>LwPU`p6r+LvO zpe1qM6#GVW)}$g`lM~^79Y)MqLi)|C6%Wr|+u||!)E_Pd?M=z zHjvEaoSHcdyIOPa5{7{Liq<37Q1=xjzq+q@@ymThyJkY17px604ENAcBmf^}g8y_f zfBC3f^mxciYSKm^nd|KIGPSfVds{*xTRNb|Z^5}_Y0cVALZ+W$0vb(%x zlyoQm*h<+e@67z@{J+ zML-KqvnXAT(6v+MVclycO2Y2x1rWRtRo1d;iW#F0#t7y6G=lnu=nNXZq(1cLBn$^! zxr3LS0#t@?IP;XVnx;OyTQuJgH$?i9sSn;GYT41~ufKL@Ik80tvUPi3;w+MiISxW? zcLT>sm?X6M`e4MM5$S3(&)$>gwQ&&JCwD_1eLIB(W}F{Lavppa_1m-fV6DZbaXw^|ZGMQX+U z*TH0(@P`B?0^4{OH>495PYEB~Z+<=Y;Cao$z$@Z;B5J{ha=C|?Lt%Q~R#_~cWhOx#&Ytm>w#=A(Lk==^KX;il)- zl;oYQUh(r&*fWpq8lSMzIgsVvkai`-X}eQ~t8tmcm9)U5fgM8Jx{GgQPNon$fjzU1 zZ)8Hof@&?`(-j0Yn3P?QFT-^!cz4AY)7!(1_LYak)6y%P&fD)UBaY%NYGh~r3XpV8 zm_MC#M^X3fq$B!?oB23LNvA1s{sQ(5jl^|DLLlk;$pX8F8;UjGZ%GB%CfQd#6`D!3ZFvlN~P%m@`5k6A1sEBLlG_;6NuD92hSbEeh++n48GCWJics`JMN|F zfn09fpJ{^z?hM$1M3DPu#>e5w&1LYJ*wQ@P7b63;wY7Vz^NSBhp;eU#C-7ltZ4DeY z1$kMG4D|OzZf`+h%_DkZ-aGr-@KTuyP~9R4Ke)PP&*2$-$akz3AcE!sFw>Z&b)N1QfEWm(@Nl%&#)r z7noJpgDf6KATOD3tv6b1O+V)w5F-I2?+oNoQO(20J>!}w_O=zKGpK4rj2(A+2uf?g zz~}2AN{%~AQ!6VgsfS>s7w{N)Cspb=vL22J9RVZu=eetP=lw}SJ@tH6$ARzV9lEVg znWciskex7Ctr${if7WW=X?u01eg)hyPzyeEIuI-0+rOEe0`hd+Xxmy@5v#87=CAhN zm>=-e2AASos;-1@*t;D-V`YFR(U+%Uk??JO0?*FWYOp7gp{yFYxTU9u*q<8!J8VI( z?3TE>GaQzj7t6KMNF`VU!?N{!w6#u$YFEaN!*=vi?&p~uuIF9S1y>;E{VPX6sZj4B z#9`jsG4IkE(0C8i6e8n5WK(bF=1wimC>U*YvDp(ol>sB*KsK*6LTozd<&uyyh~axL?uT#R3$A@df8ZgC_0uwK?qOihQ> zeSIm01lW34+4mUOk!Ry=Fx^r(a@1XN)D3;Cr>c?yMbrbE^<2w7^ zJ#TGY&PG4%p#2_AkQOU~mMnq}L(O3pJ9Ee)EG8!;{t{W3eVa8+kIQ_COK54ar-!z% z4bD@S(f*KetoZMaXa+{u%$L+Rk6@>cslT_TQGIXyc=Q6d4E;4VdI^U#arq7%8eVO@ z?zeqcd+1Jb`Hf!SdaO`d46e$7_lyiXYReYQm2CfDz;Zd8@h}Ipzi-M~SF}2Z1<)zC z4EPfmPZvSW96=+fM$KTB8jF_KD}9=_SsiC zS5Zh-hB;h~Inp6BOkPD#t~fBbw_;GRPvIet&&wX+sB~_j+Qvw`qR9&yW3)t~ALrSw zpPs+WpmxPXjh$F+ex4=Mr6633L^#94Jgmvq$tCcklk{AvKyxVZmXQe}>_s%1u?d>} z!uLi*OGzfE2Pzr!%7hc&Ha58s{vz5y^7|^~vNCiRX<0F33C3i&V`Xyix3u;h1S+8S}B&3flgr|YP7G^7-!UI zqEl)jEJ!FaoeLb_#A2H3XzHYy7Z2*FJ348q>Z+TPFh@*A&zP8agr8@Q=3?z0P;evX4!C}JXS82$e<|M%qY3?Te7cix$kRP z=vck#XdKJDMy_#z++2%|FQJdpT>Vx;bi8yxnwwOjz*ip2T;DA2EgO_i4y?ZUkSCXw zA4L!UVmUP1oZoUJm`HBzWHJIYTSDcb zt*@?MdDYr8pf+X5;>^CSjMrUO0+sFYD1c_qxfG^Dv#t5n=o^7EJ4JnmSU{B0Vxvo8 z%1C>HOCc~AAQrIfZ+`xc5gs=Dmc#A9U@5gj@eHK1QRVkM_w5q^Dt*Uu10@vZS{x;b zTcbApMTx~Ql;^A~HAigHj&Cr5lSkm6I<9Wi27OpbzcQSb8Q+TE_q%=rjwC#2bRqlz zV|)fZ1&EE9eCI+q_&t3{fg{X!toZ*Z{=0Sn=Sr@{OHt4@a`wwFGcP;-d4gewFNsVuazb|6{ zRrOGtnYYZeFNhW*@yx(9$F<_`CSol~)VJYBgE8i#V)pd((}&-SaU;OT>G}mwt67)A z=TIvvR7y~+H+)t52&^xRw42)SQw+BQgQe7}XjCZEUVbMd_<@WV2Fz`5 zIsIW5i^*~GW&<~04xCZ|6Y{;Nb4wIsd*tF*RgWw zBOo0W>3Hsw)ZZ;NfTcaA|G?`A%AoH=9iVLjUm(nsyEG9GW(u%O%k1;4N4p3r`}-vH z=Z>KIymd*z@Efb_TkBmwpK71Zje9s?)+UaU@iL%wYb;a3ZXf5>FO zw1g6XQ_uD1bv_oWk3;clR?F*364`#_xaw88dJ-{tw|`G9Ykt6x@5F2=DrVVB5@Dk@ zB}IvVl0X$#_`yVHgPKist@uaEp;qQ7xxoAz^MevF3J`8yIcVboc)H8R<;c?)e`InY z;0Bul=n_x^UP70SYJl*MYJeF~1E4?+Xos1~R1OWHYJh7=V)AzbNf-zu7?5l>HvGcF z?Hxx+4V?4uE5tvp9+rXnP^4|eUs(=ay2lR;2J8g}eK*yMUy+r20>1uz z)%{0K{eYOmCK*u@C`}A=Ap`yb*%U&Ms2PSD*?O$E2GYt7TL*(X-+|0J6)R-}5MT#*R+)omd8hro(}3?}VAER1N`a zP1%4S$1(shuN}162Fg#5&Gu3GIW1Hz_1$_1=xRO;tZLz-HrqvsfJJ|{Li@^Scqvc5 zPcDQ(tKStDUc}^W>v%DFpY*C-*8q z$dbpKC=IqfZBKeINU5=2{CWiD`9%oNqRwG#9GfmnSLJo9jCfU9G5cvW}FzE{LKd@B&@aijTkH(cbw}M zaQb!fw;wo@WY_9L%Dt;!kN_{0EZ-yga?6b0bpTugU+pw=9k2D?nA)h^Y>$;hP9EJU zjo+UKAMLA)OQ#$tBtCIM0eIEs_usPl!wsxQw|~EV_xtrv{NyD3pZLkC^AA5c5k}3X z{`uW+vVVHr>d~u4uX3`xSz9GXcG}c);x?t~uiThXmf-i3?z>=~m=5zxhYOlzD$qrg z$*XCN`Ot`COsTwD{iA)AW6or1PWC8e1DF`@*L$+F^_O>wli8 z51TZiWH#lq=3aZ>RJnUjtxR{7SHDEiJXkq2s72psJkkU1lr^1MreZqfrHAKujf7X5VGnS8>Ca`GAHM=@?Vd|t@Fa(#wK_{FR=Y2FYmsK(mxjN(9oOQ^!=z|7!CLjCI zr>+ILP_b))LHOP@yXQX+6$%XJtr1L`g56j+5;`pP>K}S4jcNHW5aXZ6dMf|=4hxj9 z)mF7)c+O^O4;lt;wz*srdQ$X6YP^vcJKo}Ql-+Q%#6+=NAfIJ{5tKK zh7Je6URs~`d}(mc(T|w7Dk!1cW|#Ubviz5Dyl?dPJ|d1XzLy4znzvjZ&E2q!czjur z3DYv4Q5e^_Pv=&*y}78E$Kpqs$21DGgk?^=Lk*t=#KS0hX^Eb4GJBQY4UCf?e@M74 zv3Z8!mBzEH0m_nSo(9XNbb^Gd!Wrxobk_aj}C>v4u)4tIW~u zf~X&NHJ$}#vps85dN zidY>hK4=cTrBLP&)UCN~*FAn~b!DoZA4Q)K5})e?Kwb3@=_C9Dea?WFA+BI}-N2rP57zYGTcKWCaA_Xk8c;{b;S|t=&4RaS4 zvNy-Y@2!dRImRJ4xKYA{EJc3hQGKoUlGz_*b)tw(7qqq3(dBq& z;v?-gpWVY**v%Do#g-(sg6}CDw(&4rtznyz@Qt^U6-bx8GEd7CegFRR>*k|#qzY(x zA+#V>9nUAjXK&Vo7O$VOQ z^s^13ST#QpU_{*L2GcLyybeHO~UYewf)G-Or(+CSBKSO1Y-#m1+m zm5GP=R3TfoHYFqO1XM1&djkA{_bzR&kw@LEU(O)|nXI3ZQ`p~gck_1$O77X^7Ay|M z-mIJ!#`xHHdn6KV@6H5KBep99>jkfg78r$&q6mf{&2EDUS1)rt8rifL=Y!$LjskIx=D|_z^ zj(7BvYBCQn5;T}kU(l8?gzwASi82$&Z;Gj0dUsKDbw!rpijJEzIw9>eb0XG6)9mtF zfwJ7;m-pN!=1=41IP|&@ZE)Fs^6-iW`P zPJuz69cjWx`h=e@cfG0*nG;1&4XsR1#%D9^BffJ9$c=c#Lr^QqV`H9mDy zZ1Y~+Zs*idmZAAz*79a0(O^8UOTvV~CFwbz})9P2I-ssQ!-}B8N9zNhyQXLq&UuH;giz+HuH6uXY=>Ii#=J8PV?H^aNQ}!im zNMk8WS+Zse3F*$>l!UqMV(g*pd-f%gEl~*BB})?74cQuoC|js3+47qizx$pu=b5Lc zuYbnN>%4u=`EKXBuJ3hCu#S2touEVA!=J2AMg!U%1nfDqrmj!Wh4u!|D@%_^Qan)8 zZ;1{;@ZVT#OE7UdeyU+5VJJo5N9*e#c2i?p-;t(F<|CzXr&|x#nx6~i2rpz@5cqT= z01=bxQ(!<|%WLCH=8!mQ&=+LEVoITU2@=P$;qnVnBoc0YEqoVoH&!y-)L3?wg7u%68<#HqS1(?M;Xyl=T??9iz$iLJo&P7FGLVXu}?j@&=GY&eg^&85S5%+&;k=FUx`U8uk zkJvQvZfPEge{J~2=DwN5x_6n+L`k9avB z&Hs;fxR<>Zc98v#$MNqpU;wd{~r3}so2_>-7Bz|p$ zELE7Y#GYumuUeB?Q4jG7exu17Okb#|b0m;gM=RhaAGBFGEdFqKLr7uvo7lUouv4T` zst;R(=oPw6?)T}D48PZIO|{FtWM@s`yS}8rYpZinkaaseOU2uTm5e9CbnA9EhxihQ zPdtyU&iS~sIJz@=k}9DKez%RAorY2MFuHwaeR}QhhpCxUOB$=zTSC36bW&Uoo`~9;)uT8p+ zYjdSVNDFU1O`{7;r&%(0X_TBuq}JFDRYQD+CT3Zhdr><*7l&SKnAkL~20b*Ey`iV6`-O3bL6VY5;hq*S4&UCmh#jqq!S5hCfs`8t=YTyeU+=B178#2 zd<5Oq3a~rtd5rk8Xy3b?mYb$=!WDc(;AoEV9v(?E`ErdXy^(xwUiK!?3)5 zBo#R}tR+b@dn#II(6k|LLYZ<@q?wt^EI&!qM@zVPufTOM$lGhK)>YWG_e+&3nX_1x z(MqC&Z~TbyIB}Vr63QftYlqk&l8>~M(@duRf?-J+Rmw&_SaB{N4zN980(P`J% zv(6{JJE$jI$)Y{PWIOcCbgZ~eM0lN)4Mu)rmtpiMjZ`(^k~&C?$+}Omo5c@AAISq<`C{5a~vgvZ+71aKuy?joX zA5tC6N;{bCQ$JISL|A$SC1kKUthYcVPJ8gfn^YidgCTX?)LQN_d|KVl+S$Xp`rl?y zElX9tAhuTE3D)n3>6aNPI7jhX$wa%%zb4A;%a-DHYD4K}voBlq{sTshTO9r>m2=;! zWhR-e$l~RXxf;ho)j~pR8d9A{UTnF0lnl&^rP?`eT_>VoOp>VxZhbOp?)f~F>#r^; zipwqETIBOj(`Uv$=SUe@b-5UvIy>f_l}~2t&vvd;P+D>M{TsKQcj69BxfSy7*M}h* zF=JfMeIk;mY)|U>`o2SK>dp^;WISBW(U@v?j8@|*`3rNeN_BB^>Jb&02ZQTc-vkHd z=GQZyGj>Hn>{mq{d%sBO-Fe)GhY9{EH^d#LYG z$~R;)%8QqqdQIKUjGxYOy20y^9qA(yX}U{KN1^Q2+kGKLnTX@3HaOFUP<5QX=>|>k z?5j~aH_9ASnhZy*jCDUO`mIE*=bgR9;^;3xh{(aL^JRp-|$ok;;#^x5rCu$acg zXNE<$e2RoAHyacVwgB?HZ&sJkTR`-0k53*y7={OyHgGUvP3?xEZjAGU2z6>i?po~& zn=M$5bDF>%D(Jgn1ZT-W&S$2gP zmO7dtqNPY7qeCpyiM+M`m#^M5XjeLCV%aTIFt71p#v~RI*_B2rx!GDeNd;?Hl;{5u zJ5|D>^2%I=~HwzfZS+}S0= zEC?iF5A;>gR_n>%o*t({_ryUrU>|*PJo3^Jaq5eYn>--3lk-cIyycG?=x(j3wULE6 zd=5>MdeXhHJZ{FUM#jvvCee}dhQhL5Zrh~c>U%yi;>BmA>docss!pg{Ut?sz`+(-D zuKUy(kT)f@V~YJzSKYshwG~nFTM`FqzCy^3-A%*Qt#=$L@Z9*^vzzWOX_L7VNqnuw9U7y4$JuN ziMn9hxGLc%1=UyG;*w3v8n*Zq16F30xtn$#KFNxGv-+v`vavyCc-C8a6b+43mAM6- z9|A(YVWLjilc$^jArTravtx215}>b8*?lvTt#-|Daev@@hy*iImD8_(P7KwhT_dW0 zxHYRq?>J{-HaxpRobtpqHZC%U&ul#xvcy|d>|3~xb0$SB>s^9G4%m)BU7kIXJ^Yn< zX>B2C=*Rj}`8%_D*ViWUUA>{A<0f@}?~S&k@{k*-I}#R>kKyt+*kPaiN{2i15%%zD znDk!9Td{hD;-~3v@19IO36$e)=_oN%b#jDTIl^yV^K!Owg<@R_DdUdgV$|q_{}uB&_pl^)vtgP$-eqrQ z4b73nrSuH!Wk+e#q`Y=IbS2DfNYzpU=G)3IM@)(~=M;=0+t!8{NEb33S`8GKm`rEk zu$1VZTdQs^!a^3M$Gu@?bY}c_r2Qp*M!fTcj#dqZC6D?i#$SshEpX>@s1Ug>)>Oux zJP>HwOh3C0_vTA_8I^oWzLoaJDC~Mu-IoQ#@(yhs9m;o75LKYvW}dTM9^Uq-D_gXe zb<=V^F16!^CTA(PZ^2jj21QaeVVakt zid028nOn%PL058koN3VacvwqtWOv~?xk)iY+QM28@M1SFwGK9I`pRP$~0 z*ACqmeYc3(#d=$3K6mprqWfMKYv~?3%t*vYL{3glB*PbELS|tTNP_+w8T$Q%6+NZ! z8+w9Sy4>=>7`^{|!1|f4X#Xii8_gy_!4}o|b4oILDIh?>7UcMI$`tw@2~e=r2L7BP zLPC&&%~bw#iU)cVCZGkI;rHj17PK`-fPzge`f~~m1wjfnSK`kp#*_pp*qGDrD41Xk zuH~ZFbpm&E5V`?lfM33cKX_NrRbJ2-KJccaC!kt$l%L5 zW(}qQWMNqZpBAuQ2)-;v7#N6^MR2|X3vS@c^0Wa{0J3nEM*jS1`PzbkSXl%Y0ATh% zzN`QTFa;nBSM~rb3z$WW|Fl9K!9c7mf-4Fz4HsWltP7X|kcG=X2Fn8GL*mOyas>mi zvIs6wz=T44Sx|Q{1t1HTNePw(%oxL$W#It^Vr3Cr^MEN=__7APz!ZQiT&@{d7BG~L zFKgHv48+PJxQGH{;rO!pV!;%EEZiXguqSpqOBoRn1pUK~{#TTRnBV>niKb>z From ddc5eb7a19f38852327bd9e4a533ae6a1ae4311f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 13:52:50 +0200 Subject: [PATCH 30/37] document template improvements for PR 718 --- CHANGELOG.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c989456f..7d7518dc 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Updated iGenomes paths [#710](https://github.com/BU-ISCIII/buisciii-tools/pull/710) - Updated the assembly template to skip busco and some paths in correspondence to refgenie system [#711](https://github.com/BU-ISCIII/buisciii-tools/pull/711) - Updated assembly template to nf-core/bacass 2.6.1 and restored BUSCO execution with matching config changes [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) +- Fixed conditional trimmed-read output and corrected MultiQC result paths in the assembly template [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) +- Added Bacass numeric sample ID compatibility to characterization, Snippy, and PlasmidID templates [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) +- Updated the Snippy IQ-TREE workflow to reuse the detected best-fit model [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) +- Enhanced the outbreak summary workbook with normalized sample IDs, updated AMRFinderPlus parsing, SNP distance tables, provenance, and consistent formatting [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) ### Modules @@ -919,4 +923,3 @@ Code contributions to the inital release: - [Alberto Lema](https://github.com/Alema91) - [Daniel Valle](https://github.com/Daniel-VM) - [Fernando Gomez](https://github.com/FGomez-Aldecoa) - From acfa2b00d4e97713de2fe1ee79318e75c9392466 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 15:05:28 +0200 Subject: [PATCH 31/37] linting --- .../RESULTS/generate_summary_outbreak.py | 34 +++++++++++++------ 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py b/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py index 5ffe6de2..2b39443e 100644 --- a/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py +++ b/buisciii/templates/characterization/RESULTS/generate_summary_outbreak.py @@ -332,14 +332,22 @@ def read_amrfinder_results(directory): # AMRFinderPlus renamed these fields in recent releases. gene_column = next( - (column for column in ("Element symbol", "Gene symbol") if column in df.columns), + ( + column + for column in ("Element symbol", "Gene symbol") + if column in df.columns + ), None, ) if gene_column: # The output can also contain virulence/stress hits; only AMR # elements belong in the resistance summary. type_column = next( - (column for column in ("Type", "Element type") if column in df.columns), + ( + column + for column in ("Type", "Element type") + if column in df.columns + ), None, ) if type_column: @@ -590,7 +598,9 @@ def read_fasta_alignment(file_path): sequences.append("".join(current_sequence).upper()) current_name = line[1:].split()[0] if not current_name: - raise ValueError(f"Empty FASTA identifier at {file_path}:{line_number}") + raise ValueError( + f"Empty FASTA identifier at {file_path}:{line_number}" + ) current_sequence = [] else: if current_name is None: @@ -669,8 +679,12 @@ def read_snippy_mapping_reference(commands_file): ] else: candidates = [ - os.path.normpath(os.path.join(os.path.dirname(commands_file), reference)), - os.path.normpath(os.path.join("../REFERENCES", os.path.basename(reference))), + os.path.normpath( + os.path.join(os.path.dirname(commands_file), reference) + ), + os.path.normpath( + os.path.join("../REFERENCES", os.path.basename(reference)) + ), ] resolved = next((path for path in candidates if os.path.isfile(path)), None) resolved_files.append(resolved or reference) @@ -702,7 +716,10 @@ def calculate_pairwise_nucleotide_differences(sequences, chunk_size=100000): independently for each pair. """ sequence_array = np.vstack( - [np.frombuffer(sequence.encode("ascii"), dtype=np.uint8) for sequence in sequences] + [ + np.frombuffer(sequence.encode("ascii"), dtype=np.uint8) + for sequence in sequences + ] ) sequence_count, alignment_length = sequence_array.shape differences = np.zeros((sequence_count, sequence_count), dtype=np.int64) @@ -1048,10 +1065,7 @@ def apply_homogeneous_formatting(wb): # Determine the real populated range rather than using template-only # formatting, which can extend hundreds of blank rows. populated_cells = [ - cell - for row in ws.iter_rows() - for cell in row - if cell.value is not None + cell for row in ws.iter_rows() for cell in row if cell.value is not None ] if populated_cells: last_row = max(cell.row for cell in populated_cells) From f8d7470c561f41780e4c204bfd52b54f9293d398 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 17:31:51 +0200 Subject: [PATCH 32/37] add Snippy SNP distance and close-pair QC tools --- .../05-snippy/evaluate_close_pairs.py | 853 ++++++++++++++++++ .../generate_snp_distance_matrices.py | 157 ++++ .../ANALYSIS01_SNIPPY/05-snippy/lablog | 48 +- 3 files changed, 1057 insertions(+), 1 deletion(-) create mode 100644 buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py create mode 100644 buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py new file mode 100644 index 00000000..f8f54b22 --- /dev/null +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py @@ -0,0 +1,853 @@ +#!/usr/bin/env python3 +"""Find closely related Snippy samples and report their differing SNPs. + +Distances are calculated directly from snippy-core's core.tab. For each pair +below the requested threshold, the script looks up supporting read counts in +each sample's snps.vcf (or snps.vcf.gz). Reference calls, which are absent +from variant-only VCFs, are measured from snps.bam with samtools mpileup. +""" + +from __future__ import annotations + +import argparse +import csv +import gzip +import itertools +import shutil +import subprocess +import sys +import tempfile +from pathlib import Path + + +VALID_BASES = frozenset("ACGT") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Calculate pairwise SNP distances from core.tab and extract VCF " + "evidence for variants in pairs with fewer than N differences." + ) + ) + parser.add_argument( + "snippy_dir", + nargs="?", + type=Path, + default=Path("."), + help="snippy-core output directory (default: current directory)", + ) + parser.add_argument( + "--core-tab", + type=Path, + help="core.tab path (default: SNIPPY_DIR/core.tab)", + ) + parser.add_argument( + "--threshold", + type=int, + default=20, + help="retain pairs with fewer than this many SNPs (default: 20)", + ) + parser.add_argument( + "--output-prefix", + type=Path, + default=Path("close_pair_snp_qc"), + help="output path prefix (default: close_pair_snp_qc)", + ) + parser.add_argument( + "--matrices-dir", + type=Path, + help="directory containing generated SNP-distance matrices (default: SNIPPY_DIR)", + ) + parser.add_argument( + "--gubbins-gff", + type=Path, + help="Gubbins recombination GFF (default: SNIPPY_DIR/gubbins.recombination_predictions.gff)", + ) + parser.add_argument( + "--missing-vcf", + choices=("error", "warn"), + default="error", + help="action when a selected sample has no VCF (default: error)", + ) + parser.add_argument( + "--reference", + type=Path, + help="reference FASTA override for mpileup (default: each sample's reference/ref.fa)", + ) + parser.add_argument( + "--samtools", + default="samtools", + help="samtools executable or path (default: samtools)", + ) + parser.add_argument( + "--no-bam-evidence", + action="store_true", + help="do not obtain evidence for reference calls from snps.bam", + ) + parser.add_argument( + "--pileup-mapq", + type=int, + default=10, + help="minimum mapping quality for mpileup evidence (default: 10)", + ) + parser.add_argument( + "--pileup-baseq", + type=int, + default=5, + help="minimum base quality for mpileup evidence (default: 5)", + ) + parser.add_argument( + "--warn-min-depth", + type=float, + default=10, + help="warn when evidence depth is below this value (default: 10)", + ) + parser.add_argument( + "--warn-min-called-fraction", + type=float, + default=0.90, + help="warn when the called allele fraction is below this value (default: 0.90)", + ) + parser.add_argument( + "--no-qc-warnings", + action="store_true", + help="suppress per-site QC warnings on standard output", + ) + parser.add_argument( + "--warn-min-strand-fraction", + type=float, + default=0.10, + help="warn if less than this fraction of called reads is on either strand (default: 0.10)", + ) + parser.add_argument( + "--warn-indel-distance", + type=int, + default=5, + help="warn for SNPs within this many bases of a raw-VCF indel/complex call (default: 5)", + ) + parser.add_argument( + "--warn-homopolymer-run", + type=int, + default=5, + help="warn when FreeBayes RUN is at least this length (default: 5)", + ) + parser.add_argument( + "--cluster-window", + type=int, + default=100, + help="window size in bases for local pairwise SNP clusters (default: 100)", + ) + parser.add_argument( + "--cluster-min-snps", + type=int, + default=3, + help="minimum pairwise SNPs in the window to flag a cluster (default: 3)", + ) + return parser.parse_args() + + +def read_core_tab(path: Path) -> tuple[list[str], list[tuple[str, int, str, list[str]]]]: + with path.open(newline="") as handle: + reader = csv.reader(handle, delimiter="\t") + try: + header = next(reader) + except StopIteration: + raise ValueError(f"Empty core table: {path}") + if len(header) < 4 or header[:3] != ["CHR", "POS", "REF"]: + raise ValueError(f"Unexpected core.tab header in {path}") + samples = header[3:] + if len(samples) != len(set(samples)): + raise ValueError("Duplicate sample names in core.tab") + sites = [] + for line_number, row in enumerate(reader, 2): + if not row: + continue + if len(row) != len(header): + raise ValueError( + f"Line {line_number} has {len(row)} columns; expected {len(header)}" + ) + sites.append((row[0], int(row[1]), row[2].upper(), [x.upper() for x in row[3:]])) + return samples, sites + + +def find_close_pairs(samples, sites, threshold): + results = [] + for left, right in itertools.combinations(range(len(samples)), 2): + differences = [] + comparable = 0 + ignored = 0 + for chrom, pos, ref, calls in sites: + a, b = calls[left], calls[right] + if a in VALID_BASES and b in VALID_BASES: + comparable += 1 + if a != b: + differences.append((chrom, pos, ref, a, b)) + else: + ignored += 1 + if len(differences) < threshold: + results.append( + { + "sample_1": samples[left], + "sample_2": samples[right], + "snp_distance": len(differences), + "comparable_core_sites": comparable, + "ignored_non_acgt_sites": ignored, + "differences": differences, + } + ) + return sorted(results, key=lambda x: (x["snp_distance"], x["sample_1"], x["sample_2"])) + + +def read_distance_matrix(path: Path) -> dict[frozenset[str], int]: + with path.open(newline="") as handle: + reader = csv.reader(handle, delimiter="\t") + try: + header = next(reader) + except StopIteration: + raise ValueError(f"Empty distance matrix: {path}") + names = header[1:] + if not names or len(names) != len(set(names)): + raise ValueError(f"Invalid or duplicate matrix sample names: {path}") + distances = {} + observed_rows = [] + for row_number, row in enumerate(reader, 2): + if len(row) != len(header): + raise ValueError(f"Wrong column count at {path}:{row_number}") + row_name = row[0] + observed_rows.append(row_name) + for column, value in enumerate(row[1:]): + if row_name != names[column]: + distances[frozenset((row_name, names[column]))] = int(value) + if observed_rows != names: + raise ValueError(f"Matrix row and column sample order differs: {path}") + return distances + + +def read_gubbins_intervals(path: Path) -> list[tuple[int, int]]: + """Read and merge 1-based inclusive recombinant intervals from Gubbins.""" + intervals = [] + with path.open() as handle: + for line_number, line in enumerate(handle, 1): + if not line.strip() or line.startswith("#"): + continue + fields = line.rstrip("\n").split("\t") + if len(fields) < 5: + raise ValueError(f"Malformed Gubbins GFF at {path}:{line_number}") + intervals.append((int(fields[3]), int(fields[4]))) + merged = [] + for start, end in sorted(intervals): + if merged and start <= merged[-1][1] + 1: + merged[-1] = (merged[-1][0], max(merged[-1][1], end)) + else: + merged.append((start, end)) + return merged + + +def position_in_intervals(position: int, intervals: list[tuple[int, int]]) -> bool: + import bisect + + starts = [start for start, _ in intervals] + index = bisect.bisect_right(starts, position) - 1 + return index >= 0 and position <= intervals[index][1] + + +def open_text(path: Path): + return gzip.open(path, "rt") if path.suffix == ".gz" else path.open() + + +def parse_list(value: str, cast=int): + if value in ("", "."): + return [] + result = [] + for item in value.split(","): + try: + result.append(cast(item)) + except ValueError: + result.append(None) + return result + + +def scalar_or_alt(value: str | None, alt_index: int): + if value in (None, "", "."): + return "" + values = value.split(",") + return values[alt_index] if alt_index < len(values) else value + + +def read_vcf(path: Path) -> dict[tuple[str, int, str], dict[str, str]]: + """Index SNP evidence by (chromosome, position, called base). + + Equal-length MNP/complex substitutions are expanded base-by-base so that + sites emitted separately in core.tab still find their source VCF record. + """ + index = {} + with open_text(path) as handle: + for line in handle: + if not line or line.startswith("#"): + continue + fields = line.rstrip("\n").split("\t") + if len(fields) < 8: + continue + chrom, pos_text, _, ref, alt_text, qual, filt, info_text = fields[:8] + pos = int(pos_text) + info = {} + for entry in info_text.split(";"): + key, sep, value = entry.partition("=") + if sep: + info[key] = value + fmt = {} + if len(fields) >= 10: + fmt = dict(zip(fields[8].split(":"), fields[9].split(":"))) + dp = fmt.get("DP", info.get("DP", "")) + ro = fmt.get("RO", info.get("RO", "")) + ao_all = fmt.get("AO", info.get("AO", "")) + for alt_index, alt in enumerate(alt_text.split(",")): + evidence = { + "vcf_record_pos": pos_text, + "vcf_ref": ref, + "vcf_alt": alt, + "qual": qual if qual != "." else "", + "filter": filt, + "genotype": fmt.get("GT", ""), + "depth": dp, + "ref_count": ro, + "alt_count": scalar_or_alt(ao_all, alt_index), + "variant_type": scalar_or_alt(info.get("TYPE"), alt_index), + "ref_forward_count": info.get("SRF", ""), + "ref_reverse_count": info.get("SRR", ""), + "alt_forward_count": scalar_or_alt(info.get("SAF"), alt_index), + "alt_reverse_count": scalar_or_alt(info.get("SAR"), alt_index), + "homopolymer_run": info.get("RUN", ""), + "nearest_indel_distance": "", + } + if len(ref) == len(alt): + for offset, (ref_base, alt_base) in enumerate(zip(ref.upper(), alt.upper())): + if ref_base != alt_base and alt_base in VALID_BASES: + index[(chrom, pos + offset, alt_base)] = evidence + elif len(ref) == 1 and len(alt) == 1 and alt.upper() in VALID_BASES: + index[(chrom, pos, alt.upper())] = evidence + return index + + +def locate_vcf(snippy_dir: Path, sample: str) -> Path | None: + for name in ("snps.vcf", "snps.vcf.gz"): + candidate = snippy_dir / sample / name + if candidate.is_file(): + return candidate + return None + + +def locate_raw_vcf(snippy_dir: Path, sample: str) -> Path | None: + for name in ("snps.raw.vcf", "snps.raw.vcf.gz"): + candidate = snippy_dir / sample / name + if candidate.is_file(): + return candidate + return None + + +def read_non_snp_positions(path: Path) -> dict[str, list[int]]: + positions: dict[str, list[int]] = {} + with open_text(path) as handle: + for line in handle: + if line.startswith("#"): + continue + fields = line.rstrip("\n").split("\t") + if len(fields) < 8: + continue + info = dict( + entry.split("=", 1) for entry in fields[7].split(";") if "=" in entry + ) + types = info.get("TYPE", "").split(",") + if any(kind and kind != "snp" for kind in types): + positions.setdefault(fields[0], []).append(int(fields[1])) + for values in positions.values(): + values.sort() + return positions + + +def nearest_distance(sorted_positions: list[int], position: int) -> int | None: + import bisect + + index = bisect.bisect_left(sorted_positions, position) + candidates = sorted_positions[max(0, index - 1):index + 1] + return min((abs(position - other) for other in candidates), default=None) + + +def parse_pileup_bases(text: str, reference: str) -> dict[str, dict[str, int]]: + """Count A/C/G/T observations in a samtools mpileup read-bases field.""" + counts = {base: {"forward": 0, "reverse": 0} for base in VALID_BASES} + i = 0 + while i < len(text): + symbol = text[i] + if symbol == "^": # start of read plus one mapping-quality character + i += 2 + elif symbol == "$": + i += 1 + elif symbol in "+-": # indel length and inserted/deleted sequence + i += 1 + start = i + while i < len(text) and text[i].isdigit(): + i += 1 + length = int(text[start:i]) if i > start else 0 + i += length + elif symbol in ".,": + if reference in counts: + strand = "forward" if symbol == "." else "reverse" + counts[reference][strand] += 1 + i += 1 + elif symbol.upper() in counts: + strand = "forward" if symbol.isupper() else "reverse" + counts[symbol.upper()][strand] += 1 + i += 1 + else: # N/n, deletion placeholders (* or #), and reference skips + i += 1 + return counts + + +def read_bam_pileup( + samtools: str, + reference: Path, + bam: Path, + positions: set[tuple[str, int]], + minimum_mapq: int, + minimum_baseq: int, +) -> dict[tuple[str, int], dict[str, object]]: + if not positions: + return {} + with tempfile.NamedTemporaryFile("w", suffix=".bed") as bed: + for chrom, pos in sorted(positions): + bed.write(f"{chrom}\t{pos - 1}\t{pos}\n") + bed.flush() + command = [ + samtools, "mpileup", + "-q", str(minimum_mapq), "-Q", str(minimum_baseq), "-f", str(reference), + "-l", bed.name, str(bam), + ] + result = subprocess.run(command, text=True, capture_output=True) + if result.returncode: + detail = result.stderr.strip() or "unknown samtools error" + raise RuntimeError(f"samtools mpileup failed for {bam}: {detail}") + index = {} + for line in result.stdout.splitlines(): + fields = line.split("\t") + if len(fields) < 5: + continue + chrom, pos_text, reference_base, depth, read_bases = fields[:5] + reference_base = reference_base.upper() + index[(chrom, int(pos_text))] = { + "depth": depth, + "counts": parse_pileup_bases(read_bases, reference_base), + } + return index + + +def evidence_for(index, bam_index, chrom, pos, core_ref, call, other_call): + record = index.get((chrom, pos, call)) + if record is not None: + return "variant_record", record + if call == core_ref and (chrom, pos) in bam_index: + pileup = bam_index[(chrom, pos)] + counts = pileup["counts"] + ref_strands = counts.get(core_ref, {"forward": 0, "reverse": 0}) + alt_strands = counts.get(other_call, {"forward": 0, "reverse": 0}) + evidence = {key: "" for key in EVIDENCE_COLUMNS} + evidence.update( + { + "vcf_ref": core_ref, + "vcf_alt": other_call, + "depth": pileup["depth"], + "ref_count": sum(ref_strands.values()), + "alt_count": sum(alt_strands.values()), + "variant_type": "reference_call", + "ref_forward_count": ref_strands["forward"], + "ref_reverse_count": ref_strands["reverse"], + "alt_forward_count": alt_strands["forward"], + "alt_reverse_count": alt_strands["reverse"], + "nearest_indel_distance": pileup.get("nearest_indel_distance", ""), + } + ) + return "reference_call_mpileup", evidence + status = "reference_call_no_bam_evidence" if call == core_ref else "variant_not_found_in_vcf" + return status, {key: "" for key in EVIDENCE_COLUMNS} + + +EVIDENCE_COLUMNS = ( + "vcf_record_pos", "vcf_ref", "vcf_alt", "qual", "filter", "genotype", + "depth", "ref_count", "alt_count", "variant_type", + "ref_forward_count", "ref_reverse_count", "alt_forward_count", + "alt_reverse_count", "homopolymer_run", "nearest_indel_distance", +) + + +def add_fractions(row: dict[str, object], prefix: str) -> None: + try: + dp = float(row[f"{prefix}_depth"]) + except (ValueError, TypeError): + dp = 0 + for count_name in ("ref", "alt"): + try: + count = float(row[f"{prefix}_{count_name}_count"]) + value = count / dp if dp else "" + except (ValueError, TypeError): + value = "" + row[f"{prefix}_{count_name}_fraction"] = value + + +def numeric(value: object) -> float | None: + try: + return float(value) + except (TypeError, ValueError): + return None + + +def collect_qc_warnings( + row: dict[str, object], + sample_number: int, + minimum_depth: float, + minimum_called_fraction: float, + minimum_strand_fraction: float, + maximum_indel_distance: int, + minimum_homopolymer_run: int, +) -> list[tuple[str, str]]: + prefix = f"sample_{sample_number}" + sample = str(row[prefix]) + location = f"{row['chrom']}:{row['pos']}" + context = ( + f"sample={sample} position={location} call={row[f'{prefix}_call']} " + f"pair={row['pair_id']}" + ) + status = row[f"{prefix}_evidence_status"] + if status not in ("variant_record", "reference_call_mpileup"): + return [("missing_evidence", f"{context} status={status}")] + + warnings = [] + depth = numeric(row[f"{prefix}_depth"]) + if depth is None: + warnings.append(("missing_depth", context)) + elif depth < minimum_depth: + warnings.append( + ("low_depth", f"{context} depth={depth:g} threshold={minimum_depth:g}") + ) + + fraction_name = "ref_fraction" if row[f"{prefix}_call"] == row["core_ref"] else "alt_fraction" + called_fraction = numeric(row[f"{prefix}_{fraction_name}"]) + if called_fraction is None: + warnings.append(("missing_called_fraction", context)) + elif called_fraction < minimum_called_fraction: + warnings.append( + ( + "low_called_fraction", + f"{context} called_fraction={called_fraction:.4f} " + f"threshold={minimum_called_fraction:.4f}", + ) + ) + + called_prefix = "ref" if row[f"{prefix}_call"] == row["core_ref"] else "alt" + forward = numeric(row[f"{prefix}_{called_prefix}_forward_count"]) + reverse = numeric(row[f"{prefix}_{called_prefix}_reverse_count"]) + if forward is not None and reverse is not None and forward + reverse > 0: + minor_fraction = min(forward, reverse) / (forward + reverse) + if minor_fraction < minimum_strand_fraction: + warnings.append( + ( + "strand_bias", + f"{context} forward={forward:g} reverse={reverse:g} " + f"minor_strand_fraction={minor_fraction:.4f} " + f"threshold={minimum_strand_fraction:.4f}", + ) + ) + + indel_distance = numeric(row[f"{prefix}_nearest_indel_distance"]) + if indel_distance is not None and indel_distance <= maximum_indel_distance: + warnings.append( + ( + "near_indel_or_complex", + f"{context} distance={indel_distance:g} threshold={maximum_indel_distance}", + ) + ) + run = numeric(row[f"{prefix}_homopolymer_run"]) + if run is not None and run >= minimum_homopolymer_run: + warnings.append( + ( + "homopolymer_context", + f"{context} run_length={run:g} threshold={minimum_homopolymer_run}", + ) + ) + return warnings + + +def main() -> int: + args = parse_args() + if args.threshold <= 0: + raise ValueError("--threshold must be greater than zero") + if args.pileup_mapq < 0 or args.pileup_baseq < 0: + raise ValueError("pileup quality thresholds cannot be negative") + if args.warn_min_depth < 0: + raise ValueError("--warn-min-depth cannot be negative") + if not 0 <= args.warn_min_called_fraction <= 1: + raise ValueError("--warn-min-called-fraction must be between 0 and 1") + if not 0 <= args.warn_min_strand_fraction <= 0.5: + raise ValueError("--warn-min-strand-fraction must be between 0 and 0.5") + if args.warn_indel_distance < 0 or args.warn_homopolymer_run < 1: + raise ValueError("indel distance must be non-negative and homopolymer run positive") + if args.cluster_window < 1 or args.cluster_min_snps < 2: + raise ValueError("cluster window must be positive and --cluster-min-snps at least 2") + snippy_dir = args.snippy_dir.resolve() + core_tab = (args.core_tab or snippy_dir / "core.tab").resolve() + samples, sites = read_core_tab(core_tab) + pairs = find_close_pairs(samples, sites, args.threshold) + + chromosomes = {chrom for chrom, *_ in sites} + gubbins_path = (args.gubbins_gff or snippy_dir / "gubbins.recombination_predictions.gff").resolve() + gubbins_intervals = [] + if gubbins_path.is_file(): + if len(chromosomes) == 1: + gubbins_intervals = read_gubbins_intervals(gubbins_path) + else: + print( + "WARNING: Gubbins coordinates cannot be mapped automatically to a " + "multi-contig core.tab; removed variant positions will be unavailable", + file=sys.stderr, + ) + else: + print(f"WARNING: Gubbins GFF not found: {gubbins_path}", file=sys.stderr) + + matrices_dir = (args.matrices_dir or snippy_dir).resolve() + matrix_specs = ( + ("core_tab_matrix_distance", "core_tab_snp_distances.tsv"), + ("phylo_aln_distance", "phylo_snp_distances.tsv"), + ("clean_core_aln_distance", "clean_core_snp_distances.tsv"), + ) + matrices = {} + for column, filename in matrix_specs: + matrix_path = matrices_dir / filename + if matrix_path.is_file(): + matrices[column] = read_distance_matrix(matrix_path) + else: + print(f"WARNING: comparison matrix not found: {matrix_path}", file=sys.stderr) + matrices[column] = {} + for pair in pairs: + key = frozenset((pair["sample_1"], pair["sample_2"])) + for column, _ in matrix_specs: + pair[column] = matrices[column].get(key, "") + phylo_distance = pair["phylo_aln_distance"] + clean_distance = pair["clean_core_aln_distance"] + pair["core_tab_matrix_matches"] = ( + "yes" if pair["core_tab_matrix_distance"] == pair["snp_distance"] else + ("" if pair["core_tab_matrix_distance"] == "" else "no") + ) + pair["phylo_matches_core_tab"] = ( + "yes" if phylo_distance == pair["snp_distance"] else + ("" if phylo_distance == "" else "no") + ) + pair["snps_removed_by_gubbins"] = ( + pair["snp_distance"] - clean_distance if clean_distance != "" else "" + ) + pair["clean_distance_valid"] = ( + "yes" if clean_distance != "" and clean_distance <= pair["snp_distance"] else + ("" if clean_distance == "" else "no") + ) + comparison_warnings = [] + if pair["core_tab_matrix_matches"] == "no": + comparison_warnings.append("core_tab_matrix_mismatch") + if pair["phylo_matches_core_tab"] == "no": + comparison_warnings.append("phylo_core_tab_mismatch") + if pair["clean_distance_valid"] == "no": + comparison_warnings.append("clean_distance_greater_than_core") + if any(pair[column] == "" for column, _ in matrix_specs): + comparison_warnings.append("comparison_matrix_missing") + pair["distance_comparison_warnings"] = ";".join(comparison_warnings) + + needed_samples = sorted({p[key] for p in pairs for key in ("sample_1", "sample_2")}) + vcf_indexes = {} + non_snp_positions = {} + for sample in needed_samples: + path = locate_vcf(snippy_dir, sample) + if path is None: + message = f"No snps.vcf or snps.vcf.gz found for sample {sample}" + if args.missing_vcf == "error": + raise FileNotFoundError(message) + print(f"WARNING: {message}", file=sys.stderr) + vcf_indexes[sample] = {} + non_snp_positions[sample] = {} + else: + vcf_indexes[sample] = read_vcf(path) + raw_path = locate_raw_vcf(snippy_dir, sample) + if raw_path is not None: + raw_index = read_vcf(raw_path) + contextual_fields = ( + "ref_forward_count", "ref_reverse_count", "alt_forward_count", + "alt_reverse_count", "homopolymer_run", + ) + for key, evidence in vcf_indexes[sample].items(): + raw_evidence = raw_index.get(key, {}) + for field in contextual_fields: + evidence[field] = raw_evidence.get(field, "") + non_snp_positions[sample] = read_non_snp_positions(raw_path) + else: + print( + f"WARNING: no snps.raw.vcf(.gz) for {sample}; strand, " + "homopolymer, and nearby-indel checks will be unavailable", + file=sys.stderr, + ) + non_snp_positions[sample] = {} + + for sample, index in vcf_indexes.items(): + sample_non_snps = non_snp_positions[sample] + for (chrom, pos, _), evidence in index.items(): + distance = nearest_distance(sample_non_snps.get(chrom, []), pos) + evidence["nearest_indel_distance"] = "" if distance is None else distance + + bam_indexes = {sample: {} for sample in needed_samples} + if not args.no_bam_evidence: + samtools = shutil.which(args.samtools) if Path(args.samtools).name == args.samtools else args.samtools + if not samtools: + raise FileNotFoundError( + f"samtools executable not found: {args.samtools!r}; load samtools, use " + "--samtools /path/to/samtools, or use --no-bam-evidence" + ) + reference_positions = {sample: set() for sample in needed_samples} + for pair in pairs: + for chrom, pos, ref, call_1, call_2 in pair["differences"]: + if call_1 == ref: + reference_positions[pair["sample_1"]].add((chrom, pos)) + if call_2 == ref: + reference_positions[pair["sample_2"]].add((chrom, pos)) + for sample, positions in reference_positions.items(): + if not positions: + continue + bam = snippy_dir / sample / "snps.bam" + if not bam.is_file(): + raise FileNotFoundError(f"BAM not found for sample {sample}: {bam}") + reference = ( + args.reference.resolve() + if args.reference + else snippy_dir / sample / "reference" / "ref.fa" + ) + if not reference.is_file(): + raise FileNotFoundError(f"Reference FASTA not found: {reference}") + bam_indexes[sample] = read_bam_pileup( + str(samtools), reference, bam, positions, + args.pileup_mapq, args.pileup_baseq, + ) + for (chrom, pos), pileup in bam_indexes[sample].items(): + distance = nearest_distance(non_snp_positions[sample].get(chrom, []), pos) + pileup["nearest_indel_distance"] = "" if distance is None else distance + + prefix = args.output_prefix + if not prefix.is_absolute(): + prefix = Path.cwd() / prefix + prefix.parent.mkdir(parents=True, exist_ok=True) + pairs_path = prefix.with_name(prefix.name + "_pairs.tsv") + variants_path = prefix.with_name(prefix.name + "_variants.tsv") + + pair_columns = [ + "pair_id", "sample_1", "sample_2", "snp_distance", + "comparable_core_sites", "ignored_non_acgt_sites", + "core_tab_matrix_distance", "core_tab_matrix_matches", + "phylo_aln_distance", "phylo_matches_core_tab", + "clean_core_aln_distance", "snps_removed_by_gubbins", + "clean_distance_valid", "distance_comparison_warnings", + ] + with pairs_path.open("w", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=pair_columns, delimiter="\t") + writer.writeheader() + for number, pair in enumerate(pairs, 1): + writer.writerow({"pair_id": f"pair_{number:04d}", **{k: pair[k] for k in pair_columns[1:]}}) + + base_columns = [ + "pair_id", "sample_1", "sample_2", "snp_distance", "chrom", "pos", + "core_ref", "sample_1_call", "sample_2_call", "pair_local_snp_count", + "pair_clustered_snp", "removed_by_gubbins", "qc_warnings", + ] + evidence_output = [] + for n in (1, 2): + evidence_output.extend( + [f"sample_{n}_evidence_status"] + + [f"sample_{n}_{key}" for key in EVIDENCE_COLUMNS] + + [f"sample_{n}_ref_fraction", f"sample_{n}_alt_fraction"] + ) + qc_warnings: dict[tuple[str, str, int, str], tuple[str, str]] = {} + with variants_path.open("w", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=base_columns + evidence_output, delimiter="\t") + writer.writeheader() + for number, pair in enumerate(pairs, 1): + for chrom, pos, ref, call_1, call_2 in pair["differences"]: + local_snp_count = sum( + other_chrom == chrom and abs(other_pos - pos) <= args.cluster_window + for other_chrom, other_pos, *_ in pair["differences"] + ) + row = { + "pair_id": f"pair_{number:04d}", + "sample_1": pair["sample_1"], "sample_2": pair["sample_2"], + "snp_distance": pair["snp_distance"], "chrom": chrom, "pos": pos, + "core_ref": ref, "sample_1_call": call_1, "sample_2_call": call_2, + "pair_local_snp_count": local_snp_count, + "pair_clustered_snp": "yes" if local_snp_count >= args.cluster_min_snps else "no", + "removed_by_gubbins": ( + "yes" if position_in_intervals(pos, gubbins_intervals) else "no" + ), + } + row_warning_types = [] + if row["removed_by_gubbins"] == "yes": + row_warning_types.append("gubbins_removed") + if local_snp_count >= args.cluster_min_snps: + row_warning_types.append("local_snp_cluster") + message = ( + f"pair={row['pair_id']} samples={pair['sample_1']},{pair['sample_2']} " + f"position={chrom}:{pos} local_snps={local_snp_count} " + f"window=+/-{args.cluster_window}" + ) + key = ("local_snp_cluster", str(row["pair_id"]), pos, chrom) + qc_warnings.setdefault(key, ("local_snp_cluster", message)) + for n, call, other_call in ((1, call_1, call_2), (2, call_2, call_1)): + sample = pair[f"sample_{n}"] + status, evidence = evidence_for( + vcf_indexes[sample], bam_indexes[sample], chrom, pos, ref, + call, other_call, + ) + row[f"sample_{n}_evidence_status"] = status + row.update({f"sample_{n}_{key}": value for key, value in evidence.items()}) + add_fractions(row, f"sample_{n}") + sample_warnings = collect_qc_warnings( + row, n, args.warn_min_depth, args.warn_min_called_fraction, + args.warn_min_strand_fraction, args.warn_indel_distance, + args.warn_homopolymer_run, + ) + for warning_type, message in sample_warnings: + row_warning_types.append(f"sample_{n}:{warning_type}") + # A sample/site can occur in several close pairs; report it once. + key = (warning_type, sample, pos, chrom) + qc_warnings.setdefault(key, (warning_type, message)) + # Preserve warning order while removing repeated labels. + row["qc_warnings"] = ";".join(dict.fromkeys(row_warning_types)) + writer.writerow(row) + + variant_count = sum(len(pair["differences"]) for pair in pairs) + print(f"Samples: {len(samples)}") + print(f"Core SNP sites: {len(sites)}") + print(f"Pairs with < {args.threshold} SNPs: {len(pairs)}") + print(f"Pair-specific variant rows: {variant_count}") + print(f"Wrote: {pairs_path}") + print(f"Wrote: {variants_path}") + if not args.no_qc_warnings: + for warning_type, message in qc_warnings.values(): + print(f"WARNING [{warning_type}] {message}") + warning_counts = {} + for warning_type, _ in qc_warnings.values(): + warning_counts[warning_type] = warning_counts.get(warning_type, 0) + 1 + if warning_counts: + summary = ", ".join(f"{key}={value}" for key, value in sorted(warning_counts.items())) + print(f"Unique QC warnings: {summary}") + else: + print("Unique QC warnings: none") + return 0 + + +if __name__ == "__main__": + try: + raise SystemExit(main()) + except (OSError, RuntimeError, ValueError) as error: + print(f"ERROR: {error}", file=sys.stderr) + raise SystemExit(2) diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py new file mode 100644 index 00000000..0c2ab915 --- /dev/null +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +"""Generate reusable SNP-distance matrices from Snippy/Gubbins outputs.""" + +from __future__ import annotations + +import argparse +import csv +import json +from pathlib import Path + +import numpy as np + + +VALID_BASES = np.frombuffer(b"ACGT", dtype=np.uint8) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Generate core.tab, phylo.aln, and clean.core.aln distance matrices." + ) + parser.add_argument( + "snippy_dir", nargs="?", type=Path, default=Path("."), + help="Snippy output directory (default: current directory)", + ) + parser.add_argument( + "--output-dir", type=Path, + help="matrix output directory (default: SNIPPY_DIR)", + ) + return parser.parse_args() + + +def read_fasta_alignment(path: Path) -> tuple[list[str], list[str]]: + names, sequences = [], [] + current_name = None + current_sequence: list[str] = [] + with path.open() as handle: + for line_number, raw_line in enumerate(handle, 1): + line = raw_line.strip() + if not line: + continue + if line.startswith(">"): + if current_name is not None: + names.append(current_name) + sequences.append("".join(current_sequence).upper()) + current_name = line[1:].split()[0] + if not current_name: + raise ValueError(f"Empty FASTA identifier at {path}:{line_number}") + current_sequence = [] + elif current_name is None: + raise ValueError(f"Sequence before first FASTA header at {path}:{line_number}") + else: + current_sequence.append(line) + if current_name is not None: + names.append(current_name) + sequences.append("".join(current_sequence).upper()) + validate_sequences(path, names, sequences) + return names, sequences + + +def read_core_tab(path: Path) -> tuple[list[str], list[str]]: + with path.open(newline="") as handle: + reader = csv.reader(handle, delimiter="\t") + try: + header = next(reader) + except StopIteration: + raise ValueError(f"Empty core table: {path}") + if header[:3] != ["CHR", "POS", "REF"] or len(header) < 4: + raise ValueError(f"Unexpected core.tab header: {path}") + names = header[3:] + sequence_parts = [[] for _ in names] + for line_number, row in enumerate(reader, 2): + if len(row) != len(header): + raise ValueError(f"Wrong column count at {path}:{line_number}") + for index, call in enumerate(row[3:]): + sequence_parts[index].append(call.upper()) + sequences = ["".join(parts) for parts in sequence_parts] + validate_sequences(path, names, sequences) + return names, sequences + + +def validate_sequences(path: Path, names: list[str], sequences: list[str]) -> None: + if not names: + raise ValueError(f"No sequences found in {path}") + if len(names) != len(set(names)): + raise ValueError(f"Duplicate sequence identifiers in {path}") + lengths = {len(sequence) for sequence in sequences} + if len(lengths) != 1: + raise ValueError(f"Sequences have different lengths in {path}: {sorted(lengths)}") + + +def calculate_distances(sequences: list[str], chunk_size: int = 100_000) -> np.ndarray: + """Count absolute differences with pairwise deletion of non-ACGT calls.""" + sequence_array = np.vstack( + [np.frombuffer(sequence.encode("ascii"), dtype=np.uint8) for sequence in sequences] + ) + sequence_count, alignment_length = sequence_array.shape + distances = np.zeros((sequence_count, sequence_count), dtype=np.int64) + for start in range(0, alignment_length, chunk_size): + block = sequence_array[:, start:start + chunk_size] + valid = np.isin(block, VALID_BASES) + for first in range(sequence_count - 1): + pair_valid = valid[first + 1:] & valid[first] + different = (block[first + 1:] != block[first]) & pair_valid + distances[first, first + 1:] += np.count_nonzero(different, axis=1) + distances += distances.T + return distances + + +def write_matrix(path: Path, names: list[str], distances: np.ndarray) -> None: + with path.open("w", newline="") as handle: + writer = csv.writer(handle, delimiter="\t", lineterminator="\n") + writer.writerow(["Sample", *names]) + for name, values in zip(names, distances): + writer.writerow([name, *(int(value) for value in values)]) + + +def main() -> int: + args = parse_args() + snippy_dir = args.snippy_dir.resolve() + output_dir = (args.output_dir or snippy_dir).resolve() + output_dir.mkdir(parents=True, exist_ok=True) + inputs = [ + ("core_tab", snippy_dir / "core.tab", "core_tab_snp_distances.tsv", read_core_tab), + ("phylo", snippy_dir / "phylo.aln", "phylo_snp_distances.tsv", read_fasta_alignment), + ("clean_core", snippy_dir / "clean.core.aln", "clean_core_snp_distances.tsv", read_fasta_alignment), + ] + metadata = { + "method": "Absolute A/C/G/T differences with pairwise deletion of other calls", + "matrices": {}, + } + for key, source, output_name, reader in inputs: + if not source.is_file(): + print(f"WARNING: input not found, skipping: {source}") + continue + names, sequences = reader(source) + matrix_path = output_dir / output_name + write_matrix(matrix_path, names, calculate_distances(sequences)) + metadata["matrices"][output_name] = { + "source": str(source), + "sample_count": len(names), + "site_count": len(sequences[0]), + } + print(f"Wrote: {matrix_path}") + metadata_path = output_dir / "snp_distance_matrices_metadata.json" + with metadata_path.open("w") as handle: + json.dump(metadata, handle, indent=2) + handle.write("\n") + print(f"Wrote: {metadata_path}") + return 0 + + +if __name__ == "__main__": + try: + raise SystemExit(main()) + except (OSError, ValueError) as error: + print(f"ERROR: {error}") + raise SystemExit(2) diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/lablog b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/lablog index 645a0a1f..f386cd96 100644 --- a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/lablog +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/lablog @@ -15,7 +15,7 @@ cat ../samples_id.txt | while read in; do echo -e "${in}\t${scratch_dir}/../02-p ## If you are going to include .fasta files in the analysis, run the following line. Run the previous line as well if you also have .fastq.gz files for your analysis. ## Bear in mind this line considers the .fasta files are inside a certain folder from REFERENCES! Change this line accordingly. -## ls ../../../REFERENCES/*/*.fasta | cut -d '/' -f6 | sed 's/.fasta//g' | while read in; do paste <(echo ${in}) <(echo ${scratch_dir}/../../../REFERENCES/*/${in}.fasta); done >> input.tab +ls ../../../REFERENCES/*/*.fasta | cut -d '/' -f6 | sed 's/.fasta//g' | while read in; do paste <(echo ${in}) <(echo ${scratch_dir}/../../../REFERENCES/*/${in}.fasta); done >> input.tab ls ../../../REFERENCES/*.fna | xargs -I %% snippy-multi ${scratch_dir}/input.tab --mincov 10 --mapqual 10 --basequal 5 --minqual 30 --ref %% --cpus 5 > commands.out @@ -55,3 +55,49 @@ echo "snp-sites -c gubbins.filtered_polymorphic_sites.fasta > clean.core.aln" >> ## Run gubbins echo "srun --chdir ${scratch_dir} --output logs/GUBBINS.%j.log --job-name GUBBINS --cpus-per-task 20 --mem 49152 --partition short_idx --time 02:00:00 bash _03_gubbins.sh &" > _03_run_gubbins.sh + +# BLOCK 6: CHECK SNP DIFFERENCES AND POSSIBLE CALLING ERRORS + +## Run after the Snippy, phylo.aln, and Gubbins/clean.core.aln steps are complete. +## The QC script calculates pairwise +## distances directly from core.tab, so the Excel distance matrix is not needed. +## It selects pairs with fewer than 20 SNP differences, reads ALT evidence from +## each sample's VCF, and uses samtools mpileup on snps.bam for REF evidence. +## Activate an environment containing samtools before running it (the Snippy +## environment should normally provide it). +# +## First generate the reusable distance matrices. These files are also consumed +## by RESULTS/generate_summary_outbreak.py, which does not recalculate distances. +echo "python3 generate_snp_distance_matrices.py ." > _04_snp_distance_matrices.sh + +## Run the SNP-distance matrix generation script. +echo "srun --chdir ${scratch_dir} --output logs/SNP_DISTANCE_MATRICES.%j.log --job-name SNP_DISTANCES --cpus-per-task 1 --mem 8192 --partition short_idx --time 02:00:00 bash _04_snp_distance_matrices.sh &" > _04_run_snp_distance_matrices.sh +# +# python3 evaluate_close_pairs.py . \ +# --threshold 20 \ +# --warn-min-depth 10 \ +# --warn-min-called-fraction 0.90 \ +# --warn-min-strand-fraction 0.10 \ +# --warn-indel-distance 5 \ +# --warn-homopolymer-run 5 \ +# --cluster-window 100 \ +# --cluster-min-snps 3 \ +# --no-qc-warnings \ +# --output-prefix close_pair_snp_qc +# +## Output files: +## core_tab_snp_distances.tsv distances calculated directly from core.tab +## phylo_snp_distances.tsv distances calculated from phylo.aln +## clean_core_snp_distances.tsv distances calculated from clean.core.aln +## snp_distance_matrices_metadata.json sources, sample counts, and site counts +## close_pair_snp_qc_pairs.tsv close pairs and their SNP distances +## close_pair_snp_qc_variants.tsv differing positions, read evidence, and a +## semicolon-separated qc_warnings column +## +## Possible qc_warnings values include low_depth, low_called_fraction, +## strand_bias, near_indel_or_complex, homopolymer_context, local_snp_cluster, +## missing_depth, missing_called_fraction, and missing_evidence. +## --no-qc-warnings suppresses the long per-position stdout list; warning totals +## are still printed, and all row-level warnings remain in the variants TSV. +## In the variants TSV, filter removed_by_gubbins=yes or qc_warnings containing +## gubbins_removed to inspect individual differences removed from clean.core.aln. From e7673d6906e11a8b6e80b297b821d014f3cca80f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 17:31:51 +0200 Subject: [PATCH 33/37] document Snippy SNP QC tools for PR 718 --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7d7518dc..2256c965 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Added Bacass numeric sample ID compatibility to characterization, Snippy, and PlasmidID templates [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) - Updated the Snippy IQ-TREE workflow to reuse the detected best-fit model [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) - Enhanced the outbreak summary workbook with normalized sample IDs, updated AMRFinderPlus parsing, SNP distance tables, provenance, and consistent formatting [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) +- Added reusable Snippy SNP-distance matrices and close-pair variant QC tooling [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) ### Modules From b648f53e83756cfc9413e0c31d279154e7ec2abe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sara=20Monz=C3=B3n?= Date: Wed, 15 Jul 2026 19:06:08 +0200 Subject: [PATCH 34/37] linting --- .../05-snippy/evaluate_close_pairs.py | 206 +++++++++++++----- .../generate_snp_distance_matrices.py | 51 +++-- 2 files changed, 192 insertions(+), 65 deletions(-) diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py index f8f54b22..ba4cdaf6 100644 --- a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/evaluate_close_pairs.py @@ -19,7 +19,6 @@ import tempfile from pathlib import Path - VALID_BASES = frozenset("ACGT") @@ -147,7 +146,9 @@ def parse_args() -> argparse.Namespace: return parser.parse_args() -def read_core_tab(path: Path) -> tuple[list[str], list[tuple[str, int, str, list[str]]]]: +def read_core_tab( + path: Path, +) -> tuple[list[str], list[tuple[str, int, str, list[str]]]]: with path.open(newline="") as handle: reader = csv.reader(handle, delimiter="\t") try: @@ -167,7 +168,9 @@ def read_core_tab(path: Path) -> tuple[list[str], list[tuple[str, int, str, list raise ValueError( f"Line {line_number} has {len(row)} columns; expected {len(header)}" ) - sites.append((row[0], int(row[1]), row[2].upper(), [x.upper() for x in row[3:]])) + sites.append( + (row[0], int(row[1]), row[2].upper(), [x.upper() for x in row[3:]]) + ) return samples, sites @@ -196,7 +199,9 @@ def find_close_pairs(samples, sites, threshold): "differences": differences, } ) - return sorted(results, key=lambda x: (x["snp_distance"], x["sample_1"], x["sample_2"])) + return sorted( + results, key=lambda x: (x["snp_distance"], x["sample_1"], x["sample_2"]) + ) def read_distance_matrix(path: Path) -> dict[frozenset[str], int]: @@ -322,7 +327,9 @@ def read_vcf(path: Path) -> dict[tuple[str, int, str], dict[str, str]]: "nearest_indel_distance": "", } if len(ref) == len(alt): - for offset, (ref_base, alt_base) in enumerate(zip(ref.upper(), alt.upper())): + for offset, (ref_base, alt_base) in enumerate( + zip(ref.upper(), alt.upper()) + ): if ref_base != alt_base and alt_base in VALID_BASES: index[(chrom, pos + offset, alt_base)] = evidence elif len(ref) == 1 and len(alt) == 1 and alt.upper() in VALID_BASES: @@ -370,7 +377,7 @@ def nearest_distance(sorted_positions: list[int], position: int) -> int | None: import bisect index = bisect.bisect_left(sorted_positions, position) - candidates = sorted_positions[max(0, index - 1):index + 1] + candidates = sorted_positions[max(0, index - 1) : index + 1] return min((abs(position - other) for other in candidates), default=None) @@ -420,9 +427,17 @@ def read_bam_pileup( bed.write(f"{chrom}\t{pos - 1}\t{pos}\n") bed.flush() command = [ - samtools, "mpileup", - "-q", str(minimum_mapq), "-Q", str(minimum_baseq), "-f", str(reference), - "-l", bed.name, str(bam), + samtools, + "mpileup", + "-q", + str(minimum_mapq), + "-Q", + str(minimum_baseq), + "-f", + str(reference), + "-l", + bed.name, + str(bam), ] result = subprocess.run(command, text=True, capture_output=True) if result.returncode: @@ -468,15 +483,31 @@ def evidence_for(index, bam_index, chrom, pos, core_ref, call, other_call): } ) return "reference_call_mpileup", evidence - status = "reference_call_no_bam_evidence" if call == core_ref else "variant_not_found_in_vcf" + status = ( + "reference_call_no_bam_evidence" + if call == core_ref + else "variant_not_found_in_vcf" + ) return status, {key: "" for key in EVIDENCE_COLUMNS} EVIDENCE_COLUMNS = ( - "vcf_record_pos", "vcf_ref", "vcf_alt", "qual", "filter", "genotype", - "depth", "ref_count", "alt_count", "variant_type", - "ref_forward_count", "ref_reverse_count", "alt_forward_count", - "alt_reverse_count", "homopolymer_run", "nearest_indel_distance", + "vcf_record_pos", + "vcf_ref", + "vcf_alt", + "qual", + "filter", + "genotype", + "depth", + "ref_count", + "alt_count", + "variant_type", + "ref_forward_count", + "ref_reverse_count", + "alt_forward_count", + "alt_reverse_count", + "homopolymer_run", + "nearest_indel_distance", ) @@ -530,7 +561,9 @@ def collect_qc_warnings( ("low_depth", f"{context} depth={depth:g} threshold={minimum_depth:g}") ) - fraction_name = "ref_fraction" if row[f"{prefix}_call"] == row["core_ref"] else "alt_fraction" + fraction_name = ( + "ref_fraction" if row[f"{prefix}_call"] == row["core_ref"] else "alt_fraction" + ) called_fraction = numeric(row[f"{prefix}_{fraction_name}"]) if called_fraction is None: warnings.append(("missing_called_fraction", context)) @@ -590,16 +623,22 @@ def main() -> int: if not 0 <= args.warn_min_strand_fraction <= 0.5: raise ValueError("--warn-min-strand-fraction must be between 0 and 0.5") if args.warn_indel_distance < 0 or args.warn_homopolymer_run < 1: - raise ValueError("indel distance must be non-negative and homopolymer run positive") + raise ValueError( + "indel distance must be non-negative and homopolymer run positive" + ) if args.cluster_window < 1 or args.cluster_min_snps < 2: - raise ValueError("cluster window must be positive and --cluster-min-snps at least 2") + raise ValueError( + "cluster window must be positive and --cluster-min-snps at least 2" + ) snippy_dir = args.snippy_dir.resolve() core_tab = (args.core_tab or snippy_dir / "core.tab").resolve() samples, sites = read_core_tab(core_tab) pairs = find_close_pairs(samples, sites, args.threshold) chromosomes = {chrom for chrom, *_ in sites} - gubbins_path = (args.gubbins_gff or snippy_dir / "gubbins.recombination_predictions.gff").resolve() + gubbins_path = ( + args.gubbins_gff or snippy_dir / "gubbins.recombination_predictions.gff" + ).resolve() gubbins_intervals = [] if gubbins_path.is_file(): if len(chromosomes) == 1: @@ -625,7 +664,9 @@ def main() -> int: if matrix_path.is_file(): matrices[column] = read_distance_matrix(matrix_path) else: - print(f"WARNING: comparison matrix not found: {matrix_path}", file=sys.stderr) + print( + f"WARNING: comparison matrix not found: {matrix_path}", file=sys.stderr + ) matrices[column] = {} for pair in pairs: key = frozenset((pair["sample_1"], pair["sample_2"])) @@ -634,19 +675,22 @@ def main() -> int: phylo_distance = pair["phylo_aln_distance"] clean_distance = pair["clean_core_aln_distance"] pair["core_tab_matrix_matches"] = ( - "yes" if pair["core_tab_matrix_distance"] == pair["snp_distance"] else - ("" if pair["core_tab_matrix_distance"] == "" else "no") + "yes" + if pair["core_tab_matrix_distance"] == pair["snp_distance"] + else ("" if pair["core_tab_matrix_distance"] == "" else "no") ) pair["phylo_matches_core_tab"] = ( - "yes" if phylo_distance == pair["snp_distance"] else - ("" if phylo_distance == "" else "no") + "yes" + if phylo_distance == pair["snp_distance"] + else ("" if phylo_distance == "" else "no") ) pair["snps_removed_by_gubbins"] = ( pair["snp_distance"] - clean_distance if clean_distance != "" else "" ) pair["clean_distance_valid"] = ( - "yes" if clean_distance != "" and clean_distance <= pair["snp_distance"] else - ("" if clean_distance == "" else "no") + "yes" + if clean_distance != "" and clean_distance <= pair["snp_distance"] + else ("" if clean_distance == "" else "no") ) comparison_warnings = [] if pair["core_tab_matrix_matches"] == "no": @@ -677,8 +721,11 @@ def main() -> int: if raw_path is not None: raw_index = read_vcf(raw_path) contextual_fields = ( - "ref_forward_count", "ref_reverse_count", "alt_forward_count", - "alt_reverse_count", "homopolymer_run", + "ref_forward_count", + "ref_reverse_count", + "alt_forward_count", + "alt_reverse_count", + "homopolymer_run", ) for key, evidence in vcf_indexes[sample].items(): raw_evidence = raw_index.get(key, {}) @@ -701,7 +748,11 @@ def main() -> int: bam_indexes = {sample: {} for sample in needed_samples} if not args.no_bam_evidence: - samtools = shutil.which(args.samtools) if Path(args.samtools).name == args.samtools else args.samtools + samtools = ( + shutil.which(args.samtools) + if Path(args.samtools).name == args.samtools + else args.samtools + ) if not samtools: raise FileNotFoundError( f"samtools executable not found: {args.samtools!r}; load samtools, use " @@ -728,11 +779,17 @@ def main() -> int: if not reference.is_file(): raise FileNotFoundError(f"Reference FASTA not found: {reference}") bam_indexes[sample] = read_bam_pileup( - str(samtools), reference, bam, positions, - args.pileup_mapq, args.pileup_baseq, + str(samtools), + reference, + bam, + positions, + args.pileup_mapq, + args.pileup_baseq, ) for (chrom, pos), pileup in bam_indexes[sample].items(): - distance = nearest_distance(non_snp_positions[sample].get(chrom, []), pos) + distance = nearest_distance( + non_snp_positions[sample].get(chrom, []), pos + ) pileup["nearest_indel_distance"] = "" if distance is None else distance prefix = args.output_prefix @@ -743,23 +800,46 @@ def main() -> int: variants_path = prefix.with_name(prefix.name + "_variants.tsv") pair_columns = [ - "pair_id", "sample_1", "sample_2", "snp_distance", - "comparable_core_sites", "ignored_non_acgt_sites", - "core_tab_matrix_distance", "core_tab_matrix_matches", - "phylo_aln_distance", "phylo_matches_core_tab", - "clean_core_aln_distance", "snps_removed_by_gubbins", - "clean_distance_valid", "distance_comparison_warnings", + "pair_id", + "sample_1", + "sample_2", + "snp_distance", + "comparable_core_sites", + "ignored_non_acgt_sites", + "core_tab_matrix_distance", + "core_tab_matrix_matches", + "phylo_aln_distance", + "phylo_matches_core_tab", + "clean_core_aln_distance", + "snps_removed_by_gubbins", + "clean_distance_valid", + "distance_comparison_warnings", ] with pairs_path.open("w", newline="") as handle: writer = csv.DictWriter(handle, fieldnames=pair_columns, delimiter="\t") writer.writeheader() for number, pair in enumerate(pairs, 1): - writer.writerow({"pair_id": f"pair_{number:04d}", **{k: pair[k] for k in pair_columns[1:]}}) + writer.writerow( + { + "pair_id": f"pair_{number:04d}", + **{k: pair[k] for k in pair_columns[1:]}, + } + ) base_columns = [ - "pair_id", "sample_1", "sample_2", "snp_distance", "chrom", "pos", - "core_ref", "sample_1_call", "sample_2_call", "pair_local_snp_count", - "pair_clustered_snp", "removed_by_gubbins", "qc_warnings", + "pair_id", + "sample_1", + "sample_2", + "snp_distance", + "chrom", + "pos", + "core_ref", + "sample_1_call", + "sample_2_call", + "pair_local_snp_count", + "pair_clustered_snp", + "removed_by_gubbins", + "qc_warnings", ] evidence_output = [] for n in (1, 2): @@ -770,7 +850,9 @@ def main() -> int: ) qc_warnings: dict[tuple[str, str, int, str], tuple[str, str]] = {} with variants_path.open("w", newline="") as handle: - writer = csv.DictWriter(handle, fieldnames=base_columns + evidence_output, delimiter="\t") + writer = csv.DictWriter( + handle, fieldnames=base_columns + evidence_output, delimiter="\t" + ) writer.writeheader() for number, pair in enumerate(pairs, 1): for chrom, pos, ref, call_1, call_2 in pair["differences"]: @@ -780,11 +862,18 @@ def main() -> int: ) row = { "pair_id": f"pair_{number:04d}", - "sample_1": pair["sample_1"], "sample_2": pair["sample_2"], - "snp_distance": pair["snp_distance"], "chrom": chrom, "pos": pos, - "core_ref": ref, "sample_1_call": call_1, "sample_2_call": call_2, + "sample_1": pair["sample_1"], + "sample_2": pair["sample_2"], + "snp_distance": pair["snp_distance"], + "chrom": chrom, + "pos": pos, + "core_ref": ref, + "sample_1_call": call_1, + "sample_2_call": call_2, "pair_local_snp_count": local_snp_count, - "pair_clustered_snp": "yes" if local_snp_count >= args.cluster_min_snps else "no", + "pair_clustered_snp": ( + "yes" if local_snp_count >= args.cluster_min_snps else "no" + ), "removed_by_gubbins": ( "yes" if position_in_intervals(pos, gubbins_intervals) else "no" ), @@ -804,15 +893,26 @@ def main() -> int: for n, call, other_call in ((1, call_1, call_2), (2, call_2, call_1)): sample = pair[f"sample_{n}"] status, evidence = evidence_for( - vcf_indexes[sample], bam_indexes[sample], chrom, pos, ref, - call, other_call, + vcf_indexes[sample], + bam_indexes[sample], + chrom, + pos, + ref, + call, + other_call, ) row[f"sample_{n}_evidence_status"] = status - row.update({f"sample_{n}_{key}": value for key, value in evidence.items()}) + row.update( + {f"sample_{n}_{key}": value for key, value in evidence.items()} + ) add_fractions(row, f"sample_{n}") sample_warnings = collect_qc_warnings( - row, n, args.warn_min_depth, args.warn_min_called_fraction, - args.warn_min_strand_fraction, args.warn_indel_distance, + row, + n, + args.warn_min_depth, + args.warn_min_called_fraction, + args.warn_min_strand_fraction, + args.warn_indel_distance, args.warn_homopolymer_run, ) for warning_type, message in sample_warnings: @@ -838,7 +938,9 @@ def main() -> int: for warning_type, _ in qc_warnings.values(): warning_counts[warning_type] = warning_counts.get(warning_type, 0) + 1 if warning_counts: - summary = ", ".join(f"{key}={value}" for key, value in sorted(warning_counts.items())) + summary = ", ".join( + f"{key}={value}" for key, value in sorted(warning_counts.items()) + ) print(f"Unique QC warnings: {summary}") else: print("Unique QC warnings: none") diff --git a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py index 0c2ab915..60fd7fc2 100644 --- a/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py +++ b/buisciii/templates/snippy/ANALYSIS/ANALYSIS01_SNIPPY/05-snippy/generate_snp_distance_matrices.py @@ -10,7 +10,6 @@ import numpy as np - VALID_BASES = np.frombuffer(b"ACGT", dtype=np.uint8) @@ -19,11 +18,15 @@ def parse_args() -> argparse.Namespace: description="Generate core.tab, phylo.aln, and clean.core.aln distance matrices." ) parser.add_argument( - "snippy_dir", nargs="?", type=Path, default=Path("."), + "snippy_dir", + nargs="?", + type=Path, + default=Path("."), help="Snippy output directory (default: current directory)", ) parser.add_argument( - "--output-dir", type=Path, + "--output-dir", + type=Path, help="matrix output directory (default: SNIPPY_DIR)", ) return parser.parse_args() @@ -47,7 +50,9 @@ def read_fasta_alignment(path: Path) -> tuple[list[str], list[str]]: raise ValueError(f"Empty FASTA identifier at {path}:{line_number}") current_sequence = [] elif current_name is None: - raise ValueError(f"Sequence before first FASTA header at {path}:{line_number}") + raise ValueError( + f"Sequence before first FASTA header at {path}:{line_number}" + ) else: current_sequence.append(line) if current_name is not None: @@ -85,23 +90,28 @@ def validate_sequences(path: Path, names: list[str], sequences: list[str]) -> No raise ValueError(f"Duplicate sequence identifiers in {path}") lengths = {len(sequence) for sequence in sequences} if len(lengths) != 1: - raise ValueError(f"Sequences have different lengths in {path}: {sorted(lengths)}") + raise ValueError( + f"Sequences have different lengths in {path}: {sorted(lengths)}" + ) def calculate_distances(sequences: list[str], chunk_size: int = 100_000) -> np.ndarray: """Count absolute differences with pairwise deletion of non-ACGT calls.""" sequence_array = np.vstack( - [np.frombuffer(sequence.encode("ascii"), dtype=np.uint8) for sequence in sequences] + [ + np.frombuffer(sequence.encode("ascii"), dtype=np.uint8) + for sequence in sequences + ] ) sequence_count, alignment_length = sequence_array.shape distances = np.zeros((sequence_count, sequence_count), dtype=np.int64) for start in range(0, alignment_length, chunk_size): - block = sequence_array[:, start:start + chunk_size] + block = sequence_array[:, start : start + chunk_size] valid = np.isin(block, VALID_BASES) for first in range(sequence_count - 1): - pair_valid = valid[first + 1:] & valid[first] - different = (block[first + 1:] != block[first]) & pair_valid - distances[first, first + 1:] += np.count_nonzero(different, axis=1) + pair_valid = valid[first + 1 :] & valid[first] + different = (block[first + 1 :] != block[first]) & pair_valid + distances[first, first + 1 :] += np.count_nonzero(different, axis=1) distances += distances.T return distances @@ -120,9 +130,24 @@ def main() -> int: output_dir = (args.output_dir or snippy_dir).resolve() output_dir.mkdir(parents=True, exist_ok=True) inputs = [ - ("core_tab", snippy_dir / "core.tab", "core_tab_snp_distances.tsv", read_core_tab), - ("phylo", snippy_dir / "phylo.aln", "phylo_snp_distances.tsv", read_fasta_alignment), - ("clean_core", snippy_dir / "clean.core.aln", "clean_core_snp_distances.tsv", read_fasta_alignment), + ( + "core_tab", + snippy_dir / "core.tab", + "core_tab_snp_distances.tsv", + read_core_tab, + ), + ( + "phylo", + snippy_dir / "phylo.aln", + "phylo_snp_distances.tsv", + read_fasta_alignment, + ), + ( + "clean_core", + snippy_dir / "clean.core.aln", + "clean_core_snp_distances.tsv", + read_fasta_alignment, + ), ] metadata = { "method": "Absolute A/C/G/T differences with pairwise deletion of other calls", From 6eb9bd6519a16be70912a185cb573d4eb022801b Mon Sep 17 00:00:00 2001 From: paupascualmas Date: Fri, 10 Jul 2026 13:46:27 +0200 Subject: [PATCH 35/37] Fix IRMA stats for flu and rsv when having missing readcounts --- .../04-irma/create_irma_stats_flu.sh | 128 +++++++++++++----- .../04-irma/create_irma_stats_rsv.sh | 46 +++++-- 2 files changed, 132 insertions(+), 42 deletions(-) diff --git a/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_flu.sh b/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_flu.sh index d36ea0f2..492a496a 100755 --- a/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_flu.sh +++ b/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_flu.sh @@ -1,43 +1,105 @@ +#!/usr/bin/env bash + +OUTPUT="irma_stats_flu.txt" + +echo -e "sample_ID\tTotalReads\tMappedReads\t%MappedReads\tFlu_type\tReads_HA\tReads_MP\tReads_NA\tReads_NP\tReads_NS\tReads_PA\tReads_PB1\tReads_PB2" > "$OUTPUT" + +while read -r in; do + [[ -z "$in" ]] && continue + + SAMPLE_ID="$in" + FILE="${in}/tables/READ_COUNTS.txt" + + if [[ -f "$FILE" ]]; then + read -r TOTAL_READS MAPPEDREADS PCTMAPPED < <( + awk -F'\t' ' + $1=="1-initial" {t=$2} + $1=="3-match" {m=$2} + END { + t += 0 + m += 0 + + if (t > 0) { + pct = (m / t) * 100 + } else { + pct = 0 + } + + printf "%d %d %.2f\n", t, m, pct + } + ' "$FILE" + ) + + FLU_TYPE=$(paste \ + <(grep '4-[A-C]' "$FILE" | cut -f1 | cut -d_ -f1 | cut -d- -f2 | head -n1) \ + <(grep '4-[A-B]_HA' "$FILE" | cut -f1 | cut -d_ -f3 | cut -d- -f2) \ + <(grep '4-[A-B]_NA' "$FILE" | cut -f1 | cut -d_ -f3) \ + | tr '\t' '_' \ + | sed 's/_*$//' + ) + + HA=$(grep -m1 '4-[A-C]_HA' "$FILE" | cut -f2) + MP=$(grep -m1 '4-[A-C]_MP' "$FILE" | cut -f2) + NA=$(grep -m1 '4-[A-C]_NA' "$FILE" | cut -f2) + NP=$(grep -m1 '4-[A-C]_NP' "$FILE" | cut -f2) + NS=$(grep -m1 '4-[A-C]_NS' "$FILE" | cut -f2) + PA=$(grep -m1 '4-[A-C]_PA' "$FILE" | cut -f2) + PB1=$(grep -m1 '4-[A-C]_PB1' "$FILE" | cut -f2) + PB2=$(grep -m1 '4-[A-C]_PB2' "$FILE" | cut -f2) + HE=$(grep -m1 '4-C_HE' "$FILE" | cut -f2) -echo -e "sample_ID\tTotalReads\tMappedReads\t%MappedReads\tFlu_type\tReads_HA\tReads_MP\tReads_NA\tReads_NP\tReads_NS\tReads_PA\tReads_PB1\tReads_PB2" > irma_stats_flu.txt - -cat ../samples_id.txt | while read in -do - SAMPLE_ID=$(echo ${in}) - TOTAL_READS=""; MAPPEDREADS=""; PCTMAPPED=""; FLU_TYPE=""; HA=""; MP=""; NA=""; NP=""; NS=""; PA=""; PB1=""; PB2=""; HE="" - - if test -f ${in}/tables/READ_COUNTS.txt; then - TOTAL_READS=$(grep '1-initial' ${in}/tables/READ_COUNTS.txt | cut -f2) - MAPPEDREADS=$(grep '3-match' ${in}/tables/READ_COUNTS.txt | cut -f2) - PCTMAPPED=$(awk "BEGIN {printf \"%.2f\", ($MAPPEDREADS/$TOTAL_READS)*100}") - FLU_TYPE=$(paste <(grep '4-[A-C]' ${in}/tables/READ_COUNTS.txt | cut -f1 | cut -d '_' -f1 | cut -d '-' -f2 | head -n1 ) <(grep '4-[A-B]_HA' ${in}/tables/READ_COUNTS.txt | cut -f1 | cut -d '_' -f3 | cut -d '-' -f2) <(grep '4-[A-B]_NA' ${in}/tables/READ_COUNTS.txt | cut -f1 | cut -d '_' -f3) | tr '\t' '_' | sed 's/_*$//') - HA=$(grep '4-[A-C]_HA' ${in}/tables/READ_COUNTS.txt | cut -f2) - MP=$(grep '4-[A-C]_MP' ${in}/tables/READ_COUNTS.txt | cut -f2) - NA=$(grep '4-[A-C]_NA' ${in}/tables/READ_COUNTS.txt | cut -f2) - NP=$(grep '4-[A-C]_NP' ${in}/tables/READ_COUNTS.txt | cut -f2) - NS=$(grep '4-[A-C]_NS' ${in}/tables/READ_COUNTS.txt | cut -f2) - PA=$(grep '4-[A-C]_PA' ${in}/tables/READ_COUNTS.txt | cut -f2) - PB1=$(grep '4-[A-C]_PB1' ${in}/tables/READ_COUNTS.txt | cut -f2) - PB2=$(grep '4-[A-C]_PB2' ${in}/tables/READ_COUNTS.txt | cut -f2) - #In case of Influenza C in samples: - HE=$(grep '4-C_HE' ${in}/tables/READ_COUNTS.txt | cut -f2) if [[ -n "$HE" ]]; then - LINE=$(paste <(echo $SAMPLE_ID) <(echo $TOTAL_READS) <(echo $MAPPEDREADS) <(echo $PCTMAPPED) <(echo $FLU_TYPE) <(echo $HA) <(echo $MP) <(echo $NA) <(echo $NP) <(echo $NS) <(echo $PA) <(echo $PB1) <(echo $PB2) <(echo $HE)) + printf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" \ + "$SAMPLE_ID" \ + "$TOTAL_READS" \ + "$MAPPEDREADS" \ + "$PCTMAPPED" \ + "$FLU_TYPE" \ + "$HA" \ + "$MP" \ + "$NA" \ + "$NP" \ + "$NS" \ + "$PA" \ + "$PB1" \ + "$PB2" \ + "$HE" \ + >> "$OUTPUT" else - LINE=$(paste <(echo $SAMPLE_ID) <(echo $TOTAL_READS) <(echo $MAPPEDREADS) <(echo $PCTMAPPED) <(echo $FLU_TYPE) <(echo $HA) <(echo $MP) <(echo $NA) <(echo $NP) <(echo $NS) <(echo $PA) <(echo $PB1) <(echo $PB2)) + printf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" \ + "$SAMPLE_ID" \ + "$TOTAL_READS" \ + "$MAPPEDREADS" \ + "$PCTMAPPED" \ + "$FLU_TYPE" \ + "$HA" \ + "$MP" \ + "$NA" \ + "$NP" \ + "$NS" \ + "$PA" \ + "$PB1" \ + "$PB2" \ + >> "$OUTPUT" fi + else echo "Sample $SAMPLE_ID doesn't have READ_COUNTS.txt file. Skipping" - TOTAL_READS=NA - MAPPEDREADS=0 - PCTMAPPED=0 - LINE=$(paste <(echo $SAMPLE_ID) <(echo $TOTAL_READS) <(echo $MAPPEDREADS) <(echo $PCTMAPPED)) + + printf "%s\tNA\t0\t0.00\t\t\t\t\t\t\t\t\t\n" \ + "$SAMPLE_ID" \ + >> "$OUTPUT" fi - echo "$LINE" >> irma_stats_flu.txt -done +done < ../samples_id.txt -ANY_C=$(grep "C" irma_stats_flu.txt) -if [[ -n "$ANY_C" ]]; then - sed -i 's/Reads_PB2/Reads_PB2\tReads_HE/g' irma_stats_flu.txt +if awk -F'\t' ' + NR > 1 && NF >= 14 && $14 != "" { + found = 1 + } + END { + exit(found ? 0 : 1) + } +' "$OUTPUT"; then + sed -i '1s/Reads_PB2$/Reads_PB2\tReads_HE/' "$OUTPUT" fi \ No newline at end of file diff --git a/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_rsv.sh b/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_rsv.sh index 781e0539..8cb40ace 100644 --- a/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_rsv.sh +++ b/buisciii/templates/IRMA/ANALYSIS/ANALYSIS01_IRMA/04-irma/create_irma_stats_rsv.sh @@ -1,11 +1,39 @@ +#!/usr/bin/env bash + echo -e "Sample_ID\tTotalReads\tMappedReads\t%MappedReads\tRSV_type" > irma_stats_rsv.txt -cat ../samples_id.txt | while read in -do -SAMPLE_ID=$(echo ${in}) -TOTAL_READS=$(grep '1-initial' ${in}/tables/READ_COUNTS.txt | cut -f2) -MAPPEDREADS=$(grep '3-match' ${in}/tables/READ_COUNTS.txt | cut -f2) -PCTMAPPED=$(awk "BEGIN {printf \"%.2f\", ($MAPPEDREADS/$TOTAL_READS)*100}") -RSV_TYPE=$(grep '4-RSV_' ${in}/tables/READ_COUNTS.txt | cut -f1 | cut -d '_' -f2) -echo -e "${SAMPLE_ID}\t${TOTAL_READS}\t${MAPPEDREADS}\t${PCTMAPPED}\t${RSV_TYPE}" >> irma_stats_rsv.txt -done +while read -r in; do + SAMPLE_ID="$in" + FILE="${in}/tables/READ_COUNTS.txt" + + if [[ -f "$FILE" ]]; then + read -r TOTAL_READS MAPPEDREADS PCTMAPPED < <( + awk -F'\t' ' + $1=="1-initial" {t=$2} + $1=="3-match" {m=$2} + END { + t+=0 + m+=0 + printf "%d %d %.2f\n", t, m, t>0 ? (m/t)*100 : 0 + } + ' "$FILE" + ) + + RSV_TYPE=$(grep '4-RSV_' "$FILE" | cut -f1 | cut -d_ -f2) + else + echo "Sample $SAMPLE_ID doesn't have READ_COUNTS.txt file. Skipping" + TOTAL_READS="NA" + MAPPEDREADS=0 + PCTMAPPED="0.00" + RSV_TYPE="" + fi + + printf "%s\t%s\t%s\t%s\t%s\n" \ + "$SAMPLE_ID" \ + "$TOTAL_READS" \ + "$MAPPEDREADS" \ + "$PCTMAPPED" \ + "$RSV_TYPE" \ + >> irma_stats_rsv.txt + +done < ../samples_id.txt \ No newline at end of file From 0e67670749bf5a985a71299e831212e6e8eaf969 Mon Sep 17 00:00:00 2001 From: paupascualmas Date: Fri, 10 Jul 2026 13:53:43 +0200 Subject: [PATCH 36/37] Updated CHANGELOG.md --- CHANGELOG.md | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2256c965..02566678 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,31 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [2.3.5]_dev - XXXX-XX-XX : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.5_dev + +### Credits + +- [Pau Pascual](https://github.com/PauPascualMas) + +### Template fixes and updates + +- Fix IRMA stats generation when 3-match is absent [#719](https://github.com/BU-ISCIII/buisciii-tools/pull/719) + +### Modules + +#### Implementation + +#### Added enhancements + +#### Fixes + +#### Changed + +#### Removed + +### Requirements + + ## [2.3.4] - 2026-07-02 : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.4 ### Credits From fba4f7257c4be5e28569ea9213390020c2ae9782 Mon Sep 17 00:00:00 2001 From: paupascualmas Date: Mon, 20 Jul 2026 13:31:27 +0200 Subject: [PATCH 37/37] Updated CHANGELOG.md --- CHANGELOG.md | 26 ++------------------------ 1 file changed, 2 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 02566678..1de04a14 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,30 +4,6 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [2.3.5]_dev - XXXX-XX-XX : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.5_dev - -### Credits - -- [Pau Pascual](https://github.com/PauPascualMas) - -### Template fixes and updates - -- Fix IRMA stats generation when 3-match is absent [#719](https://github.com/BU-ISCIII/buisciii-tools/pull/719) - -### Modules - -#### Implementation - -#### Added enhancements - -#### Fixes - -#### Changed - -#### Removed - -### Requirements - ## [2.3.4] - 2026-07-02 : https://github.com/BU-ISCIII/buisciii-tools/releases/tag/2.3.4 @@ -36,6 +12,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - [Enrique Sapena](https://github.com/ESapenaVentura) - [Victor Lopez](https://github.com/victor5lm) - [Sara Monzon](https://github.com/saramonzon) +- [Pau Pascual](https://github.com/PauPascualMas) ### Template fixes and updates @@ -48,6 +25,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Updated the Snippy IQ-TREE workflow to reuse the detected best-fit model [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) - Enhanced the outbreak summary workbook with normalized sample IDs, updated AMRFinderPlus parsing, SNP distance tables, provenance, and consistent formatting [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) - Added reusable Snippy SNP-distance matrices and close-pair variant QC tooling [#718](https://github.com/BU-ISCIII/buisciii-tools/pull/718) +- Fix IRMA stats generation when 3-match is absent [#719](https://github.com/BU-ISCIII/buisciii-tools/pull/719) ### Modules