RUST程序设计之十三 fastqc软件RUST代码实现

[package]name = "fastqc_rs"version = "0.1.0"edition = "2021"[dependencies]plotters = "0.3"rand = "0.8"clap = { version = "4.5", features = ["derive"] }anyhow = "1.0"rand_distr = "0.4" flate2 = "1.0"

use anyhow::{Context, Result};use clap::Parser;use flate2::read::MultiGzDecoder;use plotters::prelude::*;use std::fs::{self, File};use std::io::{BufRead, BufReader, BufWriter, Write};use std::path::Path;#[derive(Parser, Debug)]#[command(author, version, about = "Fast FASTQ filter and Quality Plotter")]struct Args {    /// Input R1 fastq file    #[arg(short, long)]    r1: String,    /// Input R2 fastq file    #[arg(short, long)]    r2: String,    /// Output directory    #[arg(short, long)]    out_dir: String,    /// Max allowed 'N' bases per read    #[arg(short, long, default_value_t = 30)]    max_n: usize,}struct BoxplotStats {    label: String,    q10: f64,    q1: f64,    median: f64,    q3: f64,    q90: f64,}fncalculate_quantile(hist: &[usize], total: usize, q: f64) -> f64{    if total == 0 { return 0.0; }    let target = (total as f64 * q) as usize;    let mut cumsum = 0;    for(score, &count) in hist.iter().enumerate() {        cumsum += count;        if cumsum >= target { return score as f64; }    }    0.0}fncount_char(slice: &[u8], target: u8) -> usize{    slice.iter().filter(|&&c| c == target).count()}fnopen_file<P: AsRef<Path>>(path: P) -> Result<Box<dynBufRead>> {    let f = File::open(&path).with_context(|| format!("无法打开文件: {:?}", path.as_ref()))?;    let reader: Box<dyn BufRead> = if path.as_ref().extension().map_or(false, |ext| ext == "gz") {        Box::new(BufReader::new(MultiGzDecoder::new(f)))    } else {        Box::new(BufReader::new(f))    };    Ok(reader)}fndraw_fastqc_plot(data: &[BoxplotStats], title: &str, out_path: &Path) -> Result<()> {    if data.is_empty() {        println!("警告: 绘图数据为空！");        return Ok(());    }    let root = BitMapBackend::new(out_path, (1024, 600)).into_drawing_area();    root.fill(&WHITE)?;    // 核心修复：将 X 轴放大 10 倍，使用纯 i32 坐标系，宽度分配更精准    let n_points = data.len() as i32;    let x_max = n_points * 10;    let mut chart = ChartBuilder::on(&root)        .caption(title, ("sans-serif", 20).into_font())        .margin(10)        .x_label_area_size(40)        .y_label_area_size(40)        .build_cartesian_2d(0i32..x_max, 0f64..42f64)?;    chart.configure_mesh()        .disable_x_mesh()        .disable_y_mesh()        .x_labels(n_points as usize)        // 根据放大后的坐标反推标签        .x_label_formatter(&|x| {            let idx = (*x / 10) as usize;            let safe_idx = idx.min(data.len().saturating_sub(1));            data[safe_idx].label.clone()        })        .x_desc("Position in read (bp)")        .y_desc("Quality score")        .draw()?;    // 1. 画背景色带 (红、黄、绿)    chart.draw_series(std::iter::once(Rectangle::new(        [(0, 0.0), (x_max, 20.0)], RGBColor(255, 0, 0).mix(0.2).filled(),    )))?;    chart.draw_series(std::iter::once(Rectangle::new(        [(0, 20.0), (x_max, 28.0)], RGBColor(255, 204, 0).mix(0.2).filled(),    )))?;    chart.draw_series(std::iter::once(Rectangle::new(        [(0, 28.0), (x_max, 42.0)], RGBColor(0, 255, 0).mix(0.2).filled(),    )))?;    // 2. 画斑马纹背景    for i in (0..n_points).step_by(2) {        chart.draw_series(std::iter::once(Rectangle::new(            [(i * 10, 0.0), ((i + 1) * 10, 42.0)], RGBColor(200, 200, 200).mix(0.15).filled(),        )))?;    }    // 3. 画箱线图    for(i, stat) in data.iter().enumerate() {        let x_center = (i as i32) * 10;        let box_width = 4; // 在放大 10 倍的坐标系里，宽度设 4 比较合适        // 黑须 (Q10 - Q90)        chart.draw_series(std::iter::once(PathElement::new(            vec![(x_center, stat.q10), (x_center, stat.q90)],            BLACK.stroke_width(1),        )))?;        // 黄色箱体 (Q1 - Q3)        chart.draw_series(std::iter::once(Rectangle::new(            [(x_center - box_width, stat.q1), (x_center + box_width, stat.q3)],            RGBColor(255, 255, 0).filled().stroke_width(1),        )))?;        // 红色中位数横线        chart.draw_series(std::iter::once(PathElement::new(            vec![(x_center - box_width, stat.median), (x_center + box_width, stat.median)],            RED.stroke_width(2),        )))?;    }    // 4. 画蓝色中位数折线    let medians: Vec<(i32, f64)> = data.iter().enumerate().map(|(i, s)| ((i as i32) * 10, s.median)).collect();    chart.draw_series(LineSeries::new(medians, &BLUE))?;    // 5. 画顶部红色虚线    for i in (0..x_max).step_by(20) {        chart.draw_series(std::iter::once(PathElement::new(            vec![(i, 40.0), (i + 10, 40.0)], RED.stroke_width(2),        )))?;    }    root.present()?;    Ok(())}fnmain() -> Result<()> {    let args = Args::parse();    let out_dir = Path::new(&args.out_dir);    if !out_dir.exists() {        fs::create_dir_all(out_dir).context("无法创建输出目录")?;    }    let r1_path = Path::new(&args.r1);    let r2_path = Path::new(&args.r2);    let r1_stem = r1_path.file_stem().unwrap_or_default().to_string_lossy();    let out_r1_path = out_dir.join(format!("{}.clean.fq", r1_stem));    let out_r2_path = out_dir.join(format!("{}_R2.clean.fq", r1_stem));    let plot_path = out_dir.join(format!("{}_quality.png", r1_stem));    println!("开始处理...");    let f1 = open_file(r1_path)?;    let f2 = open_file(r2_path)?;    let mut writer1 = BufWriter::new(File::create(&out_r1_path)?);    let mut writer2 = BufWriter::new(File::create(&out_r2_path)?);    let mut reader1 = BufReader::new(f1);    let mut reader2 = BufReader::new(f2);    let mut buf1 = Vec::with_capacity(1024);    let mut buf2 = Vec::with_capacity(1024);    let mut qual_stats_r1: Vec<Vec<usize>> = Vec::new();    let mut count_total = 0;    let mut count_kept = 0;    loop {        let len1 = reader1.read_until(b'\n', &mut buf1)?;        let len2 = reader2.read_until(b'\n', &mut buf2)?;        if len1 == 0 || len2 == 0 { break; }        count_total += 1;        let mut pair1_rows = vec![buf1.clone()];        let mut pair2_rows = vec![buf2.clone()];        let mut valid = true;        for line_idx in 0..3 {            buf1.clear();            buf2.clear();            reader1.read_until(b'\n', &mut buf1)?;            reader2.read_until(b'\n', &mut buf2)?;            if line_idx == 0 {                if count_char(&buf1, b'N') > args.max_n || count_char(&buf2, b'N') > args.max_n {                    valid = false;                }            }            if line_idx == 2 && valid {                let qual_str = buf1.trim_ascii_end();                if qual_stats_r1.len() < qual_str.len() {                    qual_stats_r1.resize(qual_str.len(), vec![0; 94]);                }                for(pos, &ascii) in qual_str.iter().enumerate() {                    if ascii >= 33 {                        let score = (ascii - 33) as usize;                        if score < 94 { qual_stats_r1[pos][score] += 1; }                    }                }            }            pair1_rows.push(buf1.clone());            pair2_rows.push(buf2.clone());        }        if valid {            count_kept += 1;            for row in pair1_rows { writer1.write_all(&row)?; }            for row in pair2_rows { writer2.write_all(&row)?; }        }    }    writer1.flush()?;    writer2.flush()?;    println!("过滤完成！总计: {}, 保留: {}, 丢弃: {}", count_total, count_kept, count_total - count_kept);    // 生成绘图数据    if !qual_stats_r1.is_empty() {        println!("已统计到序列长度维度: {}", qual_stats_r1.len());        let mut plot_data = Vec::new();        for(pos, hist) in qual_stats_r1.iter().enumerate() {            // 降采样：取前 10 个点，之后每 5 个取一个点            if pos < 10 || pos % 5 == 4 {                let total: usize = hist.iter().sum();                let label = if pos < 10 { format!("{}", pos + 1) } else { format!("{}-{}", pos - 3, pos + 1) };                plot_data.push(BoxplotStats {                    label,                    q10: calculate_quantile(hist, total, 0.10),                    q1: calculate_quantile(hist, total, 0.25),                    median: calculate_quantile(hist, total, 0.50),                    q3: calculate_quantile(hist, total, 0.75),                    q90: calculate_quantile(hist, total, 0.90),                });            }        }        println!("生成绘图点数量: {}", plot_data.len());        println!("第一个点的中位数质量分: {}", plot_data[0].median);        draw_fastqc_plot(&plot_data, "Quality scores across all bases (Clean R1)", &plot_path)?;        println!("质量图表已保存至: {:?}", plot_path);    }    Ok(())}