use crate::*;
use anyhow::anyhow;
use bio::io::fasta;
use itertools::Itertools;
use std::collections::HashSet;
use std::io::Write;
use std::process::Command;
use std::str;
lazy_static! {
static ref BASES: HashSet<u8> = vec![b'a', b'g', b'c', b't', b'A', b'G', b'C', b'T',]
.into_iter()
.collect();
}
pub fn pair_d(seq1: &[u8], seq2: &[u8]) -> f32 {
assert_eq!(
seq1.len(),
seq2.len(),
"Two sequences of different length ({}!={})",
seq1.len(),
seq2.len()
);
let mut comparable = 0;
let mut difference = 0;
for (base1, base2) in seq1.iter().zip(seq2) {
if BASES.contains(base1) && BASES.contains(base2) {
comparable += 1;
if base1.to_ascii_uppercase() != base2.to_ascii_uppercase() {
difference += 1;
}
}
}
assert_ne!(comparable, 0, "Comparable bases shouldn't be zero");
difference as f32 / comparable as f32
}
pub fn alignment_stat(seqs: &[&[u8]]) -> (i32, i32, i32, i32, i32, f32) {
let seq_count = seqs.len();
assert_ne!(seq_count, 0, "Need sequences");
let length = seqs[0].len();
let mut comparable = 0;
let mut difference = 0;
let mut gap = 0;
let mut ambiguous = 0;
for pos in 0..length {
let mut column = vec![];
for i in 0..seq_count {
column.push(seqs[i][pos].to_ascii_uppercase());
}
column = column.into_iter().unique().collect();
if column.clone().into_iter().all(|e| BASES.contains(&e)) {
comparable += 1;
if column.clone().into_iter().any(|e| e != column[0]) {
difference += 1;
}
} else if column.clone().into_iter().any(|e| e == b'-') {
gap += 1;
} else {
ambiguous += 1;
}
}
assert_ne!(comparable, 0, "Comparable bases shouldn't be zero");
let mut dists = vec![];
for i in 0..seq_count {
for j in i + 1..seq_count {
let dist = pair_d(seqs[i], seqs[j]);
dists.push(dist);
}
}
let mean_d = f32::trunc(dists.iter().sum::<f32>() / dists.len() as f32 * 1000.0) / 1000.0;
(
length as i32,
comparable,
difference,
gap,
ambiguous,
mean_d,
)
}
pub fn indel_intspan(seq: &[u8]) -> IntSpan {
let mut positions = vec![];
for (i, base) in seq.iter().enumerate() {
if *base == b'-' {
positions.push(i as i32 + 1);
}
}
let mut ints = IntSpan::new();
ints.add_vec(&positions);
ints
}
pub fn seq_intspan(seq: &[u8]) -> IntSpan {
IntSpan::from_pair(1, seq.len() as i32).diff(&indel_intspan(seq))
}
pub fn align_seqs(seqs: &[&[u8]], aligner: &str) -> anyhow::Result<Vec<String>> {
let mut bin = String::new();
match aligner {
"clustalw" => {
for e in &["clustalw", "clustal-w", "clustalw2"] {
if let Ok(pth) = which::which(e) {
bin = pth.to_string_lossy().to_string();
break;
}
}
}
"muscle" => {
for e in &["muscle"] {
if let Ok(pth) = which::which(e) {
bin = pth.to_string_lossy().to_string();
break;
}
}
}
"mafft" => {
for e in &["mafft"] {
if let Ok(pth) = which::which(e) {
bin = pth.to_string_lossy().to_string();
break;
}
}
}
_ => {
return Err(anyhow!("Unrecognized aligner: {}", aligner));
}
};
if bin.is_empty() {
return Err(anyhow!("Can't find the external command: {}", aligner));
}
let mut seq_in = tempfile::Builder::new()
.prefix("seq-in-")
.suffix(".fasta")
.rand_bytes(8)
.tempfile()?;
for (i, seq) in seqs.iter().enumerate() {
write!(seq_in, ">seq-{}\n{:?}\n", i, str::from_utf8(seq).unwrap())?;
}
let seq_in_path = seq_in.into_temp_path();
let seq_out = tempfile::Builder::new()
.prefix("seq-out-")
.suffix(".fasta")
.rand_bytes(8)
.tempfile()?;
let seq_out_path = seq_out.into_temp_path();
let output = match aligner {
"clustalw" => Command::new(bin)
.arg("-align")
.arg("-type=dna")
.arg("-output=fasta")
.arg("-outorder=input")
.arg("-quiet")
.arg(format!("-infile={}", seq_in_path.to_string_lossy()))
.arg(format!("-outfile={}", seq_out_path.to_string_lossy()))
.output()?,
"muscle" => Command::new(bin)
.arg("-quiet")
.arg("-in")
.arg(seq_in_path.to_string_lossy().to_string())
.arg("-out")
.arg(seq_out_path.to_string_lossy().to_string())
.output()?,
"mafft" => Command::new(bin)
.arg("-quiet")
.arg("-auto")
.arg(seq_in_path.to_string_lossy().to_string())
.arg(">")
.arg(seq_out_path.to_string_lossy().to_string())
.output()?,
_ => unreachable!(),
};
if !output.status.success() {
return Err(anyhow!("Command executed with failing error code"));
}
let mut out_seq = vec![];
let reader = reader(seq_out_path.to_string_lossy().as_ref());
let fa_in = fasta::Reader::new(reader);
for result in fa_in.records() {
let record = result.unwrap();
out_seq.push(String::from_utf8(record.seq().to_vec()).unwrap());
}
seq_in_path.close()?;
seq_out_path.close()?;
Ok(out_seq)
}
pub fn chr_to_align(ints: &IntSpan, pos: i32, chr_start: i32, strand: &str) -> anyhow::Result<i32> {
let chr_end = chr_start + ints.size() - 1;
if pos < chr_start || pos > chr_end {
return Err(anyhow!(
"[{}] out of ranges [{}, {}]",
pos,
chr_start,
chr_end
));
}
let aln_pos = match strand {
"+" => ints.at(pos - chr_start + 1),
"-" => ints.at(-(pos - chr_start + 1)),
_ => {
return Err(anyhow!("Unrecognized strand: {}", strand));
}
};
Ok(aln_pos)
}
pub fn align_to_chr(ints: &IntSpan, pos: i32, chr_start: i32, strand: &str) -> anyhow::Result<i32> {
let chr_end = chr_start + ints.size() - 1;
if pos < 1 {
return Err(anyhow!("align pos [{}] out of ranges", pos,));
}
let mut chr_pos = if ints.contains(pos) {
ints.index(pos)
} else if pos < ints.min() {
1
} else if pos > ints.max() {
ints.size()
} else {
let spans = ints.spans();
let mut cursor = pos;
for i in 0..spans.len() {
if spans[i].1 < cursor {
continue;
} else {
cursor = spans[i - 1].1;
break;
}
}
ints.index(cursor)
};
chr_pos = match strand {
"+" => chr_pos + chr_start - 1,
"-" => chr_end - chr_pos + 1,
_ => {
return Err(anyhow!("Unrecognized strand: {}", strand));
}
};
Ok(chr_pos)
}