1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

pub mod clipping;
pub mod cli;

use cli::Parser;
use clipping::ClipStat;

use log::{debug, info};
use rust_htslib::{bam, bam::Read};

/// Workflow to process an input bam file and write the pass-filter alignments
/// into a new bam file
///
/// # Arguments
/// - `in_bam`: input bam file
/// - `out_bam`: output bam file
/// - `inverse`: boolean flag to indicate whether we want to write out the failed-filter alignments only
/// - `both_end`: maximum fraction of total clipped bases relative to the read sequence length to consider as pass
/// - `left_side`: maximum fraction of of clipped bases on either side relative to the read sequence length to consider as pass
/// - `right_side`: maximum fraction of of clipped bases on either side relative to the read sequence length to consider as pass
///
/// # Examples
///
/// ```
/// use filter_clipped::run;
/// use rust_htslib::bam;
/// use rust_htslib::bam::Read;
/// fn count_bam(bam_file: String, expected_count: i32) {
///     // helper function to verify result
///     let mut bam_reader = bam::Reader::from_path(bam_file).unwrap();
///     let mut aln_count = 0;
///     for r in bam_reader.records() {
///         let _record = r.unwrap();
///         aln_count += 1;
///     }
///     assert_eq!(expected_count, aln_count);
/// }
///
/// let out_bam = "out.sam";
/// run(
///     "test/data/test.sam".to_string(),
///     out_bam.to_string(),
///     false,
///     0.1,
///     0.1,
///     0.1
///     );
/// count_bam(out_bam.to_string(), 6);
/// ```
pub fn run(
    in_bam: String,
    out_bam: String,
    inverse: bool,
    both_end: f64,
    left_side: f64,
    right_side: f64,
) -> Result<u8, String>{
    let mut out_count = 0;
    let mut in_count = 0;
    info!("Reading from alignment file: {}", in_bam);
    info!("Writing to alignment file: {}", out_bam);
    info!(
        "Thresholds: trailing clipped: {}, leading clipped: {}, total clipped: {}",
        right_side, left_side, both_end
    );
    let mut in_bam = match in_bam.eq("-") {
        true => bam::Reader::from_stdin().map_err(|e| e.to_string())?,
        _ => bam::Reader::from_path(&in_bam).map_err(|e| e.to_string())?,
    };
    let header = bam::Header::from_template(in_bam.header());

    let mut out_bam = match out_bam.eq("-") {
        true => bam::Writer::from_stdout(&header, bam::Format::Bam).map_err(|e| e.to_string())?,
        _ => bam::Writer::from_path(&out_bam, &header, bam::Format::Bam).map_err(|e| e.to_string())?,
    };

    for r in in_bam.records() {
        in_count += 1;
        let record = r.map_err(|e| e.to_string())?;
        let seq_len = record.seq().len() as f64;
        let cigar = record.cigar();

        let leading_clipped: Vec<i64> = vec![cigar.leading_softclips(), cigar.leading_hardclips()];
        let trailing_cliped: Vec<i64> =
            vec![cigar.trailing_softclips(), cigar.trailing_hardclips()];

        let clip_stat: ClipStat = ClipStat::new(leading_clipped, trailing_cliped);

        let keep: bool = clip_stat.total_fraction(seq_len)? < both_end
            && clip_stat.left_fraction(seq_len)? <= left_side
            && clip_stat.right_fraction(seq_len)? <= right_side;

        debug!("{:?} {}", clip_stat, seq_len);
        if (keep && !inverse) || (inverse && !keep) {
            out_bam.write(&record).map_err(|e| e.to_string())?;
            out_count += 1;
        }
    }
    info!(
        "Read {} alignments; Written {} alignments",
        in_count, out_count
    );
    Ok(0) // exit code 0
}

    
/// Just a wrapper function to read command line arguments and pass it to `run`
/// 
pub fn wrapper(){
    let args = cli::Command::parse();
    let result = run(
        args.in_bam, args.out_bam, args.inverse, args.both_end, args.left_side, args.right_side,
    );
    match result{
        Ok(_) => (),
        Err(err) => println!("{}", err),
    };
}

#[cfg(test)]
mod tests {
    use super::*;
    use rstest::rstest;
    use std::string::String;

    fn count_bam(bam_file: String, expected_count: i32) {
        let mut bam_reader = bam::Reader::from_path(bam_file).unwrap();
        let mut aln_count = 0;
        for r in bam_reader.records() {
            let _record = r.unwrap();
            aln_count += 1;
        }
        assert_eq!(expected_count, aln_count);
    }

    #[rstest]
    #[case(1, 0.2, 0.3, true, 0)]
    #[case(2, 0.2, 0.3, false, 9)]
    #[case(3, 0.1, 0.1, false, 6)]
    fn test_run(
        #[case] test_case: usize,
        #[case] max_both_end: f64,
        #[case] max_single_end: f64,
        #[case] inverse: bool,
        #[case] expected_count: i32,
    ) {
        let out_bam: &str = &format!("test/data/out_{}.bam", test_case);
        let result = run(
            "test/data/test.sam".to_string(),
            out_bam.to_string(),
            inverse,
            max_both_end,
            max_single_end,
            max_single_end,
        ).unwrap();
        assert_eq!(result, 0);
        count_bam(out_bam.to_string(), expected_count);
    }
}