Calpha only import for PDBReader (plus unit tests)

git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2099 5a81b35b-ba03-0410-adc8-b2c5c5119f08

modules/io/pymod/export_pdb_io.cc

@@ -53,6 +53,7 @@ void export_pdb_io()
     pdb_scope.attr("SKIP_FAULTY_RECORDS")=PDB::SKIP_FAULTY_RECORDS;
     pdb_scope.attr("WRITE_MULTIPLE_MODELS")=PDB::WRITE_MULTIPLE_MODELS;
     pdb_scope.attr("JOIN_SPREAD_ATOM_RECORDS")=PDB::JOIN_SPREAD_ATOM_RECORDS;
+    pdb_scope.attr("CALPHA_ONLY")=PDB::CALPHA_ONLY;    
   }
 
   class_<PDBReader, boost::noncopyable>("PDBReader", init<String>())

modules/io/src/mol/pdb_reader.cc

@@ -283,8 +283,18 @@ bool PDBReader::ParseAtomIdent(const StringRef& line, int line_num,
     }
     throw IOException(str(format("invalid res number on line %d") % line_num));
   }
+
   char  ins_c=line[26];  
   resnum=to_res_num(res_num.second, ins_c);
+  if (PDB::Flags() & PDB::CALPHA_ONLY) {
+    if (record_type[0]=='H' || record_type[0]=='h') {
+      return false;
+    } 
+    if (atom_name!=StringRef("CA", 2)) {
+      return false;
+    }
+    return true;
+  }
   return true;
 }
 
@@ -301,7 +311,7 @@ void PDBReader::ParseAnisou(const StringRef& line, int line_num,
   }
   double anisou[6]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};  
   for (int i=0;i<6; ++i) {
-    std::pair<bool, int> result=line.substr(29+i*7, 6).to_int();
+    std::pair<bool, int> result=line.substr(29+i*7, 6).ltrim().to_int();
     if (!result.first) {
       if (PDB::Flags() & PDB::SKIP_FAULTY_RECORDS) {
         return;
@@ -312,7 +322,8 @@ void PDBReader::ParseAnisou(const StringRef& line, int line_num,
   }
   String aname(atom_name.str());
   if (!curr_residue_.IsValid()) {
-    if (PDB::Flags() & PDB::SKIP_FAULTY_RECORDS) {
+    if (PDB::Flags() & PDB::SKIP_FAULTY_RECORDS || 
+        PDB::Flags() & PDB::CALPHA_ONLY) {
       return;
     } 
     const char* fmt_str="invalid ANISOU record for inexistent atom on line %d";
@@ -320,7 +331,8 @@ void PDBReader::ParseAnisou(const StringRef& line, int line_num,
   }
   mol::AtomHandle atom=curr_residue_.FindAtom(aname);
   if (!atom.IsValid()) {
-    if (PDB::Flags() & PDB::SKIP_FAULTY_RECORDS) {
+    if (PDB::Flags() & PDB::SKIP_FAULTY_RECORDS ||
+        PDB::Flags() & PDB::CALPHA_ONLY) {
       return;
     } 
     const char* fmt_str="invalid ANISOU record for inexistent atom on line %d";

modules/io/tests/test_io_pdb.cc

@@ -56,7 +56,8 @@ BOOST_AUTO_TEST_CASE(atom_record)
   BOOST_REQUIRE_EQUAL(ent.GetChainCount(), 2);
   BOOST_REQUIRE_EQUAL(ent.GetResidueCount(), 2);
   BOOST_REQUIRE_EQUAL(ent.GetAtomCount(), 2);
-  mol::AtomHandle a1=ent.FindChain("A").GetAtomList()[0];
+  mol::AtomHandle a1=ent.FindAtom("A", mol::ResNum(1), "N");
+  BOOST_REQUIRE(a1.IsValid());
   BOOST_CHECK_EQUAL(a1.GetName(), "N");
   BOOST_CHECK_EQUAL(a1.GetResidue().GetName(), "MET");
   BOOST_CHECK_EQUAL(a1.GetResidue().GetChain().GetName(), "A");  
@@ -66,7 +67,8 @@ BOOST_AUTO_TEST_CASE(atom_record)
   BOOST_CHECK_EQUAL(a1.GetProp().occupancy, 0.5);  
   BOOST_CHECK_EQUAL(a1.GetProp().element, "N");
   BOOST_CHECK_EQUAL(a1.GetProp().is_hetatm, false);
-  mol::AtomHandle a2=ent.FindChain(" ").GetAtomList()[0];
+  mol::AtomHandle a2=ent.FindAtom(" ", mol::ResNum(1), "CA");
+  BOOST_REQUIRE(a2.IsValid());  
   BOOST_CHECK_EQUAL(a2.GetName(), "CA");
   BOOST_CHECK_EQUAL(a2.GetResidue().GetName(), "MET");
   BOOST_CHECK_EQUAL(a2.GetResidue().GetChain().GetName(), " ");  
@@ -120,9 +122,53 @@ BOOST_AUTO_TEST_CASE(join_spread_records_off)
   BOOST_CHECK_EQUAL(res1.GetAtomCount(), 1);
   BOOST_CHECK(res1.FindAtom("N"));
   mol::ResidueHandle res2=ent.FindResidue("A", mol::ResNum(2));  
-  BOOST_CHECK(res2.IsValid());
+  BOOST_REQUIRE(res2.IsValid());
   BOOST_CHECK_EQUAL(res2.GetAtomCount(), 1);
   BOOST_CHECK(res2.FindAtom("N"));
 }
 
+BOOST_AUTO_TEST_CASE(calpha_only_import_on)
+{
+  String fname("testfiles/pdb/calpha.pdb");
+  PDBReader reader(fname);
+  PDB::PushFlags(PDB::CALPHA_ONLY);
+  mol::EntityHandle ent=mol::CreateEntity();
+  reader.Import(ent);
+  PDB::PopFlags();
+  BOOST_CHECK_EQUAL(ent.GetResidueCount(), 1);
+  BOOST_CHECK_EQUAL(ent.GetAtomCount(), 1);
+}
+
+BOOST_AUTO_TEST_CASE(calpha_only_import_off)
+{
+  String fname("testfiles/pdb/calpha.pdb");
+  PDBReader reader(fname);
+  mol::EntityHandle ent=mol::CreateEntity();
+  reader.Import(ent);
+  BOOST_CHECK_EQUAL(ent.GetResidueCount(), 2);
+  BOOST_CHECK_EQUAL(ent.GetAtomCount(), 4);
+}
+
+BOOST_AUTO_TEST_CASE(anisou_record)
+{
+  String fname("testfiles/pdb/anisou.pdb");
+  PDBReader reader(fname);
+  mol::EntityHandle ent=mol::CreateEntity();
+  reader.Import(ent);
+  BOOST_REQUIRE(ent.GetAtomCount()==1);
+  mol::AtomHandle a1=ent.FindAtom("A", mol::ResNum(7), "N");
+  BOOST_REQUIRE(a1.IsValid());
+  mol::AtomProp props=a1.GetProp();
+  BOOST_CHECK_CLOSE( 0.0100, props.anisou(0, 0), 1e-4);
+  BOOST_CHECK_CLOSE(-0.0016, props.anisou(1, 0), 1e-4);
+  BOOST_CHECK_CLOSE(-0.0026, props.anisou(2, 0), 1e-4);
+  BOOST_CHECK_CLOSE(-0.0016, props.anisou(0, 1), 1e-4);
+  BOOST_CHECK_CLOSE( 0.0110, props.anisou(1, 1), 1e-4);
+  BOOST_CHECK_CLOSE(-0.0054, props.anisou(2, 1), 1e-4);
+  BOOST_CHECK_CLOSE(-0.0026, props.anisou(0, 2), 1e-4);
+  BOOST_CHECK_CLOSE(-0.0054, props.anisou(1, 2), 1e-4);
+  BOOST_CHECK_CLOSE( 0.0120, props.anisou(2, 2), 1e-4);    
+}
+
+
 BOOST_AUTO_TEST_SUITE_END()

modules/io/tests/testfiles/pdb/anisou.pdb

+ATOM    134  N   GLU A   7      -3.391  15.554   9.400  1.00  3.46           C
+ANISOU  133  N   GLU A   7      100    110    120    -16    -26    -54       N  

modules/io/tests/testfiles/pdb/calpha.pdb

+ATOM    134  CA  GLU A   7      -3.391  15.554   9.400  1.00  3.46           C  
+ANISOU  134  CA  GLU A   7      443    511    361    -28    -37    -88       C  
+ATOM    135  C   GLU A   7      -3.799  15.410  10.865  1.00  3.42           C  
+ANISOU  135  C   GLU A   7      419    507    375    -14    -45   -103       C  
+ATOM    136  O   GLU A   7      -4.948  15.082  11.188  1.00  4.00           O
+HETATM 1469 CA    CA A1000      46.224   7.128  -8.698  1.00  7.59          CA
\ No newline at end of file