use std::{collections::HashMap, io::BufRead};
use itertools::Itertools;
use mzcore::sequence::AminoAcid;
use crate::{
NotASpecies, Species,
imgt_gene::IMGTGene,
structs::{AASequence, DataItem, Location, Region},
};
pub(crate) fn parse_dat<T: BufRead>(
reader: T,
) -> impl Iterator<Item = Result<(DataItem, Vec<String>), String>> {
reader
.lines()
.batching(|f| {
let mut data = PreDataItem::default();
let mut errors = Vec::new();
for line in f.filter_map(Result::ok) {
let (end, error) = parse_dat_line(&mut data, &line);
if !error.is_empty() {
errors.push(error);
}
if end {
return Some((data, errors));
}
}
None
})
.filter(|(pre, _)| {
pre.kw.contains(&"immunoglobulin (IG)".to_string())
&& (pre.kw.contains(&"functional".to_string())
|| pre.kw.contains(&"germline".to_string())
|| pre.kw.contains(&"productive".to_string()))
&& pre.os.is_some()
})
.map(DataItem::new)
}
fn parse_dat_line(data: &mut PreDataItem, line: &str) -> (bool, String) {
let mut error = String::new();
if line.len() < 2 {
return (false, error);
}
match &line[..2] {
"//" => return (true, error),
"ID" => data.id = line.to_string(),
"KW" => data.kw.extend(
line[5..]
.split(';')
.map(|s| s.trim().to_string())
.filter(|s| !s.is_empty()),
),
"FH" if line.starts_with("FH Key") => {
if let Some(pos) = line.find("Location") {
data.ft_key_width = pos - 5;
} else {
error = "Incorrect FH line".to_string();
}
}
"FT" => data.ft.push(line[5..].to_string()),
"OS" if data.os.is_none() => {
data.os = Species::from_imgt(line[5..].trim()).unwrap_or_else(|NotASpecies| {
error = format!("Not a species name: `{line}`");
None
});
}
" " => data.sq.extend(
line.chars()
.filter(|c| *c == 'c' || *c == 'a' || *c == 't' || *c == 'g' || *c == 'n'),
),
"DT" if line.contains("Last updated") => {
if let Ok(year) = line[12..=15].parse()
&& let Ok(day) = line[5..=6].parse()
&& let Some(month) = [
"JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV",
"DEC",
]
.iter()
.position(|i| **i == line[8..=10])
{
data.dt = Some((year, month as u8 + 1, day, line[23..=30].to_string()));
} else {
error = "Invalid DT line".to_string();
}
}
_ => (),
}
(false, error)
}
impl DataItem {
fn new(data: (PreDataItem, Vec<String>)) -> Result<(Self, Vec<String>), String> {
let (data, errors) = data;
let mut result = Self {
id: data.id[5..].split(';').next().unwrap().to_string(),
species: data.os.ok_or("No species found")?,
sequence: data.sq,
genes: Vec::new(),
regions: Vec::new(),
release: data.dt,
};
let mut current: Option<Region> = None;
let mut is_sequence = false;
for line in data.ft {
if !line.starts_with(' ') || current.is_none() {
if let Some(region) = current.take() {
result.add_region(region);
}
if let Ok(location) = line[data.ft_key_width..].parse() {
let (key, location) = (line[..data.ft_key_width].trim().to_string(), location);
current = Some(Region {
acc: result.id.clone(),
key,
location,
reported_seq: String::new(),
found_seq: Err("Not loaded".to_string()),
allele: String::new(),
functional: true,
partial: false,
shift: 0,
splice_aa: None,
});
}
continue;
}
if let Some(current) = &mut current {
Self::parse_ft_line(&line, current, &mut is_sequence)?;
}
}
if let Some(region) = current.take() {
result.add_region(region);
}
Ok((result, errors))
}
fn parse_ft_line(
line: &str,
current: &mut Region,
is_sequence: &mut bool,
) -> Result<(), String> {
let trimmed = line.trim();
let split = trimmed.split_once('=').map(|(key, tail)| (key.to_ascii_lowercase(), tail));
match split
.as_ref()
.map_or((trimmed, None), |(key, tail)| (key.trim(), Some(tail)))
{
("/translation", Some(tail)) => {
current.reported_seq += tail.trim_matches('\"');
*is_sequence = !tail.ends_with('\"');
}
("/imgt_gene" | "/allele", Some(tail)) => {
let here = tail.trim_matches('\"').to_string();
if current.allele.is_empty() || !current.allele.contains('*') && here.contains('*')
{
current.allele = here;
}
}
("/imgt_allele", Some(tail)) => {
current.allele = tail.trim_matches('\"').to_string();
}
("/codon_start", Some(tail)) => {
current.shift = tail
.trim_matches('\"')
.parse::<usize>()
.map_err(|_| format!("Not a valid codon_start: '{tail}'"))?
- 1;
}
("/splice-expectedcodon", Some(tail)) => {
if let Some(i) = tail.find(']') {
current.splice_aa = AminoAcid::try_from(tail.as_bytes()[i - 1]).ok();
}
}
("/functional", _) => {
current.functional = true;
}
("/note" | "/imgt_note", Some(s)) => {
current.functional = s.to_ascii_lowercase().contains("functional");
}
("/pseudo" | "/orf", _) => {
current.functional = false;
}
("/partial", _) => current.partial = true,
(_, None) if *is_sequence => {
current.reported_seq += trimmed.trim_end_matches('\"');
*is_sequence = !trimmed.ends_with('\"');
}
_ => (),
}
Ok(())
}
fn add_region(&mut self, mut region: Region) {
region.found_seq = self.get_sequence(®ion.location, region.shift);
if ["V-GENE", "C-GENE", "J-GENE"].contains(®ion.key.as_str()) && region.functional
&& !region.partial
&& region.allele.starts_with("IG")
{
self.genes.push(IMGTGene {
acc: region.acc,
key: region.key,
location: region.location,
allele: region.allele,
regions: HashMap::new(),
});
} else if ["V-REGION", "C-REGION", "J-REGION"].contains(®ion.key.as_str()) && region.functional
&& !region.partial
&& region.allele.starts_with("IG")
{
if let Some(existing) = self.genes.iter_mut().find(|g| g.allele == region.allele) {
existing.regions.insert(region.key.clone(), region);
} else {
self.genes.push(IMGTGene {
acc: region.acc,
key: region.key,
location: region.location,
allele: region.allele,
regions: HashMap::new(),
});
}
} else if [
"FR1-IMGT",
"FR2-IMGT",
"FR3-IMGT",
"FR4-IMGT",
"CDR1-IMGT",
"CDR2-IMGT",
"CDR3-IMGT",
"1st-CYS",
"2nd-CYS",
"CONSERVED-TRP",
"J-REGION",
"CH1",
"CH2",
"H-CH2",
"CH3",
"CH3-CHS",
"CH4",
"CH4-CHS",
"CH5",
"CH5-CHS",
"CH6",
"CH6-CHS",
"CH7",
"CH7-CHS",
"CH8",
"CH8-CHS",
"CH9",
"CH9-CHS",
"CHS",
"CL",
"C-REGION",
"H", "H1",
"H2",
"H3",
"H4",
"M",
"M1",
"M2",
]
.contains(®ion.key.as_str())
{
if region.key == "CDR3-IMGT" {
if let Some(gene) =
self.genes.iter_mut().find(|g| g.location.overlaps(®ion.location))
{
gene.regions.insert(region.key.clone(), region);
} else {
self.regions.push(region);
}
} else if let Some(gene) =
self.genes.iter_mut().find(|g| g.location.contains(®ion.location))
{
gene.regions.insert(region.key.clone(), region);
} else {
self.regions.push(region);
}
}
}
fn get_sequence(&self, slice: &Location, shift: usize) -> Result<(String, AASequence), String> {
let (inner_shift, shift) = if shift == 2 { (1, 0) } else { (0, shift) };
translate(
&match slice {
Location::Normal(range) => {
if *range.start() < inner_shift {
return Err("Shift outside of range".to_string());
}
self.sequence
.get(range.start() - inner_shift..=*range.end())
.ok_or("Normal outside of range")?
.to_string()
}
Location::SingleNormal(index) => char::from(
*self
.sequence
.as_bytes()
.get(*index)
.ok_or("Single normal outside of range")?,
)
.to_string(),
Location::Complement(range) => complement(
self.sequence
.get(*range.start()..=*range.end() + inner_shift)
.ok_or("Complement outside of range")?
.to_string(),
),
Location::SingleComplement(index) => complement(
char::from(
*self
.sequence
.as_bytes()
.get(*index)
.ok_or("Single complement outside of range")?,
)
.to_string(),
),
}[shift..],
)
.map(|(s, v)| (s.to_owned(), AASequence(v)))
}
}
#[derive(Debug, Default)]
struct PreDataItem {
id: String,
kw: Vec<String>,
ft_key_width: usize,
ft: Vec<String>,
os: Option<Species>,
sq: String,
dt: Option<(u16, u8, u8, String)>,
}
pub(crate) fn complement(s: String) -> String {
let map = HashMap::from([
(b'a', b't'),
(b't', b'a'),
(b'c', b'g'),
(b'g', b'c'),
(b'n', b'n'),
]);
String::from_utf8(
s.as_bytes()
.iter()
.map(|c| {
*map.get(c)
.unwrap_or_else(|| panic!("Invalid sequence: {} in `{s}`", char::from(*c)))
})
.rev()
.collect(),
)
.unwrap()
}
pub(crate) fn translate(s: &str) -> Result<(&str, Vec<AminoAcid>), String> {
if s.len() < 3 {
Ok((s, Vec::new()))
} else {
Ok((
s,
(0..=s.len() - 3)
.step_by(3)
.filter_map(|chunk| {
AminoAcid::from_dna(&s[chunk..chunk + 3])
.map_err(|_| format!("Not a codon {}", &s[chunk..chunk + 3]))
.transpose()
})
.collect::<Result<Vec<AminoAcid>, String>>()?,
))
}
}