ambigviz 0.1.0

Identify and plot ambiguous nucleotide bases at given positions from a BAM file.
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
159

use std::str::FromStr;

use clap::{Parser, Subcommand};

#[derive(Parser, Debug)]
#[command(version, about = "Visulise ambigous bases in a bam file.")]
pub struct Cli {
    #[clap(subcommand)]
    pub command: Command,
}

#[derive(Debug, Subcommand)]
pub enum Command {
    /// Plot Depths
    Depth {
        // Path to input bam
        #[arg(required = true, value_parser(check_input_exists))]
        input: String,

        // SAMtools region string
        #[arg(required = false)]
        region: Option<String>,

        // Output
        #[arg(short = 'o', long = "output", default_value = "ambig")]
        output: String,
    },
    /// Plot Ambigous bases
    Ambig {
        // Path to input bam
        #[arg(required = true, value_parser(check_input_exists))]
        input: String,

        // SAMtools region string
        #[arg(required = false)]
        region: Option<String>,

        // Output
        #[arg(short = 'o', long = "output", default_value = "ambig")]
        output: String,

        // Threshold for ambigous bases
        #[arg(short = 't', long = "threshold", default_value = "0.1", value_parser(check_threshold_valid))]
        threshold: f64,

        // Threshold for base quality
        #[arg(short = 'q', long = "--min-BQ", default_value = "20")]
        base_quality_threshold: u8,

        // Threshold for map quality
        #[arg(short = 'Q', long = "--min-MQ", default_value = "60")]
        map_quality_threshold: u8,

        // Threshold for total depth
        #[arg(short = 'D', long = "--depth", default_value = "100")]
        depth_threshold: u32,

        // Threshold for depth of minor allele
        #[arg(short = 'd', long = "--minor-depth", default_value = "20")]
        minor_depth_threshold: u32,

        // Threshold for strand bias
        #[arg(short = 's', long = "--strand-bias", default_value = "0.1", value_parser(check_threshold_valid))]
        strand_bias_threshold: f64,

        // Do not include indels
        #[arg(long = "no-indel")]
        no_indel: bool,

        // Do not include labels
        #[arg(long = "no-label")]
        no_label: bool,

        // Output bed file instead of plot
        #[arg(long = "bed")]
        bed: bool,
    },
}

fn check_input_exists(s: &str) -> Result<String, String> {
    if std::path::Path::new(s).exists() {
        Ok(s.to_string())
    } else {
        Err(format!("File does not exist: {}", s))
    }
}

fn check_threshold_valid(s: &str) -> Result<f64, String> {
    let threshold = f64::from_str(s).unwrap();
    if threshold >= 0.0 && threshold <= 0.5 {
        Ok(threshold)
    } else {
        Err(format!("Threshold must be between 0 and 0.5"))
    }
}

pub fn parse_region(region: &str) -> (&str, u32, Option<u32>) {
    // Check if we have just chromosome or both chromosome and region
    if region.contains(':') {
        let region_parts: Vec<&str> = region.split(':').collect();

        let chrom = region_parts[0];
        let positions = region_parts[1];

        let start: u32;
        let end: Option<u32>;
        // Check if we have a range of positions or only one
        if positions.contains('-') {
            let position_parts: Vec<u32> = positions
                .split('-')
                .map(|x| u32::from_str(x).unwrap())
                .collect();
            start = position_parts[0];
            end = Option::from(position_parts[1]);
        } else {
            start = u32::from_str(positions).unwrap();
            end = None;
        }
        (chrom, start, end)
    } else {
        // We only have a chromosome
        (region, 1, None)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_parse_region_chrom_only() {
        let region = "chr1";
        let (chrom, start, end) = parse_region(region);

        assert_eq!(chrom, "chr1");
        assert_eq!(start, 1);
        assert_eq!(end, None);
    }

    #[test]
    fn test_parse_region_chrom_and_single_position() {
        let region = "chr2:100";
        let (chrom, start, end) = parse_region(region);

        assert_eq!(chrom, "chr2");
        assert_eq!(start, 100);
        assert_eq!(end, None);
    }

    #[test]
    fn test_parse_region_chrom_and_range() {
        let region = "chr3:200-300";
        let (chrom, start, end) = parse_region(region);

        assert_eq!(chrom, "chr3");
        assert_eq!(start, 200);
        assert_eq!(end, Some(300));
    }
}