Commit f4ead3f6 authored by NOEL Philippe's avatar NOEL Philippe

Merge branch 'atoms_weights' into 'master'

Compute center of mass

See merge request pnoel/pdbparser!5
parents 21669cfe d477e7e2
Pipeline #113320 passed with stages
in 10 minutes and 23 seconds
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -3,6 +3,7 @@ use std::io::prelude::*;
use std::io::BufReader;
use std::process;
use crate::structures::atom::AtomType;
use crate::structures::structure::Structure;
/// Parse the string to return a f32. The `trim` is used to remove
......@@ -114,7 +115,8 @@ fn update_from_line(structure: &mut Structure, l: &str) {
// If the "residue" is a amino acid, continue to parse the line and add informations to the protein
// else continue to the next one line
let residue_name = &l[17..20].trim();
let atom_name = &l[12..16].trim();
let atom_name = &l[12..16];
let atom_type = parse_atom(&atom_name);
let chain = l[21..22].chars().next().unwrap();
let atom_number = parse_int(&l[6..11]);
let residue_number = parse_int(&l[22..26]);
......@@ -137,6 +139,7 @@ fn update_from_line(structure: &mut Structure, l: &str) {
residue_number as u64,
atom_name.to_string(),
atom_number as u64,
atom_type,
[x, y, z],
occupancy,
temp_factor,
......@@ -145,3 +148,130 @@ fn update_from_line(structure: &mut Structure, l: &str) {
);
}
}
/// Parse a 4 char string lenght corresponding at the column 12 to 16 in the PDB format
/// These columns contains informations on the Atom type
/// The symbol of the atom is in the 2 first char exepting for Hydrogen atoms
fn parse_atom(txt: &str) -> AtomType {
let symbol = &txt[0..2];
match symbol {
" H" => AtomType::Hydrogen,
"LI" => AtomType::Lithium,
"BE" => AtomType::Beryllium,
" B" => AtomType::Boron,
" C" => AtomType::Carbon,
" N" => AtomType::Nitrogen,
" O" => AtomType::Oxygen,
" F" => AtomType::Fluorine,
"NE" => AtomType::Neon,
"NA" => AtomType::Sodium,
"MG" => AtomType::Magnesium,
"AL" => AtomType::Aluminum,
"SI" => AtomType::Silicon,
" P" => AtomType::Phosphorus,
" S" => AtomType::Sulfur,
"CL" => AtomType::Chlorine,
"AR" => AtomType::Argon,
" K" => AtomType::Potassium,
"CA" => AtomType::Calcium,
"SC" => AtomType::Scandium,
"TI" => AtomType::Titanium,
" V" => AtomType::Vanadium,
"CR" => AtomType::Chromium,
"MN" => AtomType::Manganese,
"FE" => AtomType::Iron,
"CO" => AtomType::Cobalt,
"NI" => AtomType::Nickel,
"CU" => AtomType::Copper,
"ZN" => AtomType::Zinc,
"GA" => AtomType::Gallium,
"GE" => AtomType::Germanium,
"AS" => AtomType::Arsenic,
"SE" => AtomType::Selenium,
"BR" => AtomType::Bromine,
"KR" => AtomType::Krypton,
"RB" => AtomType::Rubidium,
"SR" => AtomType::Strontium,
" Y" => AtomType::Yttrium,
"ZR" => AtomType::Zirconium,
"NB" => AtomType::Niobium,
"MO" => AtomType::Molybdenum,
"TC" => AtomType::Technetium,
"RU" => AtomType::Ruthenium,
"RH" => AtomType::Rhodium,
"PD" => AtomType::Palladium,
"AG" => AtomType::Silver,
"CD" => AtomType::Cadmium,
"IN" => AtomType::Indium,
"SN" => AtomType::Tin,
"SB" => AtomType::Antimony,
"TE" => AtomType::Tellurium,
" I" => AtomType::Iodine,
"XE" => AtomType::Xenon,
"CS" => AtomType::Cesium,
"BA" => AtomType::Barium,
"LA" => AtomType::Lanthanum,
"CE" => AtomType::Cerium,
"PR" => AtomType::Praseodymium,
"ND" => AtomType::Neodymium,
"PM" => AtomType::Promethium,
"SM" => AtomType::Samarium,
"EU" => AtomType::Europium,
"GD" => AtomType::Gadolinium,
"TB" => AtomType::Terbium,
"DY" => AtomType::Dysprosium,
"ER" => AtomType::Erbium,
"TM" => AtomType::Thulium,
"YB" => AtomType::Ytterbium,
"LU" => AtomType::Lutetium,
"TA" => AtomType::Tantalum,
" W" => AtomType::Tungsten,
"RE" => AtomType::Rhenium,
"OS" => AtomType::Osmium,
"IR" => AtomType::Iridium,
"PT" => AtomType::Platinum,
"AU" => AtomType::Gold,
"TL" => AtomType::Thallium,
"PB" => AtomType::Lead,
"BI" => AtomType::Bismuth,
"PO" => AtomType::Polonium,
"AT" => AtomType::Astatine,
"RN" => AtomType::Radon,
"FR" => AtomType::Francium,
"RA" => AtomType::Radium,
"AC" => AtomType::Actinium,
"TH" => AtomType::Thorium,
"PA" => AtomType::Protactinium,
" U" => AtomType::Uranium,
"NP" => AtomType::Neptunium,
"PU" => AtomType::Plutonium,
"AM" => AtomType::Americium,
"CM" => AtomType::Curium,
"BK" => AtomType::Berkelium,
"CF" => AtomType::Californium,
"ES" => AtomType::Einsteinium,
"FM" => AtomType::Fermium,
"MD" => AtomType::Mendelevium,
"NO" => AtomType::Nobelium,
"LR" => AtomType::Lawrencium,
"RF" => AtomType::Rutherfordium,
"DB" => AtomType::Dubnium,
"SG" => AtomType::Seaborgium,
"BH" => AtomType::Bohrium,
"MT" => AtomType::Meitnerium,
_ => {
let next_chars: Vec<char> = txt.chars().collect();
match next_chars[1] {
_ if next_chars[1].is_digit(10) => AtomType::Hydrogen,
_ if next_chars[2].is_digit(10) => AtomType::Hydrogen,
'G' => AtomType::Mercury,
'E' => AtomType::Helium,
'O' => AtomType::Holmium,
'F' => AtomType::Hafnium,
'S' => AtomType::Hassium,
_ => AtomType::Unknown,
}
}
}
}
This diff is collapsed.
......@@ -175,6 +175,9 @@ impl GetAtom for Chain {
}
lst_atom
}
fn compute_weight(&self) -> f32 {
self.get_atom().iter().map(|x| x.get_weight()).sum()
}
}
lazy_static! {
......
......@@ -6,6 +6,7 @@ pub mod structure;
pub trait GetAtom {
fn get_atom(&self) -> Vec<&atom::Atom>;
fn compute_weight(&self) -> f32;
}
pub use atom::Atom;
......
......@@ -68,9 +68,10 @@ impl Residue {
///
/// ````
/// use lib3dmol;
/// use lib3dmol::structures::atom::AtomType;
///
/// let mut lys = lib3dmol::structures::residue::Residue::new(String::from("lysine"), 1);
/// let carbon = lib3dmol::structures::Atom::new(String::from("HT1"), 1, [0.0, 0.0, 0.0]);
/// let carbon = lib3dmol::structures::Atom::new(String::from("C"), AtomType::Carbon, 1, [0.0, 0.0, 0.0]);
///
/// lys.add_atom(carbon);
///
......@@ -89,10 +90,11 @@ impl Residue {
///
/// ```
/// use lib3dmol;
/// use lib3dmol::structures::atom::AtomType;
///
/// let mut lys = lib3dmol::structures::Residue::new(String::from("lysine"), 1);
/// let carbon = lib3dmol::structures::Atom::new(String::from("HT1"), 1, [0.0, 0.0, 0.0]);
/// let hydrogen = lib3dmol::structures::Atom::new(String::from("HT1"), 2, [0.0, 0.0, 0.0]);
/// let carbon = lib3dmol::structures::Atom::new(String::from("CA"), AtomType::Carbon, 1, [0.0, 0.0, 0.0]);
/// let hydrogen = lib3dmol::structures::Atom::new(String::from("HA1"), AtomType::Hydrogen, 2, [0.0, 0.0, 0.0]);
///
/// lys.add_atom(carbon);
/// lys.add_atom(hydrogen);
......@@ -101,7 +103,6 @@ impl Residue {
/// lys.remove_atom(2);
/// assert_eq!(lys.get_number_atom(), 1);
/// ```
pub fn remove_atom(&mut self, n: u64) {
for index in 0..self.lst_atom.len() {
if self.lst_atom[index].number == n {
......@@ -120,6 +121,10 @@ impl GetAtom for Residue {
}
lst_atom
}
fn compute_weight(&self) -> f32 {
self.get_atom().iter().map(|x| x.get_weight()).sum()
}
}
/// Enumerate all Amino acid types
......
use super::atom::Atom;
use super::atom::{Atom, AtomType};
use super::chain::{Chain, ChainTypes};
use super::residue::Residue;
use super::*;
......@@ -221,6 +221,7 @@ impl Structure {
res_number: u64,
atom_name: String,
atom_number: u64,
a_type: AtomType,
coord: [f32; 3],
occupancy: Option<f32>,
tmp_factor: Option<f32>,
......@@ -252,6 +253,7 @@ impl Structure {
atom_name,
atom_number,
coord,
a_type,
occupancy,
tmp_factor,
element,
......@@ -325,6 +327,7 @@ impl Structure {
c_res,
atom.name.clone(),
atom.number,
atom.a_type.clone(),
atom.coord,
atom.occupancy,
atom.temp_factor,
......@@ -385,8 +388,9 @@ impl Structure {
for chain in &mut self.chains {
for residue in &mut chain.lst_res {
for index in (0..residue.lst_atom.len()).rev() {
if residue.lst_atom[index].name().starts_with("H") {
residue.remove_atom(residue.lst_atom[index].number);
match residue.lst_atom[index].a_type {
AtomType::Hydrogen => residue.remove_atom(residue.lst_atom[index].number),
_ => (),
}
}
}
......@@ -407,4 +411,7 @@ impl GetAtom for Structure {
}
lst_atom
}
fn compute_weight(&self) -> f32 {
self.get_atom().iter().map(|x| x.get_weight()).sum()
}
}
use std::f32;
use crate::structures::*;
use crate::tools::tools::atom_mass;
/// Calculate the center of mass for structure implementing the GetAtom trait
///
......@@ -15,24 +14,18 @@ use crate::tools::tools::atom_mass;
///
/// let my_prot = parser::read_pdb("tests/tests_file/5jpq.pdb", "5jqp");
///
/// assert_eq!([237.98595, 241.93814, 231.15921], tools::center_of_mass(&my_prot))
/// assert_eq!([239.22723, 239.0272, 246.19086], tools::center_of_mass(&my_prot))
/// ```
pub fn center_of_mass<T: GetAtom>(atom_list: &T) -> [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.get_atom().iter() {
let mass = match ref_masses.get(&atom.name) {
Some(mass) => *mass,
None => 0.0,
};
let mass = atom.get_weight();
mass_tot += mass;
coord[0] += mass * atom.coord[0];
coord[1] += mass * atom.coord[1];
coord[2] += mass * atom.coord[2];
coord[0] = coord[0] + mass * atom.coord[0];
coord[1] = coord[1] + mass * atom.coord[1];
coord[2] = coord[2] + mass * atom.coord[2];
}
[
......@@ -50,9 +43,10 @@ pub fn center_of_mass<T: GetAtom>(atom_list: &T) -> [f32; 3] {
/// # Example
/// ```
/// use lib3dmol::{structures, tools};
/// use lib3dmol::structures::atom::{Atom, AtomType};
///
/// let atom_1 = structures::atom::Atom::new(String::from("HT1"), 1, [0.0, 0.0, 0.0]);
/// let atom_2 = structures::atom::Atom::new(String::from("HT1"), 1, [0.0, 0.0, 0.0]);
/// let atom_1 = Atom::new(String::from("HT1"), AtomType::Hydrogen, 1, [0.0, 0.0, 0.0]);
/// let atom_2 = Atom::new(String::from("HT1"), AtomType::Hydrogen, 1, [0.0, 0.0, 0.0]);
/// assert_eq!(tools::euclidian_distance(atom_1.coord, atom_2.coord), 0.0);
/// ```
pub fn euclidian_distance(coord1: [f32; 3], coord2: [f32; 3]) -> f32 {
......@@ -100,11 +94,11 @@ mod tests {
#[test]
fn test_select_atoms_around() {
use crate::structures;
use crate::structures::atom::{Atom, AtomType};
let atom_1 = structures::atom::Atom::new(String::from("HT1"), 1, [1.0, 1.0, 1.0]);
let atom_2 = structures::atom::Atom::new(String::from("C"), 2, [10.0, 20.0, 5.0]);
let atom_3 = structures::atom::Atom::new(String::from("N"), 3, [2.0, 3.0, 6.0]); // Should be a distance of 7
let atom_1 = Atom::new(String::from("HT1"), AtomType::Hydrogen, 1, [1.0, 1.0, 1.0]);
let atom_2 = Atom::new(String::from("C"), AtomType::Carbon, 2, [10.0, 20.0, 5.0]);
let atom_3 = Atom::new(String::from("N"), AtomType::Nitrogen, 3, [2.0, 3.0, 6.0]); // Should be a distance of 7
let all_atoms = vec![atom_1, atom_2, atom_3];
......
......@@ -163,4 +163,8 @@ impl GetAtom for std::vec::Vec<Atom> {
}
lst_atom
}
fn compute_weight(&self) -> f32 {
self.get_atom().iter().map(|x| x.get_weight()).sum()
}
}
......@@ -20,29 +20,45 @@ pub fn write_pdb(my_prot: &Structure, file: &str) -> io::Result<()> {
Ok(())
}
fn format_atom(atom: &Atom) -> String {
dbg!(match atom.name.len() {
0 => match atom.get_symbol().len() {
1 => format!(" {} ", atom.get_symbol()),
2 => format!(" {} ", atom.get_symbol()),
_ => format!(" "), // Unreachable because symbol length is 1 or 2
},
1 => format!(" {} ", atom.name),
x if x == 2 && atom.get_symbol().len() == 2 => format!(" {} ", atom.get_symbol()),
x if x == 2 && atom.get_symbol().len() == 1 => format!(" {} ", atom.name),
x if x == 3 && atom.get_symbol().len() == 1 => format!(" {}", atom.name),
4 => atom.name.clone(),
_ => format!(" "), // Unreachable because atom.name cannot be > 4 and < 0
})
}
fn format_pdb(atom: &Atom, res: &Residue, chain_name: &char) -> String {
format!(
"ATOM {:>5} {:<4}{:>3} {}{:>4} {:>8.3}{:>8.3}{:>8.3}{}{} {:>}{}\n",
"ATOM {:>5} {:<4} {:>3} {}{:>4} {:>8.3}{:>8.3}{:>8.3}{}{} {:>}{}\n",
atom.number,
atom.name,
dbg!(format_atom(&atom)),
res.name,
chain_name,
res.res_num,
atom.coord[0],
atom.coord[1],
atom.coord[2],
dbg!(match &atom.occupancy {
match &atom.occupancy {
Some(o) => format!("{:>6.2}", o),
None => String::from(" "),
}),
dbg!(match &atom.temp_factor {
},
match &atom.temp_factor {
Some(t) => format!("{:>6.2}", t),
None => String::from(" "),
}),
dbg!(match &atom.element {
},
match &atom.element {
Some(e) => format!("{:>2}", e),
None => String::from(" "),
}),
},
match &atom.charge {
Some(c) => c,
None => " ",
......@@ -53,11 +69,13 @@ fn format_pdb(atom: &Atom, res: &Residue, chain_name: &char) -> String {
#[test]
fn format_pdb_write() {
use crate::structures;
use crate::structures::atom::{Atom, AtomType};
let my_atom = structures::Atom::new_complete(
let my_atom = Atom::new_complete(
String::from("CA"),
1,
[0.1, 0.2, 0.3],
AtomType::Carbon,
Some(1.01),
Some(31.02),
Some(String::from("C")),
......@@ -68,4 +86,20 @@ fn format_pdb_write() {
let output =
"ATOM 1 CA ALA 1 0.100 0.200 0.300 1.01 31.02 C1+\n";
assert_eq!(output, my_pdb_string);
let my_atom = Atom::new_complete(
String::from("HG21"),
99996,
[0.1, 0.2, 0.3],
AtomType::Carbon,
Some(1.01),
None,
Some(String::from("H")),
None,
);
let my_res = structures::Residue::new(String::from("ALA"), 1);
let my_pdb_string = format_pdb(&my_atom, &my_res, &' ');
let output =
"ATOM 99996 HG21 ALA 1 0.100 0.200 0.300 1.01 H \n";
assert_eq!(output, my_pdb_string);
}
......@@ -14,6 +14,13 @@ mod test_f2 {
assert_eq!(1085, my_struct.get_atom_number());
}
#[test]
fn molecular_weight() {
use lib3dmol::structures::GetAtom;
let my_struct = parser::read_pdb("tests/tests_file/f2.pdb", "f2");
assert_eq!(7712, my_struct.compute_weight() as u64);
}
#[test]
fn get_residue_ref() {
let mut my_struct = parser::read_pdb("tests/tests_file/f2.pdb", "f2");
......@@ -36,7 +43,7 @@ mod test_f2 {
let my_struct = parser::read_pdb("tests/tests_file/f2.pdb", "f2");
assert_eq!(
[0.011766517, 20.1037, -0.007502083],
[0.00042509945, 19.999798, -0.00044671507],
tools::center_of_mass(&my_struct)
)
}
......@@ -77,20 +84,45 @@ mod test_tools {
assert_eq!(133_831, my_struct.get_atom_number());
}
#[test]
fn remove_hydrogens_f2() {
let mut my_struct = parser::read_pdb("tests/tests_file/f2.pdb", "f2");
assert_eq!(1085, my_struct.get_atom_number());
my_struct.remove_h();
assert_eq!(541, my_struct.get_atom_number());
}
#[test]
fn test_rmsd() {
use lib3dmol::structures::atom::Atom;
let lst_atom_1 = vec![Atom::new(String::from("H"), 1, [0.0, 0.0, 0.0])];
let lst_atom_2 = vec![Atom::new(String::from("H"), 1, [0.6, 0.3, 0.2])];
use lib3dmol::structures::atom::{Atom, AtomType};
let lst_atom_1 = vec![Atom::new(
String::from("H"),
AtomType::Hydrogen,
1,
[0.0, 0.0, 0.0],
)];
let lst_atom_2 = vec![Atom::new(
String::from("H"),
AtomType::Hydrogen,
1,
[0.6, 0.3, 0.2],
)];
assert_eq!(0.7, tools::rmsd(&lst_atom_1, &lst_atom_2).unwrap());
}
#[test]
fn compute_weight_5jpq() {
use lib3dmol::structures::GetAtom;
let my_struct = parser::read_pdb("tests/tests_file/5jpq.pdb", "5jpq");
assert_eq!(1294368.6, my_struct.compute_weight());
}
#[test]
fn mass_center_5jpq() {
let my_struct = parser::read_pdb("tests/tests_file/5jpq.pdb", "5jpq");
assert_eq!(
[237.98595, 241.93814, 231.15921],
[239.22723, 239.0272, 246.19086],
tools::center_of_mass(&my_struct)
)
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment