Bug fixe: generated ssdist restraints were empty

ariaec/econverter.py

@@ -686,12 +686,16 @@ class AriaEcXMLConverter(AriaXMLConverter):
                                "_hbond.tbl")
         ssdist_file = os.path.join(self.settings.infra["tbl"], self.outprefix +
                                    "_ssdist.tbl")
+        logger.info("   Dihedral restraints (%s)" % dihed_file)
         self.write_dihedral_tbl(protein.sec_struct.ss_matrix, dihed_file)
+        logger.info("   Helix bond restraints (%s)" % hb_file)
         self.write_hb_tbl(protein, hb_file,
                           hbmap=hbmap, n_hb=n_hb,
                           lr_type=self.settings.setup.config['longrange_hbtype'],
                           dminus=self.settings.setup.config['hb_dminus'],
                           dplus=self.settings.setup.config['hb_dplus'])
+        logger.info("   Secondary structure restraints (%s)" %
+                    ssdist_file)
         self.write_ssdist_tbl(protein.sec_struct.ss_matrix,
                               protein.sec_struct.ssdist,
                               ssdist_file)
@@ -721,18 +725,18 @@ class AriaEcXMLConverter(AriaXMLConverter):
         :param desc:
         :return:
         """
+        logger.info("Reading aria template file %s" % aria_template)
         if aria_template:
             template = os.path.abspath(aria_template)
 
             try:
                 t = open(template, 'r')
-            except IOError:
-                sys.exit("""Can't open "%s" file.""" % template)
+            except Exception, msg:
+                sys.exit("""Can't open "%s" file. %s""" % (template, msg))
 
             aria_project_template = t.read()
             t.close()
         else:
-            logger.info("Reading default aria template file %s" % self.settings.ARIAPROJ_TEMPLATE)
             aria_project_template = pkgr.resource_string(
                 __name__, self.settings.ARIAPROJ_TEMPLATE)
 

ariaec/ecsettings.py

@@ -92,7 +92,7 @@ class AriaEcSettings(Settings):
                                               "bbconv", "contactdef"))
         self._infra = {}
         self._scsc_min = None
-        self._ss_dist = None
+        self._ssdist = None
         self._template = None
         self.outdir = os.getcwd()
         self.name = name
@@ -123,10 +123,29 @@ class AriaEcSettings(Settings):
         if self.infra:
             self._up_infra()
 
+    @property
+    def ssdist(self):
+        if not self._ssdist:
+            if self.main.config["ss_dist_file"] and \
+                    os.path.exists(self.main.config["ss_dist_file"]):
+                self._ssdist = self.main.config["ss_dist_file"]
+            else:
+                self._ssdist = pkgr.resource_filename(__name__, self.SS_DIST)
+        return self._ssdist
+
+    @property
+    def template(self):
+        if not self._template:
+            if self.main.config["ariaproject_template"] and \
+                    os.path.exists(self.main.config["ariaproject_template"]):
+                self._template = self.main.config["ariaproject_template"]
+            else:
+                self._template = pkgr.resource_filename(__name__,
+                                                        self.ARIAPROJ_TEMPLATE)
+        return self._template
+
     @property
     def scsc_min(self):
-        # TODO: use cpickle only if we can't read scsc file in .ini conf !!!
-        #       Otherwise use another format ...
         # If scsc_min already computed
         if not self._scsc_min:
             try:

ariaec/protein.py

@@ -207,12 +207,14 @@ class SsList:
         """
         filename = os.path.abspath(ssdistpath) if ssdistpath else None
 
+        logger.info("Loading ss dist file")
         try:
             with open(filename) as f:
                 self._read_ssdist(f, filename=filename)
-        except IOError:
+        except Exception, message:
+            logger.error("%s" % message)
             logger.error("Can't load given ss dist file...")
-            logger.info("Loading default ss dist file")
+            logger.error("Loading default ss dist file")
             with pkgr.resource_stream(__name__, self.settings.SS_DIST) as f:
                 self._read_ssdist(f, filename=self.settings.SS_DIST)
 
@@ -457,8 +459,13 @@ class Protein:
         if ssdist_filename:
             # Read secondary distance matrix
             self.sec_struct.read_ssdist(ssdist_filename)
+        else:
+            logger.error("No secondary structure distance file found. Please "
+                         "check configuration file")
         if self.aa_sequence.sequence:
             # Synchronise sec structure sequence with aa sequence
+            logger.info("Align secondary structure sequence with protein "
+                        "sequence")
             self.sec_struct.seq_sublist(self.aa_sequence.sequence)
         if ssidx:
             logger.info("Using secondary structure index for amino acid "

ariaec/protmap.py

@@ -614,7 +614,7 @@ class ResAtmMap(ProteinMap):
             newmap[:] = self.copy().stack().groupby(level=0).min()
         return newmap
 
-    def contact_map(self, contactdef, scsc_min=None):
+    def contact_map(self, contactdef, scsc_min=None, def_cut=5.0):
         """
         Contact map generator from all atoms distance map. There's a contact
         with 2 residues iff dist between 2 atoms are below the given treshold
@@ -632,12 +632,13 @@ class ResAtmMap(ProteinMap):
         # Initialize contact map to a boolean matrix filled with False
         contact_map = ResAtmMap(sequence=self.sequence, mtype="contact",
                                 desc=self.desc, sym=self.sym, path=self.path)
+        def_cutoff = float(contactdef.get("default_cutoff")) if float(
+            contactdef.get("default_cutoff")) else def_cut
 
         if type(contactdef) == float:
             contact_map[:] = self.applymap(lambda x: x < contactdef)
 
-        elif sum(x is not None for x in contactdef.values()) == 1 and \
-                contactdef.get("default_cutoff"):
+        elif sum(x is not None for x in contactdef.values()) == 1 and def_cutoff:
             logger.info("Using default cutoff")
             contact_map[:] = self.applymap(lambda x: x < contactdef.get("default_cutoff"))
 
@@ -649,8 +650,12 @@ class ResAtmMap(ProteinMap):
             for pair in contactdef.keys():
                 if pair == "default_cutoff":
                     continue
-                treshold = contactdef[pair]
                 pair = tuple(pair.upper().split("_"))
+                if contactdef[pair] > def_cutoff:
+                    treshold = contactdef[pair]
+                else:
+                    logger.warning("Treshold for %s ignored (lower than "
+                                   "the default cutoff)" % str(pair))
                 logger.info(
                     "Filtering values in matrix related to %s (%s)" %
                     (str(pair), str(treshold)))

ariaec/setup.py

@@ -58,11 +58,9 @@ class AriaEcSetup:
         # ------------------------- Load sequence ---------------------------- #
         self.protein.set_aa_sequence(self.settings.setup.args.get("seq", None))
         # -------------- Load secondary structure prediction ----------------- #
-        ss_file = self.settings.main.config["ss_dist_file"] if \
-            self.settings.main.config["ss_dist_file"] else None
         self.protein.set_sec_struct(self.settings.setup.args.get("sspred",
                                                                  None),
-                                    ss_file,
+                                    ssdist_filename=self.settings.ssdist,
                                     ssidx=self.settings.setup.args.get(
                                         "ssidx", False))
         # -------------------------------------------------------------------- #
@@ -148,6 +146,7 @@ class AriaEcSetup:
         # Setting contact number limit for hbmap
         n_hb = int(len(self.protein.aa_sequence.sequence) *
                    self.settings.setup.config.get("nf_longrange_hb"))
+        logger.info("Writing tbl files ...")
         tbl_files = self.converter.write_tbl_restraints(
             self.protein, hbmap=self.hbmaps, n_hb=n_hb)
 
@@ -163,11 +162,7 @@ class AriaEcSetup:
         xmlseq_file = self.converter.write_xmlseq()
 
         # ---------------------- ARIA XML project file ----------------------- #
-        aria_template = self.settings.main.config["ariaproject_template"] if \
-            self.settings.main.config["ariaproject_template"] and \
-            os.path.exists(self.settings.main.config["ariaproject_template"])\
-            else None
-        self.converter.write_ariaproject(aria_template,
+        self.converter.write_ariaproject(self.settings.template,
                                          xmlseq_file, dist_files, tbl_files,
                                          desc="_".join(sorted(self.allresmap.keys())))
         # ------------------------------ others ------------------------------ #

version.py

 # Do not edit this file, pipeline versioning is governed by git tags
-__version__=v0.1.11-dev2
\ No newline at end of file
+__version__=v0.1.11-dev2-8-g4239f50
\ No newline at end of file