Merge branch 'master' of github.com:hng/BiomolecularStructures

docs/index.md

@@ -18,5 +18,6 @@ The BiomolecularStructures package provides several Bioinformatics-related modul
 * Kabsch - Superimposing protein structures
 * PDB - Utility functions for parsing PDB files
 * Plot - Rudimentary plotting of matrices of atomic coordinates
+* Mafft - Julia API for multisequence alignment with MAFFT
 * Modeller - 
-* Cluster -
+* Cluster -

docs/modeller.md

@@ -0,0 +1,11 @@
+# Modeller
+
+## Exported functions
+
+``
+    function gen_script(name::String)``
+
+Generates Julia scripts for MODELLER.
+name: The name of the script (minus the extension), e.g. "build_profile"
+These scripts are based on the basic example scripts from the MODELLER website.
+Scripts are generated in the current working directory. You can find all scripts that can be generated in src/MODELLER/modeller-basic-example-julia .

mkdocs.yml

@@ -8,6 +8,7 @@ pages:
 - ['webblast.md', 'WebBLAST - NCBI/EBI BLAST Web-API']
 - ['kabsch.md', 'Kabsch - Kabsch Algorithm']
 - ['mafft.md', 'MAFFT - Julia API for MAFFT commandline tool for sequence alignment']
+- ['modeller.md', 'Modeller - Functions and scripts to use MODELLER with Julia']
 - ['pdb.md', 'PDB - PDB Utilities']
 - ['plot.md', 'MatrixPlot - Plotting Utilities']
 

src/MAFFT/mafft.jl

@@ -2,7 +2,7 @@
     TODO
     - group to group alignments
 """
-module MAFFT
+module Mafft
 export mafft, mafft_from_string, mafft_from_fasta, mafft_linsi, linsi, mafft_ginsi, ginsi, mafft_einsi, einsi, mafft_fftnsi, fftnsi, mafft_fftns, fftns, mafft_nwnsi, nwnsi, mafft_nwns, nwns, print_aligned_fasta, alignment_length, to_aminoacids
 
     using FastaIO

src/MODELLER/modeller.jl

@@ -2,11 +2,24 @@
   adapted from the modeller example scripts
 "
 
-module MODELLER
-export model_single, gen_script
+module Modeller
+export model_single, gen_script, evaluate_model
     using PyCall
 
-    function model_single(alnf, know, seq)
+
+    # Generator for MODELLER julia scripts.
+    # These scripts are based on the basic example scripts from the MODELLER website.
+    # Scripts are generated in the current working directory
+    # name: The name of the script (minus the extension), e.g. "build_profile"
+    function gen_script(name::String)
+       file = Pkg.dir("BiomolecularStructures", "src/MODELLER/modeller-basic-example-julia", "$name.jl")
+       if isfile(file)
+          cp(file, "./$name.jl")
+          println("Generated $name.jl MODELLER script") 
+       end 
+    end
+
+    function model_single(alnf::String, know, seq)
 
         @pyimport modeller
         @pyimport modeller.automodel as am
@@ -22,15 +35,30 @@ export model_single, gen_script
         a[:make]()
     end
 
-    # Generator for MODELLER julia scripts.
-    # These scripts are based on the basic example scripts from the MODELLER website.
-    # Scripts are generated in the current working directory
-    # name: The name of the script (minus the extension), e.g. "build_profile"
-    function gen_script(name::String)
-       file = Pkg.dir("BiomolecularStructures", "src/MODELLER/modeller-basic-example-julia", "$name.jl")
-       if isfile(file)
-          cp(file, "./$name.jl")
-          println("Generated $name.jl MODELLER script") 
-       end 
+    function evaluate_model(pdbfile::String; outputfile::String)
+
+        # check for optional argument output and set it to the first part of the name of the PDB file
+        if !isdefined(:outputfile)
+            outputfile = string(split(a,".")[1], ".profile")
+        end
+
+        @pyimport modeller
+        @pyimport modeller.scripts as scripts
+        @pyimport _modeller
+
+        mod_lib = _modeller.mod_libdir_get()
+
+        modeller.log[:verbose]()    # request verbose output
+        env = modeller.environ()
+        env[:libs][:topology][:read](file=string(mod_lib, "/top_heav.lib")) # read topology
+        env[:libs][:parameters][:read](file=string(mod_lib, "/par.lib")) # read parameters
+
+        # read model file
+        mdl = scripts.complete_pdb(env, pdbfile)
+
+        # Assess with DOPE:
+        s = modeller.selection(mdl)   # all atom selection
+        s[:assess_dope](output="ENERGY_PROFILE NO_REPORT", file=outputfile,
+                  normalize_profile=true, smoothing_window=15)
     end
 end