xsra 0.2.29

A performant and storage-efficient CLI tool to extract sequences from an SRA archive with support for FASTA, FASTQ, and BINSEQ outputs.
Documentation 
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

use std::path::Path;

use anyhow::Result;
use log::{debug, info};
use ncbi_vdb_sys::{SegmentType, SraReader};

use crate::{
    cli::{DescribeOptions, InputOptions},
    prefetch::identify_url,
};

mod stats;
use stats::DescribeStats;

fn calculate_average_quality(qual: &[u8]) -> f64 {
    // PHRED33 has an ASCII offset of 33
    const PHRED33_OFFSET: u8 = 33;

    if qual.is_empty() {
        return 0.0;
    }

    let total_score: u32 = qual.iter().map(|byte| (byte - PHRED33_OFFSET) as u32).sum();

    total_score as f64 / qual.len() as f64
}

pub fn describe_inner(accession: &str, skip: usize, limit: usize) -> Result<DescribeStats> {
    let reader = SraReader::new(accession)?;
    let num_spots = reader.stop();

    let l_bound = skip.max(1);
    let r_bound = (l_bound + limit).min(num_spots as usize);

    let mut num_segments = 0;
    let mut segment_types = Vec::new();
    let mut segment_lengths = Vec::new();
    let mut segment_qualities = Vec::new();
    for record in reader.into_range_iter(l_bound as i64, r_bound as u64)? {
        let record = record?;
        for segment in record.into_iter() {
            // calculate mean quality
            let mean_quality = calculate_average_quality(segment.qual());

            // update segment stats
            num_segments = num_segments.max(segment.sid());

            if segment_types.len() <= segment.sid() {
                segment_types.resize(segment.sid() + 1, SegmentType::Technical);
                segment_lengths.resize(segment.sid() + 1, Vec::new());
                segment_qualities.resize(segment.sid() + 1, Vec::new());
            }
            segment_types[segment.sid()] = segment.ty();
            segment_lengths[segment.sid()].push(segment.len() as f64);
            segment_qualities[segment.sid()].push(mean_quality);
        }
    }

    let stats = DescribeStats::new(
        segment_types,
        segment_lengths,
        segment_qualities,
        r_bound - l_bound,
        l_bound,
        r_bound,
        num_spots as usize,
    );

    Ok(stats)
}

pub fn describe(input: &InputOptions, opts: &DescribeOptions) -> Result<()> {
    let accession = if !Path::new(&input.accession).exists() {
        info!(accession = input.accession.as_str(); "Identifying SRA data URL for accession");
        let runtime = tokio::runtime::Runtime::new()?;
        let url = runtime.block_on(identify_url(&input.accession, &input.options))?;
        info!(url = url.as_str(); "Streaming SRA records from URL");
        url
    } else {
        debug!(path = input.accession.as_str(); "Using local SRA file");
        input.accession.to_string()
    };
    let stats = describe_inner(&accession, opts.skip, opts.limit)?;
    stats.pprint(&mut std::io::stdout())?;

    Ok(())
}

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

    // describe_inner tests
    #[test]
    fn describe_inner_handles_invalid_sra_file() {
        // Create an empty temporary file.
        // This is not a valid SRA file.
        let temp_sra = NamedTempFile::new().unwrap();
        let path = temp_sra.path().to_str().unwrap();

        let result = describe_inner(path, 0, 10);

        assert!(
            result.is_err(),
            "describe_inner should fail when given a non-SRA file"
        );
    }
}