Skip to content

ARCHITECTURE.md

Latest commit

 

History

History
407 lines (295 loc) · 14.8 KB

ARCHITECTURE.md

File metadata and controls

407 lines (295 loc) · 14.8 KB
Raw

Architecture and Design

This document describes the technical design decisions for kmer-diff.

System Overview

┌─────────────────────────────────────────────────────────────────┐
│                           CLI Layer                              │
│  ┌─────────┐  ┌─────────────┐  ┌───────────────────┐            │
│  │  build  │  │   compare   │  │    compare-all    │            │
│  └────┬────┘  └──────┬──────┘  └─────────┬─────────┘            │
└───────┼──────────────┼───────────────────┼──────────────────────┘
        │              │                   │
┌───────▼──────────────▼───────────────────▼──────────────────────┐
│                        Core Library                              │
│  ┌──────────┐  ┌────────────┐  ┌────────────┐  ┌─────────────┐  │
│  │ FastqPar-│  │ KmerCompar-│  │ RegionMer- │  │ Reference   │  │
│  │ ser      │  │ ator       │  │ ger        │  │ Mapper      │  │
│  └──────────┘  └────────────┘  └────────────┘  └─────────────┘  │
│  ┌──────────┐  ┌────────────┐  ┌────────────┐  ┌─────────────┐  │
│  │ KmerCoun-│  │ Coverage   │  │ Distance   │  │ Output      │  │
│  │ ter      │  │ Estimator  │  │ Calculator │  │ Writer      │  │
│  └──────────┘  └────────────┘  └────────────┘  └─────────────┘  │
└─────────────────────────────────────────────────────────────────┘
        │                              │
┌───────▼──────────────────────────────▼──────────────────────────┐
│                      Data Layer                                  │
│  ┌──────────────┐  ┌──────────────────────────────────────────┐ │
│  │ KmerEncoder  │  │            KmerDatabase                   │ │
│  │              │  │  ┌────────────┐  ┌─────────────────────┐ │ │
│  │ encode()     │  │  │   Header   │  │   Sorted K-mer      │ │ │
│  │ decode()     │  │  │            │  │   Array             │ │ │
│  │ canonicalize │  │  └────────────┘  └─────────────────────┘ │ │
│  └──────────────┘  └──────────────────────────────────────────┘ │
└─────────────────────────────────────────────────────────────────┘

K-mer Encoding

Bit Representation

Each nucleotide is encoded as 2 bits:

Base Encoding
A 00
C 01
G 10
T 11

A k-mer is stored as a 64-bit unsigned integer with the first base in the most significant position:

K-mer: ACGT (k=4)
Encoding: 00 01 10 11 = 0x1B (27 decimal)

Position:  [63:62] [61:60] [59:58] [57:56] ... [7:6] [5:4] [3:2] [1:0]
           unused  unused  unused  unused      base  base  base  base
                                               k-3   k-2   k-1    k

For k=31, this uses 62 bits, leaving 2 bits unused.

Canonicalization

Since DNA is double-stranded and reads can come from either strand, we store the canonical form of each k-mer: the lexicographically smaller of the k-mer and its reverse complement.

def canonicalize(kmer: int, k: int) -> int:
    rc = reverse_complement(kmer, k)
    return min(kmer, rc)

Reverse complement algorithm:

  1. Swap all bits (XOR with all 1s) to complement bases
  2. Reverse the order of 2-bit pairs
  3. Shift right to align (if k < 32)

Maximum K-mer Size

With 64-bit integers: k_max = 32

We default to k=31 because:

  • Large enough to be unique in bacterial genomes (~5 Mb)
  • Odd length breaks palindrome symmetry (canonical form is always well-defined)
  • Leaves 2 unused bits that could be repurposed if needed

Database File Format

Design Goals

  • Fast sequential read (for merge-join comparison)
  • Compact size
  • Simple to implement
  • Self-describing (header contains parameters)

Format Specification

Offset  Size    Field           Description
------  ----    -----           -----------
0       4       magic           "KMDB" (0x4B4D4442)
4       2       version         Format version (currently 1)
6       2       k               K-mer size
8       2       min_count       Minimum count threshold used
10      6       reserved        Reserved for future use
16      8       num_kmers       Number of k-mers in database
24      8       checksum        CRC64 of body

32      N*10    body            Sorted (kmer, count) pairs
                                - kmer: uint64, little-endian
                                - count: uint16, little-endian

Why Sorted Array?

Comparison requires finding shared and unique k-mers between two samples. Options:

Approach Memory Comparison Time Simplicity
Hash table High O(n) Medium
Sorted array + merge Low O(n + m) High
Bloom filter Very low O(n) Medium (approximate)

Sorted array wins because:

  • Merge-join is optimal O(n + m) for comparing two sorted lists
  • No hash table overhead in memory
  • Simple implementation
  • Sequential disk access is cache-friendly

Count Storage

We use 16-bit counts (max 65535). For 100x coverage of a bacterial genome:

  • Expected count for single-copy k-mers: ~100
  • Maximum for high-copy repeats: rarely exceeds 10000
  • 16 bits is sufficient with large margin

Normalization Strategy

The Problem

Raw k-mer counts aren't comparable between samples:

  • Different total sequencing depth
  • Different DNA extraction efficiency
  • Different library prep

A k-mer with count 50 could be:

  • 1 copy at 50x coverage
  • 2 copies at 25x coverage
  • 0.5 copies at 100x coverage (sequencing from one haplotype)

Solution: Estimate 1x Coverage

We estimate the sequencing depth of single-copy chromosomal DNA, then express all counts as copy numbers relative to this baseline.

Method 1: Median of distribution (default)

Most k-mers in a bacterial genome are single-copy. The median of the count distribution approximates 1x coverage.

def estimate_coverage_median(counts: np.ndarray) -> float:
    return np.median(counts)

Caveats:

  • Inflated if genome has many repeats
  • Can be affected by contamination or mixed samples

Method 2: Core gene k-mers (optional, more robust)

If user provides a FASTA of known single-copy core genes:

  1. Extract all k-mers from core gene sequences
  2. Look up counts in sample's database
  3. Use median of found k-mers as 1x estimate
def estimate_coverage_core_genes(database: KmerDatabase, 
                                  core_genes_path: str) -> Tuple[float, float]:
    core_kmers = extract_kmers(core_genes_path, database.k)
    found_counts = []
    for kmer in core_kmers:
        count = database.get_count(kmer)
        if count > 0:
            found_counts.append(count)
    
    fraction_found = len(found_counts) / len(core_kmers)
    coverage = np.median(found_counts) if found_counts else 0
    
    return coverage, fraction_found

The fraction_found is a QC metric:

  • 90%: good match between reference and sample

  • 70-90%: sample may have diverged or have quality issues
  • <70%: possible wrong species or severe contamination

Comparison Algorithm

Merge-Join

Given two sorted k-mer arrays, walk through both simultaneously:

def compare(db_a: KmerDatabase, db_b: KmerDatabase, 
            cov_a: float, cov_b: float, threshold: float) -> ComparisonResult:
    
    i, j = 0, 0
    differential = []
    unique_a = []
    unique_b = []
    
    while i < len(db_a) and j < len(db_b):
        kmer_a, count_a = db_a.kmers[i], db_a.counts[i]
        kmer_b, count_b = db_b.kmers[j], db_b.counts[j]
        
        if kmer_a < kmer_b:
            unique_a.append(kmer_a)
            i += 1
        elif kmer_a > kmer_b:
            unique_b.append(kmer_b)
            j += 1
        else:  # kmer_a == kmer_b
            copy_a = count_a / cov_a
            copy_b = count_b / cov_b
            if abs(copy_a - copy_b) > threshold:
                differential.append((kmer_a, copy_a, copy_b))
            i += 1
            j += 1
    
    # Handle remaining elements
    unique_a.extend(db_a.kmers[i:])
    unique_b.extend(db_b.kmers[j:])
    
    return ComparisonResult(differential, unique_a, unique_b)

Time complexity: O(n + m) where n, m are database sizes.

Differential Threshold

We flag k-mers where estimated copy number differs by more than 0.7 (default).

Why 0.7?

  • Midpoint between 0 and 1 copy difference
  • Tolerates noise around single-copy regions
  • Catches clear duplications/deletions

The threshold is configurable because optimal value depends on:

  • Sequencing depth (higher depth = tighter threshold possible)
  • Biological question (strict vs exploratory)

Region Merging

Adjacency Definition

Two k-mers are adjacent if they could come from consecutive positions in a sequence, i.e., one's (k-1) suffix equals the other's (k-1) prefix.

K-mer A: ACGTACGTACGTACGTACGTACGTACGTACG  (k=31)
K-mer B:  CGTACGTACGTACGTACGTACGTACGTACGT

A's suffix (k-1 = 30): CGTACGTACGTACGTACGTACGTACGTACG
B's prefix (k-1 = 30): CGTACGTACGTACGTACGTACGTACGTACG
Match! → A and B are adjacent

Efficient test using bit operations:

def are_adjacent(kmer_a: int, kmer_b: int, k: int) -> bool:
    mask = (1 << (2 * (k - 1))) - 1  # Mask for k-1 bases
    suffix_a = kmer_a & mask          # Remove first base of A
    prefix_b = kmer_b >> 2            # Remove last base of B
    return suffix_a == prefix_b

Graph Construction

Build a directed graph where:

  • Nodes = differential (or unique) k-mers
  • Edges = adjacency relationships

For a clean genomic region without variants, this forms a linear chain.

Connected Components

Find connected components using union-find or DFS. Each component represents a contiguous differential region.

Linear Chain Detection

A component is a linear chain if:

  • Each node has at most 1 incoming edge and 1 outgoing edge
  • No cycles (guaranteed for genomic sequence)

Branching occurs when:

  • SNPs within the differential region create alternate paths
  • Repeat structures cause convergence/divergence

Sequence Reconstruction

For linear chains, walk from start to end, extending the sequence by one base per k-mer:

def reconstruct(chain: List[int], k: int) -> str:
    sequence = decode(chain[0], k)  # First k-mer gives k bases
    for kmer in chain[1:]:
        last_base = decode(kmer & 0b11)  # Last 2 bits = last base
        sequence += last_base
    return sequence

Handling Branching

For branching components, we could:

  1. Report all paths (combinatorial explosion)
  2. Report longest path
  3. Flag for manual inspection

We chose option 3: flag for manual inspection. Branching is rare in closely related strains and usually indicates interesting biology worth examining.

Distance Metric

Definition

Distance between samples A and B:

distance = (total bases in differential regions) / (megabases of shared genome)

Where:

  • Total differential bases = sum of lengths of all differential regions
  • Shared genome = number of shared k-mers × 1 (each k-mer represents ~1 unique base position on average)

Rationale

  • Bases, not k-mers: More interpretable (e.g., "150 bp of differences")
  • Normalized: Accounts for variable data quality and genome coverage
  • Per megabase: Familiar unit, comparable to mutation rates

Example

Sample A and B share 4.5 million k-mers. They have 3 differential regions totaling 2500 bp.

distance = 2500 / 4.5 = 556 differential bases per Mb

This could be interpreted as ~556 SNP-equivalents per megabase, which is high for very recent divergence (suggests either significant evolutionary distance or structural variation).

Memory Considerations

Build Phase

During counting, we maintain a hash table of all k-mers:

Component Size Estimate
K-mer (key) 8 bytes
Count (value) 8 bytes (Python int overhead)
Hash table overhead ~2x

For 100x coverage of 5 Mb genome:

  • Unique k-mers: ~5 million (plus some errors)
  • Memory: ~5M × 16 bytes × 2 = ~160 MB

With error k-mers (before filtering):

  • At 0.1% error rate, ~5x more unique k-mers
  • Memory: ~800 MB

This fits comfortably in modern machines. For very high coverage or larger genomes, could implement disk-based counting.

Compare Phase

Loading two databases:

  • Each: N k-mers × 10 bytes
  • For 5M k-mers: ~50 MB each
  • Total: ~100 MB

Merge-join is streaming, adds minimal overhead.

Error Handling

Sequencing Errors

At 0.1% per-base error rate, most errors create unique k-mers (seen only once). Our min_count filter (default 3) removes these.

Expected survival of error k-mers:

  • Single error: count = 1 → filtered
  • Same error twice in overlapping reads: rare
  • Same error in non-overlapping reads: extremely rare

Coverage Dropout

Some true genomic k-mers may fall below threshold due to random sampling. At 100x expected coverage:

  • Poisson probability of count < 3: ~0.01%
  • Negligible impact on results

Contamination

Cross-contamination between samples could:

  • Create false "shared" k-mers
  • Reduce apparent differences

No automated detection; users should ensure sample purity.