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

hng · Mar 26, 2015 · 21f6249 · 21f6249
2 parents e6f62d0 + 439199b
commit 21f6249
Show file tree

Hide file tree

Showing 2 changed files with 78 additions and 3 deletions.
diff --git a/docs/modeller.md b/docs/modeller.md
@@ -3,9 +3,23 @@
 ## Exported functions
 
 ``
-    function gen_script(name::String)``
+    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 .
+
+``
+build_profile(;seq_database_file::String = "", seq_database_format::String="PIR", alignment_file::String = "", alignment_format::String = "PIR", output_name::String = "build_profile", output_profile_format::String="TEXT", output_alignment_format::String="PIR")
+``
+
+**Example:**
+
+`` build_profile(seq_database_file="pdb_95.pir", alignment_file="TvLDH.ali") ``
+
+`` align2d(model_file::String, model_segment, model_align_codes::String, atom_files::String, align_file::String, align_codes::String, outputname::String)``
+
+**Example:**
+
+``align2d("1bdm.pdb", ("FIRST:A","LAST:A"), "1bdmA", "1bdm.pdb", "TvLDH.ali", "TvLDH", "TvLDH-1bdmA")``
diff --git a/src/MODELLER/modeller.jl b/src/MODELLER/modeller.jl
@@ -3,10 +3,9 @@
 "
 
 module Modeller
-export model_single, gen_script, evaluate_model
+export  gen_script, build_profile, model_single, evaluate_model, align2d
     using PyCall
 
-
     # 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
@@ -19,6 +18,55 @@ export model_single, gen_script, evaluate_model
        end 
     end
 
+    function build_profile(;seq_database_file::String = "", seq_database_format::String="PIR", alignment_file::String = "", alignment_format::String = "PIR", output_name::String = "build_profile", output_profile_format::String="TEXT", output_alignment_format::String="PIR")
+        seq_database_file_out = string(split(seq_database_file, "."), ".bin")
+
+        pyinitialize()
+        @pyimport modeller
+        @pyimport _modeller
+        mod_lib = _modeller.mod_libdir_get()
+        modeller.log[:verbose]()
+        env = modeller.environ()
+
+        #-- Prepare the input files
+
+        #-- Read in the sequence database
+        sdb = modeller.sequence_db(env)
+        sdb[:read](seq_database_file=seq_database_file, seq_database_format=seq_database_format,
+                 chains_list="ALL", minmax_db_seq_len=(30, 4000), clean_sequences=true)
+
+        #-- Write the sequence database in binary form
+        sdb[:write](seq_database_file=seq_database_file_out, seq_database_format="BINARY",
+                  chains_list="ALL")
+
+        #-- Now, read in the binary database
+        sdb[:read](seq_database_file=seq_database_file_out, seq_database_format="BINARY",
+                 chains_list="ALL")
+
+        #-- Read in the target sequence/alignment
+        aln = modeller.alignment(env)
+        aln[:append](file=alignment_file, alignment_format=alignment_format, align_codes="ALL")
+
+        #-- Convert the input sequence/alignment into
+        #   profile format
+        prf = aln[:to_profile]()
+
+        #-- Scan sequence database to pick up homologous sequences
+        prf[:build](sdb, matrix_offset=-450, rr_file=string(mod_lib,"/blosum62.sim.mat"),
+                  gap_penalties_1d=(-500, -50), n_prof_iterations=1,
+                  check_profile=false, max_aln_evalue=0.01)
+
+        #-- Write out the profile in text format
+        prf[:write](file=string(output_name, ".prf"), profile_format=output_profile_format)
+
+        #-- Convert the profile back to alignment format
+        aln = prf[:to_alignment]()
+
+        #-- Write out the alignment file
+        aln[:write](file=string(output_name, ".ali"), alignment_format=output_alignment_format)
+
+    end
+
     function model_single(alnf::String, know, seq)
 
         @pyimport modeller
@@ -61,4 +109,17 @@ export model_single, gen_script, evaluate_model
         s[:assess_dope](output="ENERGY_PROFILE NO_REPORT", file=outputfile,
                   normalize_profile=true, smoothing_window=15)
     end
+
+   function align2d(model_file::String, model_segment, model_align_codes::String, atom_files::String, align_file::String, align_codes::String, outputname::String)
+       @pyimport modeller
+
+       env = modeller.environ()
+       aln = modeller.alignment(env)
+       mdl = modeller.model(env, file=model_file, model_segment=model_segment)
+       aln[:append_model](mdl, align_codes=model_align_codes, atom_files=atom_files)
+       aln[:append](file=align_file, align_codes=align_codes)
+       aln[:align2d]()
+       aln[:write](file=string(outputname, ".ali"), alignment_format="PIR")
+       aln[:write](file=string(outputname, ".pap"), alignment_format="PAP")
+    end
 end