diff --git a/src/build.rs b/src/build.rs
index 828b46b..d0fda0d 100644
--- a/src/build.rs
+++ b/src/build.rs
@@ -632,7 +632,146 @@ fn gap_open_aligner(reference: &str, sequence: &str) -> String {
     let cigar = alignment_to_cigar(&alignment.operations);
     // println!("{:?}", cigar);
 
-    cigar
+    // A single gap-open (-5) is pricier than a single mismatch (-1), so inside
+    // a homopolymer/tandem-repeat run the optimal-scoring alignment above will
+    // sometimes fold part of a real indel into an adjacent base as a spurious
+    // substitution instead of a clean insertion/deletion (e.g. chrM:310, a
+    // lone T flanked by poly-C, gets reported as a false near-homoplasmic
+    // T>C "SNP" alongside a same-run C insertion). That substitution is
+    // score-optimal but not the minimal/correct edit - normalize it away.
+    normalize_homopolymer_indels(&cigar, reference, sequence)
+}
+
+/// Rewrites a CIGAR to eliminate a length-1 `X` sitting directly against an
+/// `I` or `D` run whenever the X's "orphaned" base (the one with no partner
+/// on the gapped side) also occurs inside that run. In that case the run can
+/// be losslessly re-anchored so the X becomes a plain match and the indel is
+/// redistributed around it: same total ref/alt content, no substitution.
+/// Left untouched when the orphaned base doesn't occur in the adjacent run -
+/// that's very likely a genuine substitution, not a homopolymer-registration
+/// artifact.
+fn normalize_homopolymer_indels(cigar: &str, ref_seq: &str, alt_seq: &str) -> String {
+    let ops = parse_cigar(cigar);
+    let ref_bytes = ref_seq.as_bytes();
+    let alt_bytes = alt_seq.as_bytes();
+
+    let mut result: Vec<(usize, char)> = Vec::new();
+    let mut ref_pos = 0usize;
+    let mut alt_pos = 0usize;
+    let mut i = 0usize;
+
+    let push_op = |result: &mut Vec<(usize, char)>, len: usize, op: char| {
+        if len == 0 {
+            return;
+        }
+        if let Some(last) = result.last_mut() {
+            if last.1 == op {
+                last.0 += len;
+                return;
+            }
+        }
+        result.push((len, op));
+    };
+
+    while i < ops.len() {
+        let (len, op) = ops[i];
+
+        // I(n) immediately followed by a single mismatch: the mismatch's ref
+        // base may really belong inside the insertion run.
+        if op == 'I' && i + 1 < ops.len() && ops[i + 1] == (1, 'X') {
+            let ins = &alt_bytes[alt_pos..alt_pos + len];
+            let ref_x = ref_bytes[ref_pos];
+            if let Some(p) = ins.iter().rposition(|&b| b == ref_x) {
+                push_op(&mut result, p, 'I');
+                push_op(&mut result, 1, '=');
+                push_op(&mut result, len - p, 'I'); // trailing insertion piece + the X's alt base
+                ref_pos += 1;
+                alt_pos += len + 1;
+                i += 2;
+                continue;
+            }
+        }
+
+        // A single mismatch immediately followed by I(n): mirror of the above.
+        if op == 'X' && len == 1 && i + 1 < ops.len() && ops[i + 1].1 == 'I' {
+            let n = ops[i + 1].0;
+            let ins = &alt_bytes[alt_pos + 1..alt_pos + 1 + n];
+            let ref_x = ref_bytes[ref_pos];
+            if let Some(p) = ins.iter().position(|&b| b == ref_x) {
+                push_op(&mut result, 1 + p, 'I'); // the X's alt base + leading insertion piece
+                push_op(&mut result, 1, '=');
+                push_op(&mut result, n - p - 1, 'I');
+                ref_pos += 1;
+                alt_pos += 1 + n;
+                i += 2;
+                continue;
+            }
+        }
+
+        // D(n) immediately followed by a single mismatch: mirror of I+X with
+        // ref/alt roles swapped.
+        if op == 'D' && i + 1 < ops.len() && ops[i + 1] == (1, 'X') {
+            let del = &ref_bytes[ref_pos..ref_pos + len];
+            let alt_x = alt_bytes[alt_pos];
+            if let Some(q) = del.iter().rposition(|&b| b == alt_x) {
+                push_op(&mut result, q, 'D');
+                push_op(&mut result, 1, '=');
+                push_op(&mut result, len - q, 'D');
+                ref_pos += len + 1;
+                alt_pos += 1;
+                i += 2;
+                continue;
+            }
+        }
+
+        // A single mismatch immediately followed by D(n): mirror of X+I.
+        if op == 'X' && len == 1 && i + 1 < ops.len() && ops[i + 1].1 == 'D' {
+            let n = ops[i + 1].0;
+            let del = &ref_bytes[ref_pos + 1..ref_pos + 1 + n];
+            let alt_x = alt_bytes[alt_pos];
+            if let Some(q) = del.iter().position(|&b| b == alt_x) {
+                push_op(&mut result, 1 + q, 'D');
+                push_op(&mut result, 1, '=');
+                push_op(&mut result, n - q - 1, 'D');
+                ref_pos += 1 + n;
+                alt_pos += 1;
+                i += 2;
+                continue;
+            }
+        }
+
+        match op {
+            '=' | 'X' => {
+                ref_pos += len;
+                alt_pos += len;
+            }
+            'I' => alt_pos += len,
+            'D' => ref_pos += len,
+            _ => {}
+        }
+        push_op(&mut result, len, op);
+        i += 1;
+    }
+
+    result
+        .iter()
+        .map(|(count, op)| format!("{}{}", count, op))
+        .collect()
+}
+
+/// Parses a run-length CIGAR string ("23=4I1X20=") into `(length, op)` pairs.
+fn parse_cigar(cigar: &str) -> Vec<(usize, char)> {
+    let mut ops = Vec::new();
+    let mut num = String::new();
+    for c in cigar.chars() {
+        if c.is_ascii_digit() {
+            num.push(c);
+        } else {
+            ops.push((num.parse().expect("cigar length"), c));
+            num.clear();
+        }
+    }
+    ops
 }
 
 pub fn generate_cigar(
@@ -861,9 +1000,128 @@ pub fn start(output: &PathBuf, k: usize, read_path: &PathBuf, reference_path: &P
     );
     let tmp_gfa_path = PathBuf::from("tmp.gfa");
     let mut graph = agg::GraphicalGenome::load_graph(&tmp_gfa_path).unwrap();
-    
+
     // generate cigar
     generate_cigar(&mut graph, &ref_header, k, maxlength, 2);
     let graph_output = output.with_extension("gfa");
     let _ = write_graph_from_graph(graph_output.to_str().unwrap(), &graph);
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Reconstructs the alt sequence a CIGAR implies, given the same ref_seq
+    /// and alt_seq the CIGAR was computed from, so tests can prove a
+    /// rewritten CIGAR still encodes exactly the same edit as the original
+    /// rather than merely "having no X".
+    fn reconstruct_alt(cigar: &str, _ref_seq: &str, alt_seq: &str) -> String {
+        let alt_bytes = alt_seq.as_bytes();
+        let mut ref_pos = 0usize;
+        let mut alt_pos = 0usize;
+        let mut out = Vec::new();
+        let mut num = String::new();
+        for c in cigar.chars() {
+            if c.is_ascii_digit() {
+                num.push(c);
+                continue;
+            }
+            let len: usize = num.parse().unwrap();
+            num.clear();
+            match c {
+                '=' | 'X' => {
+                    out.extend_from_slice(&alt_bytes[alt_pos..alt_pos + len]);
+                    ref_pos += len;
+                    alt_pos += len;
+                }
+                'I' => {
+                    out.extend_from_slice(&alt_bytes[alt_pos..alt_pos + len]);
+                    alt_pos += len;
+                }
+                'D' => {
+                    ref_pos += len;
+                }
+                _ => {}
+            }
+        }
+        String::from_utf8(out).unwrap()
+    }
+
+    fn ref_consumed(cigar: &str) -> usize {
+        let mut total = 0usize;
+        let mut num = String::new();
+        for c in cigar.chars() {
+            if c.is_ascii_digit() {
+                num.push(c);
+                continue;
+            }
+            let len: usize = num.parse().unwrap();
+            num.clear();
+            if matches!(c, '=' | 'X' | 'D') {
+                total += len;
+            }
+        }
+        total
+    }
+
+    #[test]
+    fn homopolymer_insertion_absorbs_adjacent_mismatch() {
+        // Real chrM 287-330 window from an NA12877 poly-C tract (7 C's, a
+        // lone T at what would be chrM:310, then 5 more C's). The sample's
+        // edge inserts 4 extra C's in that stretch; the raw aligner's optimal
+        // scoring path represents 3 of them as a clean "I" but folds the 4th
+        // together with the flanking T into a spurious "1X" (T>C) - a false
+        // homoplasmic SNP that doesn't exist in any real read.
+        let ref_seq = "AAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGC";
+        let alt_seq = "AAAAATTTCCACCAAACCCCCCCCCCTCCCCCCGCTTCTGGCCACAGC";
+        let raw_cigar = "23=4I1X20=";
+        assert_eq!(reconstruct_alt(raw_cigar, ref_seq, alt_seq), alt_seq);
+
+        let normalized = normalize_homopolymer_indels(raw_cigar, ref_seq, alt_seq);
+
+        assert!(
+            !normalized.contains('X'),
+            "expected the homopolymer-adjacent mismatch to be absorbed, got {}",
+            normalized
+        );
+        assert_eq!(ref_consumed(&normalized), ref_seq.len());
+        assert_eq!(reconstruct_alt(&normalized, ref_seq, alt_seq), alt_seq);
+    }
+
+    #[test]
+    fn genuine_snp_next_to_unrelated_indel_is_left_alone() {
+        // A real substitution (ref 'G', alt 'A') sitting right next to an
+        // insertion whose inserted base never equals the mismatched ref base
+        // - there is no way to re-anchor this without changing content, so
+        // it must be left as a real X.
+        let ref_seq = "ACGTACGTG";
+        let alt_seq = "ACGTACGTTTA"; // insert "TT" then mismatch G->A
+        let raw_cigar = "8=2I1X";
+        assert_eq!(reconstruct_alt(raw_cigar, ref_seq, alt_seq), alt_seq);
+
+        let normalized = normalize_homopolymer_indels(raw_cigar, ref_seq, alt_seq);
+
+        assert_eq!(normalized, raw_cigar);
+    }
+
+    #[test]
+    fn homopolymer_deletion_absorbs_adjacent_mismatch() {
+        // Exact role-swap (ref<->alt, I<->D) of homopolymer_insertion_absorbs_adjacent_mismatch:
+        // a deletion is just an insertion viewed from the other sequence, so
+        // this is the same real chrM poly-C edit, mirrored.
+        let ref_seq = "AAAAATTTCCACCAAACCCCCCCCCCTCCCCCCGCTTCTGGCCACAGC";
+        let alt_seq = "AAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGC";
+        let raw_cigar = "23=4D1X20=";
+        assert_eq!(reconstruct_alt(raw_cigar, ref_seq, alt_seq), alt_seq);
+
+        let normalized = normalize_homopolymer_indels(raw_cigar, ref_seq, alt_seq);
+
+        assert!(
+            !normalized.contains('X'),
+            "expected the homopolymer-adjacent mismatch to be absorbed, got {}",
+            normalized
+        );
+        assert_eq!(ref_consumed(&normalized), ref_seq.len());
+        assert_eq!(reconstruct_alt(&normalized, ref_seq, alt_seq), alt_seq);
+    }
+}
diff --git a/src/call.rs b/src/call.rs
index b0afb90..5a2a4e7 100644
--- a/src/call.rs
+++ b/src/call.rs
@@ -207,6 +207,18 @@ fn edge_reads(graph: &GraphicalGenome, edge: &str) -> HashSet<String> {
 fn get_graph_intervals(graph:&GraphicalGenome, length: i64) -> HashMap<&String, (i64, i64)>{
     let mut graph_intervals_dict = HashMap::new();
     for edge in graph.edges.keys() {
+        // Edges with no computed variant record were never CIGAR-processed
+        // (get_variant skips them via this same empty-cigar check) - a read
+        // reaching a bubble only through such an edge was never genotyped
+        // there, and must not be credited as "covering" it, or it later
+        // defaults to a false ref (0) call instead of the correct NaN.
+        let has_variant = graph.edges[edge]
+            .get("variants")
+            .and_then(|v| v.as_str())
+            .map_or(false, |s| !s.is_empty());
+        if !has_variant {
+            continue;
+        }
         let src = graph.edges[edge].get("src").unwrap().as_array().unwrap()[0].as_str().unwrap();
         let dst = graph.edges[edge].get("dst").unwrap().as_array().unwrap()[0].as_str().unwrap();
         let startpos = graph.anchor
@@ -240,6 +252,34 @@ fn bubble_cover_reads(graph_intervals_dict: &HashMap<&String, (i64, i64)>, pos:
     cover
 }
 
+/// Populates `read_record`/`cover_record` for one edge's variants.
+///
+/// `graph_intervals_dict` holds raw (unwrapped) anchor coordinates, including
+/// the duplicated wrap-around span past `ref_length` used to circularize the
+/// graph. `variants` carries that same raw position, while `variants_circular`
+/// carries the wrapped position used to key/name the variant. The bubble
+/// lookup must use the raw position - using the wrapped one would look up
+/// the wrong bubble (or none) for any variant on the wrap-around edge.
+fn record_variant_reads(
+    variants: &[Variant],
+    variants_circular: &[Variant],
+    graph_intervals_dict: &HashMap<&String, (i64, i64)>,
+    graph: &GraphicalGenome,
+    readlist: &[serde_json::Value],
+    read_record: &mut HashMap<String, Vec<serde_json::Value>>,
+    cover_record: &mut HashMap<String, HashSet<String>>,
+) {
+    for (raw_v, v) in variants.iter().zip(variants_circular.iter()) {
+        let key = generate_variant_name(&v.clone());
+        let cover_reads = bubble_cover_reads(&graph_intervals_dict, raw_v.pos, graph);
+        read_record
+            .entry(key.clone())
+            .or_insert_with(Vec::new)
+            .extend(readlist.to_vec());
+        cover_record.entry(key).or_insert_with(HashSet::new).extend(cover_reads);
+    }
+}
+
 pub fn get_variant(
     graph: &mut GraphicalGenome,
     k: usize,
@@ -335,19 +375,15 @@ pub fn get_variant(
                 reads.as_array().unwrap_or(&Vec::new()).to_vec()
             });
         
-        for v in &variants_circular {
-            let key = generate_variant_name(&v.clone());
-            let cover_reads = bubble_cover_reads(&graph_intervals_dict, v.pos, graph);
-            read_record
-                .entry(key.clone())
-                .or_insert_with(Vec::new)
-                .extend(readlist.clone());
-            // let entry = cover_record.entry(key).or_default();
-            cover_record.entry(key).or_insert_with(HashSet::new).extend(cover_reads.clone());
-            // for read in &cover_reads {
-            //     entry.insert(read.clone());
-            // }
-        }
+        record_variant_reads(
+            &variants,
+            &variants_circular,
+            &graph_intervals_dict,
+            graph,
+            &readlist,
+            &mut read_record,
+            &mut cover_record,
+        );
     }
     (var, coverage, read_record, cover_record)
 }
@@ -1171,4 +1207,136 @@ pub fn start(
     // write matrix
     let matrix_output = output_file.with_extension("matrix.csv");
     let _ = write_matrix_to_csv(&filtered_matrix, &permu_filtered_var, &read_set, matrix_output);
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use serde_json::json;
+
+    /// Two bubbles: a normal one spanning raw positions [0, 5), and a
+    /// wrap-around bubble spanning raw positions [8, 12) - i.e. it crosses
+    /// the circular origin for a 10bp reference.
+    fn wraparound_test_graph() -> GraphicalGenome {
+        let mut anchor = HashMap::new();
+        anchor.insert("A1".to_string(), json!({"pos": 0}));
+        anchor.insert("A2".to_string(), json!({"pos": 5}));
+        anchor.insert("A3".to_string(), json!({"pos": 8}));
+        anchor.insert("A4".to_string(), json!({"pos": 12}));
+
+        let mut edges = HashMap::new();
+        edges.insert(
+            "E1".to_string(),
+            json!({"src": ["A1"], "dst": ["A2"], "reads": ["r1", "r2"], "variants": "5="}),
+        );
+        edges.insert(
+            "E2".to_string(),
+            json!({"src": ["A1"], "dst": ["A2"], "reads": ["r3"], "variants": "2=1X2="}),
+        );
+        edges.insert(
+            "E3".to_string(),
+            json!({"src": ["A3"], "dst": ["A4"], "reads": ["r4"], "variants": "4="}),
+        );
+        edges.insert(
+            "E4".to_string(),
+            json!({"src": ["A3"], "dst": ["A4"], "reads": ["r5"], "variants": "1=1X2="}),
+        );
+
+        GraphicalGenome {
+            anchor,
+            edges,
+            outgoing: HashMap::new(),
+            incoming: HashMap::new(),
+        }
+    }
+
+    #[test]
+    fn record_variant_reads_uses_raw_position_for_wraparound_bubble() {
+        let graph = wraparound_test_graph();
+        let ref_length = 10usize;
+        let graph_intervals_dict = get_graph_intervals(&graph, ref_length as i64);
+
+        // Raw position 11 lies on the wrap-around bubble [8, 12) and
+        // circularizes down to wrapped position 1.
+        let raw_variant = Variant {
+            pos: 11,
+            ref_allele: "A".to_string(),
+            alt_allele: "G".to_string(),
+            variant_type: "SNP".to_string(),
+            allele_count: 1,
+            filter: None,
+        };
+        let variants = vec![raw_variant];
+        let variants_circular = circuliarize_variants(variants.clone(), ref_length);
+        assert_eq!(variants_circular[0].pos, 1, "sanity check: position should wrap");
+
+        let readlist = vec![json!("r4"), json!("r5")];
+        let mut read_record = HashMap::new();
+        let mut cover_record = HashMap::new();
+
+        record_variant_reads(
+            &variants,
+            &variants_circular,
+            &graph_intervals_dict,
+            &graph,
+            &readlist,
+            &mut read_record,
+            &mut cover_record,
+        );
+
+        let key = generate_variant_name(&variants_circular[0]);
+        let mut cover: Vec<String> = cover_record.get(&key).cloned().unwrap_or_default().into_iter().collect();
+        cover.sort();
+
+        // The wrap-around bubble's true reads are r4/r5. Using the wrapped
+        // position (1) instead of the raw one (11) would incorrectly match
+        // the unrelated bubble at [0, 5), pulling in r1/r2/r3 instead.
+        assert_eq!(cover, vec!["r4".to_string(), "r5".to_string()]);
+    }
+
+    /// A bubble at raw positions [0, 10) carries the real variant (edge
+    /// `E_real`, cigar `5=1X4=`), with reads r1/r2. A read that couldn't be
+    /// assembled through the fine-grained anchor chain instead crosses the
+    /// same span via a private "skip-over" edge with no computed CIGAR
+    /// (`E_private`, read r3) - mirroring what Himito's own graph does for
+    /// reads it can't cleanly anchor (see e.g. E00001050 in a real run).
+    fn skip_edge_test_graph() -> GraphicalGenome {
+        let mut anchor = HashMap::new();
+        anchor.insert("A1".to_string(), json!({"pos": 0}));
+        anchor.insert("A2".to_string(), json!({"pos": 10}));
+
+        let mut edges = HashMap::new();
+        edges.insert(
+            "E_real".to_string(),
+            json!({"src": ["A1"], "dst": ["A2"], "reads": ["r1", "r2"], "variants": "5=1X4="}),
+        );
+        edges.insert(
+            "E_private".to_string(),
+            json!({"src": ["A1"], "dst": ["A2"], "reads": ["r3"]}),
+        );
+
+        GraphicalGenome {
+            anchor,
+            edges,
+            outgoing: HashMap::new(),
+            incoming: HashMap::new(),
+        }
+    }
+
+    #[test]
+    fn bubble_cover_reads_excludes_edges_with_no_computed_variant() {
+        let graph = skip_edge_test_graph();
+        let graph_intervals_dict = get_graph_intervals(&graph, 10);
+
+        let cover = bubble_cover_reads(&graph_intervals_dict, 5, &graph);
+        let mut cover: Vec<String> = cover.into_iter().collect();
+        cover.sort();
+
+        // r3 only reaches this span via E_private, which was never CIGAR-
+        // processed (get_variant would skip it via the empty-cigar check).
+        // It must not be credited as "covering" the variant - crediting it
+        // as covered-but-not-alt would wrongly call it ref (0) in the
+        // matrix instead of the correct NaN (unassembled/unknown).
+        assert_eq!(cover, vec!["r1".to_string(), "r2".to_string()]);
+    }
 }
\ No newline at end of file
diff --git a/src/lineage.rs b/src/lineage.rs
new file mode 100644
index 0000000..a05720d
--- /dev/null
+++ b/src/lineage.rs
@@ -0,0 +1,626 @@
+use std::collections::HashMap;
+use anyhow::{bail,Context, Result};
+use rust_htslib::bcf::{Read, Reader}; // Read trait must be in scope for .records()
+use std::io::{BufWriter, Write};
+use std::fs::File;
+
+// ─── Tree data structure ──────────────────────────────────────────────────────
+
+pub struct Node {
+    pub id: usize,
+    pub label: String,
+    pub is_leaf: bool,
+    /// IDs of child nodes (empty for leaves)
+    pub children: Vec<usize>,
+    /// ID of the parent node (None for root)
+    pub parent: Option<usize>,
+    /// Branch length to parent (0.0 for root)
+    pub branch_length: f64,
+    /// Number of reads that map to this leaf (0 for internal nodes)
+    pub read_count: usize,
+}
+
+pub struct Tree {
+    pub nodes: Vec<Node>,
+    pub root: usize,
+}
+
+impl Tree {
+    pub fn leaves(&self) -> Vec<usize> {
+        self.nodes
+            .iter()
+            .filter(|n| n.is_leaf)
+            .map(|n| n.id)
+            .collect()
+    }
+
+    /// Iterative post-order traversal (children before parent).
+    pub fn post_order(&self) -> Vec<usize> {
+        let mut result = Vec::with_capacity(self.nodes.len());
+        let mut stack: Vec<(usize, bool)> = vec![(self.root, false)];
+        while let Some((id, done)) = stack.pop() {
+            if done {
+                result.push(id);
+            } else {
+                stack.push((id, true));
+                for &child in self.nodes[id].children.iter().rev() {
+                    stack.push((child, false));
+                }
+            }
+        }
+        result
+    }
+
+    /// Iterative pre-order traversal (parent before children).
+    pub fn pre_order(&self) -> Vec<usize> {
+        let mut result = Vec::with_capacity(self.nodes.len());
+        let mut stack = vec![self.root];
+        while let Some(id) = stack.pop() {
+            result.push(id);
+            for &child in self.nodes[id].children.iter().rev() {
+                stack.push(child);
+            }
+        }
+        result
+    }
+
+    /// Sum of `read_count` over all leaves in the subtree rooted at `node_id`.
+    pub fn subtree_read_count(&self, node_id: usize) -> usize {
+        let mut total = 0usize;
+        let mut stack = vec![node_id];
+        while let Some(id) = stack.pop() {
+            if self.nodes[id].is_leaf {
+                total += self.nodes[id].read_count;
+            }
+            for &child in &self.nodes[id].children {
+                stack.push(child);
+            }
+        }
+        total
+    }
+}
+
+pub struct DistMatrix {
+    /// Haplotype labels in row/column order (same order as `HaplotypeMatrix::haplotypes`)
+    pub labels: Vec<String>,
+    /// `data[i][j]` = normalised Hamming distance between haplotype i and j.
+    /// Values are in [0.0, 1.0]. Diagonal is 0.0.
+    pub data: Vec<Vec<f64>>,
+}
+
+/// Build a pairwise Hamming distance matrix from haplotype binary profiles.
+///
+/// Distance = (number of differing positions) / (total number of variant positions).
+/// If there are no variants (shouldn't happen after filtering), all distances are 0.
+pub fn hamming_distance_matrix(matrix: &HaplotypeMatrix) -> DistMatrix {
+    let n   = matrix.haplotypes.len();
+    let m   = matrix.variants.len();
+    let inv = if m > 0 { 1.0 / m as f64 } else { 0.0 };
+
+    let labels: Vec<String> = matrix.haplotypes.iter().map(|h| h.id.clone()).collect();
+    let mut data = vec![vec![0.0f64; n]; n];
+
+    for i in 0..n {
+        for j in (i + 1)..n {
+            let diff = matrix.haplotypes[i]
+                .profile
+                .iter()
+                .zip(&matrix.haplotypes[j].profile)
+                .filter(|(&a, &b)| a != b)
+                .count();
+            let d = diff as f64 * inv;
+            data[i][j] = d;
+            data[j][i] = d;
+        }
+    }
+
+    DistMatrix { labels, data }
+}
+
+
+// ─── Neighbor-Joining ─────────────────────────────────────────────────────────
+
+/// Build a rooted Neighbor-Joining tree from `dist`.
+///
+/// Rooting strategy:
+/// * If a haplotype with the all-zero profile exists (no mutations → reference
+///   ancestor), the tree is re-rooted at that haplotype's parent so the
+///   outgroup hangs directly from the root.
+/// * Otherwise the NJ root (midpoint of the final two-taxon join) is kept.
+pub fn neighbor_joining(dist: &DistMatrix, hap_matrix: &HaplotypeMatrix) -> Result<Tree> {
+    let n = dist.labels.len();
+    if n < 2 {
+        bail!("Need at least 2 haplotypes to build a tree, found {n}");
+    }
+
+    // ── Initialise leaf nodes ────────────────────────────────────────────────
+    let mut nodes: Vec<Node> = dist
+        .labels
+        .iter()
+        .enumerate()
+        .map(|(i, label)| Node {
+            id: i,
+            label: label.clone(),
+            is_leaf: true,
+            children: vec![],
+            parent: None,
+            branch_length: 0.0,
+            read_count: hap_matrix.haplotypes[i].count,
+        })
+        .collect();
+
+    // `active[pos]` = node ID; `d[pos_i][pos_j]` = distance between them.
+    let mut active: Vec<usize> = (0..n).collect();
+    let mut d: Vec<Vec<f64>> = dist.data.clone();
+
+    // ── Main NJ loop ─────────────────────────────────────────────────────────
+    while active.len() > 2 {
+        let m = active.len();
+
+        // Row sums (r_i = sum_k d[i][k], d[i][i]=0 so the diagonal is free)
+        let row_sums: Vec<f64> = (0..m)
+            .map(|i| (0..m).map(|j| d[i][j]).sum::<f64>())
+            .collect();
+
+        // Find (best_i, best_j) minimising Q[i][j] = (m-2)*d[i][j] - r_i - r_j
+        let (mut best_i, mut best_j) = (0, 1);
+        let mut best_q = f64::MAX;
+        for i in 0..m {
+            for j in (i + 1)..m {
+                let q = (m as f64 - 2.0) * d[i][j] - row_sums[i] - row_sums[j];
+                if q < best_q {
+                    best_q = q;
+                    best_i = i;
+                    best_j = j;
+                }
+            }
+        }
+
+        // Branch lengths from new internal node u to best_i and best_j
+        let d_ij    = d[best_i][best_j];
+        let delta_i = (d_ij / 2.0
+            + (row_sums[best_i] - row_sums[best_j]) / (2.0 * (m as f64 - 2.0)))
+            .max(0.0);
+        let delta_j = (d_ij - delta_i).max(0.0);
+
+        // Create internal node u
+        let u_id = nodes.len();
+        nodes[active[best_i]].branch_length = delta_i;
+        nodes[active[best_i]].parent        = Some(u_id);
+        nodes[active[best_j]].branch_length = delta_j;
+        nodes[active[best_j]].parent        = Some(u_id);
+        nodes.push(Node {
+            id:            u_id,
+            label:         format!("INT{u_id}"),
+            is_leaf:       false,
+            children:      vec![active[best_i], active[best_j]],
+            parent:        None,
+            branch_length: 0.0,
+            read_count:    0,
+        });
+
+        // New distances: d[u][k] = (d[best_i][k] + d[best_j][k] - d_ij) / 2
+        let remaining: Vec<usize> = (0..m)
+            .filter(|&k| k != best_i && k != best_j)
+            .collect();
+        let new_row: Vec<f64> = remaining
+            .iter()
+            .map(|&k| (0.5 * (d[best_i][k] + d[best_j][k] - d_ij)).max(0.0))
+            .collect();
+
+        let new_m = remaining.len() + 1;
+        let mut new_d = vec![vec![0.0f64; new_m]; new_m];
+        for (ni, &ri) in remaining.iter().enumerate() {
+            for (nj, &rj) in remaining.iter().enumerate() {
+                new_d[ni][nj] = d[ri][rj];
+            }
+            new_d[ni][new_m - 1] = new_row[ni];
+            new_d[new_m - 1][ni] = new_row[ni];
+        }
+
+        let new_active: Vec<usize> = remaining
+            .iter()
+            .map(|&ri| active[ri])
+            .chain(std::iter::once(u_id))
+            .collect();
+
+        active = new_active;
+        d = new_d;
+    }
+
+    // ── Final join: connect the last two active nodes under a root ───────────
+    let root_id  = nodes.len();
+    let d_last   = d[0][1];
+    let half     = (d_last / 2.0).max(0.0);
+
+    nodes[active[0]].branch_length = half;
+    nodes[active[0]].parent        = Some(root_id);
+    nodes[active[1]].branch_length = half;
+    nodes[active[1]].parent        = Some(root_id);
+
+    nodes.push(Node {
+        id:            root_id,
+        label:         "ROOT".to_string(),
+        is_leaf:       false,
+        children:      vec![active[0], active[1]],
+        parent:        None,
+        branch_length: 0.0,
+        read_count:    0,
+    });
+
+    let mut tree = Tree { nodes, root: root_id };
+
+    // ── Re-root at all-zero haplotype's parent (outgroup rooting) ────────────
+    let zero_profile = vec![0u8; hap_matrix.variants.len()];
+    if let Some(og_hap) = hap_matrix
+        .haplotypes
+        .iter()
+        .find(|h| h.profile == zero_profile)
+    {
+        reroot_at_outgroup_parent(&mut tree, &og_hap.id);
+        eprintln!(
+            "      Rooted at parent of all-zero haplotype ({}) — outgroup rooting",
+            og_hap.id
+        );
+    }
+
+    Ok(tree)
+}
+
+/// Re-root `tree` so that the parent of `outgroup_label` becomes the new root.
+///
+/// Algorithm: collect the path from the outgroup's parent up to the current
+/// root, then reverse all parent–child edges along that path.
+fn reroot_at_outgroup_parent(tree: &mut Tree, outgroup_label: &str) {
+    // Find the outgroup leaf
+    let og_id = match tree.nodes.iter().find(|n| n.label == outgroup_label) {
+        Some(n) => n.id,
+        None    => return,
+    };
+    // The new root will be the outgroup's direct parent
+    let new_root = match tree.nodes[og_id].parent {
+        Some(p) => p,
+        None    => return, // outgroup is already the root
+    };
+    if new_root == tree.root {
+        // Already correctly rooted (outgroup's parent IS the old root)
+        return;
+    }
+
+    // Collect path from new_root up to (and including) old root,
+    // capturing branch lengths BEFORE any modification.
+    let mut path: Vec<usize>  = vec![new_root];
+    let mut bls:  Vec<f64>    = vec![tree.nodes[new_root].branch_length];
+    let mut cur = new_root;
+    loop {
+        match tree.nodes[cur].parent {
+            Some(p) => {
+                path.push(p);
+                bls.push(tree.nodes[p].branch_length);
+                if p == tree.root {
+                    break;
+                }
+                cur = p;
+            }
+            None => break,
+        }
+    }
+
+    // Reverse edges: for each consecutive pair (child, parent) in path,
+    // make child the new parent.
+    for i in 0..path.len().saturating_sub(1) {
+        let child  = path[i];
+        let parent = path[i + 1];
+
+        tree.nodes[parent].children.retain(|&c| c != child);
+        tree.nodes[child].children.push(parent);
+        tree.nodes[parent].parent        = Some(child);
+        tree.nodes[parent].branch_length = bls[i]; // edge length stays the same
+    }
+
+    tree.nodes[new_root].parent        = None;
+    tree.nodes[new_root].branch_length = 0.0;
+    tree.root = new_root;
+
+    // If the old root has become a unary node (only one remaining child),
+    // collapse it: attach its single child directly to new_root, summing
+    // branch lengths so no information is lost.
+    let old_root = *path.last().unwrap();
+    if tree.nodes[old_root].children.len() == 1 {
+        let sole_child = tree.nodes[old_root].children[0];
+        let combined_bl =
+            tree.nodes[old_root].branch_length + tree.nodes[sole_child].branch_length;
+
+        // Replace old_root with sole_child in its parent's children
+        let old_root_parent = tree.nodes[old_root].parent.unwrap();
+        let pos = tree.nodes[old_root_parent]
+            .children
+            .iter()
+            .position(|&c| c == old_root)
+            .unwrap();
+        tree.nodes[old_root_parent].children[pos] = sole_child;
+        tree.nodes[sole_child].parent        = Some(old_root_parent);
+        tree.nodes[sole_child].branch_length = combined_bl;
+    }
+}
+
+/// One unique binary profile and the set of reads that share it.
+pub struct Haplotype {
+    pub id: String,
+    /// Binary state over filtered variants (`profile[v]` = 0 or 1)
+    pub profile: Vec<u8>,
+    pub reads: Vec<String>,
+    pub count: usize,
+}
+
+/// Haplotypes together with the (shared) variant list.
+pub struct HaplotypeMatrix {
+    pub variants: Vec<String>,
+    pub haplotypes: Vec<Haplotype>,
+}
+
+/// HF (heteroplasmic fraction) keyed by Himito variant ID (e.g. "m.13376T>C").
+pub type HfMap = HashMap<String, f64>;
+
+/// Binary matrix after filtering.
+/// `data[variant_idx][read_idx]` is 0 or 1.
+pub struct BinaryMatrix {
+    pub variants: Vec<String>,
+    pub reads: Vec<String>,
+    pub data: Vec<Vec<u8>>,
+}
+
+/// Parse a Himito VCF file and return a map from variant ID to HF value.
+///
+/// Variant IDs are reconstructed as `m.<POS><REF>><ALT>` (the format Himito
+/// uses in the matrix CSV), e.g. `m.13376T>C`.
+pub fn parse_vcf(vcf_path: &str, min_hf: f64, max_hf: f64) -> Result<HfMap> {
+    let mut reader = Reader::from_path(vcf_path)
+        .with_context(|| format!("Cannot read VCF: {vcf_path}"))?;
+    let mut map = HfMap::new();
+
+    for result in reader.records() {
+        #[allow(unused_mut)] // rec.format(b"HF") requires &mut self in rust_htslib
+        let mut rec = result?;
+
+        // Only lineage-informative if the call itself is trustworthy.
+        if !rec.has_filter("PASS".as_bytes()) {
+            continue;
+        }
+
+        // Extract position (0-based in htslib → 1-based in VCF/Himito IDs)
+        let pos = rec.pos() + 1;
+
+        // rec.alleles() borrows rec immutably; collect into owned Strings so
+        // we can release the borrow before calling rec.format() (mut borrow).
+        let (ref_allele, alt_allele) = {
+            let alleles = rec.alleles();
+            (
+                std::str::from_utf8(alleles[0]).unwrap_or("").to_owned(),
+                std::str::from_utf8(alleles[1]).unwrap_or("").to_owned(),
+            )
+        };
+
+        let vid = format!("m.{pos}{ref_allele}>{alt_allele}");
+
+        // Read the HF FORMAT field (float, first sample, first value)
+        if let Ok(hf_data) = rec.format(b"HF").float() {
+            if let Some(&hf_val) = hf_data.get(0).and_then(|s| s.first()) {
+                let hf = hf_val as f64;
+                if hf >= min_hf && hf < max_hf {
+                    map.insert(vid, hf);
+                }
+            }
+        }
+    }
+    Ok(map)
+}
+
+/// Read `matrix_path` (Himito `.matrix.csv`) and apply filters:
+///
+/// * **HF bounds** — keep rows where `min_hf ≤ HF < max_hf`.
+/// * **Prevalence** — keep rows where the variant is present in ≥ `min_presence`
+///   reads AND absent from ≥ `min_absence` reads (guarantees a bifurcation).
+///
+/// Counts are binarised: any count ≥ 1 becomes 1.
+pub fn load_and_filter_matrix(
+    matrix_path: &str,
+    hf_map: &HfMap,   // contains only variants that already passed HF filtering
+    min_presence: usize,
+    min_absence: usize,
+) -> Result<BinaryMatrix> {
+    let mut rdr = csv::Reader::from_path(matrix_path)
+        .with_context(|| format!("Cannot read matrix CSV: {matrix_path}"))?;
+
+    // Column 0 is "variant", then one column per read
+    let headers = rdr.headers()?.clone();
+    let reads: Vec<String> = headers.iter().skip(1).map(String::from).collect();
+
+    let mut variants: Vec<String> = Vec::new();
+    let mut data: Vec<Vec<u8>>   = Vec::new();
+
+    for result in rdr.records() {
+        let rec = result?;
+        let vid = rec[0].to_string();
+
+        // Keep only variants that passed HF filtering in parse_vcf.
+        // If hf_map is empty (no VCF provided), terminate the program.
+        if !hf_map.is_empty() && !hf_map.contains_key(&vid) {
+            continue;
+        }
+
+        // Binarise counts
+        let row: Vec<u8> = rec
+            .iter()
+            .skip(1)
+            .map(|v| u8::from(v.parse::<u32>().unwrap_or(0) >= 1))
+            .collect();
+
+        // Prevalence filter
+        let present = row.iter().filter(|&&b| b == 1).count();
+        let absent  = row.iter().filter(|&&b| b == 0).count();
+        if present < min_presence || absent < min_absence {
+            continue;
+        }
+
+        variants.push(vid);
+        data.push(row);
+    }
+
+    Ok(BinaryMatrix { variants, reads, data })
+}
+
+
+
+/// Collapse reads with identical binary profiles into unique haplotypes.
+///
+/// Haplotypes are sorted by ascending mutation count (number of 1-bits),
+/// then lexicographically by profile, so that the all-zero haplotype (if
+/// present) is always `H0000`.
+pub fn deduplicate(matrix: &BinaryMatrix, min_reads: usize) -> HaplotypeMatrix {
+    let n_reads    = matrix.reads.len();
+    let n_variants = matrix.variants.len();
+
+    // profile (per-read, across variants) → list of read indices
+    let mut profile_map: HashMap<Vec<u8>, Vec<usize>> = HashMap::new();
+
+    for read_idx in 0..n_reads {
+        let profile: Vec<u8> = (0..n_variants)
+            .map(|v| matrix.data[v][read_idx])
+            .collect();
+        profile_map.entry(profile).or_default().push(read_idx);
+    }
+
+    let mut haplotypes: Vec<Haplotype> = profile_map
+        .into_iter()
+        .map(|(profile, indices)| {
+            let reads = indices.iter().map(|&r| matrix.reads[r].clone()).collect();
+            let count = indices.len();
+            Haplotype {
+                id: String::new(), // assigned after sort
+                profile,
+                reads,
+                count,
+            }
+        })
+        .collect();
+    // filter out haplotypes with less than min_reads
+    haplotypes = haplotypes.into_iter().filter(|h| h.count >= min_reads).collect();
+    // Deterministic ordering: fewest mutations first, then lexicographic profile
+    haplotypes.sort_by(|a, b| {
+        let ma: usize = a.profile.iter().map(|&x| x as usize).sum();
+        let mb: usize = b.profile.iter().map(|&x| x as usize).sum();
+        ma.cmp(&mb).then_with(|| a.profile.cmp(&b.profile))
+    });
+
+    for (i, h) in haplotypes.iter_mut().enumerate() {
+        h.id = format!("H{i:04}");
+    }
+
+    HaplotypeMatrix {
+        variants: matrix.variants.clone(),
+        haplotypes,
+    }
+}
+
+
+/// Run the four-gamete test on every pair of variants.
+///
+/// Returns a list of `(variant_a, variant_b)` pairs where all four
+/// combinations {00, 01, 10, 11} are observed across haplotypes — a
+/// violation of the Infinite Sites Assumption (ISA).
+pub fn four_gamete_test(matrix: &HaplotypeMatrix) -> Vec<(String, String)> {
+    let n_var = matrix.variants.len();
+    let mut violations = Vec::new();
+
+    for i in 0..n_var {
+        for j in (i + 1)..n_var {
+            let mut seen = [false; 4]; // index = a*2 + b for states a,b ∈ {0,1}
+            for h in &matrix.haplotypes {
+                let idx = (h.profile[i] * 2 + h.profile[j]) as usize;
+                seen[idx] = true;
+                if seen.iter().all(|&s| s) {
+                    break;
+                }
+            }
+            if seen.iter().all(|&s| s) {
+                violations.push((
+                    matrix.variants[i].clone(),
+                    matrix.variants[j].clone(),
+                ));
+            }
+        }
+    }
+
+    violations
+}
+
+pub fn write_haplotype_map(hap_matrix: &HaplotypeMatrix, path: &str) -> Result<()> {
+    let mut w = BufWriter::new(
+        File::create(path).with_context(|| format!("Cannot create {path}"))?,
+    );
+    writeln!(w, "haplotype_id\tn_mutations\tn_reads\tread_name")?;
+    for hap in &hap_matrix.haplotypes {
+        let n_mut: usize = hap.profile.iter().map(|&b| b as usize).sum();
+        let readlist = hap.reads.join(",");
+        assert_eq!(hap.count, hap.reads.len());
+        writeln!(w, "{}\t{n_mut}\t{}\t{readlist}", hap.id, hap.count)?;
+        
+    }
+    Ok(())
+}
+
+
+pub fn start (matrix_file: &str, vcf_file: Option<&str>, min_hf: f64, max_hf: f64, min_presence: usize, min_absence: usize, min_reads: usize, output_prefix: &str) -> Result<()> {
+    println!("Starting lineage analysis...");
+
+    // ── Step 1: load VCF HF values, then filter the binary matrix ─────────────
+    println!("[1/6] Parsing VCF: {}", vcf_file.as_ref().unwrap_or(&""));
+    let hf_map = parse_vcf(vcf_file.as_ref().unwrap_or(&""), min_hf, max_hf)?;
+
+    eprintln!("[1/6] Loading and filtering matrix: {}", matrix_file);
+    let binary = load_and_filter_matrix(
+        &matrix_file,
+        &hf_map,
+        min_presence,
+        min_absence,
+    )?;
+    println!("{:?}", binary.variants);
+    
+    if binary.variants.is_empty() {
+        anyhow::bail!(
+            "No informative variants remain after filtering. \
+             Try adjusting --min-hf / --max-hf thresholds."
+        );
+    }
+
+    // ── Step 2: deduplicate reads into haplotypes ───────────────────────────────
+    eprintln!("[2/6] Deduplicating reads into haplotypes...");
+    let hap_matrix = deduplicate(&binary, min_reads);
+    // write how many haplotypes and how many heteroplasmic variants on each haplotype
+    eprintln!("[2/6] Found {} haplotypes across {} variants.", hap_matrix.haplotypes.len(), hap_matrix.variants.len());
+    if hap_matrix.haplotypes.is_empty() {
+        anyhow::bail!("No haplotypes remain after --min-reads filtering.");
+    }
+    // Global haplotype map (all variants)
+    let hmap_path = format!("{}.raw_haplotype_map.tsv", output_prefix);
+    write_haplotype_map(&hap_matrix, &hmap_path)?;
+
+    // ── Step 3: pairwise Hamming distance + Neighbor-Joining tree as the initial tree ─────────────
+    eprintln!("[4/5] Building Neighbor-Joining tree");
+    let dist = hamming_distance_matrix(&hap_matrix);
+    let tree = neighbor_joining(&dist, &hap_matrix)?;
+    eprintln!("      Tree has {} nodes ({} leaves)", tree.nodes.len(), tree.leaves().len());
+
+    // ── Step 4: MCMC  idenitfy the ─────────────
+
+
+
+
+
+
+
+
+    Ok(())
+
+}
\ No newline at end of file
diff --git a/src/main.rs b/src/main.rs
index 3bc7b6b..7b74ff0 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -12,6 +12,7 @@ mod methyl;
 mod correct;
 mod minorhap;
 mod callnumts;
+mod lineage;
 
 #[derive(Debug, Parser)]
 #[clap(name = "Himito")]
@@ -354,6 +355,45 @@ enum Commands {
         /// minimal allele count for variants
         #[clap(short, long, value_parser, default_value_t = 2)]
         ac_threshold: i32,
+    },
+
+    /// Group reads into haplotypes from heteroplasmic variants and test for
+    /// four-gamete (recombination) violations between them.
+    #[clap(arg_required_else_help = true)]
+    Lineage {
+        /// path for the Himito matrix CSV (<prefix>.matrix.csv)
+        #[clap(short, long, value_parser, required = true)]
+        matrix_file: PathBuf,
+
+        /// path for the Himito VCF; used to read each variant's HF and
+        /// restrict the analysis to heteroplasmic (non-fixed) sites
+        #[clap(short, long, value_parser)]
+        vcf_file: Option<PathBuf>,
+
+        /// minimal heteroplasmic frequency (inclusive) for a variant to be considered
+        #[clap(long, value_parser, default_value_t = 0.01)]
+        min_hf: f64,
+
+        /// maximal heteroplasmic frequency (exclusive) for a variant to be considered;
+        /// default of 1.0 excludes fixed/homoplasmic variants (HF == 0.95)
+        #[clap(long, value_parser, default_value_t = 0.95)]
+        max_hf: f64,
+
+        /// minimal number of reads a variant must be present in to be informative
+        #[clap(long, value_parser, default_value_t = 2)]
+        min_presence: usize,
+
+        /// minimal number of reads a variant must be absent from to be informative
+        #[clap(long, value_parser, default_value_t = 1)]
+        min_absence: usize,
+
+        /// minimal number of reads required to report a haplotype
+        #[clap(long, value_parser, default_value_t = 3)]
+        min_reads: usize,
+
+        /// output prefix; writes <prefix>.haplotype_map.tsv
+        #[clap(short, long, value_parser, required = true)]
+        output_prefix: String,
     }
 }
 
@@ -556,6 +596,34 @@ fn main() {
         } => {
             callnumts::start(&input_bam, &chromo, max_gap_threshold, &output_vcf, &reference_file, &sample_name, ac_threshold);
         }
+        Commands::Lineage {
+            matrix_file,
+            vcf_file,
+            min_hf,
+            max_hf,
+            min_presence,
+            min_absence,
+            min_reads,
+            output_prefix,
+        } => {
+            let matrix_file = matrix_file.to_str().expect("matrix-file path is not valid UTF-8");
+            let vcf_file = vcf_file
+                .as_deref()
+                .map(|p| p.to_str().expect("vcf-file path is not valid UTF-8"));
+            if let Err(e) = lineage::start(
+                matrix_file,
+                vcf_file,
+                min_hf,
+                max_hf,
+                min_presence,
+                min_absence,
+                min_reads,
+                &output_prefix,
+            ) {
+                eprintln!("Error running lineage analysis: {:#}", e);
+                std::process::exit(1);
+            }
+        }
     }
 
     
diff --git a/src/methyl.rs b/src/methyl.rs
index eae8b48..74efce8 100644
--- a/src/methyl.rs
+++ b/src/methyl.rs
@@ -765,7 +765,9 @@ fn write_methylation_to_csv<P: AsRef<Path>>(
         .collect();
 
     // construct matrix
-    let mut matrix = Array2::<f64>::zeros((refpos_list.len(), read_vec.len()));
+    // let mut matrix = Array2::<f64>::zeros((refpos_list.len(), read_vec.len()));
+    let mut matrix = Array2::<f64>::from_elem((refpos_list.len(), read_vec.len()), f64::NAN);
+    
     for ((_pos, refpos, _motif), d) in methyl.iter() {
         let row_index = ref_pos_dict.get(refpos).unwrap();
         for (read_name, &likelihood) in d.iter() {