From 558a0577658ba1eef65e3140e7cb1230dafe0b05 Mon Sep 17 00:00:00 2001 From: KatharinaHoff Date: Sun, 16 Nov 2025 15:02:42 +0100 Subject: [PATCH 1/2] modularized --- PRISTINE/PRISTINE_essentials.py | 100 ++++----------- PRISTINE/config.py | 120 ++++++++++++++++++ .../Yersinia_pseudotuberculosis_1.fasta | 2 +- .../Yersinia_pseudotuberculosis_2.fasta | 2 +- .../non-targets/Yersinia_pestis_1.fasta | 2 +- .../non-targets/Yersinia_pestis_2.fasta | 2 +- README.md | 55 +++++++- 7 files changed, 197 insertions(+), 86 deletions(-) create mode 100644 PRISTINE/config.py diff --git a/PRISTINE/PRISTINE_essentials.py b/PRISTINE/PRISTINE_essentials.py index 74ad8b3..8e6d1e3 100644 --- a/PRISTINE/PRISTINE_essentials.py +++ b/PRISTINE/PRISTINE_essentials.py @@ -25,85 +25,17 @@ import re import seaborn as sns import matplotlib.pyplot as plt -from collections import Counter -from dataclasses import dataclass, field -from typing import Optional, Dict -import time -import yaml from pathlib import Path - -@dataclass -class InputPaths: - raw_dir: Optional[str] = None - prokka_dir: Optional[str] = None - panaroo_dir: Optional[str] = None - - -@dataclass -class Primer3Params: - global_params: Dict[str, object] = field(default_factory=dict) - design_params: Dict[str, object] = field(default_factory=dict) - - -@dataclass -class SNPPrimerDesignParams: - snp_window_size: int - snp_top_n: int - min_snps: int - - -@dataclass -class ValidationConfig: - perform: str - database: str - pblat_min_identity: float - match_median_filter_tolerance: int - - -@dataclass -class Config: - input_type: str - input_paths: InputPaths - output_dir: str - max_cores: int - aligner: str - snp_avg_prop_threshold: float - primer3_config_file: Optional[str] = None - primer3: Primer3Params = field(default_factory=Primer3Params) - snp_primer_design: SNPPrimerDesignParams = field(default_factory=SNPPrimerDesignParams) - validation: ValidationConfig = field(default_factory=ValidationConfig) - - def __post_init__(self): - allowed_inputs = {"raw", "prokka", "panaroo"} - if self.input_type not in allowed_inputs: - raise ValueError(f"`input_type` must be one of {allowed_inputs}, got: {self.input_type}") - - active_dir = { - "raw": self.input_paths.raw_dir, - "prokka": self.input_paths.prokka_dir, - "panaroo": self.input_paths.panaroo_dir, - }[self.input_type] - if not active_dir: - raise ValueError(f"{self.input_type}_dir must be set in `input_paths` for input_type = {self.input_type}") - - if self.primer3_config_file and self.primer3.global_params: - print("Warning: primer3_config_file is set. Inline primer3 parameters will be ignored.") - - -class ConfigLoader: - @staticmethod - def load(path: str) -> Config: - with open(path, "r") as f: - raw = yaml.safe_load(f) - - # Parse nested dataclasses manually - raw["input_paths"] = InputPaths(**raw.get("input_paths", {})) - raw["primer3"] = Primer3Params(**raw.get("primer3", {})) - raw["snp_primer_design"] = SNPPrimerDesignParams(**raw.get("snp_primer_design", {})) - raw["validation"] = ValidationConfig(**raw.get("validation", {})) - - return Config(**raw) +# Import configuration classes from config module +from config import ( + Config, + ConfigLoader, + InputPaths, + Primer3Params, + SNPPrimerDesignParams, + ValidationConfig, +) class MSAStrategy: @@ -2319,7 +2251,19 @@ def design_primers(self, input_folder=".", output_folder="output/primers", start if __name__ == "__main__": - config = ConfigLoader.load("config.yaml") + parser = argparse.ArgumentParser( + description="PRISTINE: PRimer-based Identification Suite Targeting Informative Nucleotide Elements" + ) + parser.add_argument( + "--config", + type=str, + default="config.yaml", + help="Path to the configuration YAML file (default: config.yaml in current directory)" + ) + args = parser.parse_args() + + # Load configuration from the specified path + config = ConfigLoader.load(args.config) runner = MSAStrategy(config) runner.run() diff --git a/PRISTINE/config.py b/PRISTINE/config.py new file mode 100644 index 0000000..63e5dd7 --- /dev/null +++ b/PRISTINE/config.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +""" +Configuration classes for PRISTINE pipeline. + +This module contains all configuration-related dataclasses and the ConfigLoader +for parsing YAML configuration files. +""" + +from __future__ import annotations +from dataclasses import dataclass, field +from typing import Optional, Dict +import yaml +import os +import sys + + +@dataclass +class InputPaths: + """Container for input directory paths.""" + raw_dir: Optional[str] = None + prokka_dir: Optional[str] = None + panaroo_dir: Optional[str] = None + + +@dataclass +class Primer3Params: + """Container for Primer3 global and design parameters.""" + global_params: Dict[str, object] = field(default_factory=dict) + design_params: Dict[str, object] = field(default_factory=dict) + + +@dataclass +class SNPPrimerDesignParams: + """Configuration for SNP-aware primer design strategy.""" + snp_window_size: int + snp_top_n: int + min_snps: int + + +@dataclass +class ValidationConfig: + """Configuration for post-design validation using pBLAT.""" + perform: str + database: str + pblat_min_identity: float + match_median_filter_tolerance: int + + +@dataclass +class Config: + """ + Main configuration class for PRISTINE pipeline. + + Validates input_type and ensures required paths are set based on the + selected input type (raw/prokka/panaroo). + """ + input_type: str + input_paths: InputPaths + output_dir: str + max_cores: int + aligner: str + snp_avg_prop_threshold: float + primer3_config_file: Optional[str] = None + primer3: Primer3Params = field(default_factory=Primer3Params) + snp_primer_design: SNPPrimerDesignParams = field(default_factory=SNPPrimerDesignParams) + validation: ValidationConfig = field(default_factory=ValidationConfig) + + def __post_init__(self): + """Validate configuration after initialization.""" + allowed_inputs = {"raw", "prokka", "panaroo"} + if self.input_type not in allowed_inputs: + raise ValueError(f"`input_type` must be one of {allowed_inputs}, got: {self.input_type}") + + active_dir = { + "raw": self.input_paths.raw_dir, + "prokka": self.input_paths.prokka_dir, + "panaroo": self.input_paths.panaroo_dir, + }[self.input_type] + if not active_dir: + raise ValueError(f"{self.input_type}_dir must be set in `input_paths` for input_type = {self.input_type}") + + if self.primer3_config_file and self.primer3.global_params: + print("Warning: primer3_config_file is set. Inline primer3 parameters will be ignored.") + + +class ConfigLoader: + """Utility class for loading YAML configuration files.""" + + @staticmethod + def load(path: str) -> Config: + """ + Load and parse a YAML configuration file. + + Args: + path: Path to the YAML configuration file + + Returns: + Config object with all nested dataclasses instantiated + + Raises: + FileNotFoundError: If the configuration file is not found + SystemExit: If the configuration file is not found (with helpful error message) + """ + if not os.path.exists(path): + print(f"Error: Configuration file not found: {path}", file=sys.stderr) + print(f"Current working directory: {os.getcwd()}", file=sys.stderr) + print(f"\nPlease ensure 'config.yaml' is in the same directory where you run the container.", file=sys.stderr) + print(f"Or specify a custom path using: --config /path/to/config.yaml", file=sys.stderr) + sys.exit(1) + + with open(path, "r") as f: + raw = yaml.safe_load(f) + + # Parse nested dataclasses manually + raw["input_paths"] = InputPaths(**raw.get("input_paths", {})) + raw["primer3"] = Primer3Params(**raw.get("primer3", {})) + raw["snp_primer_design"] = SNPPrimerDesignParams(**raw.get("snp_primer_design", {})) + raw["validation"] = ValidationConfig(**raw.get("validation", {})) + + return Config(**raw) diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta index 71a8546..579e85f 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta @@ -1,4 +1,4 @@ ->Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|CP130901 Yersinia pseudotuberculosis strain SP-1303 chromosome, complete genome. [Yersinia pseudotuberculosis SP-1303 | 633.777] +>Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|CP130901 Yersinia pseudotuberculosis strain SP-1303 chromosome, complete genome. [Yersinia pseudotuberculosis SP-1303 | 633.777] gtgtcactttcgctttggcagcagtgtcttgcccgattgcaggatgagttacctgccaca gaatttagtatgtggatacgccccctacaggcggaactgagtgacaatactctggcgctt tacgcacctaatcgttttgtactggactgggtccgtgataagtacttaaacaatatcaat diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta index 85bb684..28a3fdc 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta @@ -1,4 +1,4 @@ ->Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|LT596221 Yersinia pseudotuberculosis isolate NZYP4713 genome assembly, chromosome: 1. [Yersinia pseudotuberculosis strain NZYP4713 | 633.88] +>Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|LT596221 Yersinia pseudotuberculosis isolate NZYP4713 genome assembly, chromosome: 1. [Yersinia pseudotuberculosis strain NZYP4713 | 633.88] gtgtcactttcgctttggcagcagtgtcttgcccgattgcaggatgagttacctgccaca gaatttagtatgtggatacgccccctacaggcggaactgagtgacaatactctggcgctt tacgcacctaatcgttttgtactggactgggtccgtgataagtacttaaacaatatcaat diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta index aef660d..e000d9b 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta @@ -1,4 +1,4 @@ ->Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177922 Yersinia pestis strain 50008 chromosome, complete genome. [Yersinia pestis 50008 | 632.1142] +>Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177922 Yersinia pestis strain 50008 chromosome, complete genome. [Yersinia pestis 50008 | 632.1142] ggcggactccactcgaacaagatcaataaaagaaagaaatcaggatattttgtaactcat cgttttgacacacgccatatggcctgtagcaagccaaggtgacatgcgaactgtttttta cagttgttttttactgctgattttttacgactcacattcccccacgatcccagagtagag diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta index 1cb30a5..daf6c91 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta @@ -1,4 +1,4 @@ ->Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177934 Yersinia pestis strain 96005 chromosome, complete genome. [Yersinia pestis 96005 | 632.1143] +>Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177934 Yersinia pestis strain 96005 chromosome, complete genome. [Yersinia pestis 96005 | 632.1143] ggcggactccactcgaacaagatcaataaaagaaagaaatcaggatattttgtaactcat cgttttgacacacgccatatggcctgtagcaagccaaggtgacatgcgaactgtttttta cagttgttttttactgctgattttttacgactcacattcccccacgatcccagagtagag diff --git a/README.md b/README.md index 4069679..c4839cc 100644 --- a/README.md +++ b/README.md @@ -6,9 +6,10 @@ This pipeline is distributed as a self-contained [Apptainer](https://apptainer.o To use the pipeline, you will need to have **Apptainer** installed on your system. Installation instructions are available on the [Apptainer documentation site](https://apptainer.org/docs/). -Once Apptainer is installed, no further setup is required. Simply download the pipeline container available at https://nextcloud.uni-greifswald.de/index.php/s/GGWMXpJEYkQpCxG and run the pipeline using: +Once Apptainer is installed, no further setup is required. Simply download the pipeline container available at https://nextcloud.uni-greifswald.de/index.php/s/GGWMXpJEYkQpCxG, make it executable, and run the pipeline using: ```bash +chmod +x PRISTINE.sif ./PRISTINE.sif ``` @@ -23,7 +24,12 @@ To run the pipeline, you need: ### Configuration File -The file `config.yaml` must be placed in the **same directory as the container** (`PRISTINE.sif`). It defines all global parameters, paths, and settings required for the analysis (for details, see below). +The file `config.yaml` must be placed in the **same directory where you run the container** (`PRISTINE.sif`). It defines all global parameters, paths, and settings required for the analysis (for details, see below). + +**Note:** By default, the pipeline looks for `config.yaml` in the current working directory. You can specify a different location using: +```bash +./PRISTINE.sif --config /path/to/your/config.yaml +``` ## Quick Start with Toy Dataset @@ -36,11 +42,52 @@ To quickly test the pipeline: 1. **Download the `toy_dataset/` directory** 2. **Place the Apptainer container** (`PRISTINE.sif`) inside the `toy_dataset/` folder -3. From within the `toy_dataset/` directory, run: +3. **Make the container executable** (if not already): + ```bash + chmod +x PRISTINE.sif + ``` +4. From within the `toy_dataset/` directory, run: + ```bash + ./PRISTINE.sif + ``` + +--- + +## Development: Running with Modified Source Code + +If you're developing or testing modifications to the PRISTINE source code, you can run the container with bind mounts to override the built-in Python modules without rebuilding the entire container. + +### Option 1: Using the Wrapper Script + +A convenience wrapper script is provided that automatically binds the updated modules: ```bash -./PRISTINE.sif +cd PRISTINE/toy_dataset +/path/to/run_pristine_updated.sh ``` + +The wrapper script supports all the same arguments as the main pipeline: +```bash +/path/to/run_pristine_updated.sh --config /path/to/config.yaml +``` + +### Option 2: Direct Container Execution with Bind Mounts + +For manual control, you can use Singularity/Apptainer bind mounts directly: + +```bash +cd PRISTINE/toy_dataset +singularity exec \ + --bind /path/to/PRISTINE/PRISTINE_essentials.py:/opt/PRISTINE_essentials.py \ + --bind /path/to/PRISTINE/config.py:/opt/config.py \ + /path/to/PRISTINE.sif \ + bash -c 'source /opt/conda/etc/profile.d/conda.sh && conda activate primer_design_tool_full && python3 /opt/PRISTINE_essentials.py' +``` + +Replace `/path/to/` with the actual paths to your modified source files and container. + +**Note:** This approach is intended for development and testing only. For production use, rebuild the container with the updated source code. + --- ## Pipeline overview From da1800fd0311a28431cde97c0a2da0c89c84303d Mon Sep 17 00:00:00 2001 From: KatharinaHoff Date: Sun, 16 Nov 2025 16:12:26 +0100 Subject: [PATCH 2/2] further modularization --- PRISTINE/PRISTINE_essentials.py | 334 ++---------------- PRISTINE/data_structures.py | 155 ++++++++ PRISTINE/similarity_metrics.py | 124 +++++++ .../Yersinia_pseudotuberculosis_1.fasta | 2 +- .../Yersinia_pseudotuberculosis_2.fasta | 2 +- .../non-targets/Yersinia_pestis_1.fasta | 2 +- .../non-targets/Yersinia_pestis_2.fasta | 2 +- PRISTINE/visualization.py | 94 +++++ 8 files changed, 403 insertions(+), 312 deletions(-) create mode 100644 PRISTINE/data_structures.py create mode 100644 PRISTINE/similarity_metrics.py create mode 100644 PRISTINE/visualization.py diff --git a/PRISTINE/PRISTINE_essentials.py b/PRISTINE/PRISTINE_essentials.py index 8e6d1e3..1186574 100644 --- a/PRISTINE/PRISTINE_essentials.py +++ b/PRISTINE/PRISTINE_essentials.py @@ -37,6 +37,24 @@ ValidationConfig, ) +# Import similarity metrics +from similarity_metrics import ( + manhattan_distance, + cosine_similarity, + jaccard_similarity_with_frequencies, + calculate_weighted_minhash_similarity, + add_missing_keys_with_zero, +) + +# Import data structures +from data_structures import Segment, QuasiAlignment + +# Import visualization functions +from visualization import ( + plot_informativeness_heatmap, + plot_snp_density_distribution, +) + class MSAStrategy: # --- Internal Methods --- @@ -376,80 +394,20 @@ def assess_loci_by_species(self, species_folder): return df def plot_informativeness_heatmap(self, species_folder): + """Generate a heatmap showing informative SNP proportions.""" output_dir = os.path.join(self.output_dir, species_folder, "informative_loci") csv_file = os.path.join(output_dir, "snp_summary.csv") output_file = os.path.join(output_dir, "heatmap_informativeness.png") - df = pd.read_csv(csv_file) - - prop_cols = [col for col in df.columns if col.startswith("Prop_") and col != "Avg_Prop_Informative_SNPs"] - if not prop_cols: - print("No proportion columns found in the dataset.") - return - - heatmap_data = df.set_index("Locus")[prop_cols] - - if "Avg_Prop_Informative_SNPs" in df.columns: - sorted_loci = df.sort_values("Avg_Prop_Informative_SNPs", ascending=False)["Locus"] - heatmap_data = heatmap_data.loc[sorted_loci] - - plt.figure(figsize=(12, max(6, 0.3 * len(heatmap_data)))) - sns.heatmap(heatmap_data, annot=False, cmap="YlOrRd", cbar_kws={'label': 'Proportion of Informative SNPs'}) - plt.title("Informative SNP Proportions per Locus and Non-Target Species") - plt.xlabel("Non-Target Species") - plt.ylabel("Locus") - plt.tight_layout() - plt.savefig(output_file) - plt.close() - - return output_file + return plot_informativeness_heatmap(csv_file, output_file) def plot_snp_density_lines(self, species_folder, top_n=5): - import os - import matplotlib.pyplot as plt - import seaborn as sns - import pandas as pd - + """Generate SNP position distribution plots for top-ranked loci.""" output_dir = os.path.join(self.output_dir, species_folder, "informative_loci") csv_file = os.path.join(output_dir, "snp_summary.csv") plot_dir = os.path.join(output_dir, "snp_density_plots") - os.makedirs(plot_dir, exist_ok=True) - - df = pd.read_csv(csv_file) - - if "Avg_Prop_Informative_SNPs" not in df.columns: - print("Missing ranking column 'Avg_Prop_Informative_SNPs'.") - return - - top_loci = df.sort_values("Avg_Prop_Informative_SNPs", ascending=False).head(top_n) - - prop_cols = [col for col in df.columns if col.startswith("SNP_Pos_")] - - for _, row in top_loci.iterrows(): - locus_name = row["Locus"] - plt.figure(figsize=(12, 3)) - - for col in prop_cols: - species = col.replace("SNP_Pos_", "") - if pd.isna(row[col]) or not row[col].strip(): - continue - try: - positions = list(map(int, row[col].split(','))) - except ValueError: - continue - sns.histplot(positions, bins=50, kde=False, label=species, element="step", fill=False) - - plt.title(f"SNP Position Distribution - {locus_name}") - plt.xlabel("Alignment Position") - plt.ylabel("SNP Count") - plt.legend(title="Non-Target Species", loc="upper right", fontsize="small") - plt.tight_layout() - plot_path = os.path.join(plot_dir, f"snp_density_{locus_name.replace('.','_')}.png") - plt.savefig(plot_path) - plt.close() - print(f"SNP density plots saved to: {plot_dir}") - return plot_dir + return plot_snp_density_distribution(csv_file, plot_dir, top_n) def create_consensus_sequences(self, species_folder): IUPAC_CODES = { @@ -1113,98 +1071,9 @@ def parse_species_csv(self, file_path): species_dict[target_species] = non_target_species return species_dict - - @staticmethod - def manhattan_distance(point1, point2): - """ - Calculate the Manhattan distance between two points. - - :param point1: A list or tuple representing the coordinates of the first point. - :param point2: A list or tuple representing the coordinates of the second point. - :return: The Manhattan distance between the two points. - """ - if len(point1) != len(point2): - raise ValueError("Both points must have the same number of dimensions") - - distance = sum(abs(p1 - p2) for p1, p2 in zip(point1, point2)) - return distance - - @staticmethod - def cosine_similarity(vector_a, vector_b): - """ - Calculates the cosine similarity between two vectors. - :param vector_a: A list or tuple representing the first vector. - :param vector_b: A list or tuple representing the second vector. - :return: The cosine similarity between vector_a and vector_b. - """ - # Ensure inputs are lists or tuples - if not isinstance(vector_a, (list, tuple)) or not isinstance(vector_b, (list, tuple)): - raise TypeError("Input vectors must be lists or tuples.") - - # Convert to NumPy arrays - vector_a = np.array(vector_a) - vector_b = np.array(vector_b) - - # Calculate dot product and magnitudes - dot_product = np.dot(vector_a, vector_b) - magnitude_a = np.linalg.norm(vector_a) - magnitude_b = np.linalg.norm(vector_b) - - # Avoid division by zero - if magnitude_a == 0 or magnitude_b == 0: - return 0 # Return 0 for cosine similarity if either vector is zero - - # Compute cosine similarity - return dot_product / (magnitude_a * magnitude_b) - - @staticmethod - def jaccard_similarity_with_frequencies(vector_a, vector_b): - """ - Calculates the Jaccard similarity between two lists or tuples containing frequencies of individual items (e.g., p-mers). + # Similarity metrics moved to similarity_metrics.py module - :param list_a: A list or tuple representing the first vector. - :param list_b: A list or tuple representing the second vector. - :return: The Jaccard similarity between the two collections, accounting for frequencies. - """ - # Convert lists to Counters to handle frequencies - counter_a = Counter(vector_a) - counter_b = Counter(vector_b) - - # Calculate intersection and union based on minimum and maximum frequencies - intersection = sum((min(counter_a[item], counter_b[item]) for item in counter_a if item in counter_b)) - union = sum((max(counter_a[item], counter_b[item]) for item in set(counter_a) | set(counter_b))) - - # Avoid division by zero in case both counters are empty - if union == 0: - return 0 - - # Compute Jaccard similarity with frequencies - return intersection / union - - @staticmethod - def calculate_weighted_minhash_similarity(frequency_vector_a, frequency_vector_b, num_perm=128): - """ - Calculates the Weighted MinHash Jaccard similarity between two frequency vectors. - - :param frequency_vector_a: A list or tuple representing the frequency of p-mers for the first vector. - :param frequency_vector_b: A list or tuple representing the frequency of p-mers for the second vector. - :param num_perm: The number of permutations (hash functions) for the MinHash. - :return: The Weighted MinHash Jaccard similarity between the two vectors. - """ - if not isinstance(frequency_vector_a, (list, tuple)) or not isinstance(frequency_vector_b, (list, tuple)): - raise TypeError("Input vectors must be lists or tuples.") - - if len(frequency_vector_a) != len(frequency_vector_b): - raise ValueError("The two frequency vectors must be of the same length.") - - wmg = WeightedMinHashGenerator(len(frequency_vector_a), sample_size=num_perm) - - wm_a = wmg.minhash(frequency_vector_a) - wm_b = wmg.minhash(frequency_vector_b) - - return wm_a.jaccard(wm_b) - def run_prokka(self, input_dir=".", output_dir="output/prokka"): """ Runs Prokka on all FASTA files in the specified input directory, saving results in a single output directory. @@ -1642,27 +1511,8 @@ def get_iupac_code(self, base_counts): frozenset(['A', 'C', 'G']): 'V', frozenset(['A', 'T', 'C', 'G']): 'N' } return iupac_dict[frozenset(base_counts.keys())] - - @staticmethod - def add_missing_keys_with_zero(dict1, dict2): - """ - Adds keys that are present in dict1 but not in dict2 (and vice versa) to each dictionary with values set to zero. - - :param dict1: The first dictionary. - :param dict2: The second dictionary. - :return: Two dictionaries with added keys and values set to zero where keys were missing. - """ - # Find keys unique to each dictionary - keys_in_dict1_not_in_dict2 = dict1.keys() - dict2.keys() - keys_in_dict2_not_in_dict1 = dict2.keys() - dict1.keys() - - # Add missing keys to each dictionary with values set to zero - for key in keys_in_dict1_not_in_dict2: - dict2[key] = 0 - for key in keys_in_dict2_not_in_dict1: - dict1[key] = 0 - - return dict1, dict2 + + # add_missing_keys_with_zero moved to similarity_metrics.py module def find_consensus_binding_regions(self, all_primers, window_size=10, threshold_proportion=0.5): """ @@ -2075,141 +1925,9 @@ def design_primers_from_csv(self, input_csv, output_csv): os.remove(primer3_output_file) # --- Internal Classes --- - - class Segment: - """ - A class to represent an individual segment in a quasi-alignment. It includes information about the sequence ID, species, - and the position of the segment within the sequence. - """ - def __init__(self, seq_id: str, species_name: str, gene_name: str, start: int, length: int, sequence: str, p: int = 3): - """ - Initialize a Segment object. - - :param seq_id: The ID of the sequence from which the segment comes. - :param species_name: The species of the sequence. - :param gene_name: The gene region of the sequence. - :param start: The starting position of the segment within the sequence. - :param length: The length of the segment within the sequence. - """ - self.seq_id = seq_id - self.species_name = species_name - self.gene_name = gene_name - self.start = start # 0-indexed - self.length = length - self.sequence = sequence - self.pmer_profile = self.get_pmer_composition(self, p) - - def get_segment(self, sequence): - """Returns the segment of the sequence based on start and length.""" - end = self.start + self.length - return sequence[self.start:end] - - def get_pmer_composition(self, segment, p): - """ - Takes a Segment object and computes the p-mer composition of the segment. - The output is a dictionary where keys are only the p-mers present in the segment's sequence, and values are their counts. - P-mers are stored in uppercase. - - :param segment: A Segment object containing the start and length information. - :param p: The length of the p-mer (substring). - :return: A dictionary with present p-mers as keys and their counts as values. - """ - sequence = segment.sequence.upper() # Convert the entire sequence to uppercase - pmer_dict = {} - - # Loop through the sequence and extract p-mers - for i in range(len(sequence) - p + 1): # Ensure we don't go beyond the sequence - pmer = sequence[i:i + p] - - # Add or update the count of the p-mer in the dictionary - if pmer in pmer_dict: - pmer_dict[pmer] += 1 - else: - pmer_dict[pmer] = 1 - - return pmer_dict - - def get_distance(self, other_segment, distance_func): - """ - Calculates the distance between this segment and another segment using a specified distance function. - - :param other_segment: The other Segment object to calculate the distance to. - :param distance_func: A function that takes two lists (or tuples) as input and returns the distance between them. - :return: The distance between the two p-mer profiles. - """ - if not callable(distance_func): - raise ValueError("distance_func must be a callable function") - - aligned_profile_self, aligned_profile_other = QuasiAlignmentStrategy.add_missing_keys_with_zero(self.pmer_profile, other_segment.pmer_profile) - - profile_values_self = tuple(aligned_profile_self.values()) - profile_values_other = tuple(aligned_profile_other.values()) - - # Calculate and return the distance using the specified distance function - return distance_func(profile_values_self, profile_values_other) + # Segment and QuasiAlignment classes moved to data_structures.py module - def __repr__(self): - return f"Segment(seq_id={self.seq_id}, species={self.species}, start={self.start}, end={self.end})" - - class QuasiAlignment: - """ - A class to represent a quasi-alignment that contains multiple segments. - Each segment is associated with a sequence from a particular species. - """ - def __init__(self, cluster_id: str): - """ - Initialize a QuasiAlignment object. - - :param cluster_id: A unique identifier for the quasi-alignment cluster. - """ - self.cluster_id = cluster_id - self.segments = [] - self.medoid = None - - def add_segment(self, new_segment: "QuasiAlignmentStrategy.Segment"): - """ - Adds a new Segment to the quasi-alignment and recalculates the medoid. - The medoid is the segment with the smallest average distance to all other segments. - - :param new_segment: Segment object to be added to the quasi-alignment. - """ - # Add the new segment to the list - self.segments.append(new_segment) - - # Recalculate the medoid if there is more than one segment - if len(self.segments) > 1: - # Calculate the average distance of each segment to all others - min_avg_distance = float('inf') - new_medoid = None - - for segment in self.segments: - distances = [segment.get_distance(other, QuasiAlignmentStrategy.manhattan_distance) for other in self.segments if other != segment] - avg_distance = np.mean(distances) - - # Update the medoid if a smaller average distance is found - if avg_distance < min_avg_distance: - min_avg_distance = avg_distance - new_medoid = segment - - # Set the new medoid - self.medoid = new_medoid - else: - # If only one segment, it's the medoid by default - self.medoid = new_segment - - def get_segments_by_species(self, species: str): - """ - Get all segments that belong to a specific species. - - :param species: The species to filter segments by. - :return: A list of segments that belong to the specified species. - """ - return [segment for segment in self.segments if segment.species_name == species] - # --- Public Methods --- - - def __repr__(self): - return f"QuasiAlignment(cluster_id={self.cluster_id}, segments={self.segments}, medoid={self.medoid})" def __init__(self, species_dict_file, primer_parameters_file): self.species_dict = self.parse_species_csv(species_dict_file) diff --git a/PRISTINE/data_structures.py b/PRISTINE/data_structures.py new file mode 100644 index 0000000..ecbd1dc --- /dev/null +++ b/PRISTINE/data_structures.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python3 +""" +Data structures for quasi-alignment and genomic segment representation. + +This module contains classes for representing genomic segments and their +quasi-alignments, used in the QuasiAlignmentStrategy pipeline. +""" + +from __future__ import annotations +import numpy as np +from similarity_metrics import manhattan_distance, add_missing_keys_with_zero + + +class Segment: + """ + A class to represent an individual segment in a quasi-alignment. It includes information about the sequence ID, species, + and the position of the segment within the sequence. + """ + + def __init__(self, seq_id: str, species_name: str, gene_name: str, start: int, length: int, sequence: str, p: int = 3): + """ + Initialize a Segment object. + + :param seq_id: The ID of the sequence from which the segment comes. + :param species_name: The species of the sequence. + :param gene_name: The gene region of the sequence. + :param start: The starting position of the segment within the sequence. + :param length: The length of the segment within the sequence. + :param sequence: The actual sequence string. + :param p: The length of p-mers for composition analysis (default: 3). + """ + self.seq_id = seq_id + self.species_name = species_name + self.gene_name = gene_name + self.start = start # 0-indexed + self.length = length + self.sequence = sequence + self.pmer_profile = self.get_pmer_composition(self, p) + + def get_segment(self, sequence): + """Returns the segment of the sequence based on start and length.""" + end = self.start + self.length + return sequence[self.start:end] + + def get_pmer_composition(self, segment, p): + """ + Takes a Segment object and computes the p-mer composition of the segment. + The output is a dictionary where keys are only the p-mers present in the segment's sequence, and values are their counts. + P-mers are stored in uppercase. + + :param segment: A Segment object containing the start and length information. + :param p: The length of the p-mer (substring). + :return: A dictionary with present p-mers as keys and their counts as values. + """ + sequence = segment.sequence.upper() # Convert the entire sequence to uppercase + pmer_dict = {} + + # Loop through the sequence and extract p-mers + for i in range(len(sequence) - p + 1): # Ensure we don't go beyond the sequence + pmer = sequence[i:i + p] + + # Add or update the count of the p-mer in the dictionary + if pmer in pmer_dict: + pmer_dict[pmer] += 1 + else: + pmer_dict[pmer] = 1 + + return pmer_dict + + def get_distance(self, other_segment, distance_func): + """ + Calculates the distance between this segment and another segment using a specified distance function. + + :param other_segment: The other Segment object to calculate the distance to. + :param distance_func: A function that takes two lists (or tuples) as input and returns the distance between them. + :return: The distance between the two p-mer profiles. + """ + if not callable(distance_func): + raise ValueError("distance_func must be a callable function") + + aligned_profile_self, aligned_profile_other = add_missing_keys_with_zero( + self.pmer_profile, other_segment.pmer_profile + ) + + profile_values_self = tuple(aligned_profile_self.values()) + profile_values_other = tuple(aligned_profile_other.values()) + + # Calculate and return the distance using the specified distance function + return distance_func(profile_values_self, profile_values_other) + + def __repr__(self): + return f"Segment(seq_id={self.seq_id}, species={self.species_name}, start={self.start}, length={self.length})" + + +class QuasiAlignment: + """ + A class to represent a quasi-alignment that contains multiple segments. + Each segment is associated with a sequence from a particular species. + """ + + def __init__(self, cluster_id: str): + """ + Initialize a QuasiAlignment object. + + :param cluster_id: A unique identifier for the quasi-alignment cluster. + """ + self.cluster_id = cluster_id + self.segments = [] + self.medoid = None + + def add_segment(self, new_segment: Segment): + """ + Adds a new Segment to the quasi-alignment and recalculates the medoid. + The medoid is the segment with the smallest average distance to all other segments. + + :param new_segment: Segment object to be added to the quasi-alignment. + """ + # Add the new segment to the list + self.segments.append(new_segment) + + # Recalculate the medoid if there is more than one segment + if len(self.segments) > 1: + # Calculate the average distance of each segment to all others + min_avg_distance = float('inf') + new_medoid = None + + for segment in self.segments: + distances = [ + segment.get_distance(other, manhattan_distance) + for other in self.segments if other != segment + ] + avg_distance = np.mean(distances) + + # Update the medoid if a smaller average distance is found + if avg_distance < min_avg_distance: + min_avg_distance = avg_distance + new_medoid = segment + + # Set the new medoid + self.medoid = new_medoid + else: + # If only one segment, it's the medoid by default + self.medoid = new_segment + + def get_segments_by_species(self, species: str): + """ + Get all segments that belong to a specific species. + + :param species: The species to filter segments by. + :return: A list of segments that belong to the specified species. + """ + return [segment for segment in self.segments if segment.species_name == species] + + def __repr__(self): + return f"QuasiAlignment(cluster_id={self.cluster_id}, num_segments={len(self.segments)}, medoid={self.medoid})" diff --git a/PRISTINE/similarity_metrics.py b/PRISTINE/similarity_metrics.py new file mode 100644 index 0000000..c6e04e9 --- /dev/null +++ b/PRISTINE/similarity_metrics.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python3 +""" +Similarity and distance metrics for sequence analysis. + +This module contains various distance and similarity functions used for +comparing sequence profiles, k-mer frequencies, and other vector representations. +""" + +from __future__ import annotations +import numpy as np +from collections import Counter +from datasketch import WeightedMinHashGenerator + + +def manhattan_distance(point1, point2): + """ + Calculate the Manhattan distance between two points. + + :param point1: A list or tuple representing the coordinates of the first point. + :param point2: A list or tuple representing the coordinates of the second point. + :return: The Manhattan distance between the two points. + """ + if len(point1) != len(point2): + raise ValueError("Both points must have the same number of dimensions") + + distance = sum(abs(p1 - p2) for p1, p2 in zip(point1, point2)) + return distance + + +def cosine_similarity(vector_a, vector_b): + """ + Calculates the cosine similarity between two vectors. + + :param vector_a: A list or tuple representing the first vector. + :param vector_b: A list or tuple representing the second vector. + :return: The cosine similarity between vector_a and vector_b. + """ + # Ensure inputs are lists or tuples + if not isinstance(vector_a, (list, tuple)) or not isinstance(vector_b, (list, tuple)): + raise TypeError("Input vectors must be lists or tuples.") + + # Convert to NumPy arrays + vector_a = np.array(vector_a) + vector_b = np.array(vector_b) + + # Calculate dot product and magnitudes + dot_product = np.dot(vector_a, vector_b) + magnitude_a = np.linalg.norm(vector_a) + magnitude_b = np.linalg.norm(vector_b) + + # Avoid division by zero + if magnitude_a == 0 or magnitude_b == 0: + return 0 # Return 0 for cosine similarity if either vector is zero + + # Compute cosine similarity + return dot_product / (magnitude_a * magnitude_b) + + +def jaccard_similarity_with_frequencies(vector_a, vector_b): + """ + Calculates the Jaccard similarity between two lists or tuples containing frequencies of individual items (e.g., p-mers). + + :param vector_a: A list or tuple representing the first vector. + :param vector_b: A list or tuple representing the second vector. + :return: The Jaccard similarity between the two collections, accounting for frequencies. + """ + # Convert lists to Counters to handle frequencies + counter_a = Counter(vector_a) + counter_b = Counter(vector_b) + + # Calculate intersection and union based on minimum and maximum frequencies + intersection = sum((min(counter_a[item], counter_b[item]) for item in counter_a if item in counter_b)) + union = sum((max(counter_a[item], counter_b[item]) for item in set(counter_a) | set(counter_b))) + + # Avoid division by zero in case both counters are empty + if union == 0: + return 0 + + # Compute Jaccard similarity with frequencies + return intersection / union + + +def calculate_weighted_minhash_similarity(frequency_vector_a, frequency_vector_b, num_perm=128): + """ + Calculates the Weighted MinHash Jaccard similarity between two frequency vectors. + + :param frequency_vector_a: A list or tuple representing the frequency of p-mers for the first vector. + :param frequency_vector_b: A list or tuple representing the frequency of p-mers for the second vector. + :param num_perm: The number of permutations (hash functions) for the MinHash. + :return: The Weighted MinHash Jaccard similarity between the two vectors. + """ + if not isinstance(frequency_vector_a, (list, tuple)) or not isinstance(frequency_vector_b, (list, tuple)): + raise TypeError("Input vectors must be lists or tuples.") + + if len(frequency_vector_a) != len(frequency_vector_b): + raise ValueError("The two frequency vectors must be of the same length.") + + wmg = WeightedMinHashGenerator(len(frequency_vector_a), sample_size=num_perm) + + wm_a = wmg.minhash(frequency_vector_a) + wm_b = wmg.minhash(frequency_vector_b) + + return wm_a.jaccard(wm_b) + + +def add_missing_keys_with_zero(dict1, dict2): + """ + Adds keys that are present in dict1 but not in dict2 (and vice versa) to each dictionary with values set to zero. + + :param dict1: The first dictionary. + :param dict2: The second dictionary. + :return: Two dictionaries with added keys and values set to zero where keys were missing. + """ + # Find keys unique to each dictionary + keys_in_dict1_not_in_dict2 = dict1.keys() - dict2.keys() + keys_in_dict2_not_in_dict1 = dict2.keys() - dict1.keys() + + # Add missing keys to each dictionary with values set to zero + for key in keys_in_dict1_not_in_dict2: + dict2[key] = 0 + for key in keys_in_dict2_not_in_dict1: + dict1[key] = 0 + + return dict1, dict2 diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta index 579e85f..f7a64e1 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_1.fasta @@ -1,4 +1,4 @@ ->Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|CP130901 Yersinia pseudotuberculosis strain SP-1303 chromosome, complete genome. [Yersinia pseudotuberculosis SP-1303 | 633.777] +>Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|CP130901 Yersinia pseudotuberculosis strain SP-1303 chromosome, complete genome. [Yersinia pseudotuberculosis SP-1303 | 633.777] gtgtcactttcgctttggcagcagtgtcttgcccgattgcaggatgagttacctgccaca gaatttagtatgtggatacgccccctacaggcggaactgagtgacaatactctggcgctt tacgcacctaatcgttttgtactggactgggtccgtgataagtacttaaacaatatcaat diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta index 28a3fdc..32aca81 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/Yersinia_pseudotuberculosis_2.fasta @@ -1,4 +1,4 @@ ->Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|LT596221 Yersinia pseudotuberculosis isolate NZYP4713 genome assembly, chromosome: 1. [Yersinia pseudotuberculosis strain NZYP4713 | 633.88] +>Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_Yersinia_pseudotuberculosis_accn|LT596221 Yersinia pseudotuberculosis isolate NZYP4713 genome assembly, chromosome: 1. [Yersinia pseudotuberculosis strain NZYP4713 | 633.88] gtgtcactttcgctttggcagcagtgtcttgcccgattgcaggatgagttacctgccaca gaatttagtatgtggatacgccccctacaggcggaactgagtgacaatactctggcgctt tacgcacctaatcgttttgtactggactgggtccgtgataagtacttaaacaatatcaat diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta index e000d9b..2cafc0f 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_1.fasta @@ -1,4 +1,4 @@ ->Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177922 Yersinia pestis strain 50008 chromosome, complete genome. [Yersinia pestis 50008 | 632.1142] +>Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177922 Yersinia pestis strain 50008 chromosome, complete genome. [Yersinia pestis 50008 | 632.1142] ggcggactccactcgaacaagatcaataaaagaaagaaatcaggatattttgtaactcat cgttttgacacacgccatatggcctgtagcaagccaaggtgacatgcgaactgtttttta cagttgttttttactgctgattttttacgactcacattcccccacgatcccagagtagag diff --git a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta index daf6c91..40c5334 100755 --- a/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta +++ b/PRISTINE/toy_dataset/data/raw_genomes/Yersinia_pseudotuberculosis/non-targets/Yersinia_pestis_2.fasta @@ -1,4 +1,4 @@ ->Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177934 Yersinia pestis strain 96005 chromosome, complete genome. [Yersinia pestis 96005 | 632.1143] +>Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_Yersinia_pestis_accn|CP177934 Yersinia pestis strain 96005 chromosome, complete genome. [Yersinia pestis 96005 | 632.1143] ggcggactccactcgaacaagatcaataaaagaaagaaatcaggatattttgtaactcat cgttttgacacacgccatatggcctgtagcaagccaaggtgacatgcgaactgtttttta cagttgttttttactgctgattttttacgactcacattcccccacgatcccagagtagag diff --git a/PRISTINE/visualization.py b/PRISTINE/visualization.py new file mode 100644 index 0000000..073d6bb --- /dev/null +++ b/PRISTINE/visualization.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python3 +""" +Visualization functions for PRISTINE analysis results. + +This module contains plotting functions for generating heatmaps and +distribution plots of SNP informativeness and density. +""" + +from __future__ import annotations +import os +import matplotlib.pyplot as plt +import seaborn as sns +import pandas as pd + + +def plot_informativeness_heatmap(csv_file, output_file): + """ + Generate a heatmap showing informative SNP proportions per locus and species. + + :param csv_file: Path to the SNP summary CSV file + :param output_file: Path where the heatmap PNG will be saved + :return: Path to the saved heatmap file + """ + df = pd.read_csv(csv_file) + + prop_cols = [col for col in df.columns if col.startswith("Prop_") and col != "Avg_Prop_Informative_SNPs"] + if not prop_cols: + print("No proportion columns found in the dataset.") + return None + + heatmap_data = df.set_index("Locus")[prop_cols] + + if "Avg_Prop_Informative_SNPs" in df.columns: + sorted_loci = df.sort_values("Avg_Prop_Informative_SNPs", ascending=False)["Locus"] + heatmap_data = heatmap_data.loc[sorted_loci] + + plt.figure(figsize=(12, max(6, 0.3 * len(heatmap_data)))) + sns.heatmap(heatmap_data, annot=False, cmap="YlOrRd", cbar_kws={'label': 'Proportion of Informative SNPs'}) + plt.title("Informative SNP Proportions per Locus and Non-Target Species") + plt.xlabel("Non-Target Species") + plt.ylabel("Locus") + plt.tight_layout() + plt.savefig(output_file) + plt.close() + + return output_file + + +def plot_snp_density_distribution(csv_file, output_dir, top_n=5): + """ + Generate SNP position distribution plots for top-ranked loci. + + :param csv_file: Path to the SNP summary CSV file + :param output_dir: Directory where SNP density plots will be saved + :param top_n: Number of top-ranked loci to plot (default: 5) + :return: Path to the directory containing the plots + """ + os.makedirs(output_dir, exist_ok=True) + + df = pd.read_csv(csv_file) + + if "Avg_Prop_Informative_SNPs" not in df.columns: + print("Missing ranking column 'Avg_Prop_Informative_SNPs'.") + return None + + top_loci = df.sort_values("Avg_Prop_Informative_SNPs", ascending=False).head(top_n) + + prop_cols = [col for col in df.columns if col.startswith("SNP_Pos_")] + + for _, row in top_loci.iterrows(): + locus_name = row["Locus"] + plt.figure(figsize=(12, 3)) + + for col in prop_cols: + species = col.replace("SNP_Pos_", "") + if pd.isna(row[col]) or not row[col].strip(): + continue + try: + positions = list(map(int, row[col].split(','))) + except ValueError: + continue + sns.histplot(positions, bins=50, kde=False, label=species, element="step", fill=False) + + plt.title(f"SNP Position Distribution - {locus_name}") + plt.xlabel("Alignment Position") + plt.ylabel("SNP Count") + plt.legend(title="Non-Target Species", loc="upper right", fontsize="small") + plt.tight_layout() + plot_path = os.path.join(output_dir, f"snp_density_{locus_name.replace('.', '_')}.png") + plt.savefig(plot_path) + plt.close() + + print(f"SNP density plots saved to: {output_dir}") + return output_dir