diff --git a/modules/bindings/pymod/blast.py b/modules/bindings/pymod/blast.py
index 08359de4d5aba0f473a248aca8926a6bedf6d8c2..e63007f40d5b9feb0fb0c838166f35b3142ad24b 100644
--- a/modules/bindings/pymod/blast.py
+++ b/modules/bindings/pymod/blast.py
@@ -149,21 +149,20 @@ def CreateDB(infasta, dbout, mkdb_cmd=None):
if mkdb_cmd==None:
try:
exe=settings.Locate('formatdb')
- args=[exe, '-i', infasta, '-v',str(10000),'-n',dbout]
+ args=[exe, '-i', infasta, '-n', dbout]
except:
try:
exe=settings.Locate('makeblastdb')
- args=[exe, '-in', infasta, ' -max_file_sz 10GB -out', dbout,'-dbtype','prot']
-
+ args=[exe, '-in', infasta, '-out', dbout, '-dbtype', 'prot']
except:
raise RuntimeError('could not find makeblastdb/formatdb executable')
else:
if os.path.basename(mkdb_cmd)=='makeblastdb':
exe=settings.Locate('makeblastdb',explicit_file_name=mkdb_cmd)
- args=[exe, '-in', infasta, ' -max_file_sz 10GB -out', dbout,'-dbtype','prot']
+ args=[exe, '-in', infasta, '-out', dbout, '-dbtype', 'prot']
elif os.path.basename(mkdb_cmd)=='formatdb':
exe=settings.Locate('formatdb',explicit_filename=mkdb_cmd)
- args=[exe, '-i', infasta, '-v',str(10000),'-n',dbout]
+ args=[exe, '-i', infasta, '-n', dbout]
else:
raise IOError('mkdb command must either be the path to formatdb or makeblastdb!')
@@ -186,11 +185,11 @@ def BlastVersion(blast_location=None):
if os.path.basename(blast_exe)=='blastall':
args=[blast_exe]
- pattern=re.compile(r'blastall (\d+\.\d+\.\d+)\s+arguments:\s*')
+ pattern=re.compile(r'\s*blastall (\d+\.\d+\.\d+)\s+arguments:\s*')
else:
args=[blast_exe, '-version']
- pattern=re.compile(r'Package: blast (\d+\.\d+\.\d+),\s+')
+ pattern=re.compile(r'\s*Package: blast (\d+\.\d+\.\d+),\s+')
blast_pipe=subprocess.Popen(args, stdout=subprocess.PIPE,
stderr=subprocess.PIPE)
diff --git a/modules/bindings/pymod/hhblits.py b/modules/bindings/pymod/hhblits.py
index 0a1045171c167a5e24e958d4256739443b37970b..a203a292050ec7ad663ec7f0ac555c81d3a370ae 100644
--- a/modules/bindings/pymod/hhblits.py
+++ b/modules/bindings/pymod/hhblits.py
@@ -543,6 +543,8 @@ class HHblits:
"""
if a3m_file is None:
a3m_file = '%s.a3m' % os.path.splitext(self.filename)[0]
+ else:
+ a3m_file = os.path.abspath(a3m_file)
if os.path.exists(a3m_file):
ost.LogInfo('Reusing already existing query alignment (%s)' % a3m_file)
return a3m_file
@@ -720,7 +722,7 @@ class HHblits:
self.hhblits_bin,
opt_cmd,
os.path.abspath(a3m_file),
- hhr_file,
+ os.path.abspath(hhr_file),
os.path.join(os.path.abspath(os.path.split(database)[0]),
os.path.split(database)[1]))
ost.LogInfo('searching %s' % database)
diff --git a/modules/bindings/tests/test_blast.py b/modules/bindings/tests/test_blast.py
index 48cb560685d5dd46dfb6682c6bfe5271cc7ca9c0..18311abd9342e4b8408fc4c2a2d30534ef08ed12 100644
--- a/modules/bindings/tests/test_blast.py
+++ b/modules/bindings/tests/test_blast.py
@@ -32,21 +32,23 @@ class TestBlastBindings(unittest.TestCase):
def testBlastExec(self):
- hits=blast.Blast(self.query, 'testfiles/seqdb')
-
- expected_hits=['1griA','1fhsA','3n84E','1mw9X']
- self.assertEqual(len(hits), 4)
-
- found=0
- hit_ids=[]
+ hits = blast.Blast(self.query, 'testfiles/seqdb')
+ expected_perfect_hit_ids = ["1fhsA", "3n84E", "1griA"]
+ perfect_hit_ids = []
+
for h in hits:
- hit_ids.append(h.identifier)
-
- for identifier in expected_hits:
- if identifier in hit_ids:
- found+=1
- self.assertEqual(found, 4)
+ perfect_match = False
+ for aln_patch in h.aligned_patches:
+ if aln_patch.seqid == 1.0:
+ perfect_match = True
+ if perfect_match:
+ perfect_hit_ids.append(h.identifier)
+
+ expected_perfect_hit_ids.sort()
+ perfect_hit_ids.sort()
+
+ self.assertEqual(perfect_hit_ids, expected_perfect_hit_ids)
def testBlastParseOutput(self):
@@ -54,6 +56,7 @@ class TestBlastBindings(unittest.TestCase):
parsed_out=blast.ParseBlastOutput(raw_out)
+ # numbers and ids can be specific to a certain blast version
expected_ids=['1griA','1fhsA','3n84E','1mw9X']
expected_evalues=[4.80893e-59,4.80893e-59,9.06925e-58,2.96523]
expected_scores=[534,534,523,27]
diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst
index 795d4f96e2613fcb50c66abcde51ffe91075f724..51de4cca63a8b7f794d4a671ca5cff1b779d7cdb 100644
--- a/modules/mol/alg/doc/molalg.rst
+++ b/modules/mol/alg/doc/molalg.rst
@@ -685,18 +685,20 @@ Local Distance Test scores (lDDT, DRMSD)
Dictionary-like object containing the list of interatomic distances that
originate from a single residue to be checked during a run of the Local
Distance Difference Test algorithm
- (key = pair of :class:`UniqueAtomIdentifier`, value = pair of floats)
+ (key = pair of :class:`UniqueAtomIdentifier`, value = pair of floats
+ representing min and max distance observed in the structures used to build
+ the map).
.. class:: GlobalRDMap
- Dictionary-like object containing all the :class:`~ost.mol.alg.ResidueRDMap` objects related to residues of a single structure
- (key = :class:`~ost.mol.ResNum`, value = :class:`ResidueRDMap`)
+ Dictionary-like object containing all the :class:`~ost.mol.alg.ResidueRDMap` objects related to all the residues
+ (key = :class:`~ost.mol.ResNum`, value = :class:`ResidueRDMap`).
.. function:: PrintResidueRDMap(residue_distance_list)
Prints to standard output all the distances contained in a
- :class:`~ost.mol.alg.ResidueRDMap` object
+ :class:`~ost.mol.alg.ResidueRDMap` object.
.. function:: PrintGlobalRDMap(global_distance_list)
diff --git a/modules/mol/alg/src/filter_clashes.cc b/modules/mol/alg/src/filter_clashes.cc
index e0d2bee05a80d0959b3528776d501151f2b388c0..8623d8779342b089366d6d43bf835cdc1c51377b 100644
--- a/modules/mol/alg/src/filter_clashes.cc
+++ b/modules/mol/alg/src/filter_clashes.cc
@@ -404,19 +404,28 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
}
if (other_atom.GetHashCode() > atom.GetHandle().GetHashCode()) {
Real blength = bond.GetLength();
- String bond_str = bond_string(atom,other_atom);
- std::pair<Real,Real> length_stddev = bond_table.GetParam(bond_str,residue_str);
+ String bond_str = bond_string(atom, other_atom);
+ std::pair<Real,Real> length_stddev = bond_table.GetParam(bond_str, residue_str);
Real ref_length = length_stddev.first;
Real ref_stddev = length_stddev.second;
Real min_length = ref_length - bond_tolerance*ref_stddev;
Real max_length = ref_length + bond_tolerance*ref_stddev;
Real zscore = (blength - ref_length)/ref_stddev;
if (blength < min_length || blength > max_length) {
- LOG_INFO("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "FAIL")
+ LOG_INFO("BOND:" << " " << res.GetChain() << " "
+ << res.GetName() << " " << res.GetNumber() << " "
+ << bond_str << " " << min_length << " " << max_length
+ << " " << blength << " " << zscore << " " << "FAIL")
bad_bond_count++;
- UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),atom.GetResidue().GetNumber(),atom.GetResidue().GetName(),atom.GetName());
- UniqueAtomIdentifier other_atom_ui(other_atom.GetResidue().GetChain().GetName(),other_atom.GetResidue().GetNumber(),other_atom.GetResidue().GetName(),other_atom.GetName());
- StereoChemicalBondViolation bond_v(atom_ui,other_atom_ui,blength,std::make_pair(min_length,max_length));
+ UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),
+ atom.GetResidue().GetNumber(),
+ atom.GetResidue().GetName(), atom.GetName());
+ UniqueAtomIdentifier other_atom_ui(other_atom.GetResidue().GetChain().GetName(),
+ other_atom.GetResidue().GetNumber(),
+ other_atom.GetResidue().GetName(),
+ other_atom.GetName());
+ StereoChemicalBondViolation bond_v(atom_ui, other_atom_ui, blength,
+ std::make_pair(min_length, max_length));
bond_violation_list.push_back(bond_v);
remove_sc=true;
if (always_remove_bb==true) {
@@ -432,7 +441,10 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
}
}
} else {
- LOG_VERBOSE("BOND:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << bond_str << " " << min_length << " " << max_length << " " << blength << " " << zscore << " " << "PASS")
+ LOG_VERBOSE("BOND:" << " " << res.GetChain() << " "
+ << res.GetName() << " " << res.GetNumber() << " "
+ << bond_str << " " << min_length << " " << max_length
+ << " " << blength << " " << zscore << " " << "PASS")
}
bond_count++;
running_sum_zscore_bonds+=zscore;
@@ -486,12 +498,25 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
Real max_width = ref_width + angle_tolerance*ref_stddev;
Real zscore = (awidth - ref_width)/ref_stddev;
if (awidth < min_width || awidth > max_width) {
- LOG_INFO("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "FAIL")
+ LOG_INFO("ANGLE:" << " " << res.GetChain() << " "
+ << res.GetName() << " " << res.GetNumber() << " "
+ << angle_str << " " << min_width << " " << max_width
+ << " " << awidth << " " << zscore << " " << "FAIL")
bad_angle_count++;
- UniqueAtomIdentifier atom1_ui(atom1.GetResidue().GetChain().GetName(),atom1.GetResidue().GetNumber(),atom1.GetResidue().GetName(),atom1.GetName());
- UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),atom.GetResidue().GetNumber(),atom.GetResidue().GetName(),atom.GetName());
- UniqueAtomIdentifier atom2_ui(atom2.GetResidue().GetChain().GetName(),atom2.GetResidue().GetNumber(),atom2.GetResidue().GetName(),atom2.GetName());
- StereoChemicalAngleViolation angle_v(atom1_ui,atom_ui,atom2_ui,awidth,std::make_pair(min_width,max_width));
+ UniqueAtomIdentifier atom1_ui(atom1.GetResidue().GetChain().GetName(),
+ atom1.GetResidue().GetNumber(),
+ atom1.GetResidue().GetName(),
+ atom1.GetName());
+ UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),
+ atom.GetResidue().GetNumber(),
+ atom.GetResidue().GetName(),
+ atom.GetName());
+ UniqueAtomIdentifier atom2_ui(atom2.GetResidue().GetChain().GetName(),
+ atom2.GetResidue().GetNumber(),
+ atom2.GetResidue().GetName(),
+ atom2.GetName());
+ StereoChemicalAngleViolation angle_v(atom1_ui,atom_ui,atom2_ui,awidth,
+ std::make_pair(min_width,max_width));
angle_violation_list.push_back(angle_v);
remove_sc=true;
if (always_remove_bb==true) {
@@ -502,7 +527,10 @@ std::pair<EntityView,StereoChemistryInfo> CheckStereoChemistry(const EntityView&
remove_bb=true;
}
} else {
- LOG_VERBOSE("ANGLE:" << " " << res.GetChain() << " " << res.GetName() << " " << res.GetNumber() << " " << angle_str << " " << min_width << " " << max_width << " " << awidth << " " << zscore << " " << "PASS")
+ LOG_VERBOSE("ANGLE:" << " " << res.GetChain() << " "
+ << res.GetName() << " " << res.GetNumber() << " "
+ << angle_str << " " << min_width << " " << max_width
+ << " " << awidth << " " << zscore << " " << "PASS")
}
angle_count++;
running_sum_zscore_angles+=zscore;
@@ -610,7 +638,14 @@ std::pair<EntityView,ClashingInfo> FilterClashes(const EntityView& ent, const Cl
Real threshold=distance-tolerance;
if (d<threshold*threshold) {
- LOG_INFO(atom.GetResidue().GetChain() << " " << atom.GetResidue().GetName() << " " << atom.GetResidue().GetNumber() << " " << atom.GetName() << " " << atom2.GetResidue().GetChain() << " " << atom2.GetResidue().GetName() << " " << atom2.GetResidue().GetNumber() << " " << atom2.GetName() << " " << threshold << " " << sqrt(d) << " " << sqrt(d)-threshold << " " << "FAIL")
+ LOG_INFO("CLASH:" << " " << atom.GetResidue().GetChain() << " "
+ << atom.GetResidue().GetName() << " "
+ << atom.GetResidue().GetNumber() << " " << atom.GetName()
+ << " " << atom2.GetResidue().GetChain() << " "
+ << atom2.GetResidue().GetName() << " "
+ << atom2.GetResidue().GetNumber() << " " << atom2.GetName()
+ << " " << threshold << " " << sqrt(d) << " "
+ << sqrt(d)-threshold << " " << "FAIL")
bad_distance_count++;
UniqueAtomIdentifier atom_ui(atom.GetResidue().GetChain().GetName(),atom.GetResidue().GetNumber(),atom.GetResidue().GetName(),atom.GetName());
UniqueAtomIdentifier atom2_ui(atom2.GetResidue().GetChain().GetName(),atom2.GetResidue().GetNumber(),atom2.GetResidue().GetName(),atom2.GetName());
@@ -626,7 +661,14 @@ std::pair<EntityView,ClashingInfo> FilterClashes(const EntityView& ent, const Cl
remove_bb=true;
}
} else {
- LOG_VERBOSE("CLASH:" << " " << atom.GetResidue().GetChain() << " " << atom.GetResidue().GetName() << " " << atom.GetResidue().GetNumber() << " " << atom.GetName() << " " << atom2.GetResidue().GetChain() << " " << atom2.GetResidue().GetNumber() << " " << atom2.GetResidue().GetName() << " " << atom2.GetName() << " " << threshold << " " << sqrt(d) << " " << sqrt(d)-threshold << " " << "PASS")
+ LOG_VERBOSE("CLASH:" << " " << atom.GetResidue().GetChain() << " "
+ << atom.GetResidue().GetName() << " "
+ << atom.GetResidue().GetNumber() << " " << atom.GetName()
+ << " " << atom2.GetResidue().GetChain() << " "
+ << atom2.GetResidue().GetName() << " "
+ << atom2.GetResidue().GetNumber() << " " << atom2.GetName()
+ << " " << threshold << " " << sqrt(d) << " "
+ << sqrt(d)-threshold << " " << "PASS")
}
distance_count++;
}
diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc
index a64f5bb192d5abab6777b750786858d1fd477d16..dbe122d80d0b13e250edaec00291ea268bfa2bc6 100644
--- a/modules/mol/alg/src/lddt.cc
+++ b/modules/mol/alg/src/lddt.cc
@@ -311,11 +311,14 @@ int main (int argc, char **argv)
settings.PrintParameters();
if (structural_checks) {
LOG_INFO("Log entries format:");
+ // verbosity level/format must match filter_clashes::CheckStereoChemistry
LOG_INFO("BOND INFO FORMAT: Chain Residue ResNum Bond Min Max Observed Z-score Status");
LOG_INFO("ANGLE INFO FORMAT: Chain Residue ResNum Angle Min Max Observed Z-score Status");
- LOG_INFO("CLASH INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 Observed Difference Status");
+ // verbosity level/format must match filter_clashes::FilterClashes
+ LOG_INFO("CLASH INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 Threshold Observed Difference Status");
}
- LOG_INFO("LDDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 ModelDist TargetDist Difference Tolerance Status");
+ // verbosity level/format must match local_dist_diff_test::calc_overlap1
+ LOG_VERBOSE("LDDT INFO FORMAT: Chain1 Residue1 ResNum1 Atom1 Chain2 Residue2 ResNum2 Atom2 ModelDist TargetDistMin TargetDistMax Tolerance Status");
// error if the reference structure is empty
if (glob_dist_list.size()==0) {
diff --git a/modules/mol/alg/src/local_dist_diff_test.cc b/modules/mol/alg/src/local_dist_diff_test.cc
index c52bf115666962310aef3de44b126eb5dd300e9b..1182dc375e8c44c01685cc55341ea35269c776fd 100644
--- a/modules/mol/alg/src/local_dist_diff_test.cc
+++ b/modules/mol/alg/src/local_dist_diff_test.cc
@@ -27,7 +27,8 @@ String vector_to_string(std::vector<Real> vec) {
// helper function
bool within_tolerance(Real mdl_dist, const std::pair<Real,Real>& values, Real tol)
{
-
+ // this considers min & max distance observed in all reference structures
+ // -> for single reference case values.first and values.second are the same
return (values.first-tol)<=mdl_dist && (values.second+tol)>=mdl_dist;
}
@@ -82,10 +83,15 @@ std::pair<long int, long int> calc_overlap1(const ResidueRDMap& res_distance_lis
Real tol = * tol_list_it;
if (within_tolerance(mdl_dist,values,tol)) {
if (log) {
- LOG_VERBOSE("lddt:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
- << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
- << mdl_dist << " " << values.first << " " << values.second << " " << tol << " " << "PASS")
- }
+ LOG_VERBOSE("LDDT:" << " " << av1.GetResidue().GetChain() << " "
+ << av1.GetResidue().GetName() << " "
+ << av1.GetResidue().GetNumber() << " " << av1.GetName()
+ << " " << av2.GetResidue().GetChain() << " "
+ << av2.GetResidue().GetName() << " "
+ << av2.GetResidue().GetNumber() << " " << av2.GetName()
+ << " " << mdl_dist << " " << values.first << " "
+ << values.second << " " << tol << " " << "PASS")
+ }
overlap.first+=1;
if (!only_fixed) {
overlap_list[rindex1].first+=1;
@@ -93,10 +99,15 @@ std::pair<long int, long int> calc_overlap1(const ResidueRDMap& res_distance_lis
}
} else {
if (log) {
- LOG_VERBOSE("lddt:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName()
- << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " "
- << mdl_dist << " " << values.first << " " << values.second << " " << tol << " " << "FAIL");
- }
+ LOG_VERBOSE("LDDT:" << " " << av1.GetResidue().GetChain() << " "
+ << av1.GetResidue().GetName() << " "
+ << av1.GetResidue().GetNumber() << " " << av1.GetName()
+ << " " << av2.GetResidue().GetChain() << " "
+ << av2.GetResidue().GetName() << " "
+ << av2.GetResidue().GetNumber() << " " << av2.GetName()
+ << " " << mdl_dist << " " << values.first << " "
+ << values.second << " " << tol << " " << "FAIL");
+ }
break;
}
}
diff --git a/modules/mol/alg/src/molck.cc b/modules/mol/alg/src/molck.cc
index f0b8f946d0a66f4549699db314f9c3f240c6544d..b3cfe016c342f9edbf87eca2a8b1d3524515c285 100644
--- a/modules/mol/alg/src/molck.cc
+++ b/modules/mol/alg/src/molck.cc
@@ -80,7 +80,7 @@ void ost::mol::alg::RemoveAtoms(
edi.DeleteAtom(*i);
zremoved++;
}
- std::cerr << " --> removed " << zremoved << " hydrogen atoms" << std::endl;
+ std::cerr << " --> removed " << zremoved << " atoms with zero occupancy" << std::endl;
}
if (rm_hyd_atoms) {