tools.rs 4.78 KB
Newer Older
1
use std::collections::HashMap;
Emmanuel Bresso's avatar
Emmanuel Bresso committed
2 3 4 5 6 7
use std::env;
use std::fs;
use std::fs::File;
use std::io::BufReader;
use std::io::prelude::BufRead;
use std::process;
8

9
use super::structure::Structure;
Emmanuel Bresso's avatar
Emmanuel Bresso committed
10
use super::atom::Atom;
11

12
/// Convert the all amino acid to a FASTA sequence (1 residue as 1 char)
13 14 15 16 17 18 19
/// Consult the corresponding table to have the code 1 letter <-> 3 letters
/// [Wikipedia amino acid](https://en.wikipedia.org/wiki/Amino_acid)
///
/// # Examples
/// ```
/// use pdbparser;
///
20 21
/// let my_struct = pdbparser::read_pdb("tests/tests_file/f2.pdb", "f2");
/// assert_eq!("TSPQPYSIERTIRWLTYQVANSLALVSEADKIMQTEYMKMIQNSGEITDRGEAILRLLKTNKHYEH", pdbparser::fasta_seq(&my_struct));
22
/// ```
23
pub fn fasta_seq(my_struct: &Structure) -> String {
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
    let res: HashMap<&str, char> = [
        ("ARG", 'R'),
        ("LYS", 'K'),
        ("ASN", 'N'),
        ("ASP", 'D'),
        ("GLU", 'E'),
        ("SER", 'S'),
        ("THR", 'T'),
        ("GLN", 'Q'),
        ("CYS", 'C'),
        ("HIS", 'H'),
        ("HSD", 'H'),
        ("HSP", 'H'),
        ("HSD", 'H'),
        ("SEC", 'U'),
        ("GLY", 'G'),
        ("PRO", 'P'),
        ("ALA", 'A'),
        ("VAL", 'V'),
        ("ILE", 'I'),
        ("LEU", 'L'),
        ("MET", 'M'),
        ("PHE", 'P'),
        ("TYR", 'Y'),
        ("TRP", 'W'),
    ]
    .iter()
    .cloned()
    .collect();

54 55
    //TODO: Change the with_capacity(my_struct.get_residue_number()) because get all residue (dna, lipid, etc..)
    let mut fasta = String::with_capacity(my_struct.get_residue_number() as usize);
56

57
    for chain in &my_struct.chains {
58
        for residue in &chain.lst_res {
59 60 61
            if let Some(r) = res.get(&residue.name()[..]) {
                fasta.push(*r);
            }
62 63 64
        }
    }
    fasta
65
}
Emmanuel Bresso's avatar
Emmanuel Bresso committed
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

// extract atom mass from Charmm parameter files.
fn atom_mass() -> HashMap<String, f32> {
    let mut ref_masses = HashMap::new();

    // If CHARMM_FOLDER enviroment variable is defined, then look for .parm files in this folder. Else use the ./datasets folder*.
    let key = "CHARMM_FOLDER";
    let charmm_folder = match env::var(key) {
        Ok(charmm_folder) => charmm_folder,
        _ => "./datasets".to_string(),
    };
    let paths = match fs::read_dir(charmm_folder) {
        Ok(paths) => paths,
        Err(_e) => {
            eprintln!("cannot find your charmm folder");
            process::exit(1);
        }
    };

    for path in paths {
        let test = path.unwrap().path();
        if test.is_file() && test.extension().unwrap().to_str() == Some("prm") {
            let charmmf = match File::open(test) {
                Ok(charmmf) => charmmf,
                Err(e) => {
Emmanuel Bresso's avatar
Emmanuel Bresso committed
91
                    eprintln!("Error while reading {}", e);
Emmanuel Bresso's avatar
Emmanuel Bresso committed
92 93 94 95 96 97 98 99 100 101 102 103 104
                    process::exit(1);
                }
            };

            let reader = BufReader::new(charmmf);
            for line in reader.lines() {
                // MASS  -1  H          1.00800 ! polar H
                let l = line.unwrap();
                if l.starts_with("MASS") {
                    let split: Vec<&str> = l.split_ascii_whitespace().collect();
                    let val = match split[3].parse::<f32>() {
                        Ok(val) => val,
                        Err(e) => {
Emmanuel Bresso's avatar
Emmanuel Bresso committed
105
                            eprintln!("Cannot cast {} into f32", e);
Emmanuel Bresso's avatar
Emmanuel Bresso committed
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
                            process::exit(1);
                        }
                    };
                    ref_masses.insert(split[2].to_string(), val);
                }
            }
        }
    }
    return ref_masses;
}

/// Calculate the center of mass of a selection of atoms
/// 
/// 
/// Atom masses are extracted from CHARMM Force Field .prm files.
/// If CHARMM_FOLDER environement variable is set, .prm files in this folder are used. Else it will be files in datasets folder.
/// 
/// # Example
/// ```
125 126
/// use pdbparser;
///
Emmanuel Bresso's avatar
Emmanuel Bresso committed
127
/// let my_prot = pdbparser::read_pdb("tests/tests_file/5jpq.pdb", "5jqp");
128 129
///
/// let atom_number = my_prot.get_atom_number() as usize;
Emmanuel Bresso's avatar
Emmanuel Bresso committed
130
/// let mut atoms: Vec<pdbparser::Atom> = Vec::with_capacity(atom_number);
Emmanuel Bresso's avatar
Emmanuel Bresso committed
131
/// for chain in my_prot.chains {
Emmanuel Bresso's avatar
Emmanuel Bresso committed
132 133 134 135 136 137
///     for residue in chain.lst_res {
///         for atom in residue.lst_atom {
///             atoms.push(atom);
///         }
///     }
/// }
Emmanuel Bresso's avatar
Emmanuel Bresso committed
138
/// assert_eq!([237.98595, 241.93814, 231.15921], pdbparser::center_of_mass(atoms))
Emmanuel Bresso's avatar
Emmanuel Bresso committed
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
/// ```
pub fn center_of_mass(atom_list: Vec<Atom>) -> [f32; 3] {

    let mut coord: [f32;3] = [0.0, 0.0, 0.0];
    let mut mass_tot: f32 = 0.0;

    let ref_masses = atom_mass();

    for atom in atom_list {
        let mass = match ref_masses.get(&atom.name) {
            Some(mass) => *mass,
            None => 0.0
        };

        mass_tot += mass;
        coord[0] += mass * atom.coord[0];
        coord[1] += mass * atom.coord[1];
        coord[2] += mass * atom.coord[2];
    }

    [coord[0] / mass_tot, coord[1] / mass_tot, coord[2] / mass_tot]
}