-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCholMine_backend.py
executable file
·1158 lines (1105 loc) · 50.3 KB
/
CholMine_backend.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#! /usr/bin/env python
# CLR & CHD prediction server
#v 1.2, Python 2.7
#
# Nan Liu
#10/01/13
#The server has two options: CLR and CHD. The two options cannot be used at the same time. The current version only support one binding site input, given the protein file, the ligand residue name, the ligand residue number and the ligand chain ID.
# support multiple models and the first model will be chosen.
# the bug about counting ATOM and HETATM lines has been fixed.
# fix the CLR and CHD connect record problem
# let usr define cleft residues
# fill in water in the cleft defined and generate fake ligand
# add judgement of cleft box volume, which ranges [0, 10000]
# add judgement of cleft residue numbers cutoff below 25
# add judgement of fake ligand, if file empty, return
# add judgement of miniPDB as subset of original regular PDB
# add sitemap points count judgement
# fix bug to valid PDB if pdb contains no HETATM
# add tail explanation
from __future__ import division
import argparse
import os
import sys
import numpy as np
from ASCbasePy.utils import pdb
import csv
from datetime import datetime
def parser():
parser = argparse.ArgumentParser(prog="CholMine",usage='%(prog)s [options]',description='version1.00 server',formatter_class=argparse.RawDescriptionHelpFormatter,epilog='''
Typically,
CholMine -CLR -opse -projout /usr/SimSite3D/CholMine -pfile xxxx -lnum xxx -lchain x
CholMine -CHD -opse -projout /usr/SimSite3D/CholMine -pfile xxxx -lnum xxx -lchain x
CholMine -CLR -projout /usr/SimSite3D/CholMine -pfile xxxx -keyresidue xxx''')
parser.add_argument("-CLR", "--CLR", help="CLR predictor turned on",action='store_true')
parser.add_argument("-CHD", "--CHD", help ="CHD predictor turned on", action='store_true')
parser.add_argument("-opse", "--opse", help="output PyMOL pse file",action='store_true')
parser.add_argument("-v", help = "debug, output all the resulting files", action='store_true')
parser.add_argument("-projout", "--projout", help = "usr defined pathway for the output directory, such as /usr/SimSite3D/CholMine")
parser.add_argument("-pfile", "--pfile", help='protein file name such as 2RH1')
# parser.add_argument("-lname", "--lname", help='ligand residue name such as CLR')
parser.add_argument("-lnum", "--lnum", help='ligand residue number such as 412')
parser.add_argument("-lchain", "--lchain", help='ligand chainID such as A. If the chain ID is not supplied, the default chain ID will be _.')
parser.add_argument("-keyresidue", "--keyresidue", help='keyresidue file name such as keyresidue')
return parser
def usage_args():
parser=parser()
args = parser.parse_args()
return args
def usage_help():
parser=parser()
help = parser.print_help()
return help
def directory(projpath, pdbfile):
os.system("rm -rf " + projpath + "/" + pdbfile)
os.system("mkdir " + projpath + "/" + pdbfile)
os.system("mkdir " + projpath + "/" + pdbfile + "/initialdata")
dest = projpath + "/" + pdbfile + "/"
os.system("mkdir " + dest + "dataset")
os.system("mkdir " + dest + "dataset/ligands")
os.system("mkdir " + dest + "dataset/dbase")
return dest
def prepare(dest,wholename,pdbfilefull,linput,logfile,v ):
#prepare initial protein files
input = open(dest + pdbfilefull, "r")
output = open(dest + "initialdata/" + wholename + "_p.pdb", "w")
for line in input:
line = line.strip()
if line.startswith("ATOM"):
print >>output, line
if line.startswith("CONECT"):
print >> output, line
if line.startswith("MODEL"):
line = line.split()
if float(line[1])==2:
break
input.close()
output.close()
#prepare initial ligand file
input = open(dest + pdbfilefull, "r")
output = open(dest + "initialdata/" + wholename + "_l.pdb", "w")
Lig=[]
for line in input:
line = line.strip()
if line.startswith("HETATM"):
# lname = line[17:20]
# lname = lname.strip()
lindex = line[22:26]
lindex = lindex.strip()
cid = line[21]
cid = cid.strip()
if not cid:
cid = ""
# linfo = lname + lindex + cid
linfo = cid + "_" + lindex
if cmp(linput[4:].lower(), linfo.lower())==0:
print >>output, line
Lig.append(line)
if line.startswith("CONECT"):
print >> output, line
if line.startswith("MODEL"):
line = line.split()
if float(line[1])==2:
break
input.close()
output.close()
if not Lig:
print >>logfile, "ERROR: This ligand site cannot be found in the PDB file. Please recheck the chain ID and residue number or the format of the protein file and ensure that the ligand is labeled as HETATM rather than ATOM."
print "ERROR: This ligand site cannot be found in the PDB file. Please recheck the chain ID and residue number or the format of the protein file and ensure that the ligand is labeled as HETATM rather than ATOM."
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(13)
#check whether the mini PDB is subset of the original regular PDB file
def checksubset(dest, keyresidue, pdbfilefull, v):
f = open(dest + keyresidue, "r")
count = 0
list = []
for l in f:
if l.startswith("ATOM"):
count = count +1
l = l.strip()
info = l[17:65]
ref = open(dest + pdbfilefull, "r")
for line in ref:
line = line.strip()
if (info in line):
list.append("yes")
break
ref.close()
f.close()
if len(list) < count:
print "ERROR: The mini PDB need to be a subset of the regular PDB file."
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(26)
def countResidue(dest, keyresidue,pdbfilefull, v):
f = open(dest + keyresidue, "r")
residues = []
for l in f:
if l.startswith("ATOM"):
l = l.strip()
chainID = l[21]
chainID = chainID.strip()
lindex = l[22:26]
lindex = lindex.strip()
info = chainID + lindex
inlist = 1
if len(residues) == 0:
residues.append(info)
for index in residues:
if cmp(index, info) == 0:
inlist=0
break
if inlist ==1:
residues.append(info)
f.close()
count = len(residues)
if (count > 25):
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + pdbfilefull)
sys.exit(23)
print "The recommended cleft number should be under 25. Please reduce the number of residues."
def preparekeyresidue(dest,wholename,pdbfilefull,keyresidue, linput,logfile,v ):
#prepare initial protein files
input = open(dest + pdbfilefull, "r")
output = open(dest + "initialdata/" + wholename + "_p.pdb", "w")
for line in input:
line = line.strip()
if line.startswith("ATOM"):
print >>output, line
if line.startswith("CONECT"):
print >> output, line
if line.startswith("MODEL"):
line = line.split()
if float(line[1])==2:
break
input.close()
output.close()
#prepare initial ligand file
output = open(dest + "center.out", "w")
file = open(dest + keyresidue, "r")
xlist = []
ylist = []
zlist = []
for line in file:
line = line.strip()
if line.startswith("ATOM"):
x = float(line[30:38])
y = float(line[38:46])
z = float(line[46:54])
xlist.append(x)
ylist.append(y)
zlist.append(z)
file.close()
box = []
minxyz = []
maxxyz = []
minxyz.append(min(xlist))
minxyz.append(min(ylist))
minxyz.append(min(zlist))
maxxyz.append(max(xlist))
maxxyz.append(max(ylist))
maxxyz.append(max(zlist))
box.append(minxyz)
box.append(maxxyz)
dis = [(b-a) for (a,b) in zip(*box)]
alocate = [a for (a,b) in zip(*box)]
blocate = [b for (a,b) in zip(*box)]
nx = int(dis[0]/1)
ny = int(dis[1]/1)
nz = int(dis[2]/1)
# judgement of volume
if ((nx*ny*nz)>20000):
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(21)
print "The defined cleft is too big for typical ligand binding site. Please reduce the number of residues in the cleft file."
# nwater = nx*ny*nz
if ((nx*ny*nz)==0):
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(22)
print "The defined cleft has one dimention 0. Please check the coordinates of residues in the cleft file."
f = open(dest + "initialdata/" + wholename + "_linitial.pdb", "w")
# f2 = open(initial_path + index + "twopoints_l.pdb", "w")
formatstr = 'HETATM%6d O HOH%6d %8.3f%8.3f%8.3f1.00100.00 O '
# print >> f2, formatstr % (1, 1, alocate[0], alocate[1], alocate[2])
# print >> f2, formatstr % (2, 2, blocate[0], blocate[1], blocate[2])
# f2.close()
for i in range(nx+1):
coorx = (blocate[0])-1*(i)
for j in range(ny+1):
coory = (blocate[1])-1*(j)
for k in range(nz+1):
coorz = (blocate[2])-1*(k)
print >> f, formatstr %((i+1)*(j+1)*(k+1),(i+1)*(j+1)*(k+1), coorx, coory, coorz)
print >> f, "END"
f.close()
print >> output, wholename + " " + str(nx) + " " + str(ny) + " " + str(nz)
output.close()
f = open(dest + "initialdata/" + wholename + "_linitial.pdb", "r")
fnew = open(dest + "initialdata/" + wholename + "_l.pdb", "w")
for line in f:
line = line.strip()
if line.startswith("HETATM"):
x = float(line[30:38])
y = float(line[38:46])
z = float(line[46:54])
in_list =1
ref = open(dest + "initialdata/" + wholename + "_p.pdb", "r")
for l in ref:
l = l.strip()
if l.startswith("ATOM"):
xref = float(l[30:38])
yref = float(l[38:46])
zref = float(l[46:54])
dis = ((x-xref)**2 + (y-yref)**2 + (z-zref)**2)**0.5
if (dis > 3.4):
continue
else:
# print dis
# print xref, yref, zref
in_list=0
break
if (in_list ==1):
print >> fnew, line
else:
continue
ref.close()
f.close()
fnew.close()
filepath = dest + 'initialdata/' + wholename + '_l.pdb'
if os.stat(filepath).st_size ==0:
if not v:
os.system("rm -rf " + dest + "initialdata")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
# os.system("rm -rf " + dest + "initialdata")
sys.exit(24)
print "This binding site is too small for a possible CLR binding site. Please enlarge the cleft defined."
def gen_points(dest, wholename, mol2_dbase_path, sitemap_dbase_path, logfile, v, pdbfilefull):
os.system("/soft/linux64/openeye/bin/molcharge -in " + dest + "initialdata/" + wholename + "_l.pdb -out " + dest + "initialdata/" + wholename + "_l.mol2")
os.system("mv " + dest + "initialdata/" + wholename + "_l.mol2 " + mol2_dbase_path)
if os.system("/soft/linux64/SimSite3D_v4.5/bin/gen_points -p " + dest + "initialdata/" + wholename + "_p.pdb -l " + mol2_dbase_path + wholename + "_l.mol2 --allow_small_site_maps --no_normalization --include_metals --msms_surf " + dest + wholename + "_s.csv")== 0:
os.system("cp " + dest + wholename + "_s.csv " + sitemap_dbase_path)
os.system("cp " + dest + wholename + "_s.pdb " + sitemap_dbase_path)
os.system("cp " + dest + wholename + "_a.pdb " + sitemap_dbase_path)
os.system("cp " + dest + wholename + "_rad.pdb " + sitemap_dbase_path)
os.system("cp " + dest + wholename + "_surf.face " + sitemap_dbase_path)
os.system("cp " + dest + wholename + "_surf.vert " + sitemap_dbase_path)
print >> logfile, "## The binding site shape and chemical properties have been extracted."
print >> logfile, "## ..........................................................................."
else:
print >> logfile, "ERROR: Sitemap representation not generated for crstalID: " + wholename + " because of CholMine internal problems."
print >> logfile, "..........................................................................."
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "3KDP_CLR3001D_results/")
os.system("rm -rf " + dest + "2DYR_CHD525C_results/")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + "*csv")
os.system("rm -rf " + dest + "*face")
os.system("rm -rf " + dest + "*vert")
os.system("rm -rf " + dest + "rotateddataset")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(14)
def countsitemap(dest, wholename, v, pdbfilefull):
sfile = open(dest + wholename + "_s.pdb", "r")
counthphilic = 0
counts = 0
for line in sfile:
line = line.strip()
if line.startswith("HETATM"):
counts = counts + 1
property = float(line[60:65])
if (property!=100.00):
counthphilic = counthphilic + 1
sfile.close()
if (counts<20 or (counthphilic < 6 and counts < 30)):
if not v:
os.system("rm -rf " + dest + "initialdata")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "3KDP_CLR3001D_results/")
os.system("rm -rf " + dest + "2DYR_CHD525C_results/")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + "*csv")
os.system("rm -rf " + dest + "*face")
os.system("rm -rf " + dest + "*vert")
os.system("rm -rf " + dest + "rotateddataset")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
print "the sitemap points are not enough, please include more residues."
sys.exit(24)
def CLR_searchsites(dest, sitemap_dbase_path, mol2_dbase_path, logfile, wholename, v, pdbfilefull):
if os.system("search_sitemaps --score_threshold 0.0 --lig_frag_size 0 --prot_lig_scoring --fine_tune_tier2 --dbase_sites " + sitemap_dbase_path + " --dbase_ligs " + mol2_dbase_path + " --proj_output " + dest + "3KDP_CLR3001D_results ./hardcodedfiles/3KDP_CLR3001D_s.csv")==0:
print >> logfile, "## The binding site has been compared with the CholMine cholesterol site."
print >> logfile, "## ..........................................................................."
else:
print >> logfile, "ERROR: Site comparison failed for crstalID: " + wholename + " because of CholMine internal problems."
print >> logfile, "..........................................................................."
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "3KDP_CLR3001D_results/")
os.system("rm -rf " + dest + "2DYR_CHD525C_results/")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + "*csv")
os.system("rm -rf " + dest + "*face")
os.system("rm -rf " + dest + "*vert")
os.system("rm -rf " + dest + "rotateddataset")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(15)
def CHD_searchsites(dest, sitemap_dbase_path, mol2_dbase_path, logfile, wholename,v, pdbfilefull):
if os.system("search_sitemaps --score_threshold 0.0 --lig_frag_size 0 --prot_lig_scoring --fine_tune_tier2 --dbase_sites " + sitemap_dbase_path + " --dbase_ligs " + mol2_dbase_path + " --proj_output " + dest + "2DYR_CHD525C_results ./hardcodedfiles/2DYR_CHD525C_s.csv")==0:
print >> logfile, "## The binding site has been compared with CholMine cholate site."
print >> logfile, "## ..........................................................................."
else:
print >> logfile, "Error: Site comparison failed for crstalID: " + wholename + " because of CholMine internal problems."
print >> logfile, "..........................................................................."
if not v:
os.system("rm -rf " + dest + "initialdata/")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "3KDP_CLR3001D_results/")
os.system("rm -rf " + dest + "2DYR_CHD525C_results/")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + "*csv")
os.system("rm -rf " + dest + "*face")
os.system("rm -rf " + dest + "*vert")
os.system("rm -rf " + dest + "rotateddataset")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
sys.exit(15)
def CLRmatchprint(dest):
matchprint_file = open(dest + '3KDP_CLR3001D_results/3KDP_CLR3001D_matchprint.out', 'w')
res_file = open(dest + '3KDP_CLR3001D_results/3KDP_CLR3001D_blast.out', 'r')
for line in res_file:
if(line.startswith("#")): continue
else:
toks = line.split("|")
matchprint = toks[4]
print >> matchprint_file, matchprint, toks[0], toks[1]
matchprint_file.close()
res_file.close()
def CHDmatchprint(dest):
matchprint_file = open(dest + '2DYR_CHD525C_results/2DYR_CHD525C_matchprint.out', 'w')
res_file = open(dest + '2DYR_CHD525C_results/2DYR_CHD525C_blast.out', 'r')
for line in res_file:
if(line.startswith("#")): continue
else:
toks = line.split("|")
matchprint = toks[4]
print >> matchprint_file, matchprint, toks[0], toks[1]
matchprint_file.close()
res_file.close()
def matchprint_analysis(file_path, file_path3, logfile, pdbfile, chainid, ligandindex, linput, dest, v, pdbfilefull):
matrix = []
input = open(file_path, 'r')
for line in input:
line = line.strip()
toks = line.split()
matrix.append(toks[0])
input.close()
#figure out the columns which satifies 85% conservation and put (column number-1) to columnlist
columnlist = []
columnlist2 = []
for i in range(len(matrix[0])):
sum = 0
for j in range(len(matrix)):
sum = sum + int(matrix[j][i])
if (sum < len(matrix)*0.70): # conservation 85%
continue
else:
columnlist.append(i)
number = i + 1
columnlist2.append(number)#need to plus 1 for column number starting from 1 instead of 0
#test set
input3 = open(file_path3, 'r')# open test matchprint file
matrix_test = []
proteinidtest = []
idtest = []
for line in input3:
if line == "":
break
line = line.strip()
toks = line.split()
matrix_test.append(toks[0])
proteinidtest.append(toks[1])
row_testlist = []
for j in range(len(matrix_test)):
sum = 0
for index in columnlist:
sum = sum + int(matrix_test[j][index])
if (sum < len(columnlist)*0.70):
continue
else:
row_testlist.append(j)
input3.close()
# print >> logfile, "How many sites are left after screening: ", len(row_testlist)
marker = chainid + " " + ligandindex
marker = marker.strip()
if len(row_testlist) ==0:
print >> logfile, "## Step 4: Prediction result:"
if cmp(linput[:3], "CLR")==0:
print >> logfile, "## CholMine reports that the site at the " + marker + " in " + pdbfile + " was not predicted as a cholesterol binding site after checking conserved interactions with the prototypic site."
if cmp(linput[:3], "CHD")==0:
print >> logfile, "## CholMine reports that the site at the " + marker + " in " + pdbfile + " was not predicted as a cholate binding site after checking conserved interactions with the prototypic site."
if not v:
os.system("rm -rf " + dest + "initialdata")
os.system("rm -rf " + dest + "dataset")
os.system("rm -rf " + dest + "3KDP_CLR3001D_results/")
os.system("rm -rf " + dest + "2DYR_CHD525C_results/")
os.system("rm -rf " + dest + "*.pdb")
os.system("rm -rf " + dest + "*csv")
os.system("rm -rf " + dest + "*face")
os.system("rm -rf " + dest + "*vert")
os.system("rm -rf " + dest + "rotateddataset")
os.system("rm -rf " + dest + pdbfilefull)
os.system("rm -rf " + dest + "center.out")
if cmp(linput[:3], "CLR")==0:
sys.exit(16)
if cmp(linput[:3], "CHD")==0:
sys.exit(17)
print >> logfile, "## Step 4: Prediction result:"
if cmp(linput[:3], "CLR")==0:
print >> logfile, "## Cholesterol site predicted!"
print >> logfile, "## Based on similarities in shape and conserved interactions, CholMine reports that the site at the " + marker + " in " + pdbfile+" matches known cholesterol site features."
if cmp(linput[:3], "CHD")==0:
print >> logfile, "## Cholate site predicted!"
print >> logfile, "## Based on similarities in shape and conserved interactions, CholMine reports that the site at the " + marker + " in " + pdbfile+" matches known cholate site features."
# print >> logfile, "## ID of unknown site picked by CholMine: "
for i in range(len(row_testlist)):
# print >> logfile, "## " + proteinidtest[row_testlist[i]][:-13]
idtest.append(proteinidtest[row_testlist[i]])
return idtest
def rotatedresidues(dest, input,idtest,logfile,wholename,label,linput):
from ASCbasePy.utils import pdb
os.system("mkdir " + dest + "rotateddataset")
output_path = dest + "rotateddataset/"
coormatrix = []
for line in input:
line = line.strip()
if line.startswith("HETATM"):
keycoor = []
x = float(line[30:38])
y = float(line[38:46])
z = float(line[46:54])
keycoor.append(x)
keycoor.append(y)
keycoor.append(z)
coormatrix.append(keycoor)
input.close()
for index in idtest:
input = open(dest + label + "_results/" + label + "_blast.out", "r")
for line in input:
if(line.startswith("#")): continue
else:
toks = line.split("|")
if cmp(index, toks[0])==0:
M = toks[2] + " " + toks[3]
M = [float(s) for s in M.split()]
R = np.array(M[:9]).reshape(3,3)
T = np.array(M[9:])
pdb.transform(dest + wholename + "_s.pdb", output_path + index[:-13] + "_rotateds.pdb", R.T, T)
pdb.transform(dest + "initialdata/" + wholename + "_p.pdb", output_path + index[:-13] + "_rotatedp.pdb", R.T, T)
keypointindex = []
for index in idtest:
keyeach = []
for i in range(len(coormatrix)):
dis = []
dex = []
chem = []
file = open(output_path + index[:-13] + "_rotateds.pdb", "r")
for line in file:
line = line.strip()
if line.startswith("HETATM"):
x = float(line[30:38])
y = float(line[38:46])
z = float(line[46:54])
# type = float(line[60:66])
# if type == 100.00:
distance =((x-float(coormatrix[i][0]))**2 + (y -float(coormatrix[i][1]))**2 + (z- float(coormatrix[i][2]))**2)**0.5
#print distance
if (distance > 1.5):
continue
else:
# print distance
dis.append(distance)
pointindexall = int(line[20:26])
type = float(line[60:66])
dex.append(pointindexall)
chem.append(type)
file.close()
if len(dis) == 0:
continue
else:
locate = dis.index(min(dis))
pointindex = dex[locate]
chemtype = chem[locate]
if chemtype == 100.00:
innum =1
if (len(keyeach) == 0):
keyeach.append(pointindex)
continue
for ins in keyeach:
if (ins == pointindex):
innum = 0
break
if (innum == 1):
keyeach.append(pointindex)
keypointindex.append(keyeach)
resimatrix = []
for i in range(len(idtest)):
resieach = []
f = open (dest + wholename + "_s.csv", "rb")
reader = csv.reader(f, delimiter='|')
for line in reader:
if line == []:
continue
else:
if len(line) == 3:
# print line
num = line[0]
for j in range(len(keypointindex[i])):
numrefer = keypointindex[i][j]
if cmp(num, str(numrefer)) == 0:
resi = line[1]
resieach.append(resi)
resimatrix.append(resieach)
f.close()
#newmatrix contains residues and atoms making conserved interactions
newmatrix = []
for i in range(len(idtest)):
residues = []
for j in range(len(resimatrix[i])):
test = resimatrix[i][j].split("_")
#print test
for k in range(int(len(test)/3)):
resi = test[(k+1)*3-3] + "_"+ test[(k+1)*3-2] + "_" + test[(k+1)*3-1]
in_list =1
if len(residues) == 0:
residues.append(resi)
for index in residues:
if cmp(index, resi) == 0:
in_list=0
break
if in_list ==1:
residues.append(resi)
newmatrix.append(residues)
#print newmatrix
#newmatrix2 only contains residues making conserved interactions
newmatrix2 = []
for i in range(len(idtest)):
residues2 = []
for j in range(len(resimatrix[i])):
test = resimatrix[i][j].split("_")
for k in range(int(len(test)/3)):
resi = test[(k+1)*3-3]+ "_"+ test[(k+1)*3-2]
if test[(k+1)*3-2] == "":
resi = test[(k+1)*3-3]+ "_"+ " "
in_list =1
if len(residues2) == 0:
residues2.append(resi)
for index in residues2:
if cmp(index, resi) == 0:
in_list=0
break
if in_list ==1:
residues2.append(resi)
newmatrix2.append(residues2)
#print newmatrix2
print >> logfile, "## ..........................................................................."
if cmp(linput[:3], "CLR")==0:
print >> logfile, "## The residues making conserved interactions with cholesterol:"
if cmp(linput[:3], "CHD")==0:
print >> logfile, "## The residues making conserved interactions with cholate:"
print >> logfile, "## SiteID: Residues"
#print >> finalfile, "#ID|important residues|"
for i in range(len(idtest)):
for j in range(len(newmatrix[i])):
test = newmatrix[i][j]
# all = str(test) + " " + all
print >> logfile, idtest[i][:-13] + ": " + test
print >> logfile, "## Note that in the l.pdb file (predicted ligand position) and PyMOL .pse file (if requested), the position of the steroid tail is arbitrary, as it is outside the conserved recognition motif."
logfile.close()
#read rotated pdb file
for i in range(len(idtest)):
pdb = open(output_path + idtest[i][:-13] + "_rotatedp.pdb", "r")
pdbout = open(output_path + idtest[i][:-13] + "_rotatedkeyresidue.pdb", "w")
for line in pdb:
line=line.strip()
if line.startswith("ATOM"):
resiname = line[17:20]
chainID = line[21]
resinum = int(line[22:26])
resiwhole = resiname + str(resinum) + "_" + chainID
# print resiwhole
for j in range(len(newmatrix2[i])):
test = newmatrix2[i][j]
if cmp(resiwhole, test)==0:
print >> pdbout, line
if line.startswith("CONECT"):
print >> pdbout, line
pdb.close()
pdbout.close()
def rotatedback(in_fname, out_fname, R, T):
"""
Assumptions
Assumes R is for premultiplication by coordinates and
R and T are 3x3 and 3x1 numpy arrays
apply transformation to all ATOM and HETATM records
transformation is x' =R.T*x + T
According to the equation, if the coordinates are rotated back, then x = (x'-T)*(1/(R.T))
"""
infile = open(in_fname, "r")
out = open(out_fname, "w+")
fmat = '%8.3f%8.3f%8.3f'
for line in infile:
if(line.startswith("ATOM ") or line.startswith("HETATM")):
position= [float(line[30:38]), float(line[38:46]), float(line[46:54])]
diff = position-T
position= np.dot(diff, R.I)
# position =np.dot((position-T), R.I)
x = np.array(position)[0][0]
y = np.array(position)[0][1]
z = np.array(position)[0][2]
out.write(line[0:30]+fmat%(x, y, z)+line[54:])
else:
out.write(line)
def CLR_rotatedback_proteins(dest, input,idtest,logfile,wholename,label):
for index in idtest:
input = open(dest + label + "_results/" + label + "_blast.out", "r")
output_path = dest + "rotateddataset/"
for line in input:
if(line.startswith("#")): continue
else:
toks = line.split("|")
if cmp(index, toks[0])==0:
M = toks[2] + " " + toks[3]
M = [float(s) for s in M.split()]
R = np.matrix([[M[0],M[1],M[2]],[M[3],M[4],M[5]],[M[6],M[7],M[8]]])
T = np.matrix([M[9],M[10],M[11]])
# os.system("cat " + output_path + index[:-13] + "_rotatedp.pdb ./hardcodedfiles/3KDP_CLR3001D_l.pdb > " + output_path + index[:-13] + "_addligandp.pdb" )
rotatedback(output_path + index[:-13] + "_rotatedp.pdb",output_path + index[:-13] + "_p.pdb", R.T, T)
# rotatedback(output_path + index[:-13] + "_addligandp.pdb",output_path + index[:-13] + "_p.pdb", R.T, T)
rotatedback("./hardcodedfiles/3KDP_CLR3001D_l.pdb",output_path + index[:-13] + "_l.pdb", R.T, T)
rotatedback(output_path + index[:-13] + "_rotatedkeyresidue.pdb",output_path + index[:-13] + "_keyres.pdb", R.T, T)
# rotatedback("./hardcodedfiles/conservedsitemappoints3kdp.pdb", output_path + index[:-13] + "_keypoints.pdb", R.T, T)
def CHD_rotatedback_proteins(dest, input,idtest,logfile,wholename,label):
for index in idtest:
input = open(dest + label + "_results/" + label + "_blast.out", "r")
output_path = dest + "rotateddataset/"
for line in input:
if(line.startswith("#")): continue
else:
toks = line.split("|")
if cmp(index, toks[0])==0:
M = toks[2] + " " + toks[3]
M = [float(s) for s in M.split()]
R = np.matrix([[M[0],M[1],M[2]],[M[3],M[4],M[5]],[M[6],M[7],M[8]]])
T = np.array(M[9:])
os.system("cat " + output_path + index[:-13] + "_rotatedp.pdb ./hardcodedfiles/2DYR_CHD525C_l.pdb > " + output_path + index[:-13] + "_addligandp.pdb" )
rotatedback(output_path + index[:-13] + "_rotatedp.pdb",output_path + index[:-13] + "_p.pdb", R.T, T)
# rotatedback(output_path + index[:-13] + "_addligandp.pdb",output_path + index[:-13] + "_p.pdb", R.T, T)
rotatedback("./hardcodedfiles/2DYR_CHD525C_l.pdb",output_path + index[:-13] + "_l.pdb", R.T, T)
rotatedback(output_path + index[:-13] + "_rotatedkeyresidue.pdb",output_path + index[:-13] + "_keyres.pdb", R.T, T)
# rotatedback("./hardcodedfiles/conservedsitemappoints2DYR.pdb", output_path + index[:-13] + "_keypoints.pdb", R.T, T)
def pseout(dest, logread,wholename,resn):
import __main__
__main__.pymol_argv = ['pymol','-qc'] # Pymol: quiet and no GUI
from time import sleep
import pymol
pymol.finish_launching()
pymol.cmd.do("load " + dest + wholename + ".pdb")
pymol.cmd.do("load " + dest + wholename + "_l.pdb")
pymol.cmd.do("hide ")
pymol.cmd.do("show cartoon, " + wholename)
pymol.cmd.do("cartoon loop")
pymol.cmd.do("show sticks, " + wholename + "_l")
pymol.cmd.do("zoom resn " + resn)
# pymol.cmd.do("show surface, " + wholename)
# pymol.cmd.do("set surface_cavity_mode, 1")
# pymol.cmd.do("set cartoon_oval_quality, 2")
# pymol.cmd.do("set cartoon_tube_radius, 0.075")
# pymol.cmd.do("set cartoon_loop_radius, 0.075")
pymol.cmd.do("bg_color white")
# pymol.cmd.do("load " + dest + wholename + "_keypoints.pdb")
# pymol.cmd.do("color magenta, " + wholename + "_keypoints")
# pymol.cmd.do("show spheres, " + wholename + "_keypoints")
pymol.cmd.do("load " + dest + wholename + "_keyres.pdb")
pymol.cmd.do("show sticks, " + wholename + "_keyres")
pymol.cmd.do("show surface, " + wholename + "_keyres")
pymol.cmd.do("set line_width, 2")
pymol.cmd.do("util.cbaw " + wholename + "_keyres")
# pymol.cmd.do("zoom " + wholename + "_keyres")
pymol.cmd.do("set transparency, 0.5")
for line in logread:
line = line.strip()
if line.startswith("##"):
continue
toks = line.split(":")
id = toks[0]
toks[1] = toks[1].strip()
info = toks[1].split("_")
chain = info[1]
rname = info[0][:3]
rnum = info[0][3:]
atype = info[2]
pymol.cmd.do("select a, /" + id + "_keyres//" + chain + "/" + rname + "`" + rnum+ "/"+atype)
pymol.cmd.do("show spheres, /" + id + "_keyres//" + chain + "/" + rname + "`" + rnum+ "/"+atype)
pymol.cmd.do("color green, /" + id + "_keyres//" + chain + "/" + rname + "`" +rnum+ "/"+atype)
# ref = open(dest + wholename + "_keypoints.pdb", "r")
# i =0
# for line2 in ref:
# i = i+1
# if line2.startswith("HETATM"):
# Onum = line2[24:26].strip()
# print Onum
# pymol.cmd.distance("distance"+str(i), "a", "/" + wholename + "_keypoints///HOH`" + Onum + "/O", 4.0)
# ref.close()
logread.close()
# pymol.cmd.do("hide labels, distance*")
# pymol.cmd.do("set dash_gap, 0.5")
# pymol.cmd.do("set dash_radius, 0.1")
pymol.cmd.do("set sphere_scale, 0.3")
pymol.cmd.do("set surface_quality, 1")
pymol.cmd.do("set solvent_radius, 2")
pymol.cmd.do("delete a")
pymol.cmd.save(dest + wholename + ".pse")
def validPDB(pdbfilefull, ligandindex):
#Judge whether the pdb file is valid pdb file
input = open(pdbfilefull, "r")
atom=0
hetatm = 0
for line in input:
line = line.strip()
toks = line.split()
if len(toks)==0:
continue
firstword=toks[0]
if firstword.startswith("ATOM"):
atom=atom+1
if firstword.startswith("HETATM"):
hetatm=hetatm+1
input.close()
if atom==0:
print "ERROR: Not a valid PDB file; missing protein information (ATOM lines)."
sys.exit(6)
if hetatm==0:
if cmp(ligandindex, "cleft")!=0:
print "ERROR: Not a valid PDB file; missing ligand information (HETATM lines)."
sys.exit(7)
input = open(pdbfilefull, "r")
for line in input:
line=line.strip()
if line.startswith("ATOM"):
if len(line) < 66:
print "ERROR: Not a valid PDB file; missing protein information (ATOM lines)."
sys.exit(6)
if line.startswith("HETATM"):
if len(line) < 66:
print "ERROR: Not a valid PDB file; missing ligand information (HETATM lines)."
sys.exit(7)
def validCleft(keyresidue):
input = open(keyresidue, "r")
atom=0
for line in input:
line = line.strip()
toks = line.split()
if len(toks)==0:
continue
firstword=toks[0]
if firstword.startswith("ATOM"):
atom=atom+1
input.close()
if atom==0:
print "ERROR: Not a valid Cleft PDB file; missing cleft residue information (ATOM lines)."
sys.exit(25)
input = open(keyresidue, "r")
for line in input:
line=line.strip()
if line.startswith("ATOM"):
if len(line) < 66:
print "ERROR: Not a valid Cleft PDB file; missing cleft residue information (ATOM lines)."
sys.exit(25)
def main():
if len(sys.argv) < 6:
print "ERROR: The command line has no enough arguments. Please check the usage."
print usage_help()
sys.exit(19)
args=usage_args()
# if not args.pname or not args.lname or not args.lnum:
if not args.pfile:
print "ERROR: The command line is missing an argument parameter for protein file. Please upload protein file."
print usage_help()
sys.exit(1)
if not args.lnum and not args.keyresidue:
print "ERROR: The command line is missing an argument parameter for ligand residue number or cleft pdb file. Please enter ligand residue number."
sys.exit(2)
if not args.CHD and not args.CLR:
print "ERROR: Please specify cholesterol or cholate predictor."
print usage_help()
sys.exit(3)
pdbfilefull= args.pfile
# ligandcode= args.lname
ligandindex=args.lnum
projpath = args.projout
keyresidue = args.keyresidue
if args.lnum and args.keyresidue:
print "lnum and keysidue can not be turned on at the same time."
sys.exit(20)
if args.CLR:
if args.lnum:
if not args.lchain or args.lchain ==" ":
chainid = ""
linput = "CLR_" + "" + "_" + ligandindex.strip()
else:
chainid=args.lchain
if not chainid.isalpha():
print "ERROR: Not a valid chain ID."
sys.exit(4)
linput = "CLR_" + chainid.strip() + "_" + ligandindex.strip()
if args.keyresidue:
chainid = ""
ligandindex="cleft"
linput = "CLR_" + chainid.strip() + "_" + ligandindex.strip()
if args.CHD:
if args.lnum:
if not args.lchain or args.lchain ==" ":
chainid = ""
# linput = ligandcode.strip() + ligandindex.strip()+ "A"
linput = "CHD_" + "" + "_" + ligandindex.strip()
else:
chainid=args.lchain
if not chainid.isalpha():
print "ERROR: Not a valid chain ID."
sys.exit(4)
# linput = ligandcode.strip() + ligandindex.strip() + chainid.strip()
linput = "CHD_" + chainid.strip() + "_" + ligandindex.strip()
if args.keyresidue:
chainid = ""
ligandindex="cleft"
linput = "CHD_" + "" + "_" + ligandindex.strip()
num = pdbfilefull.count(".")
if num ==0:
pdbfile = pdbfilefull
wholename = pdbfilefull + "_" + linput
elif num>=1:
for i in range(len(pdbfilefull)):
if cmp(pdbfilefull[i], ".")==0:
indexnum =i
pdbfile = pdbfilefull[0:indexnum]
wholename = pdbfilefull[0:indexnum] + "_" + linput
# Judge whether the file exist or not"
v = args.v
if os.path.isfile(pdbfilefull) == False:
print "ERROR: The PDB file doesn't exist."
sys.exit(5)
#Judge whether the pdb file is valid pdb file
print ligandindex
validPDB(pdbfilefull, ligandindex)
if args.keyresidue:
validCleft(keyresidue)
marker = chainid + " " + ligandindex
marker = marker.strip()
#build directory
dest = directory(projpath, pdbfile )
mol2_dbase_path = dest + "dataset/ligands/"
sitemap_dbase_path = dest + "dataset/dbase/"
logfile = open(dest + wholename + "_results.txt", "w")
# print >> logfile, str(len(sys.argv))
os.system("cp " + pdbfilefull + " " + dest)
if args.keyresidue:
os.system("cp " + keyresidue + " " + dest)
if args.CLR:
if (args.CLR and args.CHD):
print "ERROR: Cholesterol and cholate predictor cannot be turned on at the same time; do successive runs."
sys.exit(8)
print "CLR predictor is turned on."
print >> logfile, "## " + str(datetime.now())
print >> logfile, "## CholMine server 1.00"
print >> logfile, "## CLR predictor is turned on."
print >> logfile, "## The binding site at the " + marker + " in " + pdbfile + " is being analyzed."
if args.lnum:
print >> logfile, "## Step 1: Read protein and ligand information and generate initial files for the binding site."
print >> logfile, "## ..........................................................................."
prepare(dest,wholename,pdbfilefull,linput,logfile,v)
print >> logfile, "## Protein and ligand information have been read."
if args.keyresidue:
print >> logfile, "## Step 1: Read protein and cleft information and generate initial files for the binding site."
print >> logfile, "## ..........................................................................."
checksubset(dest, keyresidue, pdbfilefull, v)
countResidue(dest, keyresidue,pdbfilefull, v)
preparekeyresidue(dest,wholename,pdbfilefull,keyresidue, linput,logfile,v)