Added --pdbfile option to parse PDB files provided by the user.

gemme.py

@@ -259,6 +259,10 @@ def parse_command_line():
     retMet_args.add_argument('--jetfile', dest='jetfile', type=str, \
         help="If a jet file is provided, it will skip JET2 calculation and use the precalculated JET2 data in that file. Default is None",
         required=False, default=None)
+
+    retMet_args.add_argument('-p', '--pdbfile', dest='pdbfile', type=str, \
+        help="If a pdb file is provided, it will skip fake pdb file production step and use that file. Default is None",
+        required=False, default=None)
     
     retMet_args.add_argument('--normweightmode', dest='normweightmode', type=str, \
         help="It can be one of these: 'trace', 'trace+pc', 'trace+cv' or 'trace+pc+cv'. Default is 'trace'.",
@@ -277,7 +281,7 @@ def parse_command_line():
     return args
 
 
-def doit(inAli,mutFile,retMet,bFile,fFile,n,N, jetfile, normWeightMode):
+def doit(inAli,mutFile,retMet,bFile,fFile,n,N, jetfile, pdbfile, normWeightMode):
     """
     Perfect explanation for a function: typing the function call exactly!
     doit is basically the main function in disguise!
@@ -285,10 +289,18 @@ def doit(inAli,mutFile,retMet,bFile,fFile,n,N, jetfile, normWeightMode):
     """
     simple = True
 
+    
     prot,seq,nl=extractQuerySeq(inAli)
 
-    createPDB(prot,seq)
- 
+    if(pdbfile == None):
+        createPDB(prot,seq)
+    else:
+        if(os.path.exists(prot+".pdb")):
+            shutil.move(prot+".pdb", prot+"_original.pdb")
+            shutil.copy2(pdbfile, prot+".pdb")
+        else:
+            shutil.copy2(pdbfile, prot+".pdb")
+
     print("query protein: "+prot)
     
     if((jetfile) == None):
@@ -296,7 +308,7 @@ def doit(inAli,mutFile,retMet,bFile,fFile,n,N, jetfile, normWeightMode):
         #It is too much buggy and it has too many dependencies. 
         #Using it with a Docker or Singularity may be the best solution!
         print("computing conservation levels...")
-        launchJET(prot,retMet,bFile,fFile,n,N,nl)
+        launchJET(prot,retMet,bFile,fFile, pdbfile, n,N,nl)
         print("done")
     else:
         print("using previously calculated JET2 data from "+jetfile+"...")
@@ -355,7 +367,7 @@ def doit(inAli,mutFile,retMet,bFile,fFile,n,N, jetfile, normWeightMode):
 def main():
     args = parse_command_line()
     doit(args.input,args.mutations,args.retrievingMethod,args.blastFile,\
-         args.fastaFile,args.nIter,args.NSeqs, args.jetfile, args.normweightmode)
+         args.fastaFile,args.nIter,args.NSeqs, args.jetfile, args.pdbfile, args.normweightmode)
 
 if (__name__ == '__main__'):
     main()

gemmeAnal.py

@@ -47,7 +47,10 @@ def getNbSeq(filename):
 
 # Create a PDB file with dummy CA atoms based on the query sequence
 def createPDB(prot,seq):
-
+	"""
+		If there is not a real PDB file for a given sequence,
+		create a fake PDB containing only CA atoms.  
+	"""
 	d = {'C': 'CYS', 'D': 'ASP', 'S': 'SER',  'Q': 'GLN', 'K': 'LYS',
 	'I': 'ILE', 'P': 'PRO', 'T': 'THR', 'F': 'PHE', 'N': 'ASN', 
 	'G': 'GLY',  'H': 'HIS', 'L': 'LEU', 'R': 'ARG', 'W': 'TRP', 
@@ -66,7 +69,7 @@ def editConfJET(N):
 	return(reCode)
 
 # Run JET to compute TJET values
-def launchJET(prot,retMet,bFile,fFile,n,N,nl):
+def launchJET(prot, retMet, bFile, fFile, pdbfile, n, N, nl):
 	"""
 		Call JET2 and produce prot+"_jet.res" file.
 
@@ -88,6 +91,12 @@ def launchJET(prot,retMet,bFile,fFile,n,N,nl):
 	fFile: string 
         A multiple sequence alignment file obtained with psiblast.
 		It is used only if the retMet (explained above) is input.
+	pdbfile: string 
+        a Protein Data Bank file obtained with rcsb.org or any 
+		computational method like alphafold.
+		If it is None, only JET and PC scores are calculated.
+		Otherwise, CV and other structural-dynamical features also
+		can be calculated.
 	n: int
 		Number of JET2 iterations.
 	N: int
@@ -109,7 +118,13 @@ def launchJET(prot,retMet,bFile,fFile,n,N,nl):
 				shutil.copy2(bFile+".orig ", prot+"_A.psiblast")
 			else:
 				shutil.copy2(bFile+" ", prot+"_A.psiblast")
-			jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+prot+".pdb -o `pwd` -p J -r input -b "+prot+"_A.psiblast -d chain -n "+n+" > "+prot+".out"
+			if(pdbfile == None):
+				jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+\
+						prot+".pdb -o `pwd` -p J -r input -b "+prot+"_A.psiblast -d chain -n "+n+" > "+prot+".out"
+			else:
+				jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+\
+						prot+".pdb -o `pwd` -p AVJ -r input -b "+prot+"_A.psiblast -d chain -n "+n+" > "+prot+".out"
+
 		else:
 			print(N)
 			editConfJET(N)
@@ -125,9 +140,19 @@ def launchJET(prot,retMet,bFile,fFile,n,N,nl):
 			
 			print("\nRunning:\n"+grpcmd)
 			subprocess.call(grpcmd,shell=True)
-			jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+prot+".pdb -o `pwd` -p J -r input -f "+prot+"_A.fasta -d chain -n "+n+" > "+prot+".out"
+			if(pdbfile == None):
+				jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+\
+				     	prot+".pdb -o `pwd` -p J -r input -f "+prot+"_A.fasta -d chain -n "+n+" > "+prot+".out"
+			else:
+				jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+\
+				     	prot+".pdb -o `pwd` -p AVJ -r input -f "+prot+"_A.fasta -d chain -n "+n+" > "+prot+".out"
 	else:
-		jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+prot+".pdb -o `pwd` -p J -r "+retMet+" -d chain -n "+n+" > "+prot+".out"
+		if(pdbfile == None):
+			jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+\
+			    	prot+".pdb -o `pwd` -p J -r "+retMet+" -d chain -n "+n+" > "+prot+".out"
+		else:
+			jetcmd = "java -Xmx1000m -cp $JET2_PATH:$JET2_PATH/jet/extLibs/vecmath.jar jet.JET -c default.conf -i "+\
+			    	prot+".pdb -o `pwd` -p AVJ -r "+retMet+" -d chain -n "+n+" > "+prot+".out"
 	print("\nRunning:\n"+jetcmd)
 	reCode=subprocess.call(jetcmd,shell=True)
 	if os.path.isfile(prot+"/"+prot+"_jet.res"):