From b35ea9f6b416caad279c94dc5f0107d0f9a463f1 Mon Sep 17 00:00:00 2001 From: minagayid Date: Mon, 3 Aug 2026 23:08:25 +0300 Subject: [PATCH] build offline genomics toolkit --- .github/workflows/tests.yml | 22 +++ .gitignore | 6 +- README.md | 147 +++++++++---------- examples/demo_genomics.py | 28 ++-- genopedia/__init__.py | 20 +++ genopedia/__main__.py | 6 + genopedia/adapters.py | 47 ++++++ genopedia/cli.py | 109 ++++++++++++++ genopedia/core.py | 278 ++++++++++++++++++++++++++++++++++++ genopedia/io.py | 204 ++++++++++++++++++++++++++ genopedia/models.py | 117 +++++++++++++++ genopedia/reporting.py | 151 ++++++++++++++++++++ pyproject.toml | 28 ++++ requirements-dev.txt | 4 + requirements.txt | 27 +--- run_genopedia.py | 8 ++ src/genomics/__init__.py | 275 ----------------------------------- src/models/__init__.py | 70 --------- tests/test_adapters.py | 29 ++++ tests/test_cli.py | 75 ++++++++++ tests/test_genomics.py | 77 +++++----- tests/test_io.py | 57 ++++++++ tests/test_reporting.py | 28 ++++ 23 files changed, 1317 insertions(+), 496 deletions(-) create mode 100644 .github/workflows/tests.yml create mode 100644 genopedia/__init__.py create mode 100644 genopedia/__main__.py create mode 100644 genopedia/adapters.py create mode 100644 genopedia/cli.py create mode 100644 genopedia/core.py create mode 100644 genopedia/io.py create mode 100644 genopedia/models.py create mode 100644 genopedia/reporting.py create mode 100644 pyproject.toml create mode 100644 requirements-dev.txt create mode 100644 run_genopedia.py delete mode 100644 src/genomics/__init__.py delete mode 100644 src/models/__init__.py create mode 100644 tests/test_adapters.py create mode 100644 tests/test_cli.py create mode 100644 tests/test_io.py create mode 100644 tests/test_reporting.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..fe8c77a --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,22 @@ +name: tests + +on: + push: + pull_request: + +jobs: + unit: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10", "3.11", "3.12"] + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Run standard-library tests + run: python -m unittest discover -s tests -v + - name: Check compilation + run: python -m compileall -q genopedia examples tests + diff --git a/.gitignore b/.gitignore index 90eade1..f19227c 100644 --- a/.gitignore +++ b/.gitignore @@ -94,6 +94,10 @@ logs/ *.tmp *~ +# Genopedia local reports +genopedia-report.html +genopedia-comparison.html + # Data data/ *.npy @@ -107,4 +111,4 @@ data/ # Dotenv .env.local -.env.*.local \ No newline at end of file +.env.*.local diff --git a/README.md b/README.md index 1841790..73f465a 100644 --- a/README.md +++ b/README.md @@ -1,100 +1,83 @@ -# Genopedia - Genomics Machine Learning Pipeline +# Genopedia -**Genopedia** is a genomics-focused machine learning pipeline for DNA/RNA sequence analysis, variant detection, and color-coded visualization. +Genopedia is an offline-first toolkit for inspecting DNA/RNA sequences and producing reproducible local reports. It is designed to run after download on a normal Python installation without TensorFlow, PyTorch, a database, or an internet connection. -## โœ… Core Features +The current release provides: -### Genomics Pipeline -- **Nucleotide color coding**: A=๐Ÿ”ด, T=๐ŸŸก, G=๐ŸŸข, C=๐Ÿ”ต, U=๐ŸŸค -- **Data loaders**: FASTA, FASTQ, VCF -- **Sequence analysis**: GC content, motif finding, functional region detection -- **Variant detection**: SNPs, insertions/deletions, pathogenicity scoring -- **Visualization**: HTML/SVG color-coded sequences +- Streaming FASTA, FASTQ, VCF, VCF.GZ, and raw-sequence readers. +- Sequence validation, GC fraction, ambiguous-base counts, Phred summaries, and warnings. +- Overlapping motif search and transparent candidate-region heuristics. +- Sequence comparison with substitutions, insertions, deletions, and replacements. +- A deterministic dependency-free k-mer classifier. +- Self-contained HTML/SVG reports and JSON summaries. +- A research-only interpretation boundary: no diagnosis, treatment advice, or proposed DNA edits. -### ML Models -- DNA sequence classifier (basic k-mer model implemented) -- Variant pathogenicity prediction -- Gene function annotation -- Extensible architecture for deep learning models +PCR and sequencing instruments are vendor-specific. Genopedia accepts exported files by default and exposes a clean boundary for future vendor adapters; it does not pretend that generic software can control every PCR device without its model and communication protocol. -## ๐Ÿ“ฆ Project Structure +## Quick start from a fresh download + +Requirements: Python 3.10 or newer. The core runtime uses only the Python standard library. + +```powershell +python -m genopedia demo --length 120 --output report.html --json-output report.json ``` -genopedia/ -โ”œโ”€โ”€ src/ -โ”‚ โ”œโ”€โ”€ genomics/ # Core genomics pipeline -โ”‚ โ”œโ”€โ”€ models/ # ML models -โ”‚ โ”œโ”€โ”€ api/ # FastAPI backend -โ”‚ โ””โ”€โ”€ utils/ # Utilities -โ”œโ”€โ”€ notebooks/ # Jupyter notebooks -โ”œโ”€โ”€ tests/ # Unit tests -โ”œโ”€โ”€ examples/ # Demo scripts -โ”œโ”€โ”€ docs/ # Documentation -โ””โ”€โ”€ requirements.txt + +Open `report.html` locally. No server or deployment is required. + +Analyze a file: + +```powershell +python -m genopedia analyze sample.fasta --output sample-report.html --json-output sample.json +python -m genopedia analyze sample.fastq --output reads-report.html +python -m genopedia analyze variants.vcf.gz --output variants-report.html +``` + +The convenience launcher also works directly from the repository root: + +```powershell +python run_genopedia.py demo ``` -## ๐Ÿš€ Quick Start -```bash -cd ~/Desktop/genopedia -pip install -r requirements.txt +## Development -# Run a demo -python examples/demo_genomics.py +Run the built-in test suite without installing pytest: -# Run tests -python -m pytest tests/ +```powershell +python -m unittest discover -s tests -v ``` -## ๐Ÿ”ฌ Testing -Tests cover: -- Data loading/generation -- Sequence analysis -- Variant detection -- Visualization utilities +Install the package in editable mode if you want the `genopedia` command: -```bash -python -m pytest tests/test_genomics.py -v +```powershell +python -m pip install -e . +genopedia demo ``` -## ๐Ÿ“Š Visualization -Generates HTML/SVG color-coded nucleotide sequences: +Optional integrations are intentionally isolated from the core: -```python -from src.genomics import GenomicsVisualizer -visualizer = GenomicsVisualizer() -dna_sequence = "ATGCCGTAG" -html_output = visualizer.generate_color_html(dna_sequence) +```powershell +python -m pip install -e .[bio] # Biopython helpers +python -m pip install -e .[ml] # NumPy/scikit-learn extensions +python -m pip install -e .[dev] # pytest, formatting, and lint tools ``` -## ๐ŸŒ API (Planned) -FastAPI backend for: -- Sequence analysis endpoints -- Variant detection -- ML model inference -- WebSocket for real-time visualization - -## ๐Ÿ”— Cloud Dataset Access -No large files stored locally - access datasets via: -- **GRCh38**: https://www.ncbi.nlm.nih.gov/grc/human -- **1000 Genomes**: s3://1000genomes/ -- **ClinVar**: ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/ -- **Ensembl**: ftp://ftp.ensembl.org/pub/ -- **dbSNP**: ftp://ftp.ncbi.nlm.nih.gov/snp/ - -## ๐Ÿงฌ Genomics Context -Based on: -- [ChatGPT conversation](https://chatgpt.com/share/6a2aba33-10d8-83ea-992e-a1798edd9493) -- [Qwen conversation](https://chat.qwen.ai/s/d7f76f3a-5459-45be-9e8f-2cce2bd72079) -- [Kimi conversation](https://www.kimi.com/share/19eb6e89-8c42-8704-8000-0000b9440601) - -## ๐Ÿ”ฎ Future Enhancements -1. Deep learning-based variant pathogenicity prediction -2. Gene expression analysis -3. CRISPR target identification -4. Web dashboard with interactive visualization -5. Integration with real PCR machines - -## ๐Ÿ“œ License -MIT License - see [LICENSE](LICENSE) ``` - ---- -ยฉ 2026 Genopedia Project | [GitHub](#) \ No newline at end of file +## Coordinate and safety conventions + +Sequence comparisons use zero-based positions. VCF positions remain one-based and are marked as such in variant metadata. Variant interpretations are `unknown` unless supplied evidence matches a supported evidence label. Candidate promoters and start codons are heuristics, not gene annotation. + +Genopedia is intended for research and engineering workflows. Any clinical interpretation, laboratory action, or genetic intervention requires validated laboratory methods and qualified human review. + +## Project layout + +```text +genopedia/ +โ”œโ”€โ”€ genopedia/ # dependency-free package and CLI +โ”œโ”€โ”€ examples/ # runnable example +โ”œโ”€โ”€ tests/ # standard-library tests +โ”œโ”€โ”€ pyproject.toml # install metadata and optional extras +โ””โ”€โ”€ run_genopedia.py # fresh-checkout launcher +``` + +## License + +MIT. See [LICENSE](LICENSE). diff --git a/examples/demo_genomics.py b/examples/demo_genomics.py index d927151..307188e 100644 --- a/examples/demo_genomics.py +++ b/examples/demo_genomics.py @@ -1,14 +1,14 @@ -"""GenoProject - resume point.""" -from GenoProject.src.genomics import GenomicsDataLoader, SequenceAnalyzer, GenomicsVisualizer -from GenoProject.src.models import DNAClassifier - -loader = GenomicsDataLoader("data/genomics") -sequence = loader.generate_synthetic_dna(120) -analyzer = SequenceAnalyzer() -variants = analyzer.detect_variants("ATGCCGTAG", "ATGTCGTAG") -visualizer = GenomicsVisualizer() - -print("[GenoProject] sequence:", sequence) -print("[GenoProject] variants:", [v.to_dict() for v in variants]) -print("[GenoProject] html preview:") -print(visualizer.generate_color_html(sequence, 20)) +"""Run the dependency-free Genopedia demonstration from a fresh checkout.""" + +from pathlib import Path +import sys + + +# Make ``python examples/demo_genomics.py`` work before editable installation. +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + +from genopedia.cli import main # noqa: E402 + + +if __name__ == "__main__": + raise SystemExit(main(["demo", "--length", "120"])) diff --git a/genopedia/__init__.py b/genopedia/__init__.py new file mode 100644 index 0000000..b13ffa1 --- /dev/null +++ b/genopedia/__init__.py @@ -0,0 +1,20 @@ +"""Genopedia: lightweight, offline-first genomics analysis tools.""" + +from .core import ( + FunctionalRegion, + QualityReport, + SequenceAnalyzer, + SequenceRecord, + Variant, +) + +__all__ = [ + "FunctionalRegion", + "QualityReport", + "SequenceAnalyzer", + "SequenceRecord", + "Variant", +] + +__version__ = "0.2.0" + diff --git a/genopedia/__main__.py b/genopedia/__main__.py new file mode 100644 index 0000000..0b6ae7c --- /dev/null +++ b/genopedia/__main__.py @@ -0,0 +1,6 @@ +from .cli import main + + +if __name__ == "__main__": + raise SystemExit(main()) + diff --git a/genopedia/adapters.py b/genopedia/adapters.py new file mode 100644 index 0000000..8e5cd60 --- /dev/null +++ b/genopedia/adapters.py @@ -0,0 +1,47 @@ +"""Vendor-neutral local adapters for exported instrument files.""" + +from __future__ import annotations + +from pathlib import Path + +from .io import InputData, read_input + + +SUPPORTED_SUFFIXES = frozenset( + {".fa", ".fasta", ".fna", ".fq", ".fastq", ".vcf", ".gz"} +) + + +class FileDropAdapter: + """Read sequence exports placed in a local directory. + + The adapter intentionally does not assume a PCR vendor, USB protocol, or + network service. A vendor-specific adapter can implement the same + ``discover``/``read`` boundary later without changing analysis code. + """ + + def __init__(self, directory: str | Path, recursive: bool = False) -> None: + self.directory = Path(directory).resolve() + self.recursive = recursive + + def discover(self) -> list[Path]: + if not self.directory.is_dir(): + raise FileNotFoundError(f"file-drop directory does not exist: {self.directory}") + candidates = self.directory.rglob("*") if self.recursive else self.directory.iterdir() + supported: list[Path] = [] + for path in candidates: + if not path.is_file(): + continue + name = path.name.lower() + if name.endswith((".fa", ".fasta", ".fna", ".fq", ".fastq", ".vcf", ".fa.gz", ".fasta.gz", ".fq.gz", ".fastq.gz", ".vcf.gz")): + supported.append(path) + return sorted(supported, key=lambda path: path.name.lower()) + + def read(self, path: str | Path, format: str = "auto") -> InputData: + source = Path(path).resolve() + try: + source.relative_to(self.directory) + except ValueError as error: + raise ValueError("file-drop adapter cannot read outside its configured directory") from error + return read_input(source, format=format) + diff --git a/genopedia/cli.py b/genopedia/cli.py new file mode 100644 index 0000000..2dee7ee --- /dev/null +++ b/genopedia/cli.py @@ -0,0 +1,109 @@ +"""Command-line interface for Genopedia.""" + +from __future__ import annotations + +import argparse +import json +import random +from pathlib import Path +from typing import Sequence + +from .core import SequenceAnalyzer, SequenceRecord +from .io import read_input +from .reporting import render_html_report, write_html_report + + +def _synthetic_record(length: int, seed: int) -> SequenceRecord: + if length < 1: + raise ValueError("length must be at least 1") + generator = random.Random(seed) + sequence = "".join(generator.choice("ATGC") for _ in range(length)) + return SequenceRecord(identifier="synthetic", sequence=sequence, source_format="synthetic") + + +def _summary(records: list[SequenceRecord], analyzer: SequenceAnalyzer, variants: list) -> dict[str, object]: + reports = analyzer.summarize_records(records) + return { + "records": [report.to_dict() for report in reports], + "variants": [variant.to_dict() for variant in variants], + } + + +def _write_outputs( + records: list[SequenceRecord], + variants: list, + output: Path, + json_output: Path | None, +) -> None: + analyzer = SequenceAnalyzer() + reports = analyzer.summarize_records(records) + write_html_report( + output, + render_html_report(records, variants=variants, quality_reports=reports), + ) + if json_output is not None: + json_output.parent.mkdir(parents=True, exist_ok=True) + json_output.write_text( + json.dumps(_summary(records, analyzer, variants), indent=2, sort_keys=True), + encoding="utf-8", + ) + + +def build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + prog="genopedia", + description="Offline-first DNA/RNA sequence inspection and reporting.", + ) + commands = parser.add_subparsers(dest="command", required=True) + + demo = commands.add_parser("demo", help="generate a deterministic synthetic report") + demo.add_argument("--length", type=int, default=120) + demo.add_argument("--seed", type=int, default=7) + demo.add_argument("--output", type=Path, default=Path("genopedia-report.html")) + demo.add_argument("--json-output", type=Path) + + analyze = commands.add_parser("analyze", help="analyze FASTA, FASTQ, VCF, or raw sequence text") + analyze.add_argument("input", type=Path) + analyze.add_argument("--format", default="auto", choices=["auto", "fasta", "fastq", "vcf", "text"]) + analyze.add_argument("--output", type=Path, default=Path("genopedia-report.html")) + analyze.add_argument("--json-output", type=Path) + + compare = commands.add_parser("compare", help="compare a reference sequence with a sample sequence") + compare.add_argument("reference", type=Path) + compare.add_argument("sample", type=Path) + compare.add_argument("--chromosome", default="chr1") + compare.add_argument("--output", type=Path, default=Path("genopedia-comparison.html")) + compare.add_argument("--json-output", type=Path) + return parser + + +def main(argv: Sequence[str] | None = None) -> int: + parser = build_parser() + args = parser.parse_args(argv) + try: + if args.command == "demo": + records = [_synthetic_record(args.length, args.seed)] + variants = [] + elif args.command == "analyze": + data = read_input(args.input, format=args.format) + records = list(data.records) + variants = list(data.variants) + else: + reference = read_input(args.reference, format="auto") + sample = read_input(args.sample, format="auto") + if not reference.records or not sample.records: + raise ValueError("compare requires sequence files containing at least one record each") + records = [reference.records[0], sample.records[0]] + variants = SequenceAnalyzer().detect_variants( + records[0].sequence, + records[1].sequence, + chromosome=args.chromosome, + ) + _write_outputs(records, variants, args.output, args.json_output) + print(f"Report written to {args.output}") + if args.json_output: + print(f"Summary written to {args.json_output}") + return 0 + except (OSError, ValueError) as error: + parser.error(str(error)) + return 2 diff --git a/genopedia/core.py b/genopedia/core.py new file mode 100644 index 0000000..a2eecf1 --- /dev/null +++ b/genopedia/core.py @@ -0,0 +1,278 @@ +"""Core genomics data structures and deterministic analysis.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from difflib import SequenceMatcher +from typing import Any, Iterable + + +DNA_ALPHABET = frozenset("ACGTNRYKMSWBDHV") +RNA_ALPHABET = frozenset("ACGUNRYKMSWBDHV") +UNAMBIGUOUS_BASES = frozenset("ACGTU") +INTERPRETATIONS = frozenset( + {"benign", "likely_benign", "uncertain", "likely_pathogenic", "pathogenic", "unknown"} +) + + +def _json_safe(value: Any) -> Any: + if isinstance(value, dict): + return {str(key): _json_safe(item) for key, item in value.items()} + if isinstance(value, (list, tuple)): + return [_json_safe(item) for item in value] + return value + + +def normalize_sequence(sequence: str) -> str: + """Remove formatting whitespace and normalize bases to uppercase.""" + + if not isinstance(sequence, str): + raise TypeError("sequence must be a string") + return "".join(sequence.split()).upper() + + +@dataclass(frozen=True) +class SequenceRecord: + """A sequence and its optional sequencing metadata.""" + + identifier: str + sequence: str + description: str = "" + quality_scores: tuple[int, ...] = () + source_format: str = "text" + metadata: dict[str, Any] = field(default_factory=dict) + + def to_dict(self) -> dict[str, Any]: + return { + "identifier": self.identifier, + "sequence": self.sequence, + "description": self.description, + "length": len(self.sequence), + "quality_scores": list(self.quality_scores), + "source_format": self.source_format, + "metadata": _json_safe(self.metadata), + } + + +@dataclass(frozen=True) +class Variant: + """An observed sequence difference. + + Positions produced by sequence comparison are zero-based. Positions read + from VCF retain VCF's one-based coordinate system in metadata. + """ + + chromosome: str + position: int + reference: str + observed: str + quality: float | None = None + variant_type: str = "substitution" + interpretation: str = "unknown" + metadata: dict[str, Any] = field(default_factory=dict) + + def to_dict(self) -> dict[str, Any]: + return { + "chromosome": self.chromosome, + "position": self.position, + "reference": self.reference, + "observed": self.observed, + "quality": self.quality, + "variant_type": self.variant_type, + "interpretation": self.interpretation, + "metadata": _json_safe(self.metadata), + } + + +@dataclass(frozen=True) +class FunctionalRegion: + """A candidate region identified by a transparent heuristic.""" + + name: str + start: int + end: int + region_type: str + evidence: str = "heuristic" + + def to_dict(self) -> dict[str, Any]: + return { + "name": self.name, + "start": self.start, + "end": self.end, + "region_type": self.region_type, + "evidence": self.evidence, + } + + +@dataclass(frozen=True) +class QualityReport: + """Quality-control result for one sequence record.""" + + identifier: str + length: int + gc_fraction: float + ambiguous_count: int + invalid_symbols: tuple[str, ...] + mean_quality: float | None + status: str + warnings: tuple[str, ...] = () + + def to_dict(self) -> dict[str, Any]: + return { + "identifier": self.identifier, + "length": self.length, + "gc_fraction": self.gc_fraction, + "ambiguous_count": self.ambiguous_count, + "invalid_symbols": list(self.invalid_symbols), + "mean_quality": self.mean_quality, + "status": self.status, + "warnings": list(self.warnings), + } + + +class SequenceAnalyzer: + """Small, dependency-free analysis primitives for DNA and RNA.""" + + def gc_content(self, sequence: str) -> float: + normalized = normalize_sequence(sequence) + informative = [base for base in normalized if base in UNAMBIGUOUS_BASES] + if not informative: + return 0.0 + return sum(base in {"G", "C"} for base in informative) / len(informative) + + def quality_control(self, record: SequenceRecord, molecule: str = "DNA") -> QualityReport: + normalized = normalize_sequence(record.sequence) + alphabet = DNA_ALPHABET if molecule.upper() == "DNA" else RNA_ALPHABET + invalid = tuple(sorted({base for base in normalized if base not in alphabet})) + ambiguous = sum(base not in UNAMBIGUOUS_BASES for base in normalized if base in alphabet) + warnings: list[str] = [] + + if not normalized: + warnings.append("sequence is empty") + if invalid: + warnings.append("sequence contains invalid symbols") + if ambiguous: + warnings.append("sequence contains ambiguous bases") + if record.quality_scores and len(record.quality_scores) != len(normalized): + warnings.append("quality score count does not match sequence length") + mean_quality = ( + sum(record.quality_scores) / len(record.quality_scores) + if record.quality_scores + else None + ) + if mean_quality is not None and mean_quality < 20: + warnings.append("mean Phred quality is below 20") + + status = "fail" if not normalized or invalid or ( + record.quality_scores and len(record.quality_scores) != len(normalized) + ) else ("warn" if warnings else "pass") + return QualityReport( + identifier=record.identifier, + length=len(normalized), + gc_fraction=self.gc_content(normalized), + ambiguous_count=ambiguous, + invalid_symbols=invalid, + mean_quality=mean_quality, + status=status, + warnings=tuple(warnings), + ) + + def find_motifs(self, sequence: str, motif: str) -> list[int]: + sequence = normalize_sequence(sequence) + motif = normalize_sequence(motif) + if not motif: + raise ValueError("motif must not be empty") + return [ + index + for index in range(len(sequence) - len(motif) + 1) + if sequence[index : index + len(motif)] == motif + ] + + def identify_functional_regions(self, sequence: str) -> list[FunctionalRegion]: + """Return candidate regions; this is not a clinical annotation.""" + + normalized = normalize_sequence(sequence) + regions: list[FunctionalRegion] = [] + for position in self.find_motifs(normalized, "TATA"): + regions.append( + FunctionalRegion( + name=f"candidate-promoter-{position}", + start=max(0, position - 25), + end=min(len(normalized), position + 25), + region_type="promoter", + ) + ) + for position in self.find_motifs(normalized, "ATG"): + regions.append( + FunctionalRegion( + name=f"candidate-start-codon-{position}", + start=position, + end=position + 3, + region_type="start_codon", + ) + ) + return regions + + def detect_variants( + self, reference: str, sample: str, chromosome: str = "chr1" + ) -> list[Variant]: + """Compare two strings and report substitutions and indels. + + The standard-library sequence matcher keeps the core usable on a + clean Python installation. It is intended for short/medium sequence + comparisons; dedicated aligners remain optional for large genomes. + """ + + reference = normalize_sequence(reference) + sample = normalize_sequence(sample) + matcher = SequenceMatcher(None, reference, sample, autojunk=False) + variants: list[Variant] = [] + for tag, ref_start, ref_end, sample_start, sample_end in matcher.get_opcodes(): + if tag == "equal": + continue + ref_segment = reference[ref_start:ref_end] + sample_segment = sample[sample_start:sample_end] + if tag == "replace" and len(ref_segment) == len(sample_segment): + for offset, (ref_base, sample_base) in enumerate( + zip(ref_segment, sample_segment) + ): + if ref_base != sample_base: + variants.append( + Variant( + chromosome=chromosome, + position=ref_start + offset, + reference=ref_base, + observed=sample_base, + variant_type="substitution", + ) + ) + else: + if not ref_segment: + variant_type = "insertion" + elif not sample_segment: + variant_type = "deletion" + else: + variant_type = "replacement" + variants.append( + Variant( + chromosome=chromosome, + position=ref_start, + reference=ref_segment, + observed=sample_segment, + variant_type=variant_type, + ) + ) + return variants + + def interpret_variant(self, variant: Variant, evidence: dict[str, Any] | None = None) -> str: + """Apply supplied evidence, otherwise remain explicitly unknown.""" + + evidence = evidence or variant.metadata + value = str(evidence.get("interpretation", evidence.get("clinical_significance", "unknown"))) + return value if value in INTERPRETATIONS else "unknown" + + def summarize_records( + self, records: Iterable[SequenceRecord], molecule: str = "DNA" + ) -> list[QualityReport]: + return [self.quality_control(record, molecule=molecule) for record in records] + diff --git a/genopedia/io.py b/genopedia/io.py new file mode 100644 index 0000000..9b9ca90 --- /dev/null +++ b/genopedia/io.py @@ -0,0 +1,204 @@ +"""Streaming readers for common sequence and variant text formats.""" + +from __future__ import annotations + +import gzip +from dataclasses import dataclass +from pathlib import Path +from typing import Iterator, TextIO + +from .core import SequenceRecord, Variant, normalize_sequence + + +PathLike = str | Path + + +def _open_text(path: PathLike) -> TextIO: + path = Path(path) + if path.suffix.lower() == ".gz": + return gzip.open(path, "rt", encoding="utf-8", errors="replace") + return path.open("r", encoding="utf-8", errors="replace") + + +def _record_identifier(header: str, marker: str) -> tuple[str, str]: + value = header.rstrip("\n\r") + if not value.startswith(marker): + raise ValueError(f"expected {marker!r} header") + content = value[len(marker) :].strip() + if not content: + raise ValueError("sequence header has no identifier") + fields = content.split(None, 1) + return fields[0], fields[1] if len(fields) == 2 else "" + + +def read_fasta(path: PathLike) -> Iterator[SequenceRecord]: + """Yield FASTA records without loading the full file into memory.""" + + identifier: str | None = None + description = "" + chunks: list[str] = [] + with _open_text(path) as handle: + for line_number, line in enumerate(handle, start=1): + if line.startswith(">"): + if identifier is not None: + yield SequenceRecord( + identifier=identifier, + sequence=normalize_sequence("".join(chunks)), + description=description, + source_format="fasta", + ) + identifier, description = _record_identifier(line, ">") + chunks = [] + elif line.strip(): + if identifier is None: + raise ValueError(f"FASTA sequence data before header at line {line_number}") + chunks.append(line.strip()) + if identifier is not None: + yield SequenceRecord( + identifier=identifier, + sequence=normalize_sequence("".join(chunks)), + description=description, + source_format="fasta", + ) + + +def read_fastq(path: PathLike) -> Iterator[SequenceRecord]: + """Yield FASTQ records and decode Sanger/Phred+33 qualities.""" + + with _open_text(path) as handle: + while True: + header = handle.readline() + if not header: + return + sequence_line = handle.readline() + plus_line = handle.readline() + quality_line = handle.readline() + if not sequence_line or not plus_line or not quality_line: + raise ValueError("truncated FASTQ record") + identifier, description = _record_identifier(header, "@") + if not plus_line.startswith("+"): + raise ValueError(f"FASTQ record {identifier!r} is missing its '+' line") + sequence = normalize_sequence(sequence_line) + quality_text = quality_line.rstrip("\r\n") + if len(sequence) != len(quality_text): + raise ValueError(f"FASTQ record {identifier!r} has mismatched sequence and quality lengths") + quality_scores = tuple(max(0, ord(value) - 33) for value in quality_text) + yield SequenceRecord( + identifier=identifier, + sequence=sequence, + description=description, + quality_scores=quality_scores, + source_format="fastq", + ) + + +def _parse_info(value: str) -> dict[str, object]: + metadata: dict[str, object] = {} + if value in {"", "."}: + return metadata + for item in value.split(";"): + if "=" not in item: + metadata[item] = True + continue + key, item_value = item.split("=", 1) + metadata[key] = item_value + return metadata + + +def _variant_type(reference: str, observed: str) -> str: + if len(reference) == len(observed) == 1: + return "substitution" + if len(reference) < len(observed): + return "insertion" + if len(reference) > len(observed): + return "deletion" + return "replacement" + + +def read_vcf(path: PathLike) -> Iterator[Variant]: + """Yield normalized first-ALT variants while preserving VCF metadata.""" + + with _open_text(path) as handle: + for line_number, line in enumerate(handle, start=1): + if not line.strip() or line.startswith("#"): + continue + fields = line.rstrip("\r\n").split("\t") + if len(fields) < 8: + raise ValueError(f"VCF record at line {line_number} has fewer than 8 columns") + chromosome, position_text, identifier, reference, alternate_text, quality_text, filter_value, info = fields[:8] + try: + position = int(position_text) + except ValueError as error: + raise ValueError(f"invalid VCF position at line {line_number}") from error + alternatives = [value for value in alternate_text.split(",") if value and value != "."] + if not alternatives: + raise ValueError(f"VCF record at line {line_number} has no ALT allele") + quality = None if quality_text == "." else float(quality_text) + metadata = _parse_info(info) + metadata.update( + { + "id": identifier, + "filter": filter_value, + "alternatives": alternatives, + "coordinate_system": "1-based", + } + ) + yield Variant( + chromosome=chromosome, + position=position, + reference=reference.upper(), + observed=alternatives[0].upper(), + quality=quality, + variant_type=_variant_type(reference, alternatives[0]), + metadata=metadata, + ) + + +def read_raw_sequence(path: PathLike, identifier: str | None = None) -> Iterator[SequenceRecord]: + with _open_text(path) as handle: + sequence = normalize_sequence(handle.read()) + yield SequenceRecord( + identifier=identifier or Path(path).stem, + sequence=sequence, + source_format="text", + ) + + +def read_sequence_file(path: PathLike, format: str = "auto") -> Iterator[SequenceRecord]: + path = Path(path) + selected = format.lower() + if selected == "auto": + suffixes = [suffix.lower() for suffix in path.suffixes] + if ".fastq" in suffixes or ".fq" in suffixes: + selected = "fastq" + elif ".fasta" in suffixes or ".fa" in suffixes or ".fna" in suffixes: + selected = "fasta" + else: + selected = "text" + if selected in {"fasta", "fa", "fna"}: + yield from read_fasta(path) + elif selected in {"fastq", "fq"}: + yield from read_fastq(path) + elif selected in {"text", "raw"}: + yield from read_raw_sequence(path) + else: + raise ValueError(f"unsupported sequence format: {format}") + + +@dataclass(frozen=True) +class InputData: + records: tuple[SequenceRecord, ...] = () + variants: tuple[Variant, ...] = () + + +def read_input(path: PathLike, format: str = "auto") -> InputData: + path = Path(path) + selected = format.lower() + if selected == "auto" and path.suffix.lower() in {".vcf", ".gz"}: + selected = "vcf" if path.suffix.lower() == ".vcf" else ( + "vcf" if path.name.lower().endswith(".vcf.gz") else "auto" + ) + if selected == "vcf": + return InputData(variants=tuple(read_vcf(path))) + return InputData(records=tuple(read_sequence_file(path, format=selected))) + diff --git a/genopedia/models.py b/genopedia/models.py new file mode 100644 index 0000000..6873d64 --- /dev/null +++ b/genopedia/models.py @@ -0,0 +1,117 @@ +"""Optional-model-compatible primitives implemented with the standard library.""" + +from __future__ import annotations + +import json +import math +from collections import Counter +from dataclasses import dataclass +from typing import Iterable + + +@dataclass(frozen=True) +class KmerConfig: + k: int = 3 + + +def kmer_counts(sequence: str, k: int = 3) -> dict[str, int]: + if k < 1: + raise ValueError("k must be at least 1") + normalized = "".join(sequence.split()).upper() + return dict(Counter(normalized[index : index + k] for index in range(len(normalized) - k + 1))) + + +class KmerClassifier: + """A deterministic multinomial k-mer classifier with Laplace smoothing.""" + + def __init__(self, k: int = 3) -> None: + if k < 1: + raise ValueError("k must be at least 1") + self.k = k + self.classes: tuple[str, ...] = () + self._counts: dict[str, Counter[str]] = {} + self._totals: dict[str, int] = {} + self._vocabulary: tuple[str, ...] = () + self._sample_counts: Counter[str] = Counter() + + def fit(self, sequences: Iterable[str], labels: Iterable[str]) -> "KmerClassifier": + sequence_list = list(sequences) + label_list = list(labels) + if not sequence_list or not label_list or len(sequence_list) != len(label_list): + raise ValueError("sequences and labels must be non-empty and have equal lengths") + self.classes = tuple(sorted(set(label_list))) + self._counts = {label: Counter() for label in self.classes} + self._totals = {label: 0 for label in self.classes} + self._sample_counts = Counter(label_list) + vocabulary: set[str] = set() + for sequence, label in zip(sequence_list, label_list): + counts = kmer_counts(sequence, self.k) + self._counts[label].update(counts) + self._totals[label] += sum(counts.values()) + vocabulary.update(counts) + self._vocabulary = tuple(sorted(vocabulary)) + return self + + def _log_scores(self, sequence: str) -> dict[str, float]: + if not self.classes: + raise RuntimeError("classifier must be trained before prediction") + counts = kmer_counts(sequence, self.k) + vocabulary_size = max(1, len(self._vocabulary)) + total_samples = sum(self._sample_counts.values()) + scores: dict[str, float] = {} + for label in self.classes: + prior = self._sample_counts[label] / total_samples + denominator = self._totals[label] + vocabulary_size + score = math.log(prior) + for kmer, count in counts.items(): + probability = (self._counts[label].get(kmer, 0) + 1) / denominator + score += count * math.log(probability) + scores[label] = score + return scores + + def predict_proba(self, sequence: str) -> dict[str, float]: + scores = self._log_scores(sequence) + maximum = max(scores.values()) + exponentials = {label: math.exp(score - maximum) for label, score in scores.items()} + total = sum(exponentials.values()) + return {label: value / total for label, value in exponentials.items()} + + def predict(self, sequence: str) -> str: + probabilities = self.predict_proba(sequence) + return max(probabilities, key=probabilities.get) + + def to_dict(self) -> dict[str, object]: + return { + "k": self.k, + "classes": list(self.classes), + "counts": {label: dict(counts) for label, counts in self._counts.items()}, + "totals": self._totals, + "vocabulary": list(self._vocabulary), + "sample_counts": dict(self._sample_counts), + } + + def to_json(self) -> str: + return json.dumps(self.to_dict(), sort_keys=True) + + @classmethod + def from_dict(cls, payload: dict[str, object]) -> "KmerClassifier": + classifier = cls(k=int(payload["k"])) + classifier.classes = tuple(str(value) for value in payload["classes"]) + classifier._counts = { + str(label): Counter({str(kmer): int(count) for kmer, count in values.items()}) + for label, values in payload["counts"].items() + } + classifier._totals = {str(label): int(value) for label, value in payload["totals"].items()} + classifier._vocabulary = tuple(str(value) for value in payload["vocabulary"]) + classifier._sample_counts = Counter( + {str(label): int(value) for label, value in payload["sample_counts"].items()} + ) + return classifier + + @classmethod + def from_json(cls, value: str) -> "KmerClassifier": + return cls.from_dict(json.loads(value)) + + +DNAClassifier = KmerClassifier + diff --git a/genopedia/reporting.py b/genopedia/reporting.py new file mode 100644 index 0000000..a49bb11 --- /dev/null +++ b/genopedia/reporting.py @@ -0,0 +1,151 @@ +"""Self-contained HTML and SVG reporting for local/offline use.""" + +from __future__ import annotations + +import html +from pathlib import Path +from typing import Iterable + +from .core import QualityReport, SequenceRecord, Variant + + +BASE_COLORS = { + "A": ("#d62828", "Adenine"), + "T": ("#e9c46a", "Thymine"), + "G": ("#2a9d8f", "Guanine"), + "C": ("#457b9d", "Cytosine"), + "U": ("#8d5524", "Uracil"), + "N": ("#6c757d", "Ambiguous/unknown"), +} + + +def render_svg_sequence(sequence: str, width: int = 800, height: int = 80, max_bases: int = 240) -> str: + sequence = "".join(sequence.split()).upper() + visible = sequence[:max_bases] + if not visible: + return f'' + base_width = width / len(visible) + parts = [ + f'' + ] + for index, base in enumerate(visible): + color = BASE_COLORS.get(base, BASE_COLORS["N"])[0] + x = round(index * base_width, 3) + parts.append( + f'{index}: {html.escape(base)}' + ) + if len(sequence) > max_bases: + parts.append(f'โ€ฆ {len(sequence) - max_bases} more bases') + parts.append("") + return "".join(parts) + + +def _render_sequence_spans(sequence: str, max_bases: int = 1000) -> str: + sequence = "".join(sequence.split()).upper() + spans: list[str] = [] + for index, base in enumerate(sequence[:max_bases]): + safe_base = html.escape(base) + spans.append( + f'{safe_base}' + ) + if len(sequence) > max_bases: + spans.append(f'โ€ฆ {len(sequence) - max_bases} more bases') + return "".join(spans) + + +def _quality_section(reports: Iterable[QualityReport]) -> str: + rows = [] + for report in reports: + rows.append( + "" + f"{html.escape(report.identifier)}" + f"{report.length}" + f"{report.gc_fraction:.3f}" + f"{report.ambiguous_count}" + f"{html.escape(report.status)}" + f"{html.escape('; '.join(report.warnings) or 'None')}" + "" + ) + if not rows: + return "

No quality-control results.

" + return ( + '' + '' + + "".join(rows) + + "
SampleLengthGC fractionAmbiguousStatusWarnings
" + ) + + +def _variant_section(variants: Iterable[Variant]) -> str: + rows = [] + for variant in variants: + rows.append( + "" + f"{html.escape(variant.chromosome)}" + f"{variant.position}" + f"{html.escape(variant.reference or 'โˆ…')}" + f"{html.escape(variant.observed or 'โˆ…')}" + f"{html.escape(variant.variant_type)}" + f"{html.escape(variant.interpretation)}" + "" + ) + if not rows: + return "

No observed variants.

" + return ( + '' + '' + + "".join(rows) + + "
ChromosomePositionReferenceObservedTypeInterpretation
" + ) + + +def render_html_report( + records: Iterable[SequenceRecord], + variants: Iterable[Variant] = (), + quality_reports: Iterable[QualityReport] = (), + title: str = "Genopedia report", +) -> str: + records = list(records) + variants = list(variants) + quality_reports = list(quality_reports) + legend = "".join( + f'{name} ({base})' + for base, (color, name) in BASE_COLORS.items() + ) + sequence_sections = [] + for record in records: + sequence_sections.append( + f'

{html.escape(record.identifier)}

' + f'

{html.escape(record.description)}

' + f'
{_render_sequence_spans(record.sequence)}
' + f'
{render_svg_sequence(record.sequence)}
' + ) + if not sequence_sections: + sequence_sections.append("

No sequence records were supplied.

") + return f""" + +{html.escape(title)} + +

{html.escape(title)}

Offline, research-oriented sequence review. Interpretations require appropriate evidence and human review.

+

Color legend

{legend}
+

Sequences

{''.join(sequence_sections)}
+

Quality control

{_quality_section(quality_reports)}
+

Observed variants

{_variant_section(variants)}
+""" + + +def write_html_report(path: str | Path, html_text: str) -> Path: + output = Path(path) + output.parent.mkdir(parents=True, exist_ok=True) + output.write_text(html_text, encoding="utf-8") + return output + diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..721fb7d --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,28 @@ +[build-system] +requires = ["setuptools>=68"] +build-backend = "setuptools.build_meta" + +[project] +name = "genopedia" +version = "0.2.0" +description = "Offline-first DNA/RNA sequence inspection and reporting toolkit" +readme = "README.md" +requires-python = ">=3.10" +license = { file = "LICENSE" } +authors = [{ name = "Genopedia Project" }] +dependencies = [] + +[project.optional-dependencies] +bio = ["biopython>=1.85,<2"] +ml = ["numpy>=1.26", "scikit-learn>=1.4"] +dev = ["pytest>=8", "black>=24", "flake8>=7"] + +[project.scripts] +genopedia = "genopedia.cli:main" + +[tool.setuptools.packages.find] +include = ["genopedia*"] + +[tool.black] +line-length = 100 + diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 0000000..e51f8e6 --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1,4 @@ +pytest>=8.0 +black>=24.0 +flake8>=7.0 + diff --git a/requirements.txt b/requirements.txt index a6971e0..b9f4687 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,22 +1,5 @@ -biopython==1.85 -numpy==1.26.4 -pandas==2.2.2 -scikit-learn==1.6.1 -tensorflow==2.17.0 -torch==2.4.1 -matplotlib==3.9.2 -seaborn==0.13.2 -plotly==5.24.1 -fastapi==0.115.0 -uvicorn==0.30.6 -sqlalchemy==2.0.34 -neo4j==5.25.0 -jupyter==1.1.1 -scipy==1.14.1 -tqdm==4.66.5 -pytest==8.3.3 -black==24.8.0 -flake8==7.1.1 -wandb==0.18.0 -gradio==5.0.1 -flask==2.2.5 \ No newline at end of file +# Genopedia core intentionally has no third-party runtime dependencies. +# Optional integrations are declared in pyproject.toml: +# python -m pip install -e .[bio] +# python -m pip install -e .[ml] +# python -m pip install -e .[dev] diff --git a/run_genopedia.py b/run_genopedia.py new file mode 100644 index 0000000..bfadcef --- /dev/null +++ b/run_genopedia.py @@ -0,0 +1,8 @@ +"""Convenience launcher for a fresh checkout.""" + +from genopedia.cli import main + + +if __name__ == "__main__": + raise SystemExit(main()) + diff --git a/src/genomics/__init__.py b/src/genomics/__init__.py deleted file mode 100644 index 849ef90..0000000 --- a/src/genomics/__init__.py +++ /dev/null @@ -1,275 +0,0 @@ -""" -Genomics Data Pipeline (Genopedia) -================================== -Handles DNA/RNA sequence loading, preprocessing, variant detection, -and color-coded visualization. -""" - -import os -import re -import json -import hashlib -import logging -from typing import Dict, List, Tuple, Optional, Union -from dataclasses import dataclass, field -from collections import Counter -import numpy as np - -# Bioinformatics -from Bio import SeqIO -from Bio.Seq import Seq -from Bio.SeqRecord import SeqRecord -from Bio.SeqUtils import GC - -# Setup logging -logging.basicConfig(level=logging.INFO) -logger = logging.getLogger(__name__) - - -@dataclass -class NucleotideConfig: - """Configuration for nucleotide color coding.""" - A_COLOR = "#FF0000" # Red - T_COLOR = "#FFFF00" # Yellow - G_COLOR = "#00FF00" # Green - C_COLOR = "#0000FF" # Blue - U_COLOR = "#8B4513" # Brown (for RNA) - N_COLOR = "#808080" # Gray (unknown) - - COLOR_MAP = { - 'A': A_COLOR, 'a': A_COLOR, - 'T': T_COLOR, 't': T_COLOR, - 'G': G_COLOR, 'g': G_COLOR, - 'C': C_COLOR, 'c': C_COLOR, - 'U': U_COLOR, 'u': U_COLOR, - 'N': N_COLOR, 'n': N_COLOR, - } - - -@dataclass -class Variant: - """Represents a genetic variant.""" - chromosome: str - position: int - reference: str - observed: str - quality: float = 0.0 - variant_type: str = "SNP" # SNP, Insertion, Deletion - pathogenicity: str = "unknown" # benign, likely_benign, uncertain, likely_pathogenic, pathogenic - - def to_dict(self): - return { - "chromosome": self.chromosome, - "position": self.position, - "reference": self.reference, - "observed": self.observed, - "quality": self.quality, - "variant_type": self.variant_type, - "pathogenicity": self.pathogenicity - } - - -@dataclass -class FunctionalRegion: - """Represents a functional region of a gene.""" - name: str - start: int - end: int - region_type: str # exon, intron, promoter, enhancer, utr - - -class GenomicsDataLoader: - """Handles loading and preprocessing of genomic data.""" - - def __init__(self, data_dir: str = "data/genomics"): - self.data_dir = data_dir - self.sequences = {} - self.variants = [] - self.annotations = {} - logger.info(f"GenomicsDataLoader initialized with data_dir: {data_dir}") - - def load_fasta(self, filepath: str, label: str = None) -> SeqRecord: - """Load a FASTA file.""" - logger.info(f"Loading FASTA: {filepath}") - records = list(SeqIO.parse(filepath, "fasta")) - if label: - self.sequences[label] = records - logger.info(f"Loaded {len(records)} records from {filepath}") - return records - - def load_fastq(self, filepath: str) -> List[SeqRecord]: - """Load FASTQ sequencing reads.""" - logger.info(f"Loading FASTQ: {filepath}") - records = list(SeqIO.parse(filepath, "fastq")) - logger.info(f"Loaded {len(records)} reads from {filepath}") - return records - - def load_vcf(self, filepath: str) -> List[Variant]: - """Load VCF variant file.""" - logger.info(f"Loading VCF: {filepath}") - variants = [] - with open(filepath, 'r') as f: - for line in f: - if line.startswith('#') or line.startswith('##'): - continue - parts = line.strip().split('\t') - if len(parts) >= 5: - variant = Variant( - chromosome=parts[0], - position=int(parts[1]), - reference=parts[3], - observed=parts[4], - quality=float(parts[5]) if len(parts) > 5 else 0.0 - ) - variants.append(variant) - logger.info(f"Loaded {len(variants)} variants from {filepath}") - self.variants.extend(variants) - return variants - - def generate_synthetic_dna(self, length: int = 1000, label: str = "synthetic") -> str: - """Generate a synthetic DNA sequence for testing.""" - bases = ['A', 'T', 'G', 'C'] - sequence = ''.join(np.random.choice(bases, size=length)) - logger.info(f"Generated synthetic DNA sequence of length {length}") - return sequence - - def generate_synthetic_rna(self, length: int = 1000, label: str = "synthetic_rna") -> str: - """Generate a synthetic RNA sequence (with Uracil instead of Thymine).""" - bases = ['A', 'U', 'G', 'C'] - sequence = ''.join(np.random.choice(bases, size=length)) - logger.info(f"Generated synthetic RNA sequence of length {length}") - return sequence - - -class SequenceAnalyzer: - """Analyzes DNA/RNA sequences.""" - - def __init__(self): - self.config = NucleotideConfig() - - def get_color_map(self, sequence: str) -> List[str]: - """Get color mapping for a sequence.""" - return [self.config.COLOR_MAP.get(base, self.config.N_COLOR) for base in sequence] - - def calculate_gc_content(self, sequence: str) -> float: - """Calculate GC content of a sequence.""" - return GC(Seq(sequence)) - - def find_motifs(self, sequence: str, motif: str) -> List[int]: - """Find all occurrences of a motif in a sequence.""" - positions = [] - for i in range(len(sequence) - len(motif) + 1): - if sequence[i:i+len(motif)] == motif: - positions.append(i) - return positions - - def identify_functional_regions(self, sequence: str) -> List[FunctionalRegion]: - """Identify potential functional regions in a sequence.""" - regions = [] - - # Look for promoter regions (TATA box) - tata_positions = self.find_motifs(sequence, "TATA") - for pos in tata_positions: - regions.append(FunctionalRegion( - name=f"Promoter_{pos}", - start=max(0, pos - 25), - end=min(len(sequence), pos + 25), - region_type="promoter" - )) - - # Look for start codon (ATG) - start_positions = self.find_motifs(sequence, "ATG") - for pos in start_positions: - regions.append(FunctionalRegion( - name=f"Start_Codon_{pos}", - start=pos, - end=pos + 3, - region_type="exon" - )) - - return regions - - def detect_variants(self, reference: str, sample: str, chromosome: str = "chr1") -> List[Variant]: - """Detect variants between reference and sample sequences.""" - variants = [] - min_len = min(len(reference), len(sample)) - - for i in range(min_len): - if reference[i] != sample[i]: - variants.append(Variant( - chromosome=chromosome, - position=i, - reference=reference[i], - observed=sample[i], - variant_type="SNP" - )) - - # Check for insertions/deletions at the end - if len(sample) > len(reference): - variants.append(Variant( - chromosome=chromosome, - position=min_len, - reference="-", - observed=sample[min_len:], - variant_type="Insertion" - )) - elif len(reference) > len(sample): - variants.append(Variant( - chromosome=chromosome, - position=min_len, - reference=reference[min_len:], - observed="-", - variant_type="Deletion" - )) - - return variants - - def classify_pathogenicity(self, variant: Variant) -> str: - """Classify variant pathogenicity (simplified).""" - # This is a simplified classification - # In practice, this would use databases like ClinVar - if variant.observed in ['A', 'G', 'C', 'T'] and variant.reference in ['A', 'G', 'C', 'T']: - # Simple heuristic: transitions (A<->G, C<->T) are often less pathogenic - transitions = {('A', 'G'), ('G', 'A'), ('C', 'T'), ('T', 'C')} - if (variant.reference, variant.observed) in transitions: - return "likely_benign" - else: - return "uncertain" - return "unknown" - - -class GenomicsVisualizer: - """Creates color-coded visualizations of genomic data.""" - - def __init__(self): - self.config = NucleotideConfig() - - def generate_color_html(self, sequence: str, max_length: int = 100) -> str: - """Generate HTML representation of a sequence with color-coded nucleotides.""" - html = '
' - - for i, base in enumerate(sequence[:max_length]): - color = self.config.COLOR_MAP.get(base, self.config.N_COLOR) - html += f'{base}' - if (i + 1) % 10 == 0: - html += ' ' - if (i + 1) % 60 == 0: - html += '
' - - html += '
' - return html - - def generate_svg_sequence(self, sequence: str, width: int = 800, height: int = 100) -> str: - """Generate SVG visualization of a sequence.""" - svg = f'\n' - - base_width = width / min(len(sequence), 100) - - for i, base in enumerate(sequence[:100]): - color = self.config.COLOR_MAP.get(base, self.config.N_COLOR) - x = i * base_width - svg += f' \n' - - svg += '' - return svg \ No newline at end of file diff --git a/src/models/__init__.py b/src/models/__init__.py deleted file mode 100644 index 88a1d9d..0000000 --- a/src/models/__init__.py +++ /dev/null @@ -1,70 +0,0 @@ -""" -Genomics ML Models -================== -Predefined architectures for genomics tasks. -""" -from __future__ import annotations - -from dataclasses import dataclass -from typing import List, Sequence - -try: - import numpy as np -except Exception: # pragma: no cover - optional - np = None # type: ignore - - -@dataclass -class KmerConfig: - k: int = 3 - alphabet: Sequence[str] = ("A", "C", "G", "T") - - -def kmer_counts(sequence: str, k: int = 3) -> dict[str, int]: - if np is not None: - raise RuntimeError("numpy is required for counting.") - counts: dict[str, int] = {} - for i in range(len(sequence) - k + 1): - kmer = sequence[i : i + k] - counts[kmer] = counts.get(kmer, 0) + 1 - return counts - - -class DNAClassifier: - def __init__(self, k: int = 3): - self.k = k - self.probabilities: dict[str, dict[str, float]] = {} - - def fit(self, sequences: List[str], labels: List[str]): - classes = sorted({label for label in labels}) - kmer_class_counts: dict[str, dict[str, int]] = {c: {} for c in classes} - class_totals = {c: 0 for c in classes} - for sequence, label in zip(sequences, labels): - counts = kmer_counts(sequence, self.k) - for kmer, count in counts.items(): - kmer_class_counts[label][kmer] = kmer_class_counts[label].get(kmer, 0) + count - class_totals[label] += count - self.probabilities = {} - for c in classes: - vocabulary = set(kmer_class_counts[c]) - self.probabilities[c] = { - kmer: (kmer_class_counts[c].get(kmer, 0) + 1) / (class_totals[c] + len(vocabulary)) - for kmer in vocabulary - } - self.classes = classes - - def predict(self, sequence: str) -> str: - if not getattr(self, "classes", None): - raise RuntimeError("Classifier must be trained before prediction.") - counts = kmer_counts(sequence, self.k) - scores = {} - for c in self.classes: - log_prob = 0.0 - vocabulary = set(self.probabilities[c]) - for kmer, count in counts.items(): - prob = self.probabilities[c].get(kmer, 1 / (sum(self.probabilities[c].values()) + len(vocabulary))) - log_prob += count * (np.log(prob) if np is not None else __import__('math').log(prob)) - for kmer in (set(counts) - vocabulary): - log_prob += counts[kmer] * __import__('math').log(1e-9) - scores[c] = log_prob - return max(scores, key=scores.get) diff --git a/tests/test_adapters.py b/tests/test_adapters.py new file mode 100644 index 0000000..08aa965 --- /dev/null +++ b/tests/test_adapters.py @@ -0,0 +1,29 @@ +from __future__ import annotations + +import tempfile +import unittest +from pathlib import Path + +from genopedia.adapters import FileDropAdapter + + +class FileDropAdapterTests(unittest.TestCase): + def test_adapter_discovers_supported_exports_in_stable_order(self) -> None: + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + (root / "b.fastq").write_text("@b\nAT\n+\nII\n", encoding="utf-8") + (root / "a.fasta").write_text(">a\nAT\n", encoding="utf-8") + (root / "notes.txt").write_text("not a sequence", encoding="utf-8") + + discovered = FileDropAdapter(root).discover() + + self.assertEqual([path.name for path in discovered], ["a.fasta", "b.fastq"]) + + def test_adapter_reads_an_export_without_mutating_the_source(self) -> None: + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "sample.fasta" + path.write_text(">sample\nATGC\n", encoding="utf-8") + + data = FileDropAdapter(directory).read(path) + + self.assertEqual(data.records[0].sequence, "ATGC") diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 0000000..be263bb --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,75 @@ +from __future__ import annotations + +import json +import subprocess +import sys +import tempfile +import unittest +from pathlib import Path + +from genopedia.cli import main + + +class CliTests(unittest.TestCase): + def test_example_script_runs_from_repository_root(self) -> None: + repository_root = Path(__file__).resolve().parents[1] + + result = subprocess.run( + [sys.executable, "examples/demo_genomics.py"], + cwd=repository_root, + capture_output=True, + text=True, + ) + + self.assertEqual(result.returncode, 0, result.stderr) + + def test_demo_command_writes_a_report_and_json_summary(self) -> None: + with tempfile.TemporaryDirectory() as directory: + output = Path(directory) / "report.html" + summary = Path(directory) / "report.json" + + exit_code = main( + [ + "demo", + "--length", + "24", + "--output", + str(output), + "--json-output", + str(summary), + ] + ) + + self.assertEqual(exit_code, 0) + self.assertTrue(output.exists()) + self.assertTrue(summary.exists()) + payload = json.loads(summary.read_text(encoding="utf-8")) + + self.assertEqual(payload["records"][0]["length"], 24) + + def test_compare_command_writes_detected_variants(self) -> None: + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + reference = root / "reference.txt" + sample = root / "sample.txt" + output = root / "comparison.html" + summary = root / "comparison.json" + reference.write_text("ACGT", encoding="utf-8") + sample.write_text("ACGGT", encoding="utf-8") + + exit_code = main( + [ + "compare", + str(reference), + str(sample), + "--output", + str(output), + "--json-output", + str(summary), + ] + ) + + self.assertEqual(exit_code, 0) + payload = json.loads(summary.read_text(encoding="utf-8")) + + self.assertEqual(payload["variants"][0]["variant_type"], "insertion") diff --git a/tests/test_genomics.py b/tests/test_genomics.py index d16c1a6..ff10e60 100644 --- a/tests/test_genomics.py +++ b/tests/test_genomics.py @@ -1,47 +1,60 @@ -"""GenoProject tests.""" from __future__ import annotations -import pytest +import unittest -from GenoProject.src.genomics import GenomicsDataLoader, SequenceAnalyzer, GenomicsVisualizer +from genopedia import SequenceAnalyzer, SequenceRecord, Variant +from genopedia.models import KmerClassifier -@pytest.fixture -def loader(): - return GenomicsDataLoader() +class SequenceAnalyzerTests(unittest.TestCase): + def test_gc_content_is_a_fraction_and_ignores_ambiguous_bases(self) -> None: + analyzer = SequenceAnalyzer() + self.assertEqual(analyzer.gc_content("ACTG"), 0.5) + self.assertEqual(analyzer.gc_content("ACNG"), 2 / 3) + self.assertEqual(analyzer.gc_content(""), 0.0) -@pytest.fixture -def analyzer(): - return SequenceAnalyzer() + def test_quality_control_reports_invalid_symbols_without_mutating_input(self) -> None: + analyzer = SequenceAnalyzer() + record = SequenceRecord(identifier="sample-1", sequence="ACGT?N") + report = analyzer.quality_control(record) -@pytest.fixture -def visualizer(): - return GenomicsVisualizer() + self.assertEqual(record.sequence, "ACGT?N") + self.assertEqual(report.length, 6) + self.assertEqual(report.invalid_symbols, ("?",)) + self.assertEqual(report.ambiguous_count, 1) + self.assertEqual(report.status, "fail") + def test_motif_search_finds_overlapping_matches(self) -> None: + positions = SequenceAnalyzer().find_motifs("ATATAT", "ATA") -def test_synthetic_dna_generation(loader): - sequence = loader.generate_synthetic_dna(120) - assert len(sequence) == 120 - assert set(sequence).issubset({"A", "T", "G", "C"}) + self.assertEqual(positions, [0, 2]) + def test_sequence_comparison_reports_an_insertion_at_reference_position(self) -> None: + variants = SequenceAnalyzer().detect_variants("ACGT", "ACGGT", chromosome="chr7") -def test_detects_simple_variant(analyzer): - reference = "ATGCCGTAG" - sample = "ATGTCGTAG" - variants = analyzer.detect_variants(reference, sample) - assert len(variants) == 1 - assert variants[0].position == 4 - assert variants[0].reference == "C" - assert variants[0].observed == "T" + self.assertEqual(len(variants), 1) + self.assertEqual(variants[0].chromosome, "chr7") + self.assertEqual(variants[0].position, 3) + self.assertEqual(variants[0].reference, "") + self.assertEqual(variants[0].observed, "G") + self.assertEqual(variants[0].variant_type, "insertion") + def test_default_variant_interpretation_is_unknown_without_evidence(self) -> None: + variant = Variant(chromosome="chr1", position=10, reference="A", observed="G") -def test_color_map_returns_expected_colors(analyzer): - colors = analyzer.get_color_map("ATGCUN") - assert colors[0] == "#FF0000" - assert colors[1] == "#FFFF00" - assert colors[2] == "#00FF00" - assert colors[3] == "#0000FF" - assert colors[4] == "#8B4513" - assert colors[5] == "#808080" + self.assertEqual(SequenceAnalyzer().interpret_variant(variant), "unknown") + + +class KmerClassifierTests(unittest.TestCase): + def test_classifier_can_fit_and_predict_without_numpy(self) -> None: + classifier = KmerClassifier(k=2) + classifier.fit(["AAAAAA", "AAAATA", "CCCCCC", "CCCCAC"], ["A", "A", "C", "C"]) + + self.assertEqual(classifier.predict("AAAAAA"), "A") + self.assertEqual(classifier.predict("CCCCCC"), "C") + + def test_classifier_rejects_prediction_before_training(self) -> None: + with self.assertRaises(RuntimeError): + KmerClassifier(k=3).predict("ATGC") diff --git a/tests/test_io.py b/tests/test_io.py new file mode 100644 index 0000000..3284c5b --- /dev/null +++ b/tests/test_io.py @@ -0,0 +1,57 @@ +from __future__ import annotations + +import tempfile +import unittest +from pathlib import Path + +from genopedia.io import read_fasta, read_fastq, read_vcf, read_sequence_file + + +class SequenceIoTests(unittest.TestCase): + def test_fasta_reader_streams_records_and_normalizes_sequence(self) -> None: + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "sample.fa" + path.write_text(">alpha first record\nat gc\n>beta\nTTAA\n", encoding="utf-8") + + records = list(read_fasta(path)) + + self.assertEqual([record.identifier for record in records], ["alpha", "beta"]) + self.assertEqual(records[0].sequence, "ATGC") + self.assertEqual(records[0].description, "first record") + + def test_fastq_reader_decodes_phred_quality_scores(self) -> None: + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "sample.fastq" + path.write_text("@read-1\nACGT\n+\nIIII\n", encoding="utf-8") + + records = list(read_fastq(path)) + + self.assertEqual(records[0].quality_scores, (40, 40, 40, 40)) + + def test_vcf_reader_accepts_missing_quality_and_preserves_alternatives(self) -> None: + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "sample.vcf" + path.write_text( + "##fileformat=VCFv4.3\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr2\t8\t.\tA\tG,T\t.\tPASS\tDP=12\n", + encoding="utf-8", + ) + + variants = list(read_vcf(path)) + + self.assertEqual(len(variants), 1) + self.assertEqual(variants[0].position, 8) + self.assertEqual(variants[0].observed, "G") + self.assertEqual(variants[0].metadata["alternatives"], ["G", "T"]) + self.assertIsNone(variants[0].quality) + + def test_auto_reader_uses_file_extension(self) -> None: + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "sample.fa" + path.write_text(">alpha\nATGC\n", encoding="utf-8") + + records = list(read_sequence_file(path)) + + self.assertEqual(records[0].sequence, "ATGC") + diff --git a/tests/test_reporting.py b/tests/test_reporting.py new file mode 100644 index 0000000..a2dae9e --- /dev/null +++ b/tests/test_reporting.py @@ -0,0 +1,28 @@ +from __future__ import annotations + +import unittest + +from genopedia import SequenceRecord, Variant +from genopedia.reporting import render_html_report, render_svg_sequence + + +class ReportingTests(unittest.TestCase): + def test_html_report_contains_escaped_sequence_and_analysis_sections(self) -> None: + html = render_html_report( + [SequenceRecord(identifier="sample", sequence="ATGC")], + [Variant(chromosome="chr1", position=2, reference="G", observed="A")], + ) + + self.assertIn("Genopedia report", html) + self.assertIn("sample", html) + self.assertIn("Adenine", html) + self.assertIn("Observed variants", html) + self.assertIn('class="base base-A"', html) + + def test_svg_renderer_returns_empty_svg_for_empty_sequence(self) -> None: + svg = render_svg_sequence("") + + self.assertIn("", svg) + self.assertNotIn("NaN", svg) +