ciri_toolkit/cli/
merge.rs

1//! CLI implementation for `ciri-merge` shared BSJ catalog construction.
2//!
3//! The command is intentionally narrower than the main `ciri` pipeline: it
4//! consumes only completed first-pass `.out` files and writes a deterministic
5//! BED6 shared circRNA catalog. It does not read `.segments` or original
6//! BAM/SAM files, so the output can be used as the external `--circ` catalog
7//! for a later sample-specific second pass.
8
9use anyhow::{anyhow, bail, Context, Result};
10use clap::Parser;
11use std::cmp::Ordering;
12use std::collections::{BTreeMap, HashMap};
13use std::fs::File;
14use std::io::{BufRead, BufReader, BufWriter, Write};
15use std::path::Path;
16
17/// Parsed command-line arguments for `ciri-merge`.
18#[derive(Parser, Debug)]
19#[command(
20    author,
21    version,
22    about = "Merge first-pass CIRI .out files into a circRNA bed file",
23    long_about = None
24)]
25struct Args {
26    /// First-pass CIRI `.out` files to merge.
27    #[arg(short = 'i', long = "in", required = true, num_args = 1..)]
28    inputs: Vec<String>,
29
30    /// Output circRNA bed file path.
31    #[arg(long = "out", short = 'o')]
32    output: String,
33}
34
35/// One sample entry derived from a first-pass `.out` path.
36#[derive(Debug, Clone)]
37struct MergeSample {
38    sample_id: String,
39    out_path: String,
40}
41
42/// Multi-sample grouping key for one circRNA/BSJ boundary.
43#[derive(Debug, Clone, Eq, PartialEq, Hash)]
44struct CircKey {
45    chrom: String,
46    start: i64,
47    end: i64,
48    strand: String,
49}
50
51/// Stable sort key that keeps common chromosome names in genomic order.
52#[derive(Debug, Clone, Eq, PartialEq)]
53struct ChromSortKey {
54    group: u8,
55    rank: u32,
56    suffix: String,
57}
58
59impl Ord for ChromSortKey {
60    fn cmp(&self, other: &Self) -> Ordering {
61        self.group
62            .cmp(&other.group)
63            .then_with(|| self.rank.cmp(&other.rank))
64            .then_with(|| self.suffix.cmp(&other.suffix))
65    }
66}
67
68impl PartialOrd for ChromSortKey {
69    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
70        Some(self.cmp(other))
71    }
72}
73
74/// Per-sample evidence contributing to one shared circRNA boundary.
75#[derive(Debug, Clone)]
76struct SampleCircSupport {
77    source_circ_id: String,
78    junction_reads: i64,
79}
80
81/// Aggregated multi-sample support for one [`CircKey`].
82#[derive(Debug, Default)]
83struct CohortCircSupport {
84    by_sample: BTreeMap<String, SampleCircSupport>,
85}
86
87/// User-facing summary returned by the merge implementation.
88#[derive(Debug, Default)]
89struct MergeSummary {
90    samples: usize,
91    cohort_circs: usize,
92    input_rows: usize,
93    total_bsj_reads: i64,
94}
95
96/// Parses arguments, runs the merge, and prints a compact summary.
97pub fn main() -> Result<()> {
98    let args = Args::parse();
99    let summary = run_merge(&args.inputs, &args.output)?;
100    eprintln!(
101        "Merged {} samples: {} circRNAs from {} rows, {} BSJ reads",
102        summary.samples, summary.cohort_circs, summary.input_rows, summary.total_bsj_reads
103    );
104    Ok(())
105}
106
107/// Runs shared circRNA catalog construction from first-pass `.out` files.
108///
109/// This function keeps only one compact support record per `(sample, circRNA
110/// boundary)` in memory. It is therefore bounded by the number of detected
111/// circRNA sites rather than by read count or segment rows.
112fn run_merge(inputs: &[String], output_path: &str) -> Result<MergeSummary> {
113    let samples = samples_from_inputs(inputs)?;
114    if samples.is_empty() {
115        bail!("no first-pass .out inputs were provided");
116    }
117
118    let mut cohort: HashMap<CircKey, CohortCircSupport> = HashMap::new();
119    let mut input_rows = 0usize;
120    let mut total_bsj_reads = 0i64;
121    for sample in &samples {
122        let rows = read_out_file(sample, &mut cohort)?;
123        input_rows += rows.rows;
124        total_bsj_reads += rows.total_bsj_reads;
125    }
126
127    let mut entries: Vec<_> = cohort.into_iter().collect();
128    entries.sort_by(|(left, _), (right, _)| compare_circ_keys(left, right));
129
130    write_circ_bed_catalog(output_path, &entries)?;
131
132    Ok(MergeSummary {
133        samples: samples.len(),
134        cohort_circs: entries.len(),
135        input_rows,
136        total_bsj_reads,
137    })
138}
139
140/// Builds merge sample entries from CLI input paths.
141///
142/// Sample IDs are derived from the basename after stripping one trailing `.out`
143/// suffix. This keeps the common `sample1.out` and `sample1.ciri.out` cases
144/// stable while avoiding a manifest file for the first implementation.
145fn samples_from_inputs(inputs: &[String]) -> Result<Vec<MergeSample>> {
146    let mut samples = Vec::with_capacity(inputs.len());
147    let mut seen: HashMap<String, usize> = HashMap::new();
148    for input in inputs {
149        if !Path::new(input).is_file() {
150            bail!("missing first-pass .out file: {input}");
151        }
152        let mut sample_id = sample_id_from_out_path(input)?;
153        let count = seen.entry(sample_id.clone()).or_insert(0);
154        if *count > 0 {
155            sample_id = format!("{}_{}", sample_id, *count + 1);
156        }
157        *count += 1;
158        samples.push(MergeSample {
159            sample_id,
160            out_path: input.clone(),
161        });
162    }
163    Ok(samples)
164}
165
166/// Derives a stable sample ID from one `.out` path.
167fn sample_id_from_out_path(path: &str) -> Result<String> {
168    let name = Path::new(path)
169        .file_name()
170        .and_then(|name| name.to_str())
171        .ok_or_else(|| anyhow!("cannot derive sample_id from path: {path}"))?;
172    let sample = name.strip_suffix(".out").unwrap_or(name);
173    if sample.is_empty() {
174        bail!("empty sample_id derived from path: {path}");
175    }
176    Ok(sample.to_string())
177}
178
179#[derive(Debug, Default)]
180struct OutReadSummary {
181    rows: usize,
182    total_bsj_reads: i64,
183}
184
185/// Streams one CIRI `.out` file into the shared support map.
186fn read_out_file(
187    sample: &MergeSample,
188    cohort: &mut HashMap<CircKey, CohortCircSupport>,
189) -> Result<OutReadSummary> {
190    let file = File::open(&sample.out_path)
191        .with_context(|| format!("open first-pass .out {}", sample.out_path))?;
192    let reader = BufReader::new(file);
193    let mut summary = OutReadSummary::default();
194    for (line_idx, line_result) in reader.lines().enumerate() {
195        let line = line_result?;
196        if line.trim().is_empty() || line.starts_with('#') || line.starts_with("circRNA_ID\t") {
197            continue;
198        }
199        let fields: Vec<&str> = line.split('\t').collect();
200        if fields.len() < 11 {
201            bail!(
202                "{}:{} has {} columns, expected at least 11",
203                sample.out_path,
204                line_idx + 1,
205                fields.len()
206            );
207        }
208        let circ_id = fields[0].to_string();
209        let key = CircKey {
210            chrom: fields[1].to_string(),
211            start: fields[2].parse().with_context(|| {
212                format!(
213                    "parse circRNA_start at {}:{}",
214                    sample.out_path,
215                    line_idx + 1
216                )
217            })?,
218            end: fields[3].parse().with_context(|| {
219                format!("parse circRNA_end at {}:{}", sample.out_path, line_idx + 1)
220            })?,
221            strand: normalize_strand(fields[10]),
222        };
223        let junction_reads: i64 = fields[4].parse().with_context(|| {
224            format!(
225                "parse #junction_reads at {}:{}",
226                sample.out_path,
227                line_idx + 1
228            )
229        })?;
230        if junction_reads < 0 {
231            bail!(
232                "{}:{} has negative #junction_reads: {}",
233                sample.out_path,
234                line_idx + 1,
235                junction_reads
236            );
237        }
238        let support = cohort.entry(key).or_default();
239        support
240            .by_sample
241            .entry(sample.sample_id.clone())
242            .and_modify(|existing| {
243                existing.junction_reads += junction_reads;
244                if existing.source_circ_id != circ_id {
245                    existing.source_circ_id.push(',');
246                    existing.source_circ_id.push_str(&circ_id);
247                }
248            })
249            .or_insert(SampleCircSupport {
250                source_circ_id: circ_id,
251                junction_reads,
252            });
253        summary.rows += 1;
254        summary.total_bsj_reads += junction_reads;
255    }
256    Ok(summary)
257}
258
259/// Normalizes missing or unknown CIRI strand labels for stable shared keys.
260fn normalize_strand(strand: &str) -> String {
261    match strand {
262        "+" | "-" => strand.to_string(),
263        _ => "NA".to_string(),
264    }
265}
266
267/// Writes the BED6 shared circRNA catalog used by the future `--circ` second pass.
268///
269/// The first six columns follow BED6 conventions with 0-based half-open
270/// coordinates. The score uses capped total BSJ support so the file remains
271/// directly usable by BED-aware tools.
272fn write_circ_bed_catalog(path: &str, entries: &[(CircKey, CohortCircSupport)]) -> Result<()> {
273    let mut writer = BufWriter::with_capacity(1024 * 1024, File::create(path)?);
274    for (key, support) in entries {
275        let total_bsj_reads = total_bsj_reads(support);
276        writeln!(
277            writer,
278            "{}\t{}\t{}\t{}\t{}\t{}",
279            key.chrom,
280            key.start - 1,
281            key.end,
282            cohort_circ_id(key),
283            bed_score(total_bsj_reads),
284            key.strand
285        )?;
286    }
287    writer.flush()?;
288    Ok(())
289}
290
291/// Converts total BSJ read support into the standard BED score range.
292fn bed_score(total_bsj_reads: i64) -> i64 {
293    total_bsj_reads.clamp(0, 1000)
294}
295
296/// Returns the stable shared circRNA identifier for one merged BSJ boundary.
297fn cohort_circ_id(key: &CircKey) -> String {
298    format!("{}:{}|{}:{}", key.chrom, key.start, key.end, key.strand)
299}
300
301/// Sums BSJ read support across samples for one shared circRNA.
302fn total_bsj_reads(support: &CohortCircSupport) -> i64 {
303    support
304        .by_sample
305        .values()
306        .map(|support| support.junction_reads)
307        .sum()
308}
309
310/// Compares shared circRNA keys in coordinate order.
311fn compare_circ_keys(left: &CircKey, right: &CircKey) -> Ordering {
312    chrom_sort_key(&left.chrom)
313        .cmp(&chrom_sort_key(&right.chrom))
314        .then_with(|| left.start.cmp(&right.start))
315        .then_with(|| left.end.cmp(&right.end))
316        .then_with(|| left.strand.cmp(&right.strand))
317}
318
319/// Builds a human-genome-friendly chromosome key while keeping generic contigs deterministic.
320fn chrom_sort_key(chrom: &str) -> ChromSortKey {
321    let core = chrom.strip_prefix("chr").unwrap_or(chrom);
322    if let Ok(rank) = core.parse::<u32>() {
323        return ChromSortKey {
324            group: 0,
325            rank,
326            suffix: String::new(),
327        };
328    }
329    match core {
330        "X" => ChromSortKey {
331            group: 0,
332            rank: 23,
333            suffix: String::new(),
334        },
335        "Y" => ChromSortKey {
336            group: 0,
337            rank: 24,
338            suffix: String::new(),
339        },
340        "M" | "MT" => ChromSortKey {
341            group: 0,
342            rank: 25,
343            suffix: String::new(),
344        },
345        _ => ChromSortKey {
346            group: 1,
347            rank: 0,
348            suffix: chrom.to_string(),
349        },
350    }
351}
352
353#[cfg(test)]
354mod tests {
355    use super::*;
356    use tempfile::tempdir;
357
358    #[test]
359    fn merge_out_files_builds_bed_catalog() {
360        let dir = tempdir().unwrap();
361        let out1 = dir.path().join("s1.out");
362        let out2 = dir.path().join("s2.out");
363        std::fs::write(
364            &out1,
365            concat!(
366                "circRNA_ID\tchr\tcircRNA_start\tcircRNA_end\t#junction_reads\tSM_MS_SMS\t#non_junction_reads\tjunction_reads_ratio\tcircRNA_type\tgene_id\tstrand\tjunction_reads_ID\tScore\n",
367                "chr1:10|20\tchr1\t10\t20\t3\t0_0_0\t0\t1\texon\tG1\t+\tr1,r2,r3\t3\n",
368                "chr2:30|40\tchr2\t30\t40\t5\t0_0_0\t0\t1\texon\tG2\t-\tr4\t5\n",
369            ),
370        )
371        .unwrap();
372        std::fs::write(
373            &out2,
374            concat!(
375                "circRNA_ID\tchr\tcircRNA_start\tcircRNA_end\t#junction_reads\tSM_MS_SMS\t#non_junction_reads\tjunction_reads_ratio\tcircRNA_type\tgene_id\tstrand\tjunction_reads_ID\tScore\n",
376                "chr1:10|20\tchr1\t10\t20\t7\t0_0_0\t0\t1\texon\tG1\t+\tr5\t7\n",
377            ),
378        )
379        .unwrap();
380        let output = dir.path().join("cohort.out");
381        let summary = run_merge(
382            &[
383                out1.to_string_lossy().to_string(),
384                out2.to_string_lossy().to_string(),
385            ],
386            output.to_str().unwrap(),
387        )
388        .unwrap();
389        assert_eq!(summary.samples, 2);
390        assert_eq!(summary.input_rows, 3);
391        assert_eq!(summary.cohort_circs, 2);
392        assert_eq!(summary.total_bsj_reads, 15);
393
394        let circ = std::fs::read_to_string(&output).unwrap();
395        assert!(!circ.contains("chrom\tchromStart\tchromEnd\tname\tscore\tstrand"));
396        assert!(circ.contains("chr1\t9\t20\tchr1:10|20:+\t10\t+"));
397        assert!(circ.contains("chr2\t29\t40\tchr2:30|40:-\t5\t-"));
398        assert!(!output.with_extension("out.bsj_matrix.tsv").exists());
399    }
400}
ciri_toolkit/cli/merge.rs

ciri_toolkit/cli/
merge.rs
