Advertisement
Guest User

structure_fitting.py

a guest
Sep 23rd, 2016
136
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 6.44 KB | None | 0 0
  1. #!/usr/bin/python
  2. """
  3. This script pulls in homologs of proteins from PDB
  4. as determined by HHSuite. It then employs UCSF Chimera to
  5. structurally match them and get an indication of how well
  6. they score (RMSD) in order to pick the best simulation.
  7. """
  8.  
  9. from __future__ import print_function
  10. import os
  11. import subprocess
  12. import sys
  13. import argparse
  14. import traceback
  15. import warnings
  16. import glob
  17. import fnmatch
  18.  
  19. import pychimera
  20.  
  21. import chimera
  22. from chimera import openModels,\
  23.                 Molecule,\
  24.                 replyobj
  25. from chimera import runCommand as rc
  26. from MatchMaker import\
  27.          match,\
  28.          CP_BEST,\
  29.          GAP_OPEN,\
  30.          GAP_EXTEND,\
  31.          defaults,\
  32.          MATRIX,\
  33.          ITER_CUTOFF
  34.  
  35.  
  36. # If running using python interpreter and not pychimera:
  37. # os.environ['CHIMERADIR'] = '/home/wms_joe/Applications/CHIMERA1.11'
  38. # CHIMERADIR should point to the application root directory.
  39. # This can be found with:   `chimera --root`
  40.  
  41.  
  42. # Only needed if running via normal python interpreter, not pychimera
  43. # pychimera.patch_environ()
  44. # pychimera.enable_chimera()
  45.  
  46.  
  47. ##### Main code begins #####
  48.  
  49.  
  50.  
  51. # Step 1, parse the output of HHpred to get the nearest homolog.
  52.  
  53. def main():
  54.     def_cpus=12
  55.     def_dir=os.getcwd()
  56.     try:
  57.         parser = argparse.ArgumentParser(description='This script compares protein structural homologs as determined with HHpred, to PDB models using UCSF CHIMERA, to gather metrics of structural similarity.')
  58.  
  59.         parser.add_argument(
  60.             '-d',\
  61.             '--database',\
  62.             action='store',
  63.             default=None,\
  64.             help='You can specify a different HHpred database (filepath) to use if you offer this parameter. Otherwise it defaults to PDB.')
  65.  
  66.         parser.add_argument(
  67.             '-f',\
  68.             '--fasta',\
  69.             action='store',\
  70.             help='The fasta amino acid sequence that corresponds to the simulated structure and the sequence you wish to query.')
  71.  
  72.         parser.add_argument(
  73.             '-s',\
  74.             '--simulations',\
  75.             action='store',\
  76.             help='The directory where the protein structure simulations are stored. They must be in PDB format and be named logically.')
  77.  
  78.         parser.add_argument(
  79.             '-r',\
  80.             '--rmsd',\
  81.             action='store',\
  82.             default=None,\
  83.             help='The filename to store the RMSD values for each matching.')
  84.  
  85.         parser.add_argument(
  86.             '-t',\
  87.             '--threads',\
  88.             action='store',\
  89.             default=def_cpus,\
  90.             help='The number of threads that HHsearch can execute on.')
  91.  
  92.         parser.add_argument(
  93.             '-o',\
  94.             '--outdir',\
  95.             action='store',\
  96.             default=def_dir,\
  97.             help="Directory for files to be written to. Default is the current working directory.")
  98.         args = parser.parse_args()
  99.  
  100.     except:
  101.         print("An exception occured with argument parsing. Check your provided options.")
  102.         traceback.print_exc()
  103.  
  104.   # Acquire HHPred Database
  105.     if args.database is None:
  106.         db_cmd = 'find ~/Applications/HHSuite/databases/pdb70 -type f -name "pdb70_hhm.ffdata"'
  107.             db_path = subprocess.Popen(
  108.                     db_cmd,\
  109.                     shell=True,\
  110.                     stdin=subprocess.PIPE,\
  111.                     stdout=subprocess.PIPE,\
  112.                     stderr=subprocess.PIPE)
  113.  
  114.         (stdout,stderr) = db_path.communicate()
  115.             filelist = stdout.decode().split()
  116.         args.database = stdout
  117.     else:
  118.         print("No default database found and you haven't provided one directly.")
  119.  
  120.     if args.fasta is not None:
  121.         split = os.path.splitext(args.fasta)
  122.         basename = os.path.basename(split[0])
  123.     else:
  124.         print('No input fasta was provided. Check your input parameters')
  125.         sys.exit(1)
  126.  
  127.     # Run the HHsearch
  128.     if args.fasta is not None and args.database is not None:
  129.         print("\n")
  130.         print("Running " + args.fasta + " in HHsearch, against " + args.database + " on " + args.threads + " threads.")
  131.         hhresult_file = '{0}.hhr.fasta'.format(basename)
  132.         search_cmd = 'hhsearch -dbstrlen 50 -B 1 -b 1 -p 60 -Z 1 -E 1E-03 -nocons -nopred -nodssp -Ofas {3} -cpu {0} -i {1} -d {2}'.format(args.threads,args.fasta, args.database, hhresult_file)
  133.         print("\n")
  134.         print("\n")
  135.         print("Executing HHsearch with the command:")
  136.         print(search_cmd)
  137.         # hh_process = subprocess.Popen(
  138.         #           search_cmd,\
  139.         #           shell=True,\
  140.         #           stdin=subprocess.PIPE,\
  141.         #           stderr=subprocess.PIPE)
  142.         # hh_process.wait()
  143.  
  144.     else:
  145.         print("No fasta or database has been detected. Check your input parameters.")
  146.         sys.exit(1)
  147.  
  148.     # Parse HHpred output fasta and acquire the best hit
  149.  
  150.     headers = []
  151.     with open(hhresult_file) as result_fasta:
  152.         for line in result_fasta:
  153.             if line.startswith(">"):
  154.                 line = line.split("_")[0]
  155.                 headers.append(line.replace(">",""))
  156.  
  157.     top_hit = headers[1] # HHpred puts the query seq in index 0, so 1 is your top hit.
  158.     print("Your best hit was the PDB id:")
  159.     print(top_hit + "\n")
  160.     print("Full results are found in: " + hhresult_file)
  161.     result_fasta.close()
  162.  
  163.     # Acquire the protein simulations
  164.  
  165.     model_list = []
  166.     for root, dirnames, filenames in os.walk(args.simulations):
  167.         for filename in sorted(fnmatch.filter(filenames, '*.pdb')):
  168.             if filename.startswith(basename):
  169.                 model_list.append(os.path.join(root, filename))
  170.  
  171.     print("\n")
  172.     print("Found the following models:")
  173.     print("---------------------------")
  174.     for model_path in model_list:
  175.         locus_dir = os.path.dirname(os.path.abspath(model_path))
  176.         print("Found: " + os.path.basename(model_path) + " in " + locus_dir)
  177.  
  178.  
  179.     # Get reference structure from PDB
  180.     print("\n")
  181.     print("Beginning Chimera:")
  182.     print("---------------------------")
  183.     chimera.openModels.open(top_hit,type="PDB")
  184.     print("Opened reference structure: " + top_hit)
  185.  
  186.     # Open model structures
  187.     for model_path in model_list:
  188.         chimera.openModels.open(model_path,type="PDB")
  189.         print("Opened: " + os.path.basename(model_path))
  190.  
  191.     if args.rmsd is not None:
  192.         rmsd_file = args.rmsd
  193.     else:
  194.         rmsd_file = '{0}{1}.tsv'.format(locus_dir,basename)
  195.  
  196.     print("RMSD values are:")
  197.     print("---------------------------")
  198.     with open(rmsd_file, "a") as rmsd_output_file:
  199.         all_models = chimera.openModels.list(modelTypes=[chimera.Molecule])
  200.         ref = all_models[0]
  201.         sims = all_models[1:]
  202.  
  203.         for atoms1, atoms2, rmsd, fullRmsd in match(CP_BEST,[ref, sims],defaults[MATRIX],\
  204.                                                     "nw",defaults[GAP_OPEN],defaults[GAP_EXTEND]):
  205.             ref_mol = atoms2[0].molecule
  206.             sim_mol = atoms1[0].molecule
  207.  
  208.             print(ref_mol.name + "\t" + sim_mol.name + "\t" + str(rmsd))
  209.             rmsd_output_file.write(ref_mol.name + "\t" + sim_mol.name + "\t" + str(rmsd))
  210.  
  211.     rmsd_output_file.close()
  212.  
  213.     rc('save {0}{1}_session.py'.format(locus_dir,basename))
  214.     rc('stop now')
  215.  
  216. if __name__ == '__main__' :
  217.     main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement