Merge branch 'develop-1.x'

We do not ask for the Annotation script but use MC-Annotate
as default on master, to keep full backwards-compatibility.

.travis.yml

@@ -18,7 +18,7 @@ before_install:
   - "sudo apt-get install -y graphviz"
 
 install:
-  - "pip install --upgrade pip wheel" 
+  - "pip install --upgrade pip wheel"
   - "pip install --only-binary=numpy,scipy numpy>1.10.0 scipy"
   - "pip install -r requirements.txt"
   - "mkdir mctools"
@@ -28,6 +28,8 @@ install:
   - "unzip MC-Annotate.zip"
   - "mv MC-Annotate mctools/"
   - "export PATH=$PATH:$PWD/mctools/"
+  - "mkdir ~/.config/forgi/"
+  - "echo '{\"PDB_ANNOTATION_TOOL\": \"MC-Annotate\"}' >~/.config/forgi/config.json"
 
-script: nose2 test
 
+script: nose2 test

doc/download.rst

@@ -46,9 +46,20 @@ These additional dependencies include:
 * Biopython (Needed for `forgi.threedee`)
 * pylab (Needed for `forgi.visual.mplotlib`, )
 * beautifulsoap4 (Needed for fetching informations about modified residues in pdb files)
+* appdirs (Needed for `forgi.threedee`)
 
-MC-Annotate
-~~~~~~~~~~~
+Either MC-Annotate or DSSR
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-To convert PDB files into CoarseGrainRNA objects, MC-Annotate is required.
-It is available at http://major.iric.ca/MajorLabEn/MC-Tools.html
+To convert PDB files into CoarseGrainRNA objects, either MC-Annotate or DSSR is required.
+
+MC-Annotate is available at http://major.iric.ca/MajorLabEn/MC-Tools.html
+
+DSSR is available at http://home.x3dna.org/
+
+You can create a configuration file to decide which of the two programs you want to use.
+If no configuration file exists, you will be prompted to choose the program the first time you
+convert a PDB file. Your choice will be recorded in a newly generated configuration file.
+
+We use the appdirs package to define the system specific location of the
+configuration file. Under Linux it would default to ~/.config/forgi.

examples/describe_cg.py

@@ -134,9 +134,15 @@ def describe_ml_segments(cg):
             j3_familyFlat = None
             j3_familyPerp = None
             j3_Delta = None
+        loop_start = float("inf")
+        for segment in loop:
+            if cg.define_a(segment)[0]<loop_start:
+                loop_start = cg.define_a(segment)[0]
         for segment in loop:
             if segment[0] !="m":
                 continue
+            data["loop_start_after"].append(loop_start)
+            data["segment_start_after"].append(cg.define_a(segment)[0])
             data["segment"].append(segment)
             data["junction_length"].append(len(loop))
             data["segment_length"].append(cg.get_length(segment))
@@ -172,9 +178,27 @@ def describe_ml_segments(cg):
             data["j3_family3D"].append(j3_family3D)
             data["j3_familyPerp"].append(j3_familyPerp)
             data["j3_Delta_j23_j31"].append(j3_Delta)
-
+    if data:
+        data["pk_number"]=number_by(data, "loop_start_after", "is_pseudoknotted_multiloop")
+        data["loop_number"]=number_by(data, "loop_start_after", None)
+        data["reguler_multiloop_number"]=number_by(data, "loop_start_after", "is_regular_multiloop")
     return data
 
+def number_by(data, sorting_column = "loop_start_after", only_for_col = "is_pseudoknotted_multiloop"):
+    log.debug(list((key, len(data[key])) for key in data))
+    df = pd.DataFrame(data)
+    if only_for_col is not None:
+        df = df[df[only_for_col]==True]
+    sorted_vals = list(sorted(set(df[sorting_column])))
+    out_column = []
+    for i in range(len(data[sorting_column])):
+        if only_for_col is None or data[only_for_col][i]:
+            out_column.append(sorted_vals.index(data[sorting_column][i])+1)
+        else:
+            out_column.append(0)
+    log.info("number_by column is %s, len(data[%s])=%s)", out_column, sorting_column, len(data[sorting_column]))
+    return out_column
+
 parser = generateParser()
 if __name__=="__main__":
     args = parser.parse_args()

examples/input/1y26.cg

@@ -0,0 +1,25 @@
+name 1y26
+length 71
+seq CGCUUCAUAUAAUCCUAAUGAUAUGGUUUGGGAGUUUCUACCAAGAGCCUUAAACUCUUGAUUAUGAAGUG
+seq_ids 
+define s0 1 9 63 71
+define m1 10 12
+define s2 13 19 27 33
+define h0 20 26
+define m2 34 41
+define s3 42 47 55 60
+define h1 48 54
+define m0 61 62
+connect s0 m1 m0
+connect s2 h0 m1 m2
+connect s3 h1 m0 m2
+coord h1 40.3992211551000011 -21.0722649836000002 -8.7858799860799994 56.2509994507000002 -31.7719993591000005 -10.3870000839000003
+coord s3 31.5987607900999983 -17.2471501278999995 -0.4056888140200000 40.3992211551000011 -21.0722649836000002 -8.7858799860799994
+coord s2 24.2469227276999995 -28.5048147442000008 -0.4624864930680000 37.9464206031000018 -34.4588356333999997 -3.7175804801000001
+coord h0 37.9464206031000018 -34.4588356333999997 -3.7175804801000001 50.0929985045999970 -23.0610008240000006 -16.1790008544999999
+coord s0 2.4105633276299998 0.1812946855950000 3.8104106076800002 19.2708715511000008 -11.2603784822000001 4.3801360791599997
+twist s3 -0.5945111509714611 0.2819655992687657 -0.7530284803375131 0.7082218070660000 0.0975338596706000 0.6992202930510000
+twist s2 0.3326392289800000 0.2716426164470000 0.9030843993080000 -0.2406093171120000 -0.0231200541610000 -0.9703466492000000
+twist s0 0.0845426612698000 0.0751065442191000 -0.9935851978769999 -0.5560404110783073 -0.8101227737674702 0.1857959974570234
+project 1.0 1.0 1.0
+

examples/visualize_cg.py

@@ -7,6 +7,7 @@
 import subprocess as sp
 import tempfile as tf
 import time
+import os.path
 
 import forgi.threedee.model.coarse_grain as cmg
 import forgi.utilities.debug as fud
@@ -15,7 +16,7 @@
 import forgi.threedee.model.similarity as ftur
 import forgi.threedee.utilities.vector as ftuv
 import forgi.threedee.visual.pymol as cvp
-
+from forgi.utilities.stuff import make_temp_directory
 from optparse import OptionParser
 
 
@@ -25,7 +26,7 @@ def align_cgs(cgs):
 
     The points representing each coarse grain RNA molecule
     will be the virtual residues.
-    
+
     @param cgs: A list of CoarseGrainRNA structures.
     @return: Nothing, the cgs are modified in place
     '''
@@ -96,11 +97,9 @@ def main():
 
     (options, args) = parser.parse_args()
 
-    print("hi")
     if len(args) < num_args:
         parser.print_help()
         sys.exit(1)
-    print("hi1")
 
     pp = cvp.PymolPrinter()
     pp.stem_color = options.stem_color
@@ -193,38 +192,32 @@ def main():
 
         pp.add_sphere(mid1, 'green', width=2)
         pp.add_sphere(mid2, 'red', width=2)
-
-
-    with tf.NamedTemporaryFile() as f:
-        with tf.NamedTemporaryFile(suffix='.pml') as f1:
-            f.write(pp.pymol_string())
-            f.flush()
-
-            pymol_cmd = 'hide all\n'
-            pymol_cmd += 'run %s\n' % (f.name)
-            pymol_cmd += 'show cartoon, all\n'
-            pymol_cmd += 'bg white\n'
-            pymol_cmd += 'clip slab, 10000\n'
-            pymol_cmd += 'orient\n'
 
-            if options.output is not None:
-                pymol_cmd += 'ray\n'
-                pymol_cmd += 'png %s\n' % (options.output)
-                pymol_cmd += 'quit\n'
 
+    with make_temp_directory() as tmpdir:
+        # The file describing the cg-structure as cylinders
+        stru_filename = os.path.join(tmpdir, "structure")
+        with open(stru_filename, "w") as f:
+            f.write(pp.pymol_string())
+        # The file for running pymol
+        pymol_cmd = 'hide all\n'
+        pymol_cmd += 'run %s\n' % (stru_filename)
+        pymol_cmd += 'show cartoon, all\n'
+        pymol_cmd += 'bg white\n'
+        pymol_cmd += 'clip slab, 10000\n'
+        pymol_cmd += 'orient\n'
+        if options.output is not None:
+            pymol_cmd += 'ray\n'
+            pymol_cmd += 'png %s\n' % (options.output)
+            pymol_cmd += 'quit\n'
+        pml_filename = os.path.join(tmpdir, "command.pml")
+        with open(pml_filename, "w") as f1:
             f1.write(pymol_cmd)
-            f1.flush()
-
-            print("f1.name:", f1.name)
-
-            if options.batch:
-                p = sp.Popen(['pymol', '-cq', f1.name], stdout=sp.PIPE, stderr=sp.PIPE)
-            else:
-                p = sp.Popen(['pymol', f1.name], stdout=sp.PIPE, stderr=sp.PIPE)
-
-            out, err = p.communicate()
-            print("err:", err, file=sys.stderr)
+        if options.batch:
+            p = sp.Popen(['pymol', '-cq', pml_filename], stdout=sp.PIPE, stderr=sp.PIPE)
+        else:
+            p = sp.Popen(['pymol', pml_filename], stdout=sp.PIPE, stderr=sp.PIPE)
+        out, err = p.communicate()
 
 if __name__ == '__main__':
     main()
-

forgi/config.py

@@ -1,20 +1,41 @@
-from builtins import object
+import json
 import os
-import forgi
+import os.path
+import logging
+import appdirs
 
-class Configuration(object):
-    mids_method="template"
-    #mids_method="basenormals"
-    base_dir = os.path.expanduser('.')
-    #data_base_dir = os.path.expanduser('~/data/ernwin/processed')
-    #pdb_base_dir = os.path.expanduser('~/data/ernwin/pdb')
-    stem_fragment_dir = os.path.join(base_dir, 'forgi/data')
-    #lric_stats_fn = os.path.join(base_dir, 'fess/stats/temp.energy')
-    #template_residue_fn = os.path.join(base_dir, 'fess/stats/residue_template.pdb')
-    #longrange_contact_stats_fn = os.path.join(base_dir, 'fess/stats/temp.longrange.contact')
+log = logging.getLogger(__name__)
 
-    #test_input_dir = os.path.expanduser('~/data/ernwin/processed/')
-    #test_output_dir = os.path.join(base_dir, "test_output")
-    #sampling_output_dir = 'best'
-    #barnacle_dir = '/scr/plastilin/pkerp/apps/Barnacle'
-    #stem_library = dict()
+dirs = appdirs.AppDirs("forgi", "TBI")
+
+def read_config():
+    config = {}
+    for directory in [dirs.site_config_dir, dirs.user_config_dir]:
+        filename = os.path.join(directory, "config.json")
+        try:
+            with open(filename) as f:
+                conf = json.load(f)
+        except (OSError, IOError):
+            log.debug("No configuration file present at %s", filename)
+            pass
+        else:
+            log.debug("Reading configuration from %s", filename)
+            config.update(conf)
+
+    return config
+
+def set_config(key, value):
+    filename = os.path.join(dirs.user_config_dir, "config.json")
+    try:
+        with open(filename) as f:
+            config = json.load(f)
+    except (OSError, IOError):
+        config = {}
+    config[key]=value
+    try:
+        os.makedirs(dirs.user_config_dir)
+    except OSError:
+        pass
+    with open(filename, "w") as f:
+        json.dump(config, f)
+        log.info("Configuration file %s updated", filename)

forgi/graph/bulge_graph.py

@@ -35,36 +35,15 @@
 
 from logging_exceptions import log_to_exception, log_at_caller
 
+from .residue import RESID, resid_to_str, resid_from_str
+
 try:
   profile  #The @profile decorator from line_profiler (kernprof)
 except:
   def profile(x):
     return x
 
 
-RESID = col.namedtuple("complete_resid", ["chain", "resid"])
-
-def resid_to_str(resid):
-    if resid.chain is not None:
-        out="{}:{}".format(resid.chain, resid.resid[1])
-    else:
-        out=str(resid.resid[1])
-    if resid.resid[2]!=" ":
-        out+=".{}".format(resid.resid[2])
-    return out
-
-def resid_from_str(resstr):
-    if ":" in resstr:
-        chain, resid = resstr.split(":")
-    else:
-        resid=resstr
-        log.debug("No chain given in string {!r}".format(resstr))
-        chain='A'
-    idparts=resid.split(".")
-    if len(idparts)==1:
-        idparts.append(" ")
-    return RESID(chain, (' ', int(idparts[0]), idparts[1]))
-
 
 def add_bulge(bulges, bulge, context, message):
     """
@@ -1562,7 +1541,6 @@ def connection_type(self, define, connections):
                   =   ======================================================================
 
         """
-
         if define[0] == 'i':
             # interior loop, we just have to check if
             # connections[0] < connections[1]
@@ -3211,15 +3189,6 @@ def _insert_cutpoints_into_seq(self):
                 log.info("seq now has {} cutpoints".format(self.seq.count('&')))
 
 
-    def seqids_from_residue_map(self, residue_map):
-        """
-        Create the list of seq_ids from the list of MC-Annotate identifiers in the
-        residue map.
-        """
-        self.seq_ids = []
-        for i, r in enumerate(residue_map):
-            (from_chain, from_base) = ftum.parse_chain_base(r)
-            self.seq_ids += [RESID(from_chain, ftum.parse_resid(from_base))]
 
     # This function seems to be dead code, but might be useful in the future.
     # Consider adding this to whitelist.py

forgi/graph/residue.py

@@ -0,0 +1,29 @@
+from builtins import str
+import collections
+import logging
+
+log = logging.getLogger(__name__)
+
+RESID = collections.namedtuple("complete_resid", ["chain", "resid"])
+
+def resid_to_str(resid):
+    if resid.chain is not None:
+        out="{}:{}".format(resid.chain, resid.resid[1])
+    else:
+        out=str(resid.resid[1])
+    if resid.resid[2]!=" ":
+        out+=".{}".format(resid.resid[2])
+    return out
+
+def resid_from_str(resstr):
+    resstr = str(resstr) # Make sure we use future's newstring on python2
+    if ":" in resstr:
+        chain, resid = resstr.split(":")
+    else:
+        resid=resstr
+        log.debug("No chain given in string {!r}".format(resstr))
+        chain=str('A')
+    idparts=resid.split(".")
+    if len(idparts)==1:
+        idparts.append(" ")
+    return RESID(chain, (' ', int(idparts[0]), idparts[1]))