circos dev

Cargo.lock

@@ -148,6 +148,27 @@ dependencies = [
  "cfg-if",
 ]
 
+[[package]]
+name = "csv"
+version = "1.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "acdc4883a9c96732e4733212c01447ebd805833b7275a73ca3ee080fd77afdaf"
+dependencies = [
+ "csv-core",
+ "itoa",
+ "ryu",
+ "serde",
+]
+
+[[package]]
+name = "csv-core"
+version = "0.1.12"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7d02f3b0da4c6504f86e9cd789d8dbafab48c2321be74e9987593de5a894d93d"
+dependencies = [
+ "memchr",
+]
+
 [[package]]
 name = "data-url"
 version = "0.3.1"
@@ -361,6 +382,7 @@ name = "pandora_lib_graph"
 version = "0.1.0"
 dependencies = [
  "anyhow",
+ "csv",
  "env_logger",
  "log",
  "ureq",
@@ -520,6 +542,12 @@ dependencies = [
  "unicode-script",
 ]
 
+[[package]]
+name = "ryu"
+version = "1.0.20"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "28d3b2b1366ec20994f1fd18c3c594f05c5dd4bc44d8bb0c1c632c8d6829481f"
+
 [[package]]
 name = "serde"
 version = "1.0.219"

Cargo.toml

@@ -5,6 +5,7 @@ edition = "2021"
 
 [dependencies]
 anyhow = "1.0.89"
+csv = "1.3.1"
 env_logger = "0.11.8"
 log = "0.4.27"
 # typst = "0.13.0"

src/circ.rs

@@ -4,6 +4,7 @@ use std::{
     io::{BufRead, BufReader, Write},
 };
 
+use csv::ReaderBuilder;
 use log::info;
 
 use crate::cytoband::Range;
@@ -593,11 +594,16 @@ impl Circos {
         ));
                         }
                     }
+                    // Curvature is controlled by `bulge` interpreted as a fraction of `radius`.
+                    //   bulge = 0.0  → almost straight (near the rim)
+                    //   bulge = 0.3  → typical Circos-like inward arc
+                    //   bulge = 1.0  → controls at the center (very tight arc)
+                    // Values >1.0 are clamped to center; negatives give outward curves if you allow them.
                     TrackItem::Cord {
                         start,
                         end,
                         radius,
-                        bulge,
+                        bulge, // interpret as fraction of radius
                         color,
                         thickness,
                     } => {
@@ -609,30 +615,36 @@ impl Circos {
                                 }
                             }
                         };
-                        if let (Some(angle0), Some(angle1)) = (get_angle(start), get_angle(end)) {
-                            // Endpoints
-                            let x0 = track.cx + radius * angle0.cos();
-                            let y0 = track.cy + radius * angle0.sin();
-                            let x1 = track.cx + radius * angle1.cos();
-                            let y1 = track.cy + radius * angle1.sin();
-
-                            // Control point: midpoint, bulged out by `bulge`
-                            let mid_angle = 0.5 * (angle0 + angle1);
-                            let ctrl_r = radius + bulge;
-                            let cx_ctrl = track.cx + ctrl_r * mid_angle.cos();
-                            let cy_ctrl = track.cy + ctrl_r * mid_angle.sin();
-
-                            // Bounding box
-                            for &(x, y) in &[(x0, y0), (x1, y1), (cx_ctrl, cy_ctrl)] {
-                                state.min_x = state.min_x.min(x);
-                                state.max_x = state.max_x.max(x);
-                                state.min_y = state.min_y.min(y);
-                                state.max_y = state.max_y.max(y);
+
+                        if let (Some(a0), Some(a1)) = (get_angle(start), get_angle(end)) {
+                            // Endpoints on the circle
+                            let x0 = track.cx + radius * a0.cos();
+                            let y0 = track.cy + radius * a0.sin();
+                            let x1 = track.cx + radius * a1.cos();
+                            let y1 = track.cy + radius * a1.sin();
+
+                            // Map bulge fraction -> inward control radius
+                            // Allow negative for outward (optional). Clamp inward to >= 0.
+                            let r_ctrl = (radius * (1.0 - bulge)).max(0.0);
+
+                            // Control points pulled along the radial lines of each endpoint
+                            let c1x = track.cx + r_ctrl * a0.cos();
+                            let c1y = track.cy + r_ctrl * a0.sin();
+                            let c2x = track.cx + r_ctrl * a1.cos();
+                            let c2y = track.cy + r_ctrl * a1.sin();
+
+                            // Bounding box (include stroke half-width)
+                            let half_w = *thickness as f64 / 2.0;
+                            for &(x, y) in &[(x0, y0), (x1, y1), (c1x, c1y), (c2x, c2y)] {
+                                state.min_x = state.min_x.min(x - half_w);
+                                state.max_x = state.max_x.max(x + half_w);
+                                state.min_y = state.min_y.min(y - half_w);
+                                state.max_y = state.max_y.max(y + half_w);
                             }
 
                             state.svg_elements.push(format!(
-            r#"<path d="M {:.2} {:.2} Q {:.2} {:.2} {:.2} {:.2}" stroke="{}" stroke-width="{}" fill="none"/>"#,
-            x0, y0, cx_ctrl, cy_ctrl, x1, y1, color, thickness
+            r#"<path d="M {:.2} {:.2} C {:.2} {:.2} {:.2} {:.2} {:.2} {:.2}" stroke="{}" stroke-width="{}" fill="none"/>"#,
+            x0, y0, c1x, c1y, c2x, c2y, x1, y1, color, thickness
         ));
                         }
                     }
@@ -847,3 +859,67 @@ impl Circos {
         Ok(())
     }
 }
+
+// #[derive(Debug)]
+// pub struct GeneLabel {
+//     pub gene: String,
+//     pub chr: String,
+//     pub pos: u32,
+// }
+//
+// pub fn read_gene_labels(path: &str) -> anyhow::Result<Vec<GeneLabel>> {
+//     let mut rdr = ReaderBuilder::new().delimiter(b'\t').from_path(path)?;
+//     let mut records = Vec::new();
+//
+//     for result in rdr.records() {
+//         let record = result?;
+//         let gene = record[0].to_string();
+//         let chr = record[1].to_string();
+//         let pos: u32 = record[2].parse()?;
+//         records.push(GeneLabel { gene, chr, pos });
+//     }
+//     Ok(records)
+// }
+//
+// fn avg_angle_short_way(a: f64, b: f64) -> f64 {
+//     // Vector average: works across the -π/π wrap and picks the short way
+//     (a.sin() + b.sin()).atan2(a.cos() + b.cos())
+// }
+
+#[derive(Debug)]
+pub struct TranslocData {
+    pub case_id: String,
+    pub left_chr: String,
+    pub left_pos: u32,
+    pub left_gene: String,
+    pub right_chr: String,
+    pub right_pos: u32,
+    pub right_gene: String,
+}
+
+pub fn read_translocs(path: &str) -> anyhow::Result<Vec<TranslocData>> {
+    let mut rdr = ReaderBuilder::new().delimiter(b'\t').from_path(path)?;
+    let mut records = Vec::new();
+
+    for result in rdr.records() {
+        let record = result?;
+        let case_id = record[0].to_string();
+        let left_chr = record[1].to_string();
+        let left_pos: u32 = record[2].parse()?;
+        let left_gene = record[3].to_string();
+        let right_chr = record[4].to_string();
+        let right_pos: u32 = record[5].parse()?;
+        let right_gene = record[6].to_string();
+
+        records.push(TranslocData {
+            case_id,
+            left_chr,
+            left_pos,
+            left_gene,
+            right_chr,
+            right_pos,
+            right_gene,
+        });
+    }
+    Ok(records)
+}

src/lib.rs

@@ -6,12 +6,15 @@ pub mod theme;
 
 #[cfg(test)]
 mod tests {
-    use std::{f64::consts::PI, fs::File, io::Write};
+    use std::collections::HashMap;
 
     use cytoband::{read_ranges, svg_chromosome, AdditionalRect, Lollipop, RectPosition};
     use report::compile_typst_report;
 
-    use crate::circos::{classic_caps, no_caps, AngleRange, Circos, LabelMode};
+    use crate::{
+        circ::read_translocs,
+        circos::{classic_caps, no_caps, AngleRange, Circos, LabelMode},
+    };
 
     use super::*;
 
@@ -145,93 +148,226 @@ mod tests {
 
         let mut band_track = circ::Track::new(cx, cy, angle_start, angle_end, gap);
         let mut annot_track = circ::Track::new(cx, cy, angle_start, angle_end, gap);
-        for c in 1..=22 {
-            let name = format!("chr{c}");
-            let ranges = read_ranges("/data/ref/hs1/cytoBandMapped.bed", &name)?;
 
-            band_track.add_range_set(name.clone(), ranges.clone());
-            annot_track.add_range_set(name, ranges);
+        // let mut chr_to_idx = HashMap::new(:);
+        let mut chr_idx = vec![String::new(); 23];
+
+        for c in 0..23 {
+            let name = if c <= 21 {
+                format!("chr{}", c + 1)
+            } else if c == 22 {
+                "chrX".to_string()
+            } else {
+                panic!("unexpected");
+            };
+            chr_idx[c] = name;
+        }
+
+        // for c in 1..=23 {
+        //     let name = if c <= 22 { format!("chr{c}") } else if c == 23 { "chrX".to_string() };
+        //
+        //     let ranges = read_ranges("/data/ref/hs1/cytoBandMapped.bed", &name)?;
+        //
+        //     band_track.add_range_set(name.clone(), ranges.clone());
+        //     annot_track.add_range_set(name, ranges);
+        //
+        //     let set_idx = band_track.range_sets.len() - 1;
+        //     band_track.add_arcs_for_set(set_idx, r * 0.955555, r, circ::classic_caps);
+        // }
+
+        for (set_idx, name) in chr_idx.iter().enumerate() {
+            println!("{name}");
+            let ranges = read_ranges("/data/ref/hs1/cytoBandMapped.bed", name)?;
+
+            band_track.add_range_set(name.to_string(), ranges.clone());
+            annot_track.add_range_set(name.to_string(), ranges);
 
-            let set_idx = band_track.range_sets.len() - 1;
             band_track.add_arcs_for_set(set_idx, r * 0.955555, r, circ::classic_caps);
         }
 
-        let r2 = r + 10.0;
-        let r3 = r + 50.0;
+        // Create two tracks with reference genome size
+        // for c in 1..=22 {
+        //     let name = format!("chr{c}");
+        //
+        //     let ranges = read_ranges("/data/ref/hs1/cytoBandMapped.bed", &name)?;
+        //
+        //     band_track.add_range_set(name.clone(), ranges.clone());
+        //     annot_track.add_range_set(name, ranges);
+        //
+        //     let set_idx = band_track.range_sets.len() - 1;
+        //     band_track.add_arcs_for_set(set_idx, r * 0.955555, r, circ::classic_caps);
+        // }
+
+        let r2 = r + 100.0;
+        let r3 = r + 200.0;
 
-        for c in 0..22 {
-            let name = format!("chr{}", c + 1);
+        // Adding chromosome names labels
+        for (set_idx, name) in chr_idx.iter().enumerate() {
+            let ranges = read_ranges("/data/ref/hs1/cytoBandMapped.bed", name)?;
+
+            let p = ranges
+                .iter()
+                .find_map(|r| {
+                    if r.category.as_str() == "acen" {
+                        Some(r.start + r.end.saturating_sub(r.start))
+                    } else {
+                        None
+                    }
+                })
+                .unwrap_or(0);
+            println!("{p}");
 
             annot_track.add_label(
                 circ::Attach::Data {
-                    set_idx: c,
-                    start: 0,
+                    set_idx,
+                    start: p,
                     end: None,
                 },
-                // circ::Attach::Angle {
-                //     start_angle: 50f64.to_radians(),
-                //     end_angle: None,
-                // },
-                r3 + 20.0,
+                r3,
                 circ::LabelMode::Perimeter,
                 28.0,
                 Some("white".to_string()),
                 0.85,
-                name,
+                name.to_string(),
             );
             annot_track.add_arc(
                 circ::Attach::Data {
-                    set_idx: c,
-                    start: 0,
-                    end: Some(100_000),
+                    set_idx,
+                    start: p,
+                    end: Some(p),
                 },
-                r2 + 1.0,
-                r2 + 20.0,
+                r + 10.0,
+                r3 - 35.0,
                 "gpos50".to_string(),
                 circ::CapStyle::None,
             );
         }
-        // Add a label at chr2, pos 40Mb
-        band_track.add_label(
-            circ::Attach::Data {
-                set_idx: 1,
-                start: 40000000,
-                end: None,
-            },
-            r + 40.0,
-            circ::LabelMode::Perimeter,
-            22.0,
-            Some("white".to_string()),
-            0.9,
-            "Hello".to_string(),
-        );
-
-        // Add an annotation arc and label by angle
-        annot_track.add_arc(
-            circ::Attach::Angle {
-                start_angle: 25f64.to_radians(),
-                end_angle: Some(60f64.to_radians()),
-            },
-            (r3 * 0.955555) + 5.0,
-            r3 - 10.0,
-            "gpos50".to_string(),
-            circ::CapStyle::None,
-        );
-        annot_track.add_path(
-            circ::Attach::Data {
-                set_idx: 1,
-                start: 40000000,
-                end: None,
-            },
-            circ::Attach::Angle {
-                start_angle: 25f64.to_radians(),
-                end_angle: None,
-            },
-            r + 40.0,
-            r + 70.0,
-            "red".to_string(),
-            2,
-        );
+
+        let translocs = read_translocs("/data/tmp_trl.tsv")?;
+        let mut translocs_labels: HashMap<String, Vec<(String, u32)>> = HashMap::new();
+
+        for trl in translocs {
+            if let (Some(left_set_idx), Some(right_set_idx)) = (
+                chr_idx.iter().enumerate().find_map(|(i, name)| {
+                    if name == &trl.left_chr {
+                        Some(i)
+                    } else {
+                        None
+                    }
+                }),
+                chr_idx.iter().enumerate().find_map(|(i, name)| {
+                    if name == &trl.right_chr {
+                        Some(i)
+                    } else {
+                        None
+                    }
+                }),
+            ) {
+                annot_track.add_cord(
+                    circ::Attach::Data {
+                        set_idx: left_set_idx,
+                        start: trl.left_pos,
+                        end: None,
+                    },
+                    circ::Attach::Data {
+                        set_idx: right_set_idx,
+                        start: trl.right_pos,
+                        end: None,
+                    },
+                    r - 100.0,
+                    0.3,
+                    "red".to_string(),
+                    2,
+                );
+
+                translocs_labels
+                    .entry(trl.left_gene)
+                    .or_default()
+                    .push((trl.left_chr.to_string(), trl.left_pos));
+                translocs_labels
+                    .entry(trl.right_gene)
+                    .or_default()
+                    .push((trl.right_chr.to_string(), trl.right_pos));
+            }
+        }
+
+        let mut trl_labels: Vec<(String, usize, u32, usize)> = translocs_labels
+            .iter()
+            .filter_map(|(gene, positions)| {
+                let mut chrs: Vec<String> = positions.iter().map(|(c, _)| c.to_string()).collect();
+                chrs.sort();
+                chrs.dedup();
+                let chr = if chrs.len() != 1 {
+                    panic!("More than one position for a label: {gene} {positions:?}");
+                } else {
+                    chrs.first().unwrap()
+                };
+
+                let n_alt = positions.len();
+                let pos = (positions.iter().map(|(_, p)| *p as f64).sum::<f64>() / n_alt as f64)
+                    .round() as u32;
+                chr_idx
+                        .iter()
+                        .enumerate()
+                        .find_map(|(i, name)| if name == chr { Some(i) } else { None }).map(|set_idx| (gene.to_string(), set_idx, pos, n_alt))
+            })
+            .collect();
+
+        trl_labels.sort_by(|a, b| (a.1, a.2).cmp(&(b.1, b.2)));
+
+        let mut labels_pos = Vec::new();
+        let dist_thre = 10_000_000;
+        for (label, id, pos, n) in trl_labels {
+            let (_last_label, last_id, last_pos, _last_n, last_bump) = if labels_pos.is_empty() {
+                labels_pos.push((label, id, pos, n, None));
+                continue;
+            } else {
+                labels_pos.last().unwrap()
+            };
+
+            let last_ref_pos = last_bump.unwrap_or(*last_pos);
+
+            if *last_id == id && pos.saturating_sub(last_ref_pos) < dist_thre {
+                labels_pos.push((label, id, pos, n, Some(last_ref_pos + dist_thre)));
+            } else {
+                labels_pos.push((label, id, pos, n, None));
+            }
+        }
+
+        println!("{labels_pos:?}");
+
+        for (label, set_idx, pos, n, bump) in labels_pos {
+            annot_track.add_label(
+                circ::Attach::Data {
+                    set_idx,
+                    start: bump.unwrap_or(pos),
+                    end: None,
+                },
+                r2,
+                circ::LabelMode::Perimeter,
+                18.0,
+                Some("yellow".to_string()),
+                0.95,
+                format!("{label} (n={n})"),
+            );
+
+            annot_track.add_path(
+                circ::Attach::Data {
+                    set_idx,
+                    start: pos,
+                    end: None,
+                },
+                circ::Attach::Data {
+                    set_idx,
+                    start: bump.unwrap_or(pos),
+                    end: None,
+                },
+                r + 10.0,
+                r2 - 70.0,
+                "black".to_string(),
+                1
+            );
+        }
 
         circos.add_track(band_track);
         circos.add_track(annot_track);
@@ -274,7 +410,6 @@ mod tests {
                 },
                 950.0,
                 -580.0,
-
                 "red".to_string(),
                 200,
             );
@@ -298,6 +433,7 @@ mod tests {
         c.save_to_file("/data/benti.svg", 100.0)?;
         Ok(())
     }
+
     #[test]
     fn pdf() {
         compile_typst_report("LEVASSEUR").unwrap();