Implementation Plan

This document describes the implementation strategy for kmer-diff, organized into development phases.

Phase 1: Test Data Generator

Goal: Create synthetic test data with known variants to validate the tool.

Components

1.1 Reference Genome Generator

Generate a synthetic 100,000 bp bacterial genome:

50% GC content (realistic for many bacteria)
Random sequence with no extreme repetitive regions
Define 20 "core gene" regions of 500-1000 bp each
Ensure core genes are distributed across the genome

Output: reference.fasta

1.2 Core Genes Extractor

Extract the 20 core gene regions as a separate FASTA:

Each gene as a separate entry
Named systematically (e.g., core_gene_01, core_gene_02, ...)

Output: core_genes.fasta

1.3 Variant Generator

Create modified genomes with known variants:

Sample B variants:

10 SNPs at random positions (outside core genes)
1 duplicated region: 2,000 bp (appears twice in sample B, once in reference)
1 deleted region: 1,000 bp (absent in sample B, present in reference)

Sample C variants (optional):

Different set of 10 SNPs
Different structural variants

Output: sample_b_genome.fasta, variants.tsv (ground truth)

1.4 Read Simulator

Generate paired-end Illumina-like reads:

Read length: 150 bp
Insert size: 300-500 bp (normal distribution, mean 400, sd 50)
Coverage: 100x
Realistic error model:
- Position-dependent error rate (increases toward read end)
- Base error rate: ~0.1% at start, ~0.5% at end
- Mostly substitution errors
- Quality scores that correlate with error probability

Output: sample_a_R1.fastq.gz, sample_a_R2.fastq.gz, etc.

Deliverables

scripts/
  generate_test_data.py    # Main script

test_data/
  reference.fasta          # 100kb synthetic genome
  core_genes.fasta         # 20 single-copy genes
  sample_a_R1.fastq.gz     # Baseline reads
  sample_a_R2.fastq.gz
  sample_b_R1.fastq.gz     # Reads with SNPs + dup + del
  sample_b_R2.fastq.gz
  sample_c_R1.fastq.gz     # Different variant set
  sample_c_R2.fastq.gz
  variants.tsv             # Ground truth positions and types

Phase 2: Core Data Structures

Goal: Implement k-mer encoding and the .kdb file format.

Components

2.1 K-mer Encoding

Implement efficient k-mer representation:

class KmerEncoder:
    def __init__(self, k: int = 31):
        self.k = k
    
    def encode(self, sequence: str) -> int:
        """Convert k-mer string to 64-bit integer (2 bits per base)."""
        
    def decode(self, kmer_int: int) -> str:
        """Convert 64-bit integer back to k-mer string."""
        
    def canonicalize(self, kmer_int: int) -> int:
        """Return lexicographically smaller of k-mer and reverse complement."""
        
    def reverse_complement(self, kmer_int: int) -> int:
        """Compute reverse complement of encoded k-mer."""

Encoding scheme:

A = 00, C = 01, G = 10, T = 11
K-mer stored in least significant 2k bits
For k=31, uses 62 of 64 bits

2.2 K-mer Database Format

Binary .kdb file structure:

Header (32 bytes):
  - Magic number: 4 bytes ("KMDB")
  - Version: 2 bytes (1)
  - K-mer size: 2 bytes
  - Min count threshold: 2 bytes
  - Reserved: 6 bytes
  - Number of k-mers: 8 bytes
  - Checksum: 8 bytes

Body:
  - Sorted array of (kmer: uint64, count: uint16) pairs
  - Total size: num_kmers * 10 bytes

class KmerDatabase:
    def __init__(self, k: int, min_count: int):
        self.k = k
        self.min_count = min_count
        self.kmers: np.ndarray  # sorted uint64
        self.counts: np.ndarray  # uint16
    
    @classmethod
    def from_file(cls, path: str) -> 'KmerDatabase':
        """Load database from .kdb file."""
        
    def to_file(self, path: str) -> None:
        """Save database to .kdb file."""
        
    def __len__(self) -> int:
        """Number of k-mers in database."""
        
    def get_count(self, kmer: int) -> int:
        """Lookup count for a k-mer (binary search)."""

Deliverables

kmer_diff/
  encoding.py              # KmerEncoder class
  database.py              # KmerDatabase class
  
tests/
  test_encoding.py
  test_database.py

Phase 3: Build Command

Goal: Implement FASTQ parsing and k-mer counting.

Components

3.1 FASTQ Parser

Parse gzipped paired-end FASTQ files:

class FastqParser:
    def __init__(self, r1_path: str, r2_path: str):
        self.r1_path = r1_path
        self.r2_path = r2_path
    
    def __iter__(self) -> Iterator[Tuple[str, str, str, str]]:
        """Yield (seq1, qual1, seq2, qual2) for each read pair."""

Requirements:

Handle gzipped files transparently
Validate FASTQ format
Memory-efficient streaming

3.2 K-mer Counter

Count k-mers in memory using hash table:

class KmerCounter:
    def __init__(self, k: int = 31, min_count: int = 3):
        self.k = k
        self.min_count = min_count
        self.encoder = KmerEncoder(k)
        self.counts: Dict[int, int] = {}
    
    def add_sequence(self, sequence: str) -> None:
        """Extract and count all k-mers from a sequence."""
        
    def add_fastq(self, r1_path: str, r2_path: str) -> None:
        """Count k-mers from paired FASTQ files."""
        
    def to_database(self) -> KmerDatabase:
        """Convert to sorted KmerDatabase, filtering by min_count."""

3.3 CLI Entry Point

@click.command()
@click.option('--r1', required=True, help='Forward reads (R1) FASTQ file')
@click.option('--r2', required=True, help='Reverse reads (R2) FASTQ file')
@click.option('-o', '--output', required=True, help='Output .kdb file')
@click.option('--kmer-size', default=31, help='K-mer size (default: 31)')
@click.option('--min-count', default=3, help='Minimum count threshold (default: 3)')
def build(r1, r2, output, kmer_size, min_count):
    """Build k-mer database from paired-end FASTQ files."""

Deliverables

kmer_diff/
  fastq.py                 # FastqParser class
  counter.py               # KmerCounter class
  cli.py                   # CLI entry points
  
tests/
  test_fastq.py
  test_counter.py
  test_build_command.py

Phase 4: Compare Command

Goal: Implement pairwise comparison with normalization and region merging.

Components

4.1 Coverage Estimator

Estimate 1x chromosomal coverage:

class CoverageEstimator:
    def __init__(self, database: KmerDatabase):
        self.database = database
    
    def estimate_median(self) -> float:
        """Estimate 1x coverage from median k-mer count."""
        
    def estimate_from_core_genes(self, core_genes_path: str, k: int) -> Tuple[float, float]:
        """
        Estimate 1x coverage from core gene k-mers.
        Returns (coverage_estimate, fraction_found).
        """

4.2 K-mer Comparator

Compare two databases using merge-join:

class KmerComparator:
    def __init__(self, db_a: KmerDatabase, db_b: KmerDatabase, 
                 cov_a: float, cov_b: float, threshold: float = 0.7):
        self.db_a = db_a
        self.db_b = db_b
        self.cov_a = cov_a
        self.cov_b = cov_b
        self.threshold = threshold
    
    def compare(self) -> ComparisonResult:
        """
        Walk through sorted k-mer lists and categorize:
        - differential: shared k-mers with |copy_a - copy_b| > threshold
        - unique_to_a: k-mers only in A
        - unique_to_b: k-mers only in B
        """

@dataclass
class ComparisonResult:
    differential: List[Tuple[int, float, float]]  # (kmer, copy_a, copy_b)
    unique_to_a: List[int]  # kmers
    unique_to_b: List[int]  # kmers
    shared_count: int
    total_a: int
    total_b: int

4.3 Region Merger

Merge adjacent k-mers into regions:

class RegionMerger:
    def __init__(self, k: int):
        self.k = k
        self.encoder = KmerEncoder(k)
    
    def merge(self, kmers: List[int]) -> List[MergedRegion]:
        """
        Build adjacency graph and find connected components.
        Return linear chains and flag branching components.
        """
    
    def _are_adjacent(self, kmer_a: int, kmer_b: int) -> bool:
        """Check if two k-mers overlap by k-1 bases."""
        
    def _reconstruct_sequence(self, chain: List[int]) -> str:
        """Walk linear chain to reconstruct sequence."""

@dataclass
class MergedRegion:
    sequence: str
    kmer_count: int
    mean_copy_a: float
    mean_copy_b: float
    is_branching: bool

4.4 Reference Mapper (Optional)

Map regions to reference coordinates:

class ReferenceMapper:
    def __init__(self, reference_path: str):
        self.reference = self._load_reference(reference_path)
    
    def map_region(self, sequence: str) -> Optional[GenomicLocation]:
        """Find location of sequence in reference genome."""

@dataclass
class GenomicLocation:
    contig: str
    start: int
    end: int
    strand: str

4.5 Output Writer

Generate output files:

class CompareOutputWriter:
    def __init__(self, output_dir: str):
        self.output_dir = output_dir
    
    def write_summary(self, result: ComparisonResult, ...) -> None:
    def write_differential_regions(self, regions: List[MergedRegion]) -> None:
    def write_unique_regions(self, regions: List[MergedRegion], sample: str) -> None:
    def write_branching_flags(self, regions: List[MergedRegion]) -> None:
    def write_bed(self, regions: List[MergedRegion], locations: List[GenomicLocation]) -> None:

Deliverables

kmer_diff/
  coverage.py              # CoverageEstimator
  comparator.py            # KmerComparator
  merger.py                # RegionMerger
  mapper.py                # ReferenceMapper
  output.py                # Output writers
  
tests/
  test_coverage.py
  test_comparator.py
  test_merger.py
  test_mapper.py
  test_compare_command.py

Phase 5: Compare-All Command

Goal: Implement all-vs-all comparison with distance matrix.

Components

5.1 Distance Calculator

Compute pairwise distances:

class DistanceCalculator:
    def __init__(self, threshold: float = 0.7):
        self.threshold = threshold
    
    def calculate(self, db_a: KmerDatabase, db_b: KmerDatabase,
                  cov_a: float, cov_b: float) -> float:
        """
        Calculate distance as:
        total_differential_bases / shared_megabases
        """

5.2 Matrix Builder

Build and output distance matrix:

class MatrixBuilder:
    def __init__(self, databases: Dict[str, KmerDatabase],
                 coverages: Dict[str, float]):
        self.databases = databases
        self.coverages = coverages
    
    def build_matrix(self) -> pd.DataFrame:
        """Compute all pairwise distances."""
        
    def write_matrix(self, path: str) -> None:
        """Write TSV distance matrix."""

5.3 CLI Integration

@click.command()
@click.argument('databases', nargs=-1, required=True)
@click.option('-o', '--output', required=True, help='Output directory')
@click.option('--threshold', default=0.7, help='Copy number difference threshold')
@click.option('--core-genes', help='Core genes FASTA for normalization')
@click.option('--reference', help='Reference genome for coordinate mapping')
@click.option('--full-pairwise', is_flag=True, help='Generate full pairwise reports')
def compare_all(databases, output, threshold, core_genes, reference, full_pairwise):
    """Compare all samples pairwise and generate distance matrix."""

Deliverables

kmer_diff/
  distance.py              # DistanceCalculator
  matrix.py                # MatrixBuilder
  
tests/
  test_distance.py
  test_matrix.py
  test_compare_all_command.py

Phase 6: Integration and Polish

Goal: End-to-end testing, documentation, and packaging.

Components

6.1 Integration Tests

Run full pipeline on test data
Verify detected variants match ground truth
Test edge cases (empty files, single sample, etc.)

6.2 Performance Testing

Profile memory usage on realistic data
Identify bottlenecks
Document resource requirements

6.3 Documentation

Complete docstrings
Usage examples
Troubleshooting guide

6.4 Packaging

setup.py / pyproject.toml
Entry point configuration
Dependencies specification

Deliverables

tests/
  test_integration.py      # End-to-end tests
  
docs/
  usage.md                 # Detailed usage guide
  troubleshooting.md       # Common issues
  
pyproject.toml             # Package configuration

Development Order

Recommended implementation sequence:

Phase 1: Test data generator (enables testing from the start)
Phase 2: Core data structures (foundation for everything)
Phase 3: Build command (first usable feature)
Phase 4: Compare command (core functionality)
Phase 5: Compare-all command (extends compare)
Phase 6: Integration and polish

Each phase should include tests before moving to the next.

Uh oh!

FilesExpand file tree

IMPLEMENTATION.md

Latest commit

History

IMPLEMENTATION.md

File metadata and controls

Implementation Plan

Phase 1: Test Data Generator

Components

1.1 Reference Genome Generator

1.2 Core Genes Extractor

1.3 Variant Generator

1.4 Read Simulator

Deliverables

Phase 2: Core Data Structures

Components

2.1 K-mer Encoding

2.2 K-mer Database Format

Deliverables

Phase 3: Build Command

Components

3.1 FASTQ Parser

3.2 K-mer Counter

3.3 CLI Entry Point

Deliverables

Phase 4: Compare Command

Components

4.1 Coverage Estimator

4.2 K-mer Comparator

4.3 Region Merger

4.4 Reference Mapper (Optional)

4.5 Output Writer

Deliverables

Phase 5: Compare-All Command

Components

5.1 Distance Calculator

5.2 Matrix Builder

5.3 CLI Integration

Deliverables

Phase 6: Integration and Polish

Components

6.1 Integration Tests

6.2 Performance Testing

6.3 Documentation

6.4 Packaging

Deliverables

Development Order