use anyhow::{Context, Result};
use log::debug;
use plotters::prelude::*;
use plotters_svg::SVGBackend;
use std::path::Path;
const SCALE: u32 = 1;
const fn s(v: u32) -> u32 {
v * SCALE
}
const WIDTH: u32 = 1024;
const HEIGHT: u32 = 768;
const CHART_BG: RGBColor = RGBColor(230, 230, 230);
const PLOT_BG: RGBColor = RGBColor(255, 255, 255);
const GRIDLINE_COLOR: RGBColor = RGBColor(192, 168, 144);
const TITLE_COLOR: RGBColor = RGBColor(64, 64, 64);
const COVERAGE_LINE_COLOR: RGBColor = RGBColor(255, 0, 0);
const HISTOGRAM_BAR_BORDER: RGBColor = RGBColor(50, 50, 50);
const PIE_COLOR_0: RGBColor = RGBColor(0xEC, 0x62, 0x5C); const PIE_COLOR_1: RGBColor = RGBColor(0x55, 0x55, 0xF6); const PIE_COLOR_2: RGBColor = RGBColor(0x8A, 0xFC, 0x6E);
const PIE_SHADOW_COLOR: RGBColor = RGBColor(0x80, 0x80, 0x80);
const PIE_LABEL_BG: RGBColor = RGBColor(255, 255, 204);
const PIE_LABEL_BORDER: RGBColor = RGBColor(153, 153, 153);
pub fn coverage_profile_plot(
profile: &[f64; 100],
title: &str,
sample_name: &str,
output_path: &Path,
) -> Result<()> {
{
let root = BitMapBackend::new(output_path, (s(WIDTH), s(HEIGHT))).into_drawing_area();
render_coverage_profile(&root, profile, title, sample_name, SCALE as f64)?;
root.present()
.context("Failed to write coverage profile PNG")?;
}
let svg_path = output_path.with_extension("svg");
{
let root = SVGBackend::new(&svg_path, (WIDTH, HEIGHT)).into_drawing_area();
render_coverage_profile(&root, profile, title, sample_name, 1.0)?;
root.present()
.context("Failed to write coverage profile SVG")?;
}
debug!("Wrote coverage profile: {}", output_path.display());
Ok(())
}
fn render_coverage_profile<DB: DrawingBackend>(
root: &DrawingArea<DB, plotters::coord::Shift>,
profile: &[f64; 100],
title: &str,
sample_name: &str,
pxs: f64,
) -> Result<()>
where
DB::ErrorType: 'static,
{
let ps = |v: f64| (v * pxs) as u32;
root.fill(&CHART_BG).map_err(|e| anyhow::anyhow!("{}", e))?;
let y_max = profile.iter().cloned().fold(0.0f64, f64::max);
let y_ceil = nice_ceil(y_max);
use plotters::style::text_anchor::{HPos, Pos, VPos};
let title_style = ("sans-serif", ps(24.0) as f64)
.into_font()
.style(FontStyle::Bold)
.color(&TITLE_COLOR)
.pos(Pos::new(HPos::Center, VPos::Top));
let subtitle_style = ("sans-serif", ps(16.0) as f64)
.into_font()
.color(&TITLE_COLOR)
.pos(Pos::new(HPos::Center, VPos::Top));
let (w, _h) = root.dim_in_pixel();
let center_x = w as i32 / 2;
let title_y = ps(10.0) as i32;
let subtitle_y = title_y + ps(28.0) as i32;
root.draw_text(title, &title_style, (center_x, title_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
root.draw_text(sample_name, &subtitle_style, (center_x, subtitle_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let chart_top = ps(55.0);
let chart_area = root.clone().margin(chart_top, ps(20.0), ps(20.0), ps(20.0));
let mut chart = ChartBuilder::on(&chart_area)
.x_label_area_size(ps(40.0))
.y_label_area_size(ps(70.0))
.build_cartesian_2d(0.0..100.0f64, 0.0..y_ceil)
.map_err(|e| anyhow::anyhow!("{}", e))?;
chart
.plotting_area()
.fill(&PLOT_BG)
.map_err(|e| anyhow::anyhow!("{}", e))?;
chart
.configure_mesh()
.x_desc("Transcript position")
.y_desc("Counts")
.x_labels(21) .label_style(("sans-serif", ps(15.0) as f64).into_font().color(&BLACK))
.axis_desc_style(("sans-serif", ps(16.0) as f64).into_font().color(&BLACK))
.x_label_formatter(&|v| format!("{}", *v as i32))
.y_label_formatter(&|v| format_y_label(*v))
.axis_style(BLACK.stroke_width(1))
.light_line_style(GRIDLINE_COLOR.mix(0.0)) .bold_line_style(GRIDLINE_COLOR.stroke_width(1))
.set_all_tick_mark_size(ps(4.0))
.draw()
.map_err(|e| anyhow::anyhow!("{}", e))?;
let data: Vec<(f64, f64)> = profile
.iter()
.enumerate()
.map(|(i, &v)| (i as f64, v))
.collect();
chart
.draw_series(LineSeries::new(
data,
COVERAGE_LINE_COLOR.stroke_width(ps(2.0)),
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
Ok(())
}
pub fn coverage_histogram_plot(
histogram: &[u64; 51],
sample_name: &str,
output_path: &Path,
) -> Result<()> {
{
let root = BitMapBackend::new(output_path, (s(WIDTH), s(HEIGHT))).into_drawing_area();
render_coverage_histogram(&root, histogram, sample_name, SCALE as f64)?;
root.present()
.context("Failed to write coverage histogram PNG")?;
}
let svg_path = output_path.with_extension("svg");
{
let root = SVGBackend::new(&svg_path, (WIDTH, HEIGHT)).into_drawing_area();
render_coverage_histogram(&root, histogram, sample_name, 1.0)?;
root.present()
.context("Failed to write coverage histogram SVG")?;
}
debug!("Wrote coverage histogram: {}", output_path.display());
Ok(())
}
fn render_coverage_histogram<DB: DrawingBackend>(
root: &DrawingArea<DB, plotters::coord::Shift>,
histogram: &[u64; 51],
sample_name: &str,
pxs: f64,
) -> Result<()>
where
DB::ErrorType: 'static,
{
let ps = |v: f64| (v * pxs) as u32;
root.fill(&CHART_BG).map_err(|e| anyhow::anyhow!("{}", e))?;
let y_max = *histogram.iter().max().unwrap_or(&1) as f64;
let y_ceil = nice_ceil(y_max);
let title = "Coverage Histogram (0-50X)";
use plotters::style::text_anchor::{HPos, Pos, VPos};
let title_style = ("sans-serif", ps(24.0) as f64)
.into_font()
.style(FontStyle::Bold)
.color(&TITLE_COLOR)
.pos(Pos::new(HPos::Center, VPos::Top));
let subtitle_style = ("sans-serif", ps(16.0) as f64)
.into_font()
.color(&TITLE_COLOR)
.pos(Pos::new(HPos::Center, VPos::Top));
let (w, _h) = root.dim_in_pixel();
let center_x = w as i32 / 2;
let title_y = ps(10.0) as i32;
let subtitle_y = title_y + ps(28.0) as i32;
root.draw_text(title, &title_style, (center_x, title_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
root.draw_text(sample_name, &subtitle_style, (center_x, subtitle_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let chart_area = root.clone().margin(ps(55.0), ps(20.0), ps(20.0), ps(20.0));
let mut chart = ChartBuilder::on(&chart_area)
.x_label_area_size(ps(40.0))
.y_label_area_size(ps(70.0))
.build_cartesian_2d((-0.5f64)..50.5f64, 0.0..y_ceil)
.map_err(|e| anyhow::anyhow!("{}", e))?;
chart
.plotting_area()
.fill(&PLOT_BG)
.map_err(|e| anyhow::anyhow!("{}", e))?;
chart
.configure_mesh()
.x_desc("Coverage (X)")
.y_desc("Number of transcripts")
.x_labels(11) .label_style(("sans-serif", ps(15.0) as f64).into_font().color(&BLACK))
.axis_desc_style(("sans-serif", ps(16.0) as f64).into_font().color(&BLACK))
.x_label_formatter(&|v| {
let iv = v.round() as i32;
if (0..=50).contains(&iv) && (v - iv as f64).abs() < 0.01 {
format!("{iv}")
} else {
String::new()
}
})
.y_label_formatter(&|v| format_y_label(*v))
.axis_style(BLACK.stroke_width(1))
.light_line_style(GRIDLINE_COLOR.mix(0.0)) .bold_line_style(GRIDLINE_COLOR.stroke_width(1))
.set_all_tick_mark_size(ps(4.0))
.draw()
.map_err(|e| anyhow::anyhow!("{}", e))?;
let bar_width = 0.45; chart
.draw_series(histogram.iter().enumerate().map(|(i, &count)| {
let x0 = i as f64 - bar_width;
let x1 = i as f64 + bar_width;
Rectangle::new(
[(x0, 0.0), (x1, count as f64)],
ShapeStyle {
color: RGBColor(100, 100, 250).mix(0.59),
filled: true,
stroke_width: ps(1.0),
},
)
}))
.map_err(|e| anyhow::anyhow!("{}", e))?;
chart
.draw_series(histogram.iter().enumerate().map(|(i, &count)| {
let x0 = i as f64 - bar_width;
let x1 = i as f64 + bar_width;
Rectangle::new(
[(x0, 0.0), (x1, count as f64)],
ShapeStyle {
color: HISTOGRAM_BAR_BORDER.to_rgba(),
filled: false,
stroke_width: ps(1.0),
},
)
}))
.map_err(|e| anyhow::anyhow!("{}", e))?;
Ok(())
}
pub struct PieSlice<'a> {
pub label: &'a str,
pub value: f64,
pub color: RGBColor,
}
pub fn reads_genomic_origin_plot(
exonic: u64,
intronic: u64,
intergenic: u64,
sample_name: &str,
output_path: &Path,
) -> Result<()> {
let total = exonic + intronic + intergenic;
if total == 0 {
return Ok(());
}
let slices = vec![
PieSlice {
label: "Exonic",
value: exonic as f64,
color: PIE_COLOR_0,
},
PieSlice {
label: "Intergenic",
value: intergenic as f64,
color: PIE_COLOR_1,
},
PieSlice {
label: "Intronic",
value: intronic as f64,
color: PIE_COLOR_2,
},
];
pie_chart("Reads Genomic Origin", sample_name, &slices, output_path)
}
pub fn junction_analysis_plot(
known: u64,
partly_known: u64,
novel: u64,
sample_name: &str,
output_path: &Path,
) -> Result<()> {
let total = known + partly_known + novel;
if total == 0 {
return Ok(());
}
let slices = vec![
PieSlice {
label: "Known",
value: known as f64,
color: PIE_COLOR_0,
},
PieSlice {
label: "Partly known",
value: partly_known as f64,
color: PIE_COLOR_1,
},
PieSlice {
label: "Novel",
value: novel as f64,
color: PIE_COLOR_2,
},
];
pie_chart("Junction Analysis", sample_name, &slices, output_path)
}
fn pie_chart(
title: &str,
sample_name: &str,
slices: &[PieSlice<'_>],
output_path: &Path,
) -> Result<()> {
{
let root = BitMapBackend::new(output_path, (s(WIDTH), s(HEIGHT))).into_drawing_area();
render_pie_chart(&root, title, sample_name, slices, SCALE as f64)?;
root.present().context("Failed to write pie chart PNG")?;
}
let svg_path = output_path.with_extension("svg");
{
let root = SVGBackend::new(&svg_path, (WIDTH, HEIGHT)).into_drawing_area();
render_pie_chart(&root, title, sample_name, slices, 1.0)?;
root.present().context("Failed to write pie chart SVG")?;
}
debug!("Wrote pie chart: {}", output_path.display());
Ok(())
}
fn render_pie_chart<DB: DrawingBackend>(
root: &DrawingArea<DB, plotters::coord::Shift>,
title: &str,
sample_name: &str,
slices: &[PieSlice<'_>],
pxs: f64,
) -> Result<()>
where
DB::ErrorType: 'static,
{
let ps = |v: f64| (v * pxs) as i32;
let psu = |v: f64| (v * pxs) as u32;
root.fill(&CHART_BG).map_err(|e| anyhow::anyhow!("{}", e))?;
let (w, h) = root.dim_in_pixel();
let w = w as i32;
let h = h as i32;
use plotters::style::text_anchor::{HPos, Pos, VPos};
let title_style = ("sans-serif", ps(24.0) as f64)
.into_font()
.style(FontStyle::Bold)
.color(&TITLE_COLOR)
.pos(Pos::new(HPos::Center, VPos::Top));
let subtitle_style = ("sans-serif", ps(16.0) as f64)
.into_font()
.color(&TITLE_COLOR)
.pos(Pos::new(HPos::Center, VPos::Top));
let center_x = w / 2;
let title_y = ps(12.0);
let subtitle_y = title_y + ps(28.0);
root.draw_text(title, &title_style, (center_x, title_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
root.draw_text(sample_name, &subtitle_style, (center_x, subtitle_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let frame_margin = ps(8.0);
let frame_top = subtitle_y + ps(12.0);
let frame_bottom = h - frame_margin;
let frame_left = frame_margin;
let frame_right = w - frame_margin;
root.draw(&Rectangle::new(
[(frame_left, frame_top), (frame_right, frame_bottom)],
ShapeStyle {
color: CHART_BG.to_rgba(),
filled: true,
stroke_width: 0,
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let inner_pad = ps(1.0);
root.draw(&Rectangle::new(
[
(frame_left + inner_pad, frame_top + inner_pad),
(frame_right - inner_pad, frame_bottom - inner_pad),
],
ShapeStyle {
color: PLOT_BG.to_rgba(),
filled: true,
stroke_width: 0,
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
root.draw(&Rectangle::new(
[(frame_left, frame_top), (frame_right, frame_bottom)],
ShapeStyle {
color: BLACK.to_rgba(),
filled: false,
stroke_width: 1,
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let plot_margin = ps(30.0);
let plot_top = frame_top + ps(10.0);
let plot_w = frame_right - frame_left - 2 * plot_margin;
let plot_h = frame_bottom - plot_top - plot_margin;
let cx = frame_left + (frame_right - frame_left) / 2;
let cy = plot_top + plot_h / 2;
let radius = (plot_w.min(plot_h) as f64 * 0.38) as i32;
let total: f64 = slices.iter().map(|s| s.value).sum();
if total <= 0.0 {
return Ok(());
}
let shadow_offset = ps(4.0);
let n_shadow_pts = 360;
let shadow_points: Vec<(i32, i32)> = (0..=n_shadow_pts)
.map(|i| {
let angle = 2.0 * std::f64::consts::PI * i as f64 / n_shadow_pts as f64;
(
cx + shadow_offset + (radius as f64 * angle.cos()) as i32,
cy + shadow_offset + (radius as f64 * angle.sin()) as i32,
)
})
.collect();
root.draw(&Polygon::new(
shadow_points,
ShapeStyle {
color: PIE_SHADOW_COLOR.mix(0.5),
filled: true,
stroke_width: 0,
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let mut start_angle = std::f64::consts::FRAC_PI_2;
let mut label_info: Vec<(f64, String, RGBColor)> = Vec::new();
for slice in slices {
let pct = slice.value / total * 100.0;
if pct <= 0.0 {
continue;
}
let sweep = slice.value / total * 2.0 * std::f64::consts::PI;
let n_segments = ((sweep / 0.02) as usize).max(4);
let mut points = Vec::with_capacity(n_segments + 3);
points.push((cx, cy));
for j in 0..=n_segments {
let angle = start_angle - j as f64 * sweep / n_segments as f64;
points.push((
cx + (radius as f64 * angle.cos()) as i32,
cy - (radius as f64 * angle.sin()) as i32,
));
}
points.push((cx, cy));
root.draw(&Polygon::new(points, slice.color.filled()))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let mid_angle = start_angle - sweep / 2.0;
let label_text = format!(
"{} - {}%",
slice.label,
trim_trailing_zeros(&format!("{:.2}", pct))
);
label_info.push((mid_angle, label_text, slice.color));
start_angle -= sweep;
}
for (mid_angle, label_text, _color) in &label_info {
let edge_x = cx + (radius as f64 * mid_angle.cos()) as i32;
let edge_y = cy - (radius as f64 * mid_angle.sin()) as i32;
let label_radius = radius as f64 * 1.15;
let label_x = cx + (label_radius * mid_angle.cos()) as i32;
let label_y = cy - (label_radius * mid_angle.sin()) as i32;
root.draw(&PathElement::new(
vec![(edge_x, edge_y), (label_x, label_y)],
ShapeStyle {
color: RGBColor(80, 80, 80).to_rgba(),
filled: false,
stroke_width: psu(1.0),
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let font_size = ps(14.0) as f64;
let text_w = estimate_text_width(label_text, font_size) as i32;
let text_h = ps(14.0);
let pad_x = ps(3.0);
let pad_y = ps(2.0);
let (box_x, box_y) = if mid_angle.cos() >= 0.0 {
(label_x + ps(3.0), label_y - text_h / 2 - pad_y)
} else {
(
label_x - text_w - 2 * pad_x - ps(3.0),
label_y - text_h / 2 - pad_y,
)
};
root.draw(&Rectangle::new(
[
(box_x, box_y),
(box_x + text_w + 2 * pad_x, box_y + text_h + 2 * pad_y),
],
ShapeStyle {
color: PIE_LABEL_BG.to_rgba(),
filled: true,
stroke_width: psu(1.0),
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
root.draw(&Rectangle::new(
[
(box_x, box_y),
(box_x + text_w + 2 * pad_x, box_y + text_h + 2 * pad_y),
],
ShapeStyle {
color: PIE_LABEL_BORDER.to_rgba(),
filled: false,
stroke_width: psu(1.0),
},
))
.map_err(|e| anyhow::anyhow!("{}", e))?;
let text_style = ("sans-serif", font_size)
.into_font()
.color(&BLACK)
.into_text_style(root);
root.draw_text(label_text, &text_style, (box_x + pad_x, box_y + pad_y))
.map_err(|e| anyhow::anyhow!("{}", e))?;
}
Ok(())
}
fn format_y_label(val: f64) -> String {
if val == 0.0 {
return "0".to_string();
}
if val.abs() < 1.0 {
let s = format!("{:.2}", val);
return s.trim_end_matches('0').trim_end_matches('.').to_string();
}
let n = val.round() as i64;
let s = n.to_string();
let bytes = s.as_bytes();
let len = bytes.len();
let negative = bytes[0] == b'-';
let start = if negative { 1 } else { 0 };
let digit_len = len - start;
if digit_len <= 3 {
return s;
}
let mut result = String::with_capacity(len + (digit_len - 1) / 3);
for (i, &b) in bytes.iter().enumerate() {
if i > start && (len - i).is_multiple_of(3) {
result.push(',');
}
result.push(b as char);
}
result
}
fn nice_ceil(val: f64) -> f64 {
if val <= 0.0 {
return 1.0;
}
let magnitude = 10.0f64.powf(val.log10().floor());
let normalized = val / magnitude;
let nice = if normalized <= 1.0 {
1.0
} else if normalized <= 1.5 {
1.5
} else if normalized <= 2.0 {
2.0
} else if normalized <= 2.5 {
2.5
} else if normalized <= 3.0 {
3.0
} else if normalized <= 4.0 {
4.0
} else if normalized <= 5.0 {
5.0
} else if normalized <= 7.5 {
7.5
} else {
10.0
};
nice * magnitude
}
fn estimate_text_width(text: &str, font_size: f64) -> f64 {
text.len() as f64 * font_size * 0.6
}
fn trim_trailing_zeros(s: &str) -> String {
if !s.contains('.') {
return s.to_string();
}
let trimmed = s.trim_end_matches('0');
if trimmed.ends_with('.') {
format!("{trimmed}0")
} else {
trimmed.to_string()
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_format_y_label() {
assert_eq!(format_y_label(0.0), "0");
assert_eq!(format_y_label(500.0), "500");
assert_eq!(format_y_label(1000.0), "1,000");
assert_eq!(format_y_label(675000.0), "675,000");
assert_eq!(format_y_label(1234567.0), "1,234,567");
}
#[test]
fn test_nice_ceil() {
assert_eq!(nice_ceil(675000.0), 750000.0);
assert_eq!(nice_ceil(87500.0), 100000.0);
assert_eq!(nice_ceil(210000.0), 250000.0);
assert_eq!(nice_ceil(0.2), 0.2); assert_eq!(nice_ceil(0.25), 0.25); assert!((nice_ceil(0.26) - 0.3).abs() < 1e-15); }
#[test]
fn test_trim_trailing_zeros() {
assert_eq!(trim_trailing_zeros("76.00"), "76.0");
assert_eq!(trim_trailing_zeros("18.38"), "18.38");
assert_eq!(trim_trailing_zeros("5.20"), "5.2");
assert_eq!(trim_trailing_zeros("100.00"), "100.0");
}
}