homopolymer filter

helps a tiny bit with recall/precision - genotyping concordance
essentially the same

Author: ACEnglish

experiments/bedtest.sh

@@ -7,10 +7,10 @@ create() {
         --input test_rs/test2.vcf.gz \
         --bam /Users/english/code/kanpig/experiments/test_rs/NA24385.chr20.bam \
         --reference /Users/english/code/references/grch38/GRCh38_1kg_mainchrs.fa \
-        --sizemin 20 \
+        --sizemin 50 \
         --sizesim 0.95 --seqsim 0.90 --threads 4 \
-        --maxpaths 10000 --mapq 20 --hapsim 0.98 \
-        --chunksize 1000 --prune --try-exact \
+        --maxpaths 10000 --mapq 5 --hapsim 0.98 \
+        --chunksize 100 --maxhom 0 \
         -o test_rs/hc.vcf --bed $bed 
     # --bed /Users/english/code/kanpig/test/GRCh38_HG002-T2TQ100-V1.0_stvar.benchmark.bed \
     # --bam /Users/english/code/kanpig/experiments/test_rs/GIABHG002.bam \

src/kplib/bamparser.rs

@@ -124,7 +124,12 @@ impl BamParser {
             };
 
             let n_hap = Haplotype::new(
-                seq_to_kmer(&sequence, self.params.kmer, p.indel == Svtype::Del),
+                seq_to_kmer(
+                    &sequence,
+                    self.params.kmer,
+                    p.indel == Svtype::Del,
+                    self.params.maxhom,
+                ),
                 p.size,
                 1,
                 1,

src/kplib/cli.rs

@@ -109,6 +109,10 @@ pub struct KDParams {
     /// Don't require alignments to span vargraph region
     #[arg(long, default_value_t = true)]
     pub spanoff: bool,
+
+    /// Maximum homopolymer length to kmerize (off=0)
+    #[arg(long, default_value_t = 5)]
+    pub maxhom: usize,
 }
 
 impl ArgParser {

src/kplib/cluster.rs

@@ -43,6 +43,9 @@ pub fn cluster_haplotypes(
             hap_b.push(hap);
         }
     }
+    // Averaging Haplotypes could be done here.
+    // However, 'low-quality' long-reads are now ~97% (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10070092/)
+    // Which is reasonably within the thresholds of size/seqsim's search space
     hap_a.sort();
     hap_b.sort();
 

src/kplib/haplotype.rs

@@ -25,7 +25,7 @@ impl Haplotype {
 
     // Create an empty haplotype
     pub fn blank(kmer: u8, coverage: u64) -> Haplotype {
-        let mk = seq_to_kmer(&[], kmer, false);
+        let mk = seq_to_kmer(&[], kmer, false, 0);
         Haplotype {
             size: 0,
             n: 0,

src/kplib/kmeans.rs

@@ -67,6 +67,7 @@ pub fn kmeans(data: &[Point], k: usize) -> Vec<Cluster> {
         .map(|centroid| Cluster::new(centroid.clone()))
         .collect();
 
+    let mut ntries = 10;
     loop {
         // Assign each point to the nearest cluster
         for (idx, point) in data.iter().enumerate() {
@@ -93,7 +94,7 @@ pub fn kmeans(data: &[Point], k: usize) -> Vec<Cluster> {
 
         // Check for convergence
         let new_centroids: Vec<Point> = clusters.iter().map(|c| c.centroid.clone()).collect();
-        if new_centroids == old_centroids {
+        if new_centroids == old_centroids || ntries == 0 {
             break;
         } else {
             centroids = new_centroids;
@@ -102,6 +103,7 @@ pub fn kmeans(data: &[Point], k: usize) -> Vec<Cluster> {
                 .map(|centroid| Cluster::new(centroid.clone()))
                 .collect();
         }
+        ntries -= 1;
     }
 
     clusters

src/kplib/kmer.rs

@@ -10,7 +10,12 @@ fn encode_nuc(nuc: u8) -> u64 {
 }
 
 /// Count kmers in a sequence
-pub fn seq_to_kmer(sequence: &[u8], kmer: u8, negative: bool) -> Vec<f32> {
+pub fn seq_to_kmer(sequence: &[u8], kmer: u8, negative: bool, maxhom: usize) -> Vec<f32> {
+    let sequence = if maxhom != 0 {
+        compress_homopolymer(sequence, maxhom)
+    } else {
+        sequence.to_vec()
+    };
     let ukmer = kmer as usize;
     let mut kcounts = vec![0f32; 1 << (2 * ukmer)];
     let cnt = if negative { -1.0 } else { 1.0 };
@@ -47,3 +52,27 @@ pub fn seq_to_kmer(sequence: &[u8], kmer: u8, negative: bool) -> Vec<f32> {
 
     kcounts
 }
+
+pub fn compress_homopolymer(vector: &[u8], maxspan: usize) -> Vec<u8> {
+    let mut result = Vec::new();
+    let mut count = 0;
+    let mut prev_byte = None;
+
+    for byte in vector {
+        match prev_byte {
+            Some(prev) if prev == byte => {
+                count += 1;
+                if count < maxspan {
+                    result.push(*byte);
+                }
+            }
+            _ => {
+                count = 1;
+                result.push(*byte);
+            }
+        }
+        prev_byte = Some(byte);
+    }
+
+    result
+}

src/kplib/vargraph.rs

@@ -21,11 +21,11 @@ pub struct VarNode {
 }
 
 impl VarNode {
-    pub fn new(entry: vcf::Record, kmer: u8) -> Self {
+    pub fn new(entry: vcf::Record, kmer: u8, maxhom: usize) -> Self {
         // Want to make a hash for these names for debugging, I think.
         let name = "".to_string();
         let (start, end) = entry.boundaries();
-        let (kfeat, size) = entry.to_kfeat(kmer);
+        let (kfeat, size) = entry.to_kfeat(kmer, maxhom);
         Self {
             name,
             start,
@@ -70,7 +70,7 @@ pub struct Variants {
 /// The graph has an upstream 'src' node that point to every variant node
 /// The graph has a dnstream 'snk' node that is pointed to by every variant node and 'src'
 impl Variants {
-    pub fn new(mut variants: Vec<vcf::Record>, kmer: u8) -> Self {
+    pub fn new(mut variants: Vec<vcf::Record>, kmer: u8, maxhom: usize) -> Self {
         if variants.is_empty() {
             panic!("Cannot create a graph from no variants");
         }
@@ -84,7 +84,7 @@ impl Variants {
         node_indices.append(
             &mut variants
                 .drain(..) // drain lets us move the entry without a clone
-                .map(|entry| graph.add_node(VarNode::new(entry, kmer)))
+                .map(|entry| graph.add_node(VarNode::new(entry, kmer, maxhom)))
                 .collect(),
         );
 

src/kplib/vcf_traits.rs

@@ -32,7 +32,7 @@ impl FromStr for Svtype {
 
 /// Convert vcf::Record to kfeat
 pub trait KdpVcf {
-    fn to_kfeat(&self, kmer: u8) -> (Vec<f32>, i64);
+    fn to_kfeat(&self, kmer: u8, maxhom: usize) -> (Vec<f32>, i64);
     fn boundaries(&self) -> (u64, u64);
     fn size(&self) -> u64;
     fn is_filtered(&self) -> bool;
@@ -41,7 +41,7 @@ pub trait KdpVcf {
 
 impl KdpVcf for vcf::Record {
     /// Convert variant sequence to Kfeat
-    fn to_kfeat(&self, kmer: u8) -> (Vec<f32>, i64) {
+    fn to_kfeat(&self, kmer: u8, maxhom: usize) -> (Vec<f32>, i64) {
         let ref_seq = self.reference_bases().to_string();
         let alt_seq = self
             .alternate_bases()
@@ -51,8 +51,8 @@ impl KdpVcf for vcf::Record {
 
         let size = alt_seq.len() as i64 - ref_seq.len() as i64;
 
-        let m_ref = seq_to_kmer(&ref_seq.as_bytes()[1..], kmer, false);
-        let m_alt = seq_to_kmer(&alt_seq.as_bytes()[1..], kmer, false);
+        let m_ref = seq_to_kmer(&ref_seq.as_bytes()[1..], kmer, false, maxhom);
+        let m_alt = seq_to_kmer(&alt_seq.as_bytes()[1..], kmer, false, maxhom);
 
         let m_ret: Vec<_> = m_alt
             .iter()

src/main.rs

@@ -67,7 +67,7 @@ fn main() {
         thread::spawn(move || {
             let mut m_bam = BamParser::new(m_args.io.bam, m_args.io.reference, m_args.kd.clone());
             for chunk in receiver.into_iter().flatten() {
-                let mut m_graph = Variants::new(chunk, m_args.kd.kmer);
+                let mut m_graph = Variants::new(chunk, m_args.kd.kmer, m_args.kd.maxhom);
                 let (haps, coverage) = m_bam.find_haps(&m_graph.chrom, m_graph.start, m_graph.end);
                 let (h1, h2) = cluster_haplotypes(haps, coverage, &m_args.kd);
                 let p1 = m_graph.apply_coverage(&h1, &m_args.kd);