seqtk_rs/
nc_comp.rs

1use crate::bed::BedMap;
2use crate::dna::SeqComp;
3use crate::io_utils::{FaReader, FqReader, Output};
4use crate::record::RecordType;
5
6/// Parses FASTQ file and compute the statistic w/o masked sequences.
7/// Outputs the results to [`std::io::stdout()`].
8///The output columns are:
9/// - `#seq`: Number of reads
10/// - `#bases`: Number of bases
11/// - `#A`, `#C`, `#G`, `#T`: Number of each nucleotide
12/// - `#2`: Number of `R`, `Y`, `S`, `W`, `K`, `M`
13/// - `#3`: Number of `B`, `D`, `H`, `V`
14/// - `#4`: Number of `N`
15/// - `CG`: Number of `CG` on the template strand
16/// - `GC`: Number of `GC` on the template strand
17///
18/// # Arguments
19///
20/// * `path` - FASTQ path
21/// * `exclude_masked` - If true, masked sequences (e.g., lowercases or char other then IUPAC code) will be excluded from the output.
22///
23/// # Errors
24///
25/// Return an error if the operation cannot be completed.
26pub fn calc_fq_comp_wo_bed(path: &str, exclude_masked: bool) -> Result<(), std::io::Error> {
27    let fq_iter = FqReader::new(path)?;
28    let mut output = Output::new();
29    if exclude_masked {
30        for record in fq_iter.records() {
31            match record {
32                Ok(read) => {
33                    let result = cal_unmasked_seq(&read);
34                    print(&mut output, read.id(), read.seq().len(), &result)?;
35                }
36                Err(e) => eprintln!("Error read fASTQ: {}", e),
37            }
38        }
39    } else {
40        for record in fq_iter.records() {
41            match record {
42                Ok(read) => {
43                    let result = cal_all_seq(&read);
44                    print(&mut output, read.id(), read.seq().len(), &result)?;
45                }
46                Err(e) => eprintln!("Error read fASTQ: {}", e),
47            }
48        }
49    }
50    Ok(())
51}
52/// Parses FASTA file and compute the statistic w/o masked sequences.
53/// Outputs the results to [`std::io::stdout()`].
54///The output columns are:
55/// - `#seq`: Number of reads
56/// - `#bases`: Number of bases
57/// - `#A`, `#C`, `#G`, `#T`: Number of each nucleotide
58/// - `#2`: Number of `R`, `Y`, `S`, `W`, `K`, `M`
59/// - `#3`: Number of `B`, `D`, `H`, `V`
60/// - `#4`: Number of `N`
61/// - `CG`: Number of `CG` on the template strand
62/// - `GC`: Number of `GC` on the template strand
63///
64/// # Arguments
65///
66/// * `path` - FASTA path
67/// * `exclude_masked` - If true, masked sequences (e.g., lowercases or char other then IUPAC code) will be excluded from the output.
68///
69/// # Errors
70///
71/// Return an error if the operation cannot be completed.
72pub fn calc_fa_comp_wo_bed(path: &str, exclude_masked: bool) -> Result<(), std::io::Error> {
73    let fa_iter = FaReader::new(path)?;
74    let mut output = Output::new();
75    if exclude_masked {
76        for record in fa_iter.records() {
77            match record {
78                Ok(read) => {
79                    let result = cal_unmasked_seq(&read);
80                    print(&mut output, read.id(), read.seq().len(), &result)?;
81                }
82                Err(e) => eprintln!("Error read fASTA: {}", e),
83            }
84        }
85    } else {
86        for record in fa_iter.records() {
87            match record {
88                Ok(read) => {
89                    let result = cal_all_seq(&read);
90                    print(&mut output, read.id(), read.seq().len(), &result)?;
91                }
92                Err(e) => eprintln!("Error read fASTA: {}", e),
93            }
94        }
95    }
96    Ok(())
97}
98/// Parses FASTA file and compute the statistic w/o masked sequences with a BED file.
99/// Outputs the results to [`std::io::stdout()`].
100///The output columns are:
101/// - `#seq`: Number of reads
102/// - `#bases`: Number of bases
103/// - `#A`, `#C`, `#G`, `#T`: Number of each nucleotide
104/// - `#2`: Number of `R`, `Y`, `S`, `W`, `K`, `M`
105/// - `#3`: Number of `B`, `D`, `H`, `V`
106/// - `#4`: Number of `N`
107/// - `CG`: Number of `CG` on the template strand
108/// - `GC`: Number of `GC` on the template strand
109///
110/// # Arguments
111///
112/// * `path` - FASTA path
113/// * `bed` - BED path
114/// * `exclude_masked` - If true, masked sequences (e.g., lowercases or char other then IUPAC code) will be excluded from the output.
115///
116/// # Errors
117///
118/// Return an error if the operation cannot be completed.
119pub fn calc_fq_comp_with_bed(
120    path: &str,
121    bed: &str,
122    exclude_masked: bool,
123) -> Result<(), std::io::Error> {
124    let fq_iter = FqReader::new(path)?;
125    let mut output = Output::new();
126    if exclude_masked {
127        let bedmap = BedMap::from(bed)?;
128        for record in fq_iter.records() {
129            match record {
130                Ok(read) => {
131                    if let Some(result) = cal_unmasked_seq_with_bed(&read, &bedmap) {
132                        print(&mut output, read.id(), read.seq().len(), &result)?;
133                    }
134                }
135                Err(e) => eprintln!("Error read fASTQ: {}", e),
136            }
137        }
138    } else {
139        let bedmap = BedMap::from(bed)?;
140        for record in fq_iter.records() {
141            match record {
142                Ok(read) => {
143                    if let Some(result) = cal_all_seq_with_bed(&read, &bedmap) {
144                        print(&mut output, read.id(), read.seq().len(), &result)?;
145                    }
146                }
147                Err(e) => eprintln!("Error read fASTQ: {}", e),
148            }
149        }
150    }
151    Ok(())
152}
153/// Parses FASTQ file and compute the statistic w/o masked sequences with a BED file.
154/// Outputs the results to [`std::io::stdout()`].
155///The output columns are:
156/// - `#seq`: Number of reads
157/// - `#bases`: Number of bases
158/// - `#A`, `#C`, `#G`, `#T`: Number of each nucleotide
159/// - `#2`: Number of `R`, `Y`, `S`, `W`, `K`, `M`
160/// - `#3`: Number of `B`, `D`, `H`, `V`
161/// - `#4`: Number of `N`
162/// - `CG`: Number of `CG` on the template strand
163/// - `GC`: Number of `GC` on the template strand
164///
165/// # Arguments
166///
167/// * `path` - FASTQ path
168/// * `bed` - BED path
169/// * `exclude_masked` - If true, masked sequences (e.g., lowercases or char other then IUPAC code) will be excluded from the output.
170///
171/// # Errors
172///
173/// Return an error if the operation cannot be completed.
174pub fn calc_fa_comp_with_bed(
175    path: &str,
176    bed: &str,
177    exclude_masked: bool,
178) -> Result<(), std::io::Error> {
179    let fa_iter = FaReader::new(path)?;
180    let mut output = Output::new();
181    if exclude_masked {
182        let bedmap = BedMap::from(bed)?;
183        for record in fa_iter.records() {
184            match record {
185                Ok(read) => {
186                    if let Some(result) = cal_unmasked_seq_with_bed(&read, &bedmap) {
187                        print(&mut output, read.id(), read.seq().len(), &result)?;
188                    }
189                }
190                Err(e) => eprintln!("Error read fASTA: {}", e),
191            }
192        }
193    } else {
194        let bedmap = BedMap::from(bed)?;
195        for record in fa_iter.records() {
196            match record {
197                Ok(read) => {
198                    if let Some(result) = cal_all_seq_with_bed(&read, &bedmap) {
199                        print(&mut output, read.id(), read.seq().len(), &result)?;
200                    }
201                }
202                Err(e) => eprintln!("Error read fASTA: {}", e),
203            }
204        }
205    }
206    Ok(())
207}
208fn cal_all_seq<T: RecordType>(read: &T) -> [usize; 9] {
209    let mut count: [usize; 23] = [0; 23];
210    SeqComp::count_all_nc(&mut count, read.seq(), 0, read.seq().len());
211    SeqComp::get_all_result(&count)
212}
213fn cal_unmasked_seq<T: RecordType>(read: &T) -> [usize; 9] {
214    let mut count: [usize; 23] = [0; 23];
215    SeqComp::count_unmasked_nc(&mut count, read.seq(), 0, read.seq().len());
216    SeqComp::get_unmasked_result(&count)
217}
218fn cal_all_seq_with_bed<T: RecordType>(read: &T, bedmap: &BedMap) -> Option<[usize; 9]> {
219    let mut count: [usize; 23] = [0; 23];
220    if let Some(bedvec) = bedmap.get(read.id()) {
221        bedvec.iter().for_each(|pos| {
222            SeqComp::count_all_nc(&mut count, read.seq(), pos.0, pos.1);
223        });
224        let result = SeqComp::get_all_result(&count);
225        return Some(result);
226    }
227    None
228}
229fn cal_unmasked_seq_with_bed<T: RecordType>(read: &T, bedmap: &BedMap) -> Option<[usize; 9]> {
230    let mut count: [usize; 23] = [0; 23];
231    if let Some(bedvec) = bedmap.get(read.id()) {
232        bedvec.iter().for_each(|pos| {
233            SeqComp::count_unmasked_nc(&mut count, read.seq(), pos.0, pos.1);
234        });
235        let result = SeqComp::get_unmasked_result(&count);
236        return Some(result);
237    }
238    None
239}
240fn print(
241    output: &mut Output,
242    id: &str,
243    size: usize,
244    count: &[usize; 9],
245) -> Result<(), std::io::Error> {
246    output.write(format!(
247        "{}\t{}\t{}\n",
248        id,
249        size,
250        count
251            .iter()
252            .map(|v| v.to_string())
253            .collect::<Vec<_>>()
254            .join("\t")
255    ))?;
256    Ok(())
257}
seqtk_rs/nc_comp.rs

seqtk_rs/
nc_comp.rs
