SCHWED-776: added unit tests for full modelling pipeline (full) with oligomers.

modelling/tests/CMakeLists.txt

@@ -15,8 +15,14 @@ set(MODELLING_TEST_DATA
   data/1crn_sc.pdb
   data/1ITX-A-11_sidechain.pdb
   data/1mcg.pdb
+  data/2aoh-1_cut.pdb
+  data/2aoh-1_cut_A.fasta
+  data/2aoh-1_cut_B.fasta
   data/2dbs.pdb
   data/4R6K-A-13_sidechain.pdb
+  data/5d52-1_cut.pdb
+  data/5d52-1_cut_A.fasta
+  data/5d52-1_cut_B.fasta
   data/cbeta.fasta
   data/cbeta.pdb
   data/compounds.chemlib

modelling/tests/data/2aoh-1_cut_A.fasta

+>trg
+PQITLWKRPLVTIKIGGQLKEALLDTGADDTVIEEMSLPGRWKPKMIGGIGGFIKVRQYDQIIIEIAGHKAIGTVLVGPTPANIIGRNLLTQIGATLNF
+>A
+PQITL------TIKIGGQLKEALLDTGADDTVIEEMSLPGRWKPKMIGGIGGFIKVRQYDQIIIEIAGHKAIGTVLVGPTPANIIGRNLLTQIGATLNF

modelling/tests/data/2aoh-1_cut_B.fasta

+>trg
+PQITLWKRPLVTIKIGGQLKEALLDTGADDTVIEEMSLPGRWKPKMIGGIGGFIKVRQYDQIIIEIAGHKAIGTVLVGPTPANIIGRNLLTQIGATLNF
+>B
+PQITLWKRPLVTIKIG------LLDTGADDTVIEEMSLPGRWKPKMIGGIGGFIKVRQYDQIIIEIAGHKAIGTVLVGPTPANIIGRNLLTQIGATLNF

modelling/tests/data/5d52-1_cut_A.fasta

+>trg
+GIVEQCCTSICSLYQLENYCN
+>A
+GIVEQ------SLYQLENYCN

modelling/tests/data/5d52-1_cut_B.fasta

+>trg
+FVNQHLCGSHLVEALYLVCGERGFFYTPKA
+>B
+FVNQHLCGSHLVEALY------GFFYTPKA

modelling/tests/test_close_gaps.py

@@ -4,8 +4,7 @@ Unit tests for modelling components to close gaps.
 import unittest
 import ost
 from ost import io, seq, geom
-from promod3 import modelling
-from promod3 import loop
+from promod3 import loop, modelling
 
 # setting up an OST LogSink to capture messages
 class _FetchLog(ost.LogSink):

modelling/tests/test_pipeline.py

@@ -4,8 +4,9 @@ We go from a fake cut in crambin to reconstruct that same protein.
 """
 import unittest
 import math
-from promod3 import loop, modelling
+import ost
 from ost import io, mol, seq
+from promod3 import loop, modelling
 
 ################################################################
 # Code to generate 1crn_cut and 1crn.fasta:
@@ -60,8 +61,65 @@ from ost import io, mol, seq
 # modelling.MinimizeModelEnergy(mhandle)
 # io.SavePDB(mhandle.model, 'data/1crn_final.pdb')
 ################################################################
+# Code to generate oligomers
+# REQ: 2aoh-1.pdb and 5d52-1.pdb from SMTL (or here in data)
+################################################################
+# import os
+# from ost import io,mol,seq
+# # choose protein
+# out_path = "oligos"
+# tpl_id = '5d52-1' # hetero
+# #tpl_id = '2aoh-1' # homo (with useless chain C)
+# # load protein 
+# prot = io.LoadPDB(os.path.join(out_path, tpl_id + ".pdb"))
+# # get only amino acids
+# prot = mol.CreateEntityFromView(prot.Select("peptide=true"), True)
+# # get chain names and seqres
+# nameA = prot.chains[0].name
+# nameB = prot.chains[1].name
+# seqresA = ''.join([r.one_letter_code for r in prot.chains[0].residues])
+# seqresB = ''.join([r.one_letter_code for r in prot.chains[1].residues])
+# # choose what to cut (rnum=X1..X2 removed)
+# rnumA1 = 6
+# rnumA2 = 11
+# rnumB1 = 17
+# rnumB2 = 22
+# # cut out stuff for fake pdb file
+# myquery = "(cname=%s and (rnum < %d or rnum > %d))" \
+#           " or (cname= %s and (rnum < %d or rnum > %d))" \
+#           % (nameA, rnumA1, rnumA2, nameB, rnumB1, rnumB2)
+# prot_cut = mol.CreateEntityFromView(prot.Select(myquery), True)
+# io.SavePDB(prot_cut, os.path.join(out_path, tpl_id + "_cut.pdb"))
+# # create alignments
+# seqresA_cut = seqresA[:rnumA1-1] + '-'*(rnumA2-rnumA1+1) + seqresA[rnumA2:]
+# seqresB_cut = seqresB[:rnumB1-1] + '-'*(rnumB2-rnumB1+1) + seqresB[rnumB2:]
+# alnA = seq.CreateAlignment(seq.CreateSequence('trg', seqresA),
+#                            seq.CreateSequence(nameA, seqresA_cut))
+# alnB = seq.CreateAlignment(seq.CreateSequence('trg', seqresB),
+#                            seq.CreateSequence(nameB, seqresB_cut))
+# # dump and report
+# print 'ALIGNMENT-A\n' + str(alnA)
+# print 'ALIGNMENT-B\n' + str(alnB)
+# io.SaveAlignment(alnA, os.path.join(out_path, tpl_id + "_cut_A.fasta"))
+# io.SaveAlignment(alnB, os.path.join(out_path, tpl_id + "_cut_B.fasta"))
+################################################################
+
+# setting up an OST LogSink to capture messages
+class _FetchLog(ost.LogSink):
+    def __init__(self):
+        ost.LogSink.__init__(self)
+        self.messages = dict()
+
+    def LogMessage(self, message, severity):
+        levels = ['ERROR', 'WARNING', 'SCRIPT', 'INFO', 'VERBOSE', 'DEBUG',
+                  'TRACE']
+        level = levels[severity]
+        if not level in self.messages.keys():
+            self.messages[level] = list()
+        self.messages[level].append(message.strip())
 
 class PipelineTests(unittest.TestCase):
+
     @classmethod
     def setUpClass(cls):
         # load dbs etc here for all tests
@@ -69,6 +127,9 @@ class PipelineTests(unittest.TestCase):
         cls.structure_db = loop.LoadStructureDB()
         cls.torsion_sampler = loop.LoadTorsionSamplerCoil()
 
+    #######################################################################
+    # HELPERs
+    #######################################################################
     def compare(self, model, filename, delta=0.01):
         '''Compare model with whatever is stored in filename.'''
         model_ref = io.LoadPDB(filename)
@@ -98,6 +159,17 @@ class PipelineTests(unittest.TestCase):
         self.raw_model = mhandle.model.Copy()
         return mhandle
 
+    def getRawModelOligo(self, file_prefix):
+        '''Get raw model for tests with oligomers.'''
+        tpl = io.LoadPDB(file_prefix + '.pdb')
+        alns = seq.AlignmentList()
+        alns.append(io.LoadAlignment(file_prefix + "_A.fasta"))
+        alns.append(io.LoadAlignment(file_prefix + "_B.fasta"))
+        alns[0].AttachView(1, tpl.Select("cname=A").CreateFullView())
+        alns[1].AttachView(1, tpl.Select("cname=B").CreateFullView())
+        mhandle = modelling.BuildRawModel(alns)
+        return mhandle
+
     def getMockModel(self, model):
         '''Get ModellingHandle for given model (assume no gaps).'''
         mhandle = modelling.ModellingHandle()
@@ -109,6 +181,35 @@ class PipelineTests(unittest.TestCase):
             mhandle.seqres.AddSequence(seq.CreateSequence(ch.name, myseq))
         return mhandle
 
+    def checkFinalModel(self, mhandle, exp_gaps=0, num_chains=0):
+        '''Check if model has exp_gaps and that CheckFinalModel reports it.'''
+        # setup
+        self.assertEqual(len(mhandle.gaps), exp_gaps)
+        log = _FetchLog()
+        ost.PushLogSink(log)
+        ost.PushVerbosityLevel(1)
+        # check
+        log.messages['WARNING'] = list()
+        modelling.CheckFinalModel(mhandle)
+        if exp_gaps == 0:
+            self.assertEqual(len(log.messages['WARNING']), 0)
+        else:
+            self.assertEqual(len(log.messages['WARNING']), 1)
+            self.assertEqual(log.messages['WARNING'][-1],
+                             "Failed to close %d gap(s). " \
+                             "Returning incomplete model!" % exp_gaps)
+            # fake remove gaps to get sequence mismatch warnings
+            mhandle.gaps = modelling.StructuralGapList()
+            log.messages['WARNING'] = list()
+            modelling.CheckFinalModel(mhandle)
+            # 1 warning per chain expected
+            exps = "Sequence mismatch in chain"
+            self.assertGreaterEqual(len(log.messages['WARNING']), num_chains)
+            for i in range(num_chains):
+                self.assertTrue(log.messages['WARNING'][i].startswith(exps))
+
+    #######################################################################
+
     def testLoopReconstruction(self):
         '''Check loop reconstruction.'''
         mhandle = self.getRawModel()
@@ -157,6 +258,43 @@ class PipelineTests(unittest.TestCase):
         # check
         self.assertEqual(mhandle.model, final_model)
         self.checkFullModel(mhandle)
+        self.checkFinalModel(mhandle)
+
+    def testHomoOligomer(self):
+        '''Check that BuildFromRawModel works for a homo oligomer.'''
+        # get raw model
+        mhandle = self.getRawModelOligo("data/2aoh-1_cut")
+        # build final model
+        final_model = modelling.BuildFromRawModel(mhandle)
+        # check
+        self.assertEqual(mhandle.model, final_model)
+        self.assertEqual(mhandle.model.residue_count, 198)
+        self.assertEqual(mhandle.model.atom_count, 1508)
+        self.checkFinalModel(mhandle)
+
+    def testHeteroOligomer(self):
+        '''Check that BuildFromRawModel works for a homo oligomer.'''
+        # get raw model
+        mhandle = self.getRawModelOligo("data/5d52-1_cut")
+        # build final model
+        final_model = modelling.BuildFromRawModel(mhandle)
+        # check
+        self.assertEqual(mhandle.model, final_model)
+        self.assertEqual(mhandle.model.residue_count, 51)
+        self.assertEqual(mhandle.model.atom_count, 403)
+        self.checkFinalModel(mhandle)
+
+    def testCheckFinalModel(self):
+        '''Check that we report unclosed models.'''
+        # single chain
+        mhandle = self.getRawModel()
+        self.checkFinalModel(mhandle, exp_gaps=1, num_chains=1)
+        # homo-dimer
+        mhandle = self.getRawModelOligo("data/2aoh-1_cut")
+        self.checkFinalModel(mhandle, exp_gaps=2, num_chains=2)
+        # hetero-dimer
+        mhandle = self.getRawModelOligo("data/5d52-1_cut")
+        self.checkFinalModel(mhandle, exp_gaps=2, num_chains=2)
 
 if __name__ == "__main__":
     from ost import testutils