haddock_restraints/core/commands/
ti.rs

1use crate::*;
2use std::error::Error;
3
4/// Analyzes the true interface of a protein structure and generates Ambiguous Interaction Restraints (AIRs).
5///
6/// This function reads a PDB file, identifies the true interface between chains based on a distance cutoff,
7/// creates interactors for each chain involved in the interface, and generates AIRs.
8///
9/// # Arguments
10///
11/// * `input_file` - A string slice that holds the path to the input PDB file.
12/// * `cutoff` - A reference to a f64 value specifying the distance cutoff (in Angstroms) for determining interfaces.
13///
14/// # Returns
15///
16/// A `Result<String, Box<dyn Error>>` which is Ok(String) if the function completes successfully, or an Error if something goes wrong.
17///
18pub fn true_interface(
19    pdb: pdbtbx::PDB,
20    cutoff: &f64,
21    pml: &Option<String>,
22    pdb_path: &str,
23) -> Result<String, Box<dyn Error>> {
24    let true_interface = get_true_interface(&pdb, *cutoff);
25    let chains_in_contact = get_chains_in_contact(&pdb, *cutoff);
26
27    // Sort the true_interface by chain id
28    let mut true_interface: Vec<_> = true_interface.iter().collect();
29    true_interface.sort_by(|a, b| a.0.cmp(b.0));
30
31    // NOTE: Here the IDs of the interactors are their position in the PDB file; this is so that
32    // we can handle any order of chains.
33    let mut interactors: Vec<Interactor> = Vec::new();
34    for (chain_id, residues) in true_interface.iter() {
35        // Get what is the position of this chain in the PDB, this will be its ID
36        let target_id = pdb
37            .chains()
38            .position(|chain| chain.id() == *chain_id)
39            .unwrap();
40
41        let mut interactor = Interactor::new(target_id as u16);
42        interactor.set_chain(chain_id);
43        interactor.set_active(residues.iter().map(|&residue| residue as i16).collect());
44        interactor.set_pdb(pdb.clone());
45        interactor.remove_buried_residues();
46
47        // Assign the targets
48        for (chain_i, chain_j) in chains_in_contact.iter() {
49            let target_chain = if chain_i == *chain_id {
50                chain_j
51            } else {
52                chain_i
53            };
54            if let Some(target_index) = pdb.chains().position(|chain| chain.id() == *target_chain) {
55                interactor.add_target(target_index as u16);
56            }
57        }
58
59        interactors.push(interactor);
60    }
61
62    // Make the restraints
63    let air = Air::new(interactors);
64    let tbl = air.gen_tbl().unwrap();
65
66    println!("{}", tbl);
67
68    if let Some(output_f) = pml {
69        air.gen_pml(output_f, Some(pdb_path))
70    };
71
72    Ok(tbl)
73}
74
75/// Generates Unambiguous Topological Interactions (TIs) from a protein structure.
76///
77/// This function reads a PDB file, identifies the closest residue pairs based on a specified distance cutoff,
78/// and creates unambiguous interactors for each residue pair.
79///
80/// # Arguments
81///
82/// * input_file - A string slice that holds the path to the input PDB file.
83/// * cutoff - A reference to a f64 value specifying the distance cutoff (in Angstroms) for determining interactions.
84///
85/// # Returns
86///
87/// A Result<String, Box<dyn Error>> containing the generated TBL (Topological Restraints List) if successful.
88///
89pub fn unambig_ti(
90    pdb: pdbtbx::PDB,
91    cutoff: &f64,
92    pml: &Option<String>,
93    pdb_path: &str,
94) -> Result<String, Box<dyn Error>> {
95    let pairs = get_closest_residue_pairs(&pdb, *cutoff);
96
97    let mut interactors: Vec<Interactor> = Vec::new();
98    let mut counter = 0;
99    pairs.iter().for_each(|g| {
100        let mut interactor_i = Interactor::new(counter);
101        counter += 1;
102        let mut interactor_j = Interactor::new(counter);
103        interactor_j.add_target(counter - 1);
104        interactor_i.add_target(counter);
105        counter += 1;
106
107        interactor_i.set_chain(g.chain_i.as_str());
108        interactor_i.set_active(vec![g.res_i as i16]);
109        interactor_i.set_active_atoms(vec![g.atom_i.clone()]);
110        interactor_i.set_target_distance(g.distance);
111
112        interactor_j.set_chain(g.chain_j.as_str());
113        interactor_j.set_passive(vec![g.res_j as i16]);
114        interactor_j.set_passive_atoms(vec![g.atom_j.clone()]);
115
116        interactors.push(interactor_i);
117        interactors.push(interactor_j);
118    });
119
120    // Make the restraints
121    let air = Air::new(interactors);
122    let tbl = air.gen_tbl().unwrap();
123
124    println!("{}", tbl);
125
126    if let Some(output_f) = pml {
127        air.gen_pml(output_f, Some(pdb_path))
128    };
129
130    Ok(tbl)
131}
132
133/// Lists the interface residues for each chain in a protein structure.
134///
135/// This function analyzes a PDB file to identify the interface residues between chains
136/// based on a specified distance cutoff, and prints the results.
137///
138/// # Arguments
139///
140/// * `input_file` - A string slice that holds the path to the input PDB file.
141/// * `cutoff` - A reference to a f64 value specifying the distance cutoff (in Angstroms)
142///   for determining interface residues.
143///
144/// # Returns
145///
146/// A `Result<(), Box<dyn Error>>` which is Ok(()) if the function completes successfully,
147/// or an Error if something goes wrong.
148///
149pub fn list_interface(pdb: pdbtbx::PDB, cutoff: &f64) -> Result<(), Box<dyn Error>> {
150    let true_interface = get_true_interface(&pdb, *cutoff);
151
152    for (chain_id, residues) in true_interface.iter() {
153        let mut sorted_res = residues.iter().collect::<Vec<_>>();
154        sorted_res.sort();
155
156        println!("Chain {}: {:?}", chain_id, sorted_res);
157    }
158
159    Ok(())
160}
161#[cfg(test)]
162mod tests {
163
164    use super::*;
165    use std::io::{BufReader, Cursor};
166
167    #[test]
168    fn test_true_interface() {
169        let expected_tbl = r#"assign ( resid 933 and segid A )
170       (
171        ( resid 46 and segid B )
172     or
173        ( resid 47 and segid B )
174       ) 2.0 2.0 0.0
175
176assign ( resid 950 and segid A )
177       (
178        ( resid 46 and segid B )
179     or
180        ( resid 47 and segid B )
181       ) 2.0 2.0 0.0
182
183assign ( resid 46 and segid B )
184       (
185        ( resid 933 and segid A )
186     or
187        ( resid 950 and segid A )
188       ) 2.0 2.0 0.0
189
190assign ( resid 47 and segid B )
191       (
192        ( resid 933 and segid A )
193     or
194        ( resid 950 and segid A )
195       ) 2.0 2.0 0.0
196
197"#;
198
199        let content = std::fs::read_to_string("tests/data/complex.pdb").unwrap();
200        let mut opts = pdbtbx::ReadOptions::new();
201        opts.set_format(pdbtbx::Format::Pdb)
202            .set_level(pdbtbx::StrictnessLevel::Loose);
203        let cursor = Cursor::new(content.into_bytes());
204        let reader = BufReader::new(cursor);
205        let (pdb, _) = opts.read_raw(reader).unwrap();
206
207        let opt: Option<String> = None;
208        match true_interface(pdb, &3.0, &opt, "tests/data/complex.pdb") {
209            Ok(tbl) => assert_eq!(tbl, expected_tbl),
210            Err(_e) => (),
211        };
212    }
213    #[test]
214    fn test_true_interface_ba() {
215        let expected_tbl = r#"assign ( resid 933 and segid A )
216       (
217        ( resid 46 and segid B )
218     or
219        ( resid 47 and segid B )
220       ) 2.0 2.0 0.0
221
222assign ( resid 950 and segid A )
223       (
224        ( resid 46 and segid B )
225     or
226        ( resid 47 and segid B )
227       ) 2.0 2.0 0.0
228
229assign ( resid 46 and segid B )
230       (
231        ( resid 933 and segid A )
232     or
233        ( resid 950 and segid A )
234       ) 2.0 2.0 0.0
235
236assign ( resid 47 and segid B )
237       (
238        ( resid 933 and segid A )
239     or
240        ( resid 950 and segid A )
241       ) 2.0 2.0 0.0
242
243"#;
244
245        let content = std::fs::read_to_string("tests/data/complex_BA.pdb").unwrap();
246        let mut opts = pdbtbx::ReadOptions::new();
247        opts.set_format(pdbtbx::Format::Pdb)
248            .set_level(pdbtbx::StrictnessLevel::Loose);
249        let cursor = Cursor::new(content.into_bytes());
250        let reader = BufReader::new(cursor);
251        let (pdb, _) = opts.read_raw(reader).unwrap();
252
253        let opt: Option<String> = None;
254        match true_interface(pdb, &3.0, &opt, "tests/data/complex_BA.pdb") {
255            Ok(tbl) => assert_eq!(tbl, expected_tbl),
256            Err(_e) => (),
257        };
258    }
259
260    #[test]
261    fn test_unambigti() {
262        let content = std::fs::read_to_string("tests/data/two_res.pdb").unwrap();
263        let mut opts = pdbtbx::ReadOptions::new();
264        opts.set_format(pdbtbx::Format::Pdb)
265            .set_level(pdbtbx::StrictnessLevel::Loose);
266        let cursor = Cursor::new(content.into_bytes());
267        let reader = BufReader::new(cursor);
268        let (pdb, _) = opts.read_raw(reader).unwrap();
269
270        let opt: Option<String> = None;
271        let expected_tbl = "assign ( resid 2 and segid A and name CA ) ( resid 10 and segid B and name CA ) 9.1 2.0 0.0\n\n";
272
273        match unambig_ti(pdb, &5.0, &opt, "tests/data/two_res.pdb") {
274            Ok(tbl) => assert_eq!(tbl, expected_tbl),
275            Err(_e) => (),
276        }
277    }
278}
haddock_restraints/core/commands/ti.rs

haddock_restraints/core/commands/
ti.rs
