lib3dmol/tools/
tools.rs

1use std::collections::HashMap;
2use std::env;
3use std::fs;
4use std::fs::File;
5use std::io::prelude::BufRead;
6use std::io::BufReader;
7use std::process;
8
9use crate::structures::atom::Atom;
10use crate::structures::structure::Structure;
11use crate::structures::*;
12
13/// Convert the all amino acid to a FASTA sequence (1 residue as 1 char)
14/// Consult the corresponding table to have the code 1 letter <-> 3 letters
15/// [Wikipedia amino acid](https://en.wikipedia.org/wiki/Amino_acid)
16///
17/// # Examples
18/// ```
19/// use lib3dmol::{parser, tools};
20///
21/// let my_struct = parser::read_pdb("tests/tests_file/f2.pdb", "f2");
22/// assert_eq!("TSPQPYSIERTIRWLTYQVANSLALVSEADKIMQTEYMKMIQNSGEITDRGEAILRLLKTNKHYEH", tools::fasta_seq(&my_struct));
23/// ```
24pub fn fasta_seq(my_struct: &Structure) -> String {
25    let res: HashMap<&str, char> = [
26        ("ARG", 'R'),
27        ("LYS", 'K'),
28        ("ASN", 'N'),
29        ("ASP", 'D'),
30        ("GLU", 'E'),
31        ("SER", 'S'),
32        ("THR", 'T'),
33        ("GLN", 'Q'),
34        ("CYS", 'C'),
35        ("HIS", 'H'),
36        ("HSD", 'H'),
37        ("HSP", 'H'),
38        ("HSD", 'H'),
39        ("SEC", 'U'),
40        ("GLY", 'G'),
41        ("PRO", 'P'),
42        ("ALA", 'A'),
43        ("VAL", 'V'),
44        ("ILE", 'I'),
45        ("LEU", 'L'),
46        ("MET", 'M'),
47        ("PHE", 'P'),
48        ("TYR", 'Y'),
49        ("TRP", 'W'),
50    ]
51    .iter()
52    .cloned()
53    .collect();
54
55    //TODO: Change the with_capacity(my_struct.get_residue_number()) because get all residue (dna, lipid, etc..)
56    let mut fasta = String::with_capacity(my_struct.get_residue_number() as usize);
57
58    for chain in &my_struct.chains {
59        for residue in &chain.lst_res {
60            if let Some(r) = res.get(&residue.name()[..]) {
61                fasta.push(*r);
62            }
63        }
64    }
65    fasta
66}
67
68/// Compute RMSD distance between 2 structures wich have the GetAtom trait
69///
70/// Be careful, you should check if atoms are the same
71///
72/// # Example:
73/// ```
74/// use lib3dmol::{parser, tools};
75///
76/// // f2 and f2_adn are the same protein in files
77/// let my_first_prot = parser::read_pdb("tests/tests_file/f2.pdb", "f2");
78/// let my_second_prot = parser::read_pdb("tests/tests_file/f2_adn.pdb", "f2_adn");
79///
80/// let my_second_prot = my_second_prot.select_atoms("chain A").unwrap();
81///
82/// let rmsd = tools::rmsd(&my_first_prot, &my_second_prot);
83///
84/// assert_eq!(Ok(0.), rmsd)
85/// ```
86pub fn rmsd<T: GetAtom>(first: &T, second: &T) -> Result<f32, &'static str> {
87    let atoms_first = first.get_atom();
88    let atoms_second = second.get_atom();
89
90    // RMSD can be compute only if structures have same number of atoms
91    if atoms_first.len() != atoms_second.len() {
92        return Err("Different number of atoms");
93    }
94    let mut total = 0.0;
95    for index in 0..atoms_first.len() {
96        let coord_first = &atoms_first[index].coord;
97        let coord_second = &atoms_second[index].coord;
98        total += (coord_first[0] - coord_second[0]).powi(2)
99            + (coord_first[1] - coord_second[1]).powi(2)
100            + (coord_first[2] - coord_second[2]).powi(2)
101    }
102    Ok((total / atoms_first.len() as f32).sqrt())
103}
104
105/// Extract atom mass from Charmm parameter files.
106/// Charmm files must have the "prm" extension
107///
108pub fn atom_mass() -> HashMap<String, f32> {
109    let mut ref_masses: HashMap<String, f32> = HashMap::new();
110
111    // If CHARMM_FOLDER enviroment variable is defined, then look for .parm files in this folder. Else use the ./datasets folder*.
112    let key = "CHARMM_FOLDER";
113    let charmm_folder = match env::var(key) {
114        Ok(charmm_folder) => charmm_folder,
115        _ => "./datasets".to_string(),
116    };
117    let paths = match fs::read_dir(charmm_folder) {
118        Ok(paths) => paths,
119        Err(_e) => {
120            eprintln!("cannot find your charmm folder");
121            process::exit(1);
122        }
123    };
124
125    for path in paths {
126        let test = path.unwrap().path();
127        if test.is_file() && test.extension().unwrap().to_str() == Some("prm") {
128            let charmmf = match File::open(test) {
129                Ok(charmmf) => charmmf,
130                Err(e) => {
131                    eprintln!("Error while reading {}", e);
132                    process::exit(1);
133                }
134            };
135
136            let reader = BufReader::new(charmmf);
137            for line in reader.lines() {
138                // MASS  -1  H          1.00800 ! polar H
139                let l = line.unwrap();
140                if l.starts_with("MASS") {
141                    let split: Vec<&str> = l.split_ascii_whitespace().collect();
142                    let val = match split[3].parse::<f32>() {
143                        Ok(val) => val,
144                        Err(e) => {
145                            eprintln!("Cannot cast {} into f32", e);
146                            process::exit(1);
147                        }
148                    };
149                    ref_masses.insert(split[2].to_string(), val);
150                }
151            }
152        }
153    }
154    return ref_masses;
155}
156
157/// Implement the GetAtom trait on vector of Atoms
158impl GetAtom for std::vec::Vec<Atom> {
159    fn get_atom(&self) -> Vec<&Atom> {
160        let mut lst_atom: Vec<&Atom> = Vec::with_capacity(self.len());
161        for atom in self.iter() {
162            lst_atom.push(&atom)
163        }
164        lst_atom
165    }
166
167    fn compute_weight(&self) -> f32 {
168        self.get_atom().iter().map(|x| x.get_weight()).sum()
169    }
170}
lib3dmol/tools/tools.rs

lib3dmol/tools/
tools.rs
