-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathFlaGs_linux.py
1950 lines (1781 loc) · 72 KB
/
FlaGs_linux.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
__author__ = "Chayan Kumar Saha, Gemma C. Atkinson"
__copyright__ = "GNU General Public License v3.0"
__email__ = "chayan.sust7@gmail.com"
from Bio import SeqIO
from Bio import Entrez
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_protein, IUPAC
import math, re
import argparse
import ftplib
import socket
import random
import time
from random import randint
import colorsys
import os, sys, os.path, math
import gzip
import getopt
from collections import OrderedDict
import subprocess
from tkinter import *
master = Tk()
os.environ['DISPLAY'] = ':0' #linux
master = Tk()
usage= ''' Description: Identify flanking genes and cluster them based on similarity and visualize the structure; Requirement= Python3, BioPython; tkinter ; Optional Requirement= ETE3. '''
parser = argparse.ArgumentParser(description=usage)
parser.add_argument("-a", "--assemblyList", help=" Protein Accession with assembly Identifier eg. GCF_000001765.3 in a text file separated by newline. ")
parser.add_argument("-p", "--proteinList", help=" Protein Accession eg. XP_ or WP_047256880.1 in a text file separated by newline. ")
parser.add_argument("-l", "--localGenomeList", help=" Genome File name and Protein Accession ")
parser.add_argument("-ld", "--localGenomeDirectory", help=" Path for Local Files, Default directory is './' which is the same directory where the script is located or running from. ")
parser.add_argument("-r", "--redundant", help=" To search all assembly type -r A or -r a but for selected number of assembly eg.,5 for each query use -r 5. ")
parser.add_argument("-e", "--ethreshold", help=" E value threshold. Default = 1e-10 ")
parser.add_argument("-n", "--number", help=" Number of Jackhmmer iterations. Default = 3")
parser.add_argument("-g", "--gene", help=" Number of genes for looking up or downstream. Default = 4 ")
parser.add_argument("-t", "--tree", action="store_true", help=" If you want to see flanking genes along with phylogenetic tree, requires ETE3 installation. By default it will not produce. ")
parser.add_argument("-ts", "--tshape", help=" Size of triangle shapes that represent flanking genes, this option only works when -t is used. Default = 12 ")
parser.add_argument("-tf", "--tfontsize", help=" Size of font inside triangles that represent flanking genes, this option only works when -t is used. Default = 4 ")
parser.add_argument("-to", "--tree_order", action="store_true", help=" Generate Output with Tree, and then use the tree order to generate other view. ")
parser.add_argument("-u", "--user_email", required=True, action="append", metavar="RECIPIENT",default=[], dest="recipients", help=" User Email Address (at least one required) ")
parser.add_argument("-api", "--api_key", help=" NCBI API Key, To get this key kindly check https://ncbiinsights.ncbi.nlm.nih.gov/2017/11/02/new-api-keys-for-the-e-utilities/ ")
parser.add_argument("-o", "--out_prefix", required= True, help=" Any Keyword to define your output eg. MyQuery ")
parser.add_argument("-c", "--cpu", help="Maximum number of parallel CPU workers to use for multithreads. ")
parser.add_argument("-k", "--keep", action="store_true", help=" If you want to keep the intermediate files eg. gff3 use [-k]. By default it will remove. ")
parser.add_argument("-v", "--version", action="version", version='%(prog)s 1.2.6')
parser.add_argument("-vb", "--verbose", action="store_true", help=" Use this option to see the work progress for each query as stdout. ")
args = parser.parse_args()
parser.parse_args()
print('\nStarting FlaGs version 1.2.6 \nPlease only run one instance of FlaGs at a time to avoid making more queries than NCBI’s limit.')
print('For more information, please check https://ncbiinsights.ncbi.nlm.nih.gov/2017/11/02/new-api-keys-for-the-e-utilities/ \n')
print('Checking for RefSeq and Genbank summary files and downloading if needed ... \n')
Entrez.tool = 'FlaGs'
ncbi_time=0.4
timeout = 10
socket.setdefaulttimeout(timeout)
def checkBioPython(): #Checking Biopython Version
import Bio
return (Bio.__version__)
from tkinter.font import Font #Font for postscript to-scale output
myFont12 = Font(size=12)
myFont7 = Font(size=7)
Entrez.email = args.recipients[0] #User email
Entrez.max_tries = 5
Entrez.sleep_between_tries = 60
if not args.localGenomeList:
if args.api_key:
Entrez.api_key = args.api_key #Valid API-key allows 10 queries per seconds, which makes the tool run faster
else:
if args.api_key:
print('Since FlaGs will use Local Data api_key is not necessary, Thanks!')
sys.exit()
#Color generator
def random_color(h=None):
"""Generates a random color in RGB format."""
if not h:
c = random.random()
d = 0.5
e = 0.5
return _hls2hex(c, d, e)
def _hls2hex(c, d, e):
return '#%02x%02x%02x' %tuple(map(lambda f: int(f*255),colorsys.hls_to_rgb(c, d, e)))
def outliner (item):
if item =='#ffffff':
return '#bebebe'
elif item =='#f2f2f2':
return '#008000'
elif item =='#f2f2f3':
return '#000080'
else:
return item
if args.redundant: #Search flanking genes in limited or all available GCFs for each query
if args.redundant.isdigit():
if int(args.redundant)==0:
print('Please use -r option correctly, kindly check help or manual and try again. Thanks!')
sys.exit()
else:
pass
if not args.redundant.isdigit():
if args.redundant.lower()=='a':
pass
else:
print('Please use -r option correctly, kindly check help or manual and try again. Thanks!')
sys.exit()
if args.cpu:
if int(args.cpu)>0:
core=int(args.cpu)
else:
print('Please use number eg, 1,2...')
sys.exit()
if args.assemblyList:
if args.redundant:
print('"-r" option is only works with "-p", please try again with proper command')
sys.exit()
if args.localGenomeList:
if args.redundant:
print('"-r" option is only works with "-p", please try again with proper command')
sys.exit()
if args.tree:
if args.tshape:
if int(args.tshape)>0:
size=int(args.tshape)
else:
print("Kindly input the size of triangles, recommended 12. Not applicable for 0 and negative values")
else:
size=12
else:
if args.tshape:
print("Kindly make sure that you are using -t to make this -ts argument working")
sys.exit()
if args.tree:
if args.tfontsize:
if int(args.tfontsize)>0:
fsize=str(args.tfontsize)
else:
print("Kindly input the font Size required inside triangles, recommended 4. Not applicable for 0 and negative values")
else:
fsize=str(4)
else:
if args.tfontsize:
print("Kindly make sure that you are using -t to make this -tf argument working")
sys.exit()
if args.tree_order:
if not args.tree:
print("Kindly make sure that you are using -t to make this -to argument working")
sys.exit()
if args.ethreshold:
evthresh=args.ethreshold
else:
evthresh="1e-10"
if args.number:
iters=args.number
else:
iters="3"
if args.gene:
if int(args.gene)>0:
s= str(int(args.gene)+1)
else:
print('Please insert positive values, starting from 1')
sys.exit()
else:
s= "5"
if args.localGenomeList:
if args.localGenomeDirectory:
if os.path.isdir(args.localGenomeDirectory):
if args.localGenomeDirectory[-1]=='/':
localDirIn=args.localGenomeDirectory
print('Local Data path : ', localDirIn, '\n')
else:
localDirIn=args.localGenomeDirectory+'/'
print('Local Data path : ', localDirIn, '\n')
else:
print('No directory Found as : '+ args.localGenomeDirectory)
sys.exit()
else:
localDirIn='./'
else:
if args.localGenomeDirectory:
print('Please use -l flag to make -ld flag working')
sys.exit()
if not args.localGenomeList:
localDir='./'
else:
localDir=localDirIn
def checkChar(item): #removing characters
import re
items=item.replace('\t','').replace(' ','')
return re.sub("[a-zA-Z0-9_.]","",items)
queryList=[] #Formatting input as a list
#(queryList)
#eg 1. [['WP_019504790.1', 'GCF_000332195.1'], ['WP_028108719.1', 'GCF_000422645.1']]
#eg 2. [['WP_019504790.1'], ['WP_028108719.1']]
if args.localGenomeList:
with open (args.localGenomeList, 'r') as gList:
for line in gList:
if checkChar(line.rstrip().replace(' ',''))=='':
Line=line.rstrip().replace(' ','').split('\t')
if len(Line)<2:
print('Check Input file, Incorrect Format.')
sys.exit()
else:
newFormat=Line[1]+'\t'+Line[0]
queryList.append(newFormat.split('\t'))
else:
print('The submitted query might include characters not found in NCBI protein accessions eg. > , # , ! etc. Please provide correct format, Thanks!')
sys.exit()
else:
if args.proteinList and not args.assemblyList:
with open (args.proteinList, 'r') as pList:
for line in pList:
if checkChar(line.rstrip().replace(' ',''))=='':
Line=line.rstrip().replace(' ','').split('\t')
if len(Line)>1:
print('Check Input file, Incorrect Format.')
sys.exit()
else:
queryList.append(Line)
else:
print('The submitted query might include characters not found in NCBI protein accessions eg. > , # , ! etc. Please provide correct format, Thanks!')
sys.exit()
elif args.assemblyList and not args.proteinList :
with open (args.assemblyList, 'r') as apList:
for line in apList:
if checkChar(line.rstrip().replace(' ',''))=='':
Line=line.rstrip().replace(' ','').split('\t')
if len(Line)<2:
print('Check Input file, Incorrect Format.')
sys.exit()
else:
newFormat=Line[1]+'\t'+Line[0]
queryList.append(newFormat.split('\t'))
else:
print('The submitted query might include characters not found in NCBI protein accessions eg. > , # , ! etc. Please provide correct format, Thanks!')
sys.exit()
else:
print('Incorrect Input!')
sys.exit()
def accession_from_xp(accession_nr):
"""
:param accession_nr: NCBI protein accession
:return: Bioproject number of all species for that protein which is used to grab Assembly number
"""
i=1
retry=True
while (retry) and i<6: #Retry 5 times with changed socket timeout
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
variableIsquare=i**2
changedtimeout = 10*variableIsquare #10-250 seconds (when i =1 , changedtimeout= 10*1 =10 | when i=2, changedtimeout= 10 * 2square =40 )
socket.setdefaulttimeout(changedtimeout)
time.sleep(ncbi_time)
handle = Entrez.efetch(db="protein", id=accession_nr, rettype="gbwithparts", retmode="text")
if handle:
retry=False
record = SeqIO.read(handle, "genbank")
bioproj=record.dbxrefs
handle.close()
bio=[]
for item in bioproj:
if item.split(':')[0]=='BioProject':
bio.append(item.split(':')[1])
if bio:
return set(bio)
else:
return {'NAI'}
else:
i+=1
retry=True
except Exception as e:
retry=True
i+=1
if not i<6:
print("\t\tQuery {}, not found in database. \n" "\t\tContinuing with the next protein in the list ... \n".format(accession_nr))
return False
def accession_from_wp(accession_nr):
"""
:param accession_nr: NCBI protein accession
:return: Set of assembly number of all species for particular protein
"""
i=1
retry=True
while (retry) and i<6: #Retry 5 times with changed socket timeout
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
variableIsquare=i**2
changedtimeout = 10*variableIsquare #10-250 seconds (when i =1 , changedtimeout= 10*1 =10 | when i=2, changedtimeout= 10 * 2square =40 )
socket.setdefaulttimeout(changedtimeout)
time.sleep(ncbi_time)
handle = Entrez.efetch(db="protein", id=accession_nr, rettype="ipg", retmode="xml")
if handle:
retry=False
if float(checkBioPython())<=1.72:
record = list(Entrez.parse(handle))
handle.close()
assembly = re.findall("GC._\d*\.\d", str(record))
if assembly and len(assembly)>0:
return (set(assembly))
else:
return {'NAI'}
else:
record = Entrez.read(handle)
handle.close()
assembly = re.findall("GC._\d*\.\d", str(record['IPGReport']))
if assembly and len(assembly)>0:
return (set(assembly))
else:
return {'NAI'}
else:
i+=1
retry=True
except Exception as e:
retry=True
i+=1
if not i<6:
print("\t\tQuery {}, not found in database. \n" "\t\tContinuing with the next protein in the list ... \n".format(accession_nr))
return False
def seq_from_wp(accession_nr):
"""
:param accession_nr: NCBI protein accession
:return: Protein Sequence
"""
if accession_nr[-1]!='*':
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
time.sleep(ncbi_time)
handle = Entrez.efetch(db="protein", id=accession_nr, rettype="gbwithparts", retmode="text")
except Exception as e:
print(str(e), ", error in entrez-fetch protein accession, {}, not found in database. \n" "Continuing with the next protein in the list. \nError in function: {}".format(accession_nr, seq_from_wp.__name__))
return False
record = SeqIO.read(handle, "genbank")
handle.close()
return record.description.split('[')[0]+'\t'+record.seq
else:
return accession_nr[:-1]+'\t'+'--'
def remBadChar(item): #removing characters from species name
import re
return re.sub("[^a-zA-Z0-9]"," ",item).replace(" ","_")
def identicalProtID(accnr): #searching for identical proteins
i=1
retry=True
while (retry) and i<6: #Retry 5 times with changed socket timeout
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
variableIsquare=i**2
changedtimeout = 10*variableIsquare #10-250 seconds (when i =1 , changedtimeout= 10*1 =10 | when i=2, changedtimeout= 10 * 2square =40 )
socket.setdefaulttimeout(changedtimeout)
time.sleep(ncbi_time)
epost_1 = Entrez.read(Entrez.epost(db="protein", id=accnr))
webenv = epost_1["WebEnv"]
query_key = epost_1["QueryKey"]
iden_prots = Entrez.efetch(db="protein", rettype='ipg', retmode='text', webenv=epost_1["WebEnv"], query_key=epost_1["QueryKey"])
iAccSet=set()
sAccSet=set()
inrAccSet=set()
snrAccSet=set()
for item in iden_prots:
if item[0:2]!='Id':
if re.search("GC._\d*\.\d", item):
itemLine=item.rstrip().split('\t')
iAccession=itemLine[6]
iAssembly=itemLine[-1]
if iAccession!=accnr and iAccession[2]=='_':
iAccSet.add(iAccession)
if iAccession!=accnr and iAccession[2]!='_':
inrAccSet.add(iAccession)
if iAccession==accnr and iAccession[2]=='_':
sAccSet.add(iAccession)
if iAccession==accnr and iAccession[2]!='_':
snrAccSet.add(iAccession)
if len(iAccSet)>0 and len(sAccSet)==0:
for ids in random.sample(iAccSet,1):
return ids
elif len(sAccSet)!=0:
for ids in random.sample(sAccSet,1):
return ids
else:
if len(inrAccSet)>0 and len(snrAccSet)==0 and len(iAccSet)==0:
for ids in random.sample(inrAccSet,1):
return ids
else:
return accnr
retry=False
except:
retry=True
i+=1
if not i<6:
return accnr
def identicalProtID_WP(accnr): #searching for identical proteins
i=1
retry=True
while (retry) and i<6: #Retry 5 times with changed socket timeout
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
variableIsquare=i**2
changedtimeout = 10*variableIsquare #10-250 seconds (when i =1 , changedtimeout= 10*1 =10 | when i=2, changedtimeout= 10 * 2square =40 )
socket.setdefaulttimeout(changedtimeout)
time.sleep(ncbi_time)
epost_1 = Entrez.read(Entrez.epost(db="protein", id=accnr))
webenv = epost_1["WebEnv"]
query_key = epost_1["QueryKey"]
iden_prots = Entrez.efetch(db="protein", rettype='ipg', retmode='text', webenv=epost_1["WebEnv"], query_key=epost_1["QueryKey"])
iAccSet=set()
for item in iden_prots:
if item[0:2]!='Id':
if re.search("GC._\d*\.\d", item):
itemLine=item.rstrip().split('\t')
iAccession=itemLine[6]
iAssembly=itemLine[-1]
if iAccession!=accnr and iAccession[:3]=='WP_':
iAccSet.add(iAccession)
if len(iAccSet)>0:
for ids in random.sample(iAccSet,1):
return ids
else:
return accnr
retry=False
except:
retry=True
i+=1
if not i<6:
return accnr
def identicalProtID_WP_Sp(accnr): #searching for identical proteins with same assembly for 'NP_417570.1' > WP_000785722.1|GCF_000005845.2
i=1
retry=True
while (retry) and i<6: #Retry 5 times with changed socket timeout
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
variableIsquare=i**2
changedtimeout = 10*variableIsquare #10-250 seconds (when i =1 , changedtimeout= 10*1 =10 | when i=2, changedtimeout= 10 * 2square =40 )
socket.setdefaulttimeout(changedtimeout)
time.sleep(ncbi_time)
epost_1 = Entrez.read(Entrez.epost(db="protein", id=accnr))
webenv = epost_1["WebEnv"]
query_key = epost_1["QueryKey"]
iden_prots = Entrez.efetch(db="protein", rettype='ipg', retmode='text', webenv=epost_1["WebEnv"], query_key=epost_1["QueryKey"])
iAccSetSpecial=set()
iAccNRSetSpecial=set()
iAssemblyList=[]
iAssemblyListNR=[]
for item in iden_prots:
if item[0:2]!='Id':
if re.search("GC._\d*\.\d", item):
itemLine=item.rstrip().split('\t')
iAccession=itemLine[6]
iAssembly=itemLine[-1]
if iAccession==accnr and iAccession[2]=='_':
iAssemblyList.append(iAssembly)
if iAccession==accnr and iAccession[-2]=='.' and iAssembly[0]=='G':
iAssemblyListNR.append(iAssembly)
if iAccession!=accnr and iAccession[-2]=='.' and iAssembly[0]=='G':
assNR=iAccession+'|'+iAssembly
iAccNRSetSpecial.add(assNR)
if iAccession!=accnr and iAccession[:3]=='WP_':
if iAssemblyList:
if iAssembly==iAssemblyList[0]:
assWp=iAccession+'|'+iAssembly
iAccSetSpecial.add(assWp)
if iAccSetSpecial:
for ids in random.sample(iAccSetSpecial,1):
return ids
else:
if iAssemblyList:
for ids in random.sample(iAssemblyList,1):
return accnr+'|'+ids
else:
if iAssemblyListNR:
for nrids in random.sample(iAssemblyListNR,1):
return accnr+'|'+nrids
else:
if iAccNRSetSpecial:
for ids in random.sample(iAccNRSetSpecial,1):
return ids
else:
return '#'
retry=False
except:
retry=True
i+=1
if not i<6:
return '#'
def identicalProtID_redundant(accnr): #searching for identical proteins with same assembly for 'NP_417570.1' > WP_000785722.1|GCF_000005845.2
i=1
retry=True
while (retry) and i<6: #Retry 5 times with changed socket timeout
#Entrez.email = "_@gmail.com" # If you do >3 entrez searches on NCBI per second, your ip will be
# blocked, warning is sent to this email.
try:
variableIsquare=i**2
changedtimeout = 10*variableIsquare #10-250 seconds (when i =1 , changedtimeout= 10*1 =10 | when i=2, changedtimeout= 10 * 2square =40 )
socket.setdefaulttimeout(changedtimeout)
time.sleep(ncbi_time)
epost_1 = Entrez.read(Entrez.epost(db="protein", id=accnr))
webenv = epost_1["WebEnv"]
query_key = epost_1["QueryKey"]
iden_prots = Entrez.efetch(db="protein", rettype='ipg', retmode='text', webenv=epost_1["WebEnv"], query_key=epost_1["QueryKey"])
iAssemblyset=set()
for item in iden_prots:
if item[0:2]!='Id':
if re.search("GC._\d*\.\d", item):
itemLine=item.rstrip().split('\t')
iAccession=itemLine[6]
iAssembly=itemLine[-1]
if iAccession==accnr:
iAssemblyset.add(iAssembly)
if iAssemblyset and len(iAssemblyset)>0:
return iAssemblyset
else:
return '#'
retry=False
except:
retry=True
i+=1
if not i<6:
return '#'
def sortGCFvsGCA(gcagcfSet):
if gcagcfSet!='NAI':
Aset=set()
Fset=set()
for items in gcagcfSet:
if items[2]=='A':
Aset.add(items)
if items[2]=='F':
Fset.add(items)
if len(Fset)>0:
return Fset
elif len(Fset)==0 and len(Aset)>0:
return Aset
else:
return gcagcfSet
else:
gcagcfSet
def des_check(item):
if item:
return item
else:
return 'notFound'
def normalize_strand(item1, item2): #Strand direction change
if item1=='+':
return item2
else:
if item2=='+':
return '-'
else:
return '+'
def up(item):
if item=='+':
return 'Upstream '
else:
return 'Downstream '
def down(item):
if item=='+':
return 'Downstream '
else:
return 'Upstream '
def ups(item):
if item=='+':
return '-'
else:
return '+'
def downs(item):
if item=='+':
return '+'
else:
return '-'
def lcheck(item):
if 1 in item:
return 1
else:
return 0
def postscriptSize(item):
if int(item)<1000:
return(0)
else:
return(int(item)/1000)
def getSpeciesFromGCF(faa,item):
if item!='':
if args.redundant:
return remBadChar(item)+'_'+remBadChar(faa)
else:
return remBadChar(item)
else:
return 'Nothing'
def spLocal(faa,acc): #getting species name using assembly number or accession
if faa in speciesNameFromOnlineDict:
return speciesNameFromOnlineDict[faa]
else:
faaFile=faa+'.faa.gz'
fastaSeq = gzip.open(localDir+faaFile, "rt")
for record in SeqIO.parse(fastaSeq, "fasta"):
if record.id==acc:
if args.redundant:
return remBadChar(record.description.split('[')[-1][:-1])+'_'+remBadChar(faa)
else:
return remBadChar(record.description.split('[')[-1][:-1])
def desLocal(faa,acc):
faaFile=faa+'.faa.gz'
fastaSeq = gzip.open(localDir+faaFile, "rt")
for record in SeqIO.parse(fastaSeq, "fasta"):
if record.id==acc:
return record.description.split('[')[0]
def seqLocal(faa,acc):
faaFile=faa+'.faa.gz'
fastaSeq = gzip.open(localDir+faaFile, "rt")
for record in SeqIO.parse(fastaSeq, "fasta"):
if record.id==acc:
return str(record.seq)
def localNone(item):
if item==None:
return '--'
else:
return item
def seqFasLocal(faa,acc): #making fasta file from accession
if spLocal(faa,acc):
faaFile=faa+'.faa.gz'
fastaSeq = gzip.open(localDir+faaFile, "rt")
for record in SeqIO.parse(fastaSeq, "fasta"):
if record.id==acc.split('#')[0]:
record.id=acc+'|'+spLocal(faa,acc)#.replace(':','_').replace('[','_').replace(']','_')
record.description=''
return record.format("fasta")
else:
faaFile=faa+'.faa.gz'
fastaSeq = gzip.open(localDir+faaFile, "rt")
for record in SeqIO.parse(fastaSeq, "fasta"):
if record.id==acc.split('#')[0]:
if args.redundant:
record.id=acc+'|'+remBadChar(record.description.split('[')[-1][:-1])+'_'+remBadChar(faa)#.replace(':','_').replace('[','_').replace(']','_')
else:
record.id=acc+'|'+remBadChar(record.description.split('[')[-1][:-1])
record.description=''
return record.format("fasta")
def redundantCreate(setDict,nums):
if nums=='A' or nums=='a':
newList=random.sample(setDict,len(setDict))
else:
if len(setDict)>int(nums):
newList=random.sample(setDict,int(nums))
else:
newList=random.sample(setDict,len(setDict))
return newList
#Download assembly summary report from NCBI Refseq and genBank
if not args.localGenomeList:
refDb='./refSeq.db'
genDb='./genBank.db'
if os.path.isfile(refDb):
refDbSize=os.path.getsize(refDb)
else:
refDbSize='0'
if os.path.isfile(genDb):
genDbSize=os.path.getsize(genDb)
else:
genDbSize='0'
ftp = ftplib.FTP('ftp.ncbi.nih.gov', 'anonymous', 'anonymous@ftp.ncbi.nih.gov')
ftp.cwd("/genomes/refseq") # move to refseq directory
filenames = ftp.nlst() # get file/directory names within the directory
if 'assembly_summary_refseq.txt' in filenames:
ftp.sendcmd("TYPE i")
if int(ftp.size('assembly_summary_refseq.txt'))!=int(refDbSize):#check if the previously downloaded db exists and if that's updated to recent one
ftp.retrbinary('RETR ' + 'assembly_summary_refseq.txt', open('refSeq.db', 'wb').write) # get the assembly summary from refseq
else:
pass
ftp_gen = ftplib.FTP('ftp.ncbi.nih.gov', 'anonymous', 'anonymous@ftp.ncbi.nih.gov')
ftp_gen.cwd("/genomes/genbank") # move to refseq directory
filenames = ftp_gen.nlst() # get file/directory names within the directory
if 'assembly_summary_genbank.txt' in filenames:
ftp_gen.sendcmd("TYPE i")
if int(ftp_gen.size('assembly_summary_genbank.txt'))!=int(genDbSize):#check if the previously downloaded db exists and if that's updated to recent one
ftp_gen.retrbinary('RETR ' + 'assembly_summary_genbank.txt', open('genBank.db', 'wb').write) # get the assembly summary from refseq
else:
pass
assemblyName={}
bioDict={} #bioproject as keys and assemble number (eg.GCF_000001765.1) as value
accnr_list_dict={} #create a dictionary accessionNumber is a key and Organism name and ftp Gff3 download Link as value
with open('refSeq.db', 'r') as fileIn:
for line in fileIn:
if line[0]!='#':
Line=line.rstrip().split('\t')
accnr_list_dict[Line[0]]= Line[7]+'\t'+Line[19]
bioDict[Line[1]]=Line[0]
assemblyName[Line[0]]=Line[0]
ftp_gen = ftplib.FTP('ftp.ncbi.nih.gov', 'anonymous', 'anonymous@ftp.ncbi.nih.gov')
ftp_gen.cwd("/genomes/genbank") # move to refseq directory
assemblyName_GCA={}
bioDict_gen={}
accnr_list_dict_gen={} #create a dictionary accessionNumber is a key and Organism name and ftp Gff3 download Link as value
with open('genBank.db', 'r') as fileIn:
for line in fileIn:
if line[0]!='#':
Line=line.rstrip().split('\t')
if len(Line)>19:
if Line[18]=='identical':
if Line[17] in accnr_list_dict:
bioDict_gen[Line[1]]=Line[0]
accnr_list_dict_gen[Line[0]]= accnr_list_dict[Line[17]]
assemblyName_GCA[Line[0]]=Line[17]
else:
accnr_list_dict_gen[Line[0]]=Line[7]+'\t'+Line[19]
else:
accnr_list_dict_gen[Line[0]]=Line[7]+'\t'+Line[19]
bioDict.update(bioDict_gen)
accnr_list_dict.update(accnr_list_dict_gen)
assemblyName.update(assemblyName_GCA)
print ('\n'+ '>> Database Downloaded. Cross-checking of the accession list in progress ...'+ '\n')
q=0
ne=0
queryDict={} #protein Id as query and a set of assembly number as value [either All or Species of interest]
#{'WP_019504790.1#1': {'GCF_000332195.1'}, 'WP_028108719.1#2': {'GCF_000422645.1'}, 'WP_087820443.1#3': {'GCF_900185565.1'}}
#{'WP_019504790.1#1': 'GCF_000332195.1', 'WP_028108719.1#2': 'GCF_000422645.1', 'WP_087820443.1#3': 'GCF_900185565.1'} local
if not args.localGenomeList:
with open (args.out_prefix+'_NameError.txt', 'w') as fbad:
for query in queryList:
q+=1
accession_from_wp_out=''
accession_from_xp_out=''
identicalProtID_Out=''
accession_from_wp_ID_out=''
accession_from_xp_ID_out=''
accession_from_wp_IDSame_out=''
exceptionalWP_out=''
accession_from_wp_exceptional=''
special_out = ''
assembly_from_identical=''
if args.verbose:
print('\t Checking Query '+ query[0] +' ....'+ '('+str(q)+'/'+str(len(queryList))+')')
if len(query)<2:
if query[0][:2]=='WP' and query[0][-2]=='.': #WP Accession full WP_000785722.1
accession_from_wp_out=accession_from_wp(query[0])
if accession_from_wp_out:
queryDict[query[0]+'#'+str(q)]=sortGCFvsGCA(accession_from_wp_out)
else:
ne+=1
print(query[0], file= fbad)
elif query[0][:2]=='XP' and query[0][-2]=='.': #XP Accession full XP_003256407.1
accession_from_xp_out=accession_from_xp(query[0])
if accession_from_xp_out:
assemList=[]
for bioprojs in accession_from_xp_out:
if bioprojs in bioDict:
assemList.append(bioDict[bioprojs])
if assemList:
queryDict[query[0]+'#'+str(q)]=sortGCFvsGCA(set(assemList))
else:
ne+=1
print(query[0], file= fbad)
else: #other accessions can be XP_003256407 , WP_000785722 too
identicalProtID_Out=identicalProtID(query[0])
if identicalProtID_Out!=query[0]: #can be anything XP_ or WP_ or YP_ or NP_
if identicalProtID_Out[:-3]!='XP_': #Not XPs
accession_from_wp_ID_out=accession_from_wp(identicalProtID_Out)
if accession_from_wp_ID_out:
asset=set()
for elements in accession_from_wp_ID_out:
asset.add(elements)
if len(asset)>0:
queryDict[identicalProtID_Out+'#'+str(q)+'.'+query[0]]=sortGCFvsGCA(asset)
else:
ne+=1
print(query[0], file= fbad)
if identicalProtID_Out[:-3]=='XP_': #if XPs
accession_from_xp_ID_out=accession_from_xp(identicalProtID_Out)
if accession_from_xp_ID_out:
assemList=[]
for bioprojs in accession_from_xp_ID_out:
if bioprojs in bioDict:
assemList.append(bioDict[bioprojs])
if assemList:
queryDict[identicalProtID_Out+'#'+str(q)+'.'+query[0]]=sortGCFvsGCA(set(assemList))
else:
ne+=1
print(query[0], file= fbad)
elif identicalProtID_Out==query[0]: #excluding XP Wp pre #YP NP
exceptionalWP_out = identicalProtID_WP(identicalProtID_Out)
special_out = identicalProtID_WP_Sp(identicalProtID_Out) #list Query GCF
if not args.redundant:
if special_out!='#':
asset=set()
asset.add(special_out.split('|')[1])
queryDict[query[0]+'#'+str(q)]=asset
#asset=set()
#asset.add(special_out.split('|')[1])
#if query[0][2]=='_':
# queryDict[query[0]+'#'+str(q)]=asset
#else:
# if query[0]!=special_out.split('|')[0]:
# queryDict[special_out.split('|')[0]+'#'+str(q)+'.'+query[0]]=asset
# else:
# queryDict[query[0]+'#'+str(q)]=asset
else:
if exceptionalWP_out[:3]=='WP_':
accession_from_wp_exceptional=accession_from_wp(exceptionalWP_out)
if accession_from_wp_exceptional:
asset=set()
for elements in accession_from_wp_exceptional:
asset.add(elements)
if len(asset)>0:
if query[0]!=exceptionalWP_out:
queryDict[exceptionalWP_out+'#'+str(q)+'.'+query[0]]=sortGCFvsGCA(asset)
else:
queryDict[exceptionalWP_out+'#'+str(q)]=sortGCFvsGCA(asset)
else:
ne+=1
print(query[0], file= fbad)
if args.redundant:
if exceptionalWP_out[1:3]=='P_':
accession_from_wp_exceptional=accession_from_wp(exceptionalWP_out)
if accession_from_wp_exceptional:
asset=set()
for elements in accession_from_wp_exceptional:
asset.add(elements)
if len(asset)>0:
if query[0]!=exceptionalWP_out:
queryDict[exceptionalWP_out+'#'+str(q)+'.'+query[0]]=sortGCFvsGCA(asset)
else:
queryDict[exceptionalWP_out+'#'+str(q)]=sortGCFvsGCA(asset)
else:
ne+=1
print(query[0], file= fbad)
elif exceptionalWP_out[2]!='_':
assembly_from_identical=identicalProtID_redundant(identicalProtID_Out)
if assembly_from_identical!='#':
asset=set()
for elements in assembly_from_identical:
asset.add(elements)
if len(asset)>0:
#queryDict[identicalProtID_Out+'#'+str(q)+'.'+query[0]]=sortGCFvsGCA(asset)
if query[0]!=exceptionalWP_out:
queryDict[exceptionalWP_out+'#'+str(q)+'.'+query[0]]=sortGCFvsGCA(asset)
else:
queryDict[exceptionalWP_out+'#'+str(q)]=sortGCFvsGCA(asset)
else:
ne+=1
print(query[0], file= fbad)
else:
if query[0][:3]=='XP_' and query[0][-2]=='.': #XP Accession
asset=set()
accession_from_xp_out=accession_from_xp(query[0])
if accession_from_xp_out:
for bioprojs in accession_from_xp_out:
if bioprojs in bioDict:
if bioDict[bioprojs]==query[1]:
asset.add(query[1])
if len(asset)>0:
queryDict[query[0]+'#'+str(q)]=asset
else:
ne+=1
print(query[0], file= fbad)
elif query[0][:3]!='XP_' and query[0][-2]=='.': #not XP Accession
asset=set()
accession_from_wp_out=accession_from_wp(query[0])
if accession_from_wp_out:
for elements in accession_from_wp_out:
if query[1]==elements:
asset.add(query[1])
if len(asset)>0:
queryDict[query[0]+'#'+str(q)]=asset
else:
ne+=1
print(query[0], file= fbad)
else:
for query in queryList:
q+=1
queryDict[query[0]+'#'+str(q)]=query[1]
nai=0
NqueryDict={} #{'WP_019504790.1#1': ['GCF_000332195.1'], 'WP_028108719.1#2': ['GCF_000422645.1'], 'WP_087820443.1#3': ['GCF_900185565.1']}
if args.localGenomeList:
with open (args.out_prefix+'_Insufficient_Info_In_DB.txt', 'w') as fNai:
for query in queryDict:
assemblyIdlist=[]
if queryDict[query]!={'NAI'}:
assemblyId=queryDict[query]
faa_gz=localDir+queryDict[query]+'.faa.gz'
if os.path.isfile(faa_gz):
with gzip.open(localDir+assemblyId+'.faa.gz', 'rb') as faaIn:
for line in faaIn:
if line.decode('utf-8')[0]=='>':
Line=line.decode('utf-8').rstrip()
if '>'+query.split('#')[0]==Line.split(' ')[0]:
gff_gz=localDir+assemblyId+'.gff.gz'
if os.path.isfile(gff_gz):
with gzip.open(localDir+assemblyId+'.gff.gz', 'rb') as lgffIn: #Download and read gff.gz
cds_c=0
name_c=0
prot_c=0
cset=set()
nset=set()
pset=set()
for line in lgffIn:
if line.decode('utf-8')[0]!='#':
Line=line.decode('utf-8').rstrip().split('\t')
if Line[2]=='CDS':
cds_c=1
cset.add(cds_c)
if Line[8].split(';')[3][:5]=='Name=': #eliminates pseudo gene as they don't have 'Name='
name_c=1
nset.add(name_c)
if Line[8].split(';')[3].split('=')[1]==query.split('#')[0]:
assemblyIdlist.append(assemblyId)
NqueryDict[query]=list(set(assemblyIdlist))
prot_c=1
pset.add(prot_c)
else:
prot_c=0
pset.add(prot_c)
else:
name_c=0
nset.add(name_c)
else:
cds_c=0
cset.add(cds_c)
if lcheck(pset)>0:
pass
elif lcheck(pset)==0:
if lcheck(cset)>0 and lcheck(nset)>0:
print(query.split('#')[0],' did not match with supplement local GFF File')