files/DockQ.py

   1 #!/usr/bin/env python
   2
   3 import Bio.PDB
   4 import warnings
   5 from Bio import BiopythonWarning
   6 warnings.simplefilter('ignore', BiopythonWarning)
   7 import sys
   8 import os
   9 import re
  10 import tempfile
  11 import numpy as np
  12 from Bio.SVDSuperimposer import SVDSuperimposer
  13 from math import sqrt
  14 from argparse import ArgumentParser
  15 import itertools
  16 import subprocess
  17
  18 def parse_fnat(fnat_out):
  19     fnat=-1;
  20     nat_correct=-1
  21     nat_total=-1
  22     fnonnat=-1
  23     nonnat_count=-1
  24     model_total=-1
  25     inter=[]
  26     for line in fnat_out.split("\n"):
  27 #        print line
  28         line=line.rstrip('\n')
  29         match=re.search(r'NATIVE: (\d+)(\w) (\d+)(\w)',line)
  30         if(re.search(r'^Fnat',line)):
  31             list=line.split(' ')
  32             fnat=float(list[3])
  33             nat_correct=int(list[1])
  34             nat_total=int(list[2])
  35         elif(re.search(r'^Fnonnat',line)):
  36             list=line.split(' ')
  37             fnonnat=float(list[3])
  38             nonnat_count=int(list[1])
  39             model_total=int(list[2])
  40         elif(match):
  41             #print line
  42             res1=match.group(1)
  43             chain1=match.group(2)
  44             res2=match.group(3)
  45             chain2=match.group(4)
  46            # print res1 + ' ' + chain1 + ' ' + res2 + ' ' + chain2
  47             inter.append(res1 + chain1)
  48             inter.append(res2 + chain2)
  49     return (fnat,nat_correct,nat_total,fnonnat,nonnat_count,model_total,inter)
  50
  51 def capri_class(fnat,iRMS,LRMS,capri_peptide=False):
  52
  53
  54     if capri_peptide:
  55
  56         if(fnat < 0.2 or (LRMS > 5.0 and iRMS > 2.0)):
  57             return 'Incorrect'
  58         elif((fnat >= 0.2 and fnat < 0.5) and (LRMS <= 5.0 or iRMS <= 2.0) or (fnat >= 0.5 and LRMS > 2.0 and iRMS > 1.0)):
  59             return 'Acceptable'
  60         elif((fnat >= 0.5 and fnat < 0.8) and (LRMS <= 2.0 or iRMS <= 1.0) or (fnat >= 0.8 and LRMS > 1.0 and iRMS > 0.5)):
  61             return 'Medium'
  62         elif(fnat >= 0.8 and (LRMS <= 1.0 or iRMS <= 0.5)):
  63             return 'High'
  64         else:
  65             return 'Undef'
  66     else:
  67
  68         if(fnat < 0.1 or (LRMS > 10.0 and iRMS > 4.0)):
  69             return 'Incorrect'
  70         elif((fnat >= 0.1 and fnat < 0.3) and (LRMS <= 10.0 or iRMS <= 4.0) or (fnat >= 0.3 and LRMS > 5.0 and iRMS > 2.0)):
  71             return 'Acceptable'
  72         elif((fnat >= 0.3 and fnat < 0.5) and (LRMS <= 5.0 or iRMS <= 2.0) or (fnat >= 0.5 and LRMS > 1.0 and iRMS > 1.0)):
  73             return 'Medium'
  74         elif(fnat >= 0.5 and (LRMS <= 1.0 or iRMS <= 1.0)):
  75             return 'High'
  76         else:
  77             return 'Undef'
  78
  79
  80 def capri_class_DockQ(DockQ,capri_peptide=False):
  81
  82     if capri_peptide:
  83         return 'Undef for capri_peptides'
  84
  85     (c1,c2,c3)=(0.23,0.49,0.80)
  86     if(DockQ < c1):
  87         return 'Incorrect'
  88     elif(DockQ >= c1 and DockQ < c2):
  89         return 'Acceptable'
  90     elif(DockQ >= c2 and DockQ < c3):
  91         return 'Medium'
  92     elif(DockQ >= c3):
  93         return 'High'
  94     else:
  95         return 'Undef'
  96
  97
  98 def calc_DockQ(model,native,use_CA_only=False,capri_peptide=False):
  99
 100     exec_path=os.path.dirname(os.path.abspath(sys.argv[0]))
 101     atom_for_sup=['CA','C','N','O']
 102     if(use_CA_only):
 103         atom_for_sup=['CA']
 104
 105     cmd_fnat=exec_path + '/fnat ' + model + ' ' + native + ' 5 -all'
 106     cmd_interface=exec_path + '/fnat ' + model + ' ' + native + ' 10 -all'
 107
 108     if capri_peptide:
 109         cmd_fnat=exec_path + '/fnat ' + model + ' ' + native + ' 4 -all'
 110         cmd_interface=exec_path + '/fnat ' + model + ' ' + native + ' 8 -cb'
 111
 112     fnat_out = os.popen(cmd_fnat).read()
 113
 114     #fnat_out = subprocess.getoutput(cmd_fnat)
 115     #print(fnat_out)
 116     #    sys.exit()
 117     (fnat,nat_correct,nat_total,fnonnat,nonnat_count,model_total,interface5A)=parse_fnat(fnat_out)
 118     assert fnat!=-1, "Error running cmd: %s\n" % (cmd_fnat)
 119     inter_out = os.popen(cmd_interface).read()
 120 #   inter_out = subprocess.getoutput(cmd_interface)
 121
 122     (fnat_bb,nat_correct_bb,nat_total_bb,fnonnat_bb,nonnat_count_bb,model_total_bb,interface)=parse_fnat(inter_out)
 123     assert fnat_bb!=-1, "Error running cmd: %s\n" % (cmd_interface)
 124
 125     #print fnat
 126     #Use same interface as for fnat for iRMS
 127     #interface=interface5A
 128
 129
 130     # Start the parser
 131     pdb_parser = Bio.PDB.PDBParser(QUIET = True)
 132
 133     # Get the structures
 134     ref_structure = pdb_parser.get_structure("reference", native)
 135     sample_structure = pdb_parser.get_structure("model", model)
 136
 137     # Use the first model in the pdb-files for alignment
 138     # Change the number 0 if you want to align to another structure
 139     ref_model    = ref_structure[0]
 140     sample_model = sample_structure[0]
 141
 142     # Make a list of the atoms (in the structures) you wish to align.
 143     # In this case we use CA atoms whose index is in the specified range
 144     ref_atoms = []
 145     sample_atoms = []
 146
 147     common_interface=[]
 148
 149     chain_res={}
 150
 151
 152     #find atoms common in both sample and native
 153     atoms_def_sample=[]
 154     atoms_def_in_both=[]
 155     #first read in sample
 156     for sample_chain in sample_model:
 157 #        print sample_chain
 158         chain=sample_chain.id
 159 #        print chain
 160         for sample_res in sample_chain:
 161            # print sample_res
 162             if sample_res.get_id()[0] != ' ': #Skip hetatm.
 163                 continue
 164             resname=sample_res.get_id()[1]
 165             key=str(resname) + chain
 166             for a in atom_for_sup:
 167                 atom_key=key + '.' + a
 168                 if a in sample_res:
 169                     if atom_key in atoms_def_sample:
 170                         print(atom_key + ' already added (MODEL)!!!')
 171                     atoms_def_sample.append(atom_key)
 172
 173     #then read in native also present in sample
 174     for ref_chain in ref_model:
 175         chain=ref_chain.id
 176         for ref_res in ref_chain:
 177             #print ref_res
 178             if ref_res.get_id()[0] != ' ': #Skip hetatm.
 179 #                print ref_res.get_id()
 180                 continue
 181             resname=ref_res.get_id()[1]
 182             key=str(resname) + chain
 183             for a in atom_for_sup:
 184                 atom_key=key + '.' + a
 185                 if a in ref_res and atom_key in atoms_def_sample:
 186                     if atom_key in atoms_def_in_both:
 187                         print(atom_key + ' already added (Native)!!!')
 188                     atoms_def_in_both.append(atom_key)
 189
 190
 191 #    print atoms_def_in_both
 192     for sample_chain in sample_model:
 193         chain=sample_chain.id
 194         if chain not in list(chain_res.keys()):
 195             chain_res[chain]=[]
 196         for sample_res in sample_chain:
 197             if sample_res.get_id()[0] != ' ': #Skip hetatm.
 198                 continue
 199             resname=sample_res.get_id()[1]
 200             key=str(resname) + chain
 201             chain_res[chain].append(key)
 202             if key in interface:
 203                 for a in atom_for_sup:
 204                     atom_key=key + '.' + a
 205                     if a in sample_res and atom_key in atoms_def_in_both:
 206                         sample_atoms.append(sample_res[a])
 207                 common_interface.append(key)
 208
 209     #print inter_pairs
 210
 211     chain_ref={}
 212     common_residues=[]
 213
 214
 215
 216     # Iterate of all chains in the model in order to find all residues
 217     for ref_chain in ref_model:
 218         # Iterate of all residues in each model in order to find proper atoms
 219         #  print dir(ref_chain)
 220         chain=ref_chain.id
 221         if chain not in list(chain_ref.keys()):
 222             chain_ref[chain]=[]
 223         for ref_res in ref_chain:
 224             if ref_res.get_id()[0] != ' ': #Skip hetatm.
 225                 continue
 226             resname=ref_res.get_id()[1]
 227             key=str(resname) + chain
 228
 229             #print ref_res
 230             #      print key
 231             # print chain_res.values()
 232             if key in chain_res[chain]: # if key is present in sample
 233                 #print key
 234                 for a in atom_for_sup:
 235                     atom_key=key + '.' + a
 236                     if a in ref_res and atom_key in atoms_def_in_both:
 237                         chain_ref[chain].append(ref_res[a])
 238                         common_residues.append(key)
 239                       #chain_sample.append((ref_res['CA'])
 240             if key in common_interface:
 241               # Check if residue number ( .get_id() ) is in the list
 242               # Append CA atom to list
 243                 #print key
 244                 for a in atom_for_sup:
 245                     atom_key=key + '.' + a
 246                     #print atom_key
 247                     if a in ref_res and atom_key in atoms_def_in_both:
 248                         ref_atoms.append(ref_res[a])
 249
 250
 251
 252     #get the ones that are present in native
 253     chain_sample={}
 254     for sample_chain in sample_model:
 255         chain=sample_chain.id
 256         if chain not in list(chain_sample.keys()):
 257             chain_sample[chain]=[]
 258         for sample_res in sample_chain:
 259             if sample_res.get_id()[0] != ' ': #Skip hetatm.
 260                 continue
 261             resname=sample_res.get_id()[1]
 262             key=str(resname) + chain
 263             if key in common_residues:
 264                 for a in atom_for_sup:
 265                     atom_key=key + '.' + a
 266                     if a in sample_res and atom_key in atoms_def_in_both:
 267                         chain_sample[chain].append(sample_res[a])
 268
 269         #if key in common_residues:
 270         #     print key
 271         #sample_atoms.append(sample_res['CA'])
 272         #common_interface.append(key)
 273
 274
 275     assert len(ref_atoms)!=0, "length of native is zero"
 276     assert len(sample_atoms)!=0, "length of model is zero"
 277     assert len(ref_atoms)==len(sample_atoms), "Different number of atoms in native and model %d %d\n" % (len(ref_atoms),len(sample_atoms))
 278
 279     super_imposer = Bio.PDB.Superimposer()
 280     super_imposer.set_atoms(ref_atoms, sample_atoms)
 281     super_imposer.apply(sample_model.get_atoms())
 282
 283     # Print RMSD:
 284     irms=super_imposer.rms
 285
 286     (chain1,chain2)=list(chain_sample.keys())
 287
 288     ligand_chain=chain1
 289     receptor_chain=chain2
 290     len1=len(chain_res[chain1])
 291     len2=len(chain_res[chain2])
 292
 293     assert len1!=0, "%s chain has zero length!\n" % chain1
 294     assert len2!=0, "%s chain has zero length!\n" % chain2
 295
 296     class1='ligand'
 297     class2='receptor'
 298     if(len(chain_sample[chain1]) > len(chain_sample[chain2])):
 299         receptor_chain=chain1
 300         ligand_chain=chain2
 301         class1='receptor'
 302         class2='ligand'
 303
 304
 305
 306     #print len1
 307     #print len2
 308     #print chain_sample.keys()
 309
 310     #Set to align on receptor
 311     assert len(chain_ref[receptor_chain])==len(chain_sample[receptor_chain]), "Different number of atoms in native and model receptor (chain %c) %d %d\n" % (receptor_chain,len(chain_ref[receptor_chain]),len(chain_sample[receptor_chain]))
 312
 313     super_imposer.set_atoms(chain_ref[receptor_chain], chain_sample[receptor_chain])
 314     super_imposer.apply(sample_model.get_atoms())
 315     receptor_chain_rms=super_imposer.rms
 316     #print receptor_chain_rms
 317     #print dir(super_imposer)
 318     #print chain1_rms
 319
 320     #Grep out the transformed ligand coords
 321
 322     #print ligand_chain
 323
 324     #print chain_ref[ligand_chain]
 325     #print chain_sample[ligand_chain]
 326     #l1=len(chain_ref[ligand_chain])
 327     #l2=len(chain_sample[ligand_chain])
 328
 329
 330
 331
 332     assert len(chain_ref[ligand_chain])!=0 or len(chain_sample[ligand_chain])!=0, "Zero number of equivalent atoms in native and model ligand (chain %s) %d %d.\nCheck that the residue numbers in model and native is consistent\n" % (ligand_chain,len(chain_ref[ligand_chain]),len(chain_sample[ligand_chain]))
 333
 334
 335     assert len(chain_ref[ligand_chain])==len(chain_sample[ligand_chain]), "Different number of atoms in native and model ligand (chain %c) %d %d\n" % (ligand_chain,len(chain_ref[ligand_chain]),len(chain_sample[ligand_chain]))
 336
 337     coord1=np.array([atom.coord for atom in chain_ref[ligand_chain]])
 338     coord2=np.array([atom.coord for atom in chain_sample[ligand_chain]])
 339
 340     #coord1=np.array([atom.coord for atom in chain_ref[receptor_chain]])
 341     #coord2=np.array([atom.coord for atom in chain_sample[receptor_chain]])
 342
 343     #print len(coord1)
 344     #print len(coord2)
 345
 346     sup=SVDSuperimposer()
 347     Lrms = sup._rms(coord1,coord2) #using the private _rms function which does not superimpose
 348
 349
 350     #super_imposer.set_atoms(chain_ref[ligand_chain], chain_sample[ligand_chain])
 351     #super_imposer.apply(sample_model.get_atoms())
 352     #coord1=np.array([atom.coord for atom in chain_ref[receptor_chain]])
 353     #coord2=np.array([atom.coord for atom in chain_sample[receptor_chain]])
 354     #Rrms= sup._rms(coord1,coord2)
 355     #should give same result as above line
 356     #diff = coord1-coord2
 357     #l = len(diff) #number of atoms
 358     #from math import sqrt
 359     #print sqrt(sum(sum(diff*diff))/l)
 360     #print np.sqrt(np.sum(diff**2)/l)
 361     DockQ=(float(fnat) + 1/(1+(irms/1.5)*(irms/1.5)) + 1/(1+(Lrms/8.5)*(Lrms/8.5)))/3
 362     info={}
 363     info['DockQ']=DockQ
 364     info['irms']=irms
 365     info['Lrms']=Lrms
 366     info['fnat']=fnat
 367     info['nat_correct']=nat_correct
 368     info['nat_total']=nat_total
 369
 370     info['fnonnat']=fnonnat
 371     info['nonnat_count']=nonnat_count
 372     info['model_total']=model_total
 373
 374     info['chain1']=chain1
 375     info['chain2']=chain2
 376     info['len1']=len1
 377     info['len2']=len2
 378     info['class1']=class1
 379     info['class2']=class2
 380
 381     return info
 382
 383 def get_pdb_chains(pdb):
 384     pdb_parser = Bio.PDB.PDBParser(QUIET = True)
 385     pdb_struct = pdb_parser.get_structure("reference", pdb)[0]
 386     chain=[]
 387     for c in pdb_struct:
 388         chain.append(c.id)
 389     return chain
 390 #ATOM   3312  CA
 391 #ATOM   3315  CB  ALA H 126     -21.802  31.674  73.597  1.00 58.05           C
 392
 393 def make_two_chain_pdb(pdb,group1,group2): #renumber from 1
 394     pdb_parser = Bio.PDB.PDBParser(QUIET = True)
 395     pdb_struct = pdb_parser.get_structure("reference", pdb)[0]
 396     for c in pdb_struct:
 397         if c.id in group1:
 398             c.id='A'
 399         if c.id in group2:
 400             c.id='B'
 401     (code,outfile)=tempfile.mkstemp()
 402     io=Bio.PDB.PDBIO()
 403     io.set_structure(pdb_struct)
 404     io.save(outfile)
 405     exec_path=os.path.dirname(os.path.abspath(sys.argv[0]))
 406     cmd=exec_path + '/scripts/renumber_pdb.pl ' + outfile
 407     os.system(cmd)
 408     os.remove(outfile)
 409     return outfile +'.renum'
 410
 411 def change_chain(pdb_string,chain):
 412     new_str=[];
 413     for line in pdb_string:
 414         s=list(line);
 415         s[21]=chain;
 416         new_str.append("".join(s));
 417     return "\n".join(new_str);
 418
 419 def make_two_chain_pdb_perm(pdb,group1,group2): #not yet ready
 420     pdb_chains={}
 421     f=open(pdb);
 422     for line in f.readlines():
 423         if line[0:4] == "ATOM":
 424             #       s=list(line);
 425             #print line
 426             chain=line[21]
 427             atom=line[13:16]
 428             resnum=int(line[22:26])
 429             #        print atom + ':' + str(resnum) +':'
 430             if chain not in pdb_chains:
 431                 pdb_chains[chain]=[]
 432             pdb_chains[chain].append(line)
 433  #       print line
 434  #       s[21]='B'
 435  #       print "".join(s)
 436 #        print chain
 437
 438
 439     f.close()
 440     #sys.exit()
 441     (code,outfile)=tempfile.mkstemp()
 442     f=open(outfile,'w')
 443     for c in group1:
 444      #   print pdb_chains[c]
 445         f.write(change_chain(pdb_chains[c],"A"))
 446     f.write("TER\n");
 447     for c in group2:
 448         f.write(change_chain(pdb_chains[c],"B"))
 449     f.close();
 450     #print outfile
 451     exec_path=os.path.dirname(os.path.abspath(sys.argv[0]))
 452     cmd=exec_path + '/scripts/renumber_pdb.pl ' + outfile
 453     os.system(cmd)
 454     os.remove(outfile)
 455     return outfile +'.renum'
 456
 457 def main():
 458
 459     parser=ArgumentParser(description="DockQ - Quality measure for protein-protein docking models")
 460     parser.add_argument('model',metavar='<model>',type=str,nargs=1,help='path to model file')
 461     parser.add_argument('native',metavar='<native>',type=str,nargs=1,help='path to native file')
 462     parser.add_argument('-capri_peptide',default=False,action='store_true',help='use version for capri_peptide (DockQ cannot not be trusted for this setting)')
 463     parser.add_argument('-short',default=False,action='store_true',help='short output')
 464     parser.add_argument('-verbose',default=False,action='store_true',help='talk a lot!')
 465     parser.add_argument('-quiet',default=False,action='store_true',help='keep quiet!')
 466     parser.add_argument('-useCA',default=False,action='store_true',help='use CA instead of backbone')
 467     parser.add_argument('-skip_check',default=False,action='store_true',help='skip initial check fo speed up on two chain examples')
 468     parser.add_argument('-no_needle',default=False,action='store_true',help='do not use global alignment to fix residue numbering between native and model during chain permutation (use only in case needle is not installed, and the residues between the chains are identical')
 469     parser.add_argument('-perm1',default=False,action='store_true',help='use all chain1 permutations to find maximum DockQ (number of comparisons is n! = 24, if combined with -perm2 there will be n!*m! combinations')
 470     parser.add_argument('-perm2',default=False,action='store_true',help='use all chain2 permutations to find maximum DockQ (number of comparisons is n! = 24, if combined with -perm1 there will be n!*m! combinations')
 471 #    parser.add_argument('-comb',default=False,action='store_true',help='use all cyclicchain permutations to find maximum DockQ (number of comparisons is n!*m! = 24*24 = 576 for two tetramers interacting')
 472     parser.add_argument('-model_chain1',metavar='model_chain1', type=str,nargs='+', help='pdb chain order to group together partner 1')
 473     parser.add_argument('-model_chain2',metavar='model_chain2', type=str,nargs='+', help='pdb chain order to group together partner 2 (complement to partner 1 if undef)')
 474     parser.add_argument('-native_chain1',metavar='native_chain1', type=str,nargs='+', help='pdb chain order to group together from native partner 1')
 475     parser.add_argument('-native_chain2',metavar='native_chain2', type=str,nargs='+', help='pdb chain order to group together from native partner 2 (complement to partner 1 if undef)')
 476
 477
 478     args = parser.parse_args()
 479     #bio_ver=1.64
 480     bio_ver=1.61
 481     if(float(Bio.__version__) < bio_ver):
 482         print("Biopython version (%s) is too old need at least >=%.2f" % (Bio.__version__,bio_ver))
 483         sys.exit()
 484
 485 #    if(len(sys.argv)!=3):
 486 #        print "Usage: ./Dock.py <model> <native>"
 487 #        sys.exit()
 488
 489 #    print args
 490 #    print args.model[0]
 491 #    sys.exit()
 492 #    model=sys.argv[1]
 493 #    native=sys.argv[2]
 494
 495     exec_path=os.path.dirname(os.path.abspath(sys.argv[0]))
 496     fix_numbering=exec_path + '/scripts/fix_numbering.pl'
 497     model=args.model[0]
 498     model_in=model
 499     native=args.native[0]
 500     native_in=native
 501     use_CA_only=args.useCA
 502     capri_peptide=args.capri_peptide
 503
 504     model_chains=[]
 505     native_chains=[]
 506     best_info=''
 507     if(not args.skip_check):
 508         model_chains=get_pdb_chains(model)
 509         native_chains=get_pdb_chains(native)
 510     files_to_clean=[]
 511
 512 #    print native_chains
 513     if((len(model_chains) > 2 or len(native_chains) > 2) and
 514        (args.model_chain1 == None and args.native_chain1 == None)):
 515         print("Multi-chain model need sets of chains to group\nuse -native_chain1 and/or -model_chain1 if you want a different mapping than 1-1")
 516         print("Model chains  : " + str(model_chains))
 517         print("Native chains : " + str(native_chains))
 518         sys.exit()
 519     if not args.skip_check and (len(model_chains) < 2 or len(native_chains)< 2):
 520         print("Need at least two chains in the two inputs\n");
 521         sys.exit()
 522
 523     if len(model_chains) > 2 or len(native_chains)> 2:
 524         group1=model_chains[0]
 525         group2=model_chains[1]
 526         nat_group1=native_chains[0]
 527         nat_group2=native_chains[1]
 528         if(args.model_chain1 != None):
 529             group1=args.model_chain1
 530             nat_group1=group1
 531             if(args.model_chain2 != None):
 532                 group2=args.model_chain2
 533             else:
 534                 #will use the complement from group1
 535                 group2=[]
 536                 for c in model_chains:
 537                     if c not in group1:
 538                         group2.append(c)
 539             nat_group1=group1
 540             nat_group2=group2
 541
 542
 543         if(args.native_chain1 != None):
 544             nat_group1=args.native_chain1
 545             if(args.native_chain2 != None):
 546                 nat_group2=args.native_chain2
 547             else:
 548                 #will use the complement from group1
 549                 nat_group2=[]
 550                 for c in native_chains:
 551                     if c not in nat_group1:
 552                         nat_group2.append(c)
 553
 554         if(args.model_chain1 == None):
 555             group1=nat_group1
 556             group2=nat_group2
 557
 558         #print group1
 559         #print group2
 560
 561         #print "native"
 562         #print nat_group1
 563         #print nat_group2
 564         if(args.verbose):
 565             print('Merging ' + ''.join(group1) + ' -> ' + ''.join(nat_group1) + ' to chain A')
 566             print('Merging ' + ''.join(group2) + ' -> ' + ''.join(nat_group2) + ' to chain B')
 567         native=make_two_chain_pdb_perm(native,nat_group1,nat_group2)
 568         files_to_clean.append(native)
 569         pe=0
 570         if args.perm1 or args.perm2:
 571             best_DockQ=-1;
 572             best_g1=[]
 573             best_g2=[]
 574
 575             iter_perm1=itertools.combinations(group1,len(group1))
 576             iter_perm2=itertools.combinations(group2,len(group2))
 577             if args.perm1:
 578                 iter_perm1=itertools.permutations(group1)
 579             if args.perm2:
 580                 iter_perm2=itertools.permutations(group2)
 581
 582             combos1=[];
 583             combos2=[];
 584             for g1 in iter_perm1:#_temp:
 585                 combos1.append(g1)
 586             for g2 in iter_perm2:
 587                 combos2.append(g2)
 588
 589             for g1 in combos1:
 590                 for g2 in combos2:
 591                     pe=pe+1
 592                    # print str(g1)+' '+str(g2)
 593 #            print pe
 594 #            print group1
 595 #            print group2
 596             pe_tot=pe
 597             pe=1
 598             #sys.exit()
 599             if args.verbose:
 600                 print('Starting chain order permutation search (number of permutations: ' + str(pe_tot) + ')')
 601             for g1 in combos1:
 602                 for g2 in combos2:
 603                 #g2=group2
 604                     model_renum=make_two_chain_pdb_perm(model_in,g1,g2)
 605                     model_fixed=model_renum
 606                     if not args.no_needle:
 607                         fix_numbering_cmd=fix_numbering + ' ' + model_renum + ' ' + native + ' > /dev/null'
 608                         model_fixed=model_renum + ".fixed"
 609                         #                print fix_numbering_cmd
 610                         os.system(fix_numbering_cmd)
 611                         os.remove(model_renum)
 612                         if not os.path.exists(model_fixed):
 613                             print('If you are sure the residues are identical you can use the options -no_needle')
 614                             sys.exit()
 615                     test_dict=calc_DockQ(model_fixed,native,use_CA_only)
 616                     os.remove(model_fixed)
 617                     if not args.quiet:
 618                         print(str(pe)+'/'+str(pe_tot) + ' ' + ''.join(g1) + ' -> ' + ''.join(g2) + ' ' + str(test_dict['DockQ']))
 619                     if(test_dict['DockQ'] > best_DockQ):
 620                         best_DockQ=test_dict['DockQ'];
 621                         info=test_dict
 622                         best_g1=g1
 623                         best_g2=g2
 624                         best_info='Best score ( ' + str(best_DockQ) +' ) found for model -> native, chain1:' + ''.join(best_g1) + ' -> ' + ''.join(nat_group1) + ' chain2:' + ''.join(best_g2) + ' -> ' + ''.join(nat_group2)
 625
 626                         if args.verbose:
 627                             print(best_info)
 628                         if not args.quiet:
 629                             print("Current best " + str(best_DockQ))
 630                     pe=pe+1
 631             if not args.quiet:
 632                 print(best_info)
 633 #            print 'Best score ( ' + str(best_DockQ) +' ) found for ' + str(best_g1) + ' ' + str(best_g2)
 634         else:
 635             model_renum=make_two_chain_pdb_perm(model,group1,group2)
 636             model_fixed=model_renum
 637             if not args.no_needle:
 638                 fix_numbering_cmd=fix_numbering + ' ' + model_renum + ' ' + native + ' > /dev/null'
 639                 model_fixed=model_renum + ".fixed"
 640 #                print fix_numbering_cmd
 641                 os.system(fix_numbering_cmd)
 642                 os.remove(model_renum)
 643                 if not os.path.exists(model_fixed):
 644                     print('If you are sure the residues are identical you can use the options -no_needle')
 645                     sys.exit()
 646             info=calc_DockQ(model_fixed,native,use_CA_only)
 647
 648             #os.system('cp ' + native + ' native_multichain.pdb')
 649             #os.system('cp ' + model_fixed + ' .')
 650             os.remove(model_fixed)
 651 #            files_to_clean.append(model)
 652 #            files_to_clean.append(model_fixed)
 653    #    sys.exit()
 654
 655  #   print native
 656  #   print model
 657     else:
 658         info=calc_DockQ(model,native,use_CA_only=use_CA_only,capri_peptide=capri_peptide) #False):
 659 #        info=calc_DockQ(model,native,use_CA_only=)
 660
 661     irms=info['irms']
 662     Lrms=info['Lrms']
 663     fnat=info['fnat']
 664     DockQ=info['DockQ']
 665     fnonnat=info['fnonnat']
 666
 667     if(args.short):
 668         if capri_peptide:
 669             print(("DockQ-capri_peptide %.3f Fnat %.3f iRMS %.3f LRMS %.3f Fnonnat %.3f %s %s %s" % (DockQ,fnat,irms,Lrms,fnonnat,model_in,native_in,best_info)))
 670         else:
 671             print(("DockQ %.3f Fnat %.3f iRMS %.3f LRMS %.3f Fnonnat %.3f %s %s %s" % (DockQ,fnat,irms,Lrms,fnonnat,model_in,native_in,best_info)))
 672
 673     else:
 674         if capri_peptide:
 675             print('****************************************************************')
 676             print('*                DockQ-CAPRI peptide                           *')
 677             print('*   Do not trust any thing you read....                        *')
 678             print('*   OBS THE DEFINITION OF Fnat and iRMS are different for      *')
 679             print('*   peptides in CAPRI                                          *')
 680             print('*                                                              *')
 681             print('*   For the record:                                            *')
 682             print('*   Definition of contact <4A all heavy atoms (Fnat)           *')
 683             print('*   Definition of interface <8A CB (iRMS)                      *')
 684             print('*   For comments, please email: bjorn.wallner@.liu.se          *')
 685             print('****************************************************************')
 686         else:
 687             print('****************************************************************')
 688             print('*                       DockQ                                  *')
 689             print('*   Scoring function for protein-protein docking models        *')
 690             print('*   Statistics on CAPRI data:                                  *')
 691             print('*    0.00 <= DockQ <  0.23 - Incorrect                         *')
 692             print('*    0.23 <= DockQ <  0.49 - Acceptable quality                *')
 693             print('*    0.49 <= DockQ <  0.80 - Medium quality                    *')
 694             print('*            DockQ >= 0.80 - High quality                      *')
 695             print('*   Reference: Sankar Basu and Bjorn Wallner, DockQ: A quality *')
 696             print('*   measure for protein-protein docking models, submitted      *')
 697             print('*                                                              *')
 698             print('*   For the record:                                            *')
 699             print('*   Definition of contact <5A (Fnat)                           *')
 700             print('*   Definition of interface <10A all heavy atoms (iRMS)        *')
 701             print('*   For comments, please email: bjorn.wallner@.liu.se          *')
 702             print('*                                                              *')
 703             print('****************************************************************')
 704         print(("Model  : %s" % model_in))
 705         print(("Native : %s" % native_in))
 706         if len(best_info):
 707             print(best_info)
 708         print('Number of equivalent residues in chain ' + info['chain1'] + ' ' + str(info['len1']) + ' (' + info['class1'] + ')')
 709         print('Number of equivalent residues in chain ' + info['chain2'] + ' ' + str(info['len2']) + ' (' + info['class2'] + ')')
 710         print(("Fnat %.3f %d correct of %d native contacts" % (info['fnat'],info['nat_correct'],info['nat_total'])))
 711         print(("Fnonnat %.3f %d non-native of %d model contacts" % (info['fnonnat'],info['nonnat_count'],info['model_total'])))
 712         print(("iRMS %.3f" % irms))
 713         print(("LRMS %.3f" % Lrms))
 714        # print 'CAPRI ' + capri_class(fnat,irms,Lrms,capri_peptide=capri_peptide)
 715         peptide_suffix=''
 716         if capri_peptide:
 717             peptide_suffix='_peptide'
 718         print(('CAPRI{} {}'.format(peptide_suffix,capri_class(fnat,irms,Lrms,capri_peptide=capri_peptide))))
 719         print('DockQ_CAPRI ' + capri_class_DockQ(DockQ,capri_peptide=capri_peptide))
 720         peptide_disclaimer=''
 721         if capri_peptide:
 722             peptide_disclaimer='DockQ not reoptimized for CAPRI peptide evaluation'
 723         print(("DockQ {:.3f} {}".format(DockQ,peptide_disclaimer)))
 724
 725     for f in files_to_clean:
 726         os.remove(f)
 727
 728 if __name__ == '__main__':
 729   main()