-
Notifications
You must be signed in to change notification settings - Fork 0
/
phyloutils.py
executable file
·1982 lines (1745 loc) · 89 KB
/
phyloutils.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
# Copyright (C) 2016 Shengwei Hou
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import os
import sys
import argparse
import numpy as np
from Bio import SeqIO
from Bio import AlignIO
from collections import Counter
from ftplib import FTP
import ftputil
import subprocess
import gzip
import urllib2
def main_usage(parser):
""" display usage for main parser
"""
print >>sys.stderr, parser.format_help()
def subparser_usage(argv, parser):
""" display usage for subparser
"""
cmd = argv[1]
found = 0
for action in parser._actions:
if isinstance(action, argparse._SubParsersAction):
for choice, subparser in action.choices.items():
if cmd == choice:
print >>sys.stderr, subparser.format_help()
found = 1
if not found:
if cmd in ("-v", "--version", "-h", "--help"):
args = parser.parse_args()
else:
print >>sys.stderr, "\n\nERROR:%s is not a valid command!!!\n\n" % cmd
main_usage(parser)
def display_help(argv, parser):
""" display help information
"""
if len(argv) == 1:
main_usage(parser)
sys.exit(1)
elif len(argv) == 2:
subparser_usage(argv, parser)
sys.exit(1)
else:
pass
class SubCommandError(Exception):
def __init__(self, value):
self.value = value
def __str__(self):
return repr(self.value)
def html2text(strText):
"""
This is a function to convert html file into plain text. Please credit
the original code:
http://stackoverflow.com/questions/14694482/converting-html-to-text-with-python
"""
str1 = strText
int2 = str1.lower().find("<body")
if int2 > 0:
str1 = str1[int2:]
int2 = str1.lower().find("</body>")
if int2 > 0:
str1 = str1[:int2]
list1 = ['<br>', '<tr', '<td', '</p>', 'span>', 'li>', '</h', 'div>']
list2 = [chr(13), chr(13), chr(9), chr(13), chr(13),
chr(13), chr(13), chr(13)]
bolFlag1 = True
bolFlag2 = True
strReturn = ""
for int1 in range(len(str1)):
str2 = str1[int1]
for int2 in range(len(list1)):
if str1[int1:int1+len(list1[int2])].lower() == list1[int2]:
strReturn = strReturn + list2[int2]
if str1[int1:int1+7].lower() == '<script' or str1[int1:int1+9].lower() == '<noscript':
bolFlag1 = False
if str1[int1:int1+6].lower() == '<style':
bolFlag1 = False
if str1[int1:int1+7].lower() == '</style':
bolFlag1 = True
if str1[int1:int1+9].lower() == '</script>' or str1[int1:int1+11].lower() == '</noscript>':
bolFlag1 = True
if str2 == '<':
bolFlag2 = False
if bolFlag1 and bolFlag2 and (ord(str2) != 10):
strReturn = strReturn + str2
if str2 == '>':
bolFlag2 = True
if bolFlag1 and bolFlag2:
strReturn = strReturn.replace(chr(32)+chr(13), chr(13))
strReturn = strReturn.replace(chr(9)+chr(13), chr(13))
strReturn = strReturn.replace(chr(13)+chr(32), chr(13))
strReturn = strReturn.replace(chr(13)+chr(9), chr(13))
strReturn = strReturn.replace(chr(13)+chr(13), chr(13))
strReturn = strReturn.replace(chr(13), '\n')
return strReturn
def merge_multiple_fasta_to_single_record(input_mfa, delim="N"*100):
""" given a mfa file, read the contents of this mfa.
If only one fasta, return this fasta sequence;
If more than one, connect these sequences using delim then return.
"""
# handle gzip compressed file
if input_mfa.endswith(".gz"):
ih = gzip.open(input_mfa, "r")
else:
ih = open(input_mfa, "r")
# merge fasta records in mfa
single_seq = ""
fasta_records = SeqIO.parse(ih, "fasta")
for fasta in fasta_records:
# if single_seq already assigned, first add delim, then add seq
if single_seq:
single_seq += delim
single_seq += str(fasta.seq)
# if single_seq is not assigned, then assign fasta seq to it
else:
single_seq = str(fasta.seq)
return single_seq
def mergeMultipleFastaRecords(args):
""" given a mfa file, call merge_multiple_fasta_to_single_record,
the
"""
input_mfa = args.input_mfa
delim = args.numberOfNs*'N'
header = args.header
with open(header+".fa", "w") as oh:
genome_seq = merge_multiple_fasta_to_single_record(input_mfa, delim)
oh.write(">"+header+"\n")
oh.write(str(genome_seq)+"\n")
def ftp_download_assemblies(assembly2ftp, outdir):
""" assembly2ftp is a dictionary contains assemblyID and it's ftp address
"""
# login NCBI ftp site
host = ftputil.FTPHost('ftp.ncbi.nlm.nih.gov', 'anonymous',
'housw2010@gmail.com')
# download all the files in the directory
for assembly, ftp_path in assembly2ftp.items():
# create a local folder
local_folder = os.path.join(outdir, assembly)
subprocess.check_call(['mkdir', '-p', local_folder])
# move to the ftp_path folder on host
server_folder = ftp_path.split("ftp.ncbi.nlm.nih.gov")[-1]
host.chdir(server_folder)
file_list = host.listdir(host.curdir)
for remote_file in file_list:
if host.path.isfile(remote_file) and \
(not host.path.islink(remote_file)) and \
(remote_file.endswith('.gz')):
print "[getGenomesFromGenbank]: Downloading %s ..." % \
remote_file
local_file = os.path.join(local_folder, remote_file)
try:
# arguments: remote filename, local filename, callback
# ftputil.error.FTPOSError: Debugging info: ftputil 3.2,
# Python 2.7.12 (linux2)
host.download(remote_file, local_file)
except Exception as e:
print "[getGenomesFromGenbank]:ERROR: CANNOT DOWNLOAD %s"%remote_file
print "[getGenomesFromGenbank]:ERROR: MESSAGE IS: %s"%e
def download_assemblies(assembly2ftp, outdir):
""" assembly2url is a dictionary contains assemblyID and it's NCBI ftp file address,
this is an alternative function of ftp_download_assemblies, in case that
doesn't work
"""
# download all the files in the directory
for assembly, ftp_path in assembly2ftp.items():
# create a local folder
local_folder = os.path.join(outdir, assembly)
remote_folder = os.path.split(ftp_path)[-1]
subprocess.check_call(['mkdir', '-p', local_folder])
try:
ftp_folder = urllib2.urlopen(ftp_path)
except ValueError as e:
print "[getGenomesFromGenbank]: ERROR: can NOT download %s using address: %s"%(assembly, ftp_path)
print "[getGenomesFromGenbank]: The error message is : %s"%e
continue
for f in ftp_folder:
# here convert the html content to text, and split into lines
for line in html2text(f).split("\n"):
line = line.strip().split()
if len(line) <=1:
continue
remote_file = line[-1]
if remote_file.startswith(remote_folder) and remote_file.endswith(".gz"):
print "[getGenomesFromGenbank]: Downloading %s ..."%remote_file
local_file = os.path.join(local_folder, remote_file)
try:
success = False
tries = 0
cmd="wget %s -O %s"%(os.path.join(ftp_path, remote_file), local_file)
while not success:
tries += 1
p = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = p.communicate()
if int(p.returncode) == 0:
success = True
if tries >=10:
print "ERROR: %s was NOT downloaded after 10 tries!"%remote_file
break
except Exception as e:
print "[getGenomesFromGenbank]:ERROR: CANNOT DOWNLOAD %s"%remote_file
print "[getGenomesFromGenbank]:ERROR: MESSAGE IS: %s"%e
def getGenomesFromGenbank(args):
"""download genome sequences from NCBI, given an input assembly report file
"""
input_file = args.input_file
type = args.type
if type == "report":
# parse assembly report to get assembly ID and ftp path
assembly2ftp = {} # {assembly: ftp}
with open(input_file, "r") as ih:
header = ih.readline().strip().split('\t')
assemby_idx = header.index('Assembly')
ftp_idx = header.index('GenBank FTP')
for line in ih:
line = line.strip().split("\t")
assembly = line[assemby_idx].strip()
ftp = line[ftp_idx].strip()
assert assembly not in assembly2ftp, \
"Repeating genbank assembly ID %s has been found!"%assembly
print "[getGenomesFromGenbank]: %s ==> %s"%(assembly, ftp)
assembly2ftp.update({assembly:ftp})
ftp_download_assemblies(assembly2ftp, args.outdir)
else:
# read assembly summary from genbank, store ftp path for each assembly
ID2path = {}
assembly_summary_genbank = "ftp://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt"
response = urllib2.urlopen(assembly_summary_genbank)
for line in response:
# consume comments
if line.startswith("#"):
continue
sline = line.strip().split("\t")
ID = sline[0].strip()
path = sline[-1].strip()
if path:
ID2path[ID] = path
# read input genbank ID file, get IDs
assembly2ftp = {} # {genbank_ID: ftp}
with open(input_file, "r") as ih:
for line in ih:
if line.startswith("GenBank Assembly ID") or line.startswith('#'):
continue
line = line.strip().split("\t")
assembly = line[0].strip()
ftp = ID2path[assembly]
print "[getGenomesFromGenbank]: %s ==> %s"%(assembly, ftp)
assembly2ftp[assembly] = ftp
download_assemblies(assembly2ftp, args.outdir)
def generateShortNameTable(args):
""" generate a tab-deliminated lookUpTable of accession numbers and
specie names
"""
input_file = args.input_file
outdir = args.outdir
suffix = args.suffix
prefix = args.prefix
def _clean_name(sciName):
""" replace characters to underscores or empty
"""
sciName = sciName.replace(" ", "_").replace("'", "").replace("-", "_")
sciName = sciName.replace("/", "_").replace("=", "_")
sciName = sciName.replace("[", "").replace("]", "")
sciName = sciName.replace("(", "_").replace(")", "_")
sciName = sciName.replace("___", "_" ).replace("__", "_")
return sciName
# save all the names
names = set()
if prefix:
output_file = prefix+"_shortNames.tab"
else:
basename = os.path.basename(input_file)
filestem = os.path.splitext(basename)[0]
output_file = filestem+"_shortNames.tab"
# read in genome assembly report file
with open(input_file, "r") as ih, open(output_file, "w") as oh:
oh.write("#filename\tshortname\tgenus\tstrain\n")
header = ih.readline().strip().split("\t")
genus_pos = header.index("Organism/Name")
strain_pos = header.index("Strain")
accession_pos = header.index("Assembly")
Genbank_FTP_pos = header.index("GenBank FTP")
for line in ih:
line = line.strip().split("\t")
genus_full = line[genus_pos].strip().split(" ")[0].lstrip("[").rstrip("]")
genus_short = genus_full[:5]
accession = line[accession_pos].strip()
assemblyID = line[Genbank_FTP_pos].strip().split("_")[-1] if "ASM" in line[-1].strip().split("_")[-1] else accession
strain = line[strain_pos].strip().replace(" ", "")
if strain:
sciName = genus_short+"_"+strain+"_"+assemblyID
else:
sciName = genus_short+"_"+assemblyID
sciName = _clean_name(sciName)
accession += "."+suffix
oh.write(accession+"\t"+sciName+"\t"+genus_full+"\t"+strain+"\n")
def mergeGenbankGenomes(args):
""" given the genbank genome records downloaded using getGenomesFromGenbank,
merge fasta records for each assembly into one fasta record.
"""
delim = args.numberOfNs*'N'
input_folder = args.input_GenbankGenomes
output_folder = args.outdir
subprocess.check_call(['mkdir', '-p', output_folder])
# input folder has a lot of subfolders, each subfolder contains its genome information
for folder in os.listdir(input_folder):
# take the fodler name as its genome name
genome_name = folder
print "[mergeGenbankGenomes]: Merging genome sequence: "+genome_name
with open(os.path.join(output_folder, genome_name+".mfa"), "w") as oh:
for f in os.listdir(os.path.join(input_folder, folder)):
if f.endswith("_genomic.fna.gz") and \
("_cds_from_genomic" not in f) and \
("_rna_from_genomic" not in f):
gz_fna = os.path.join(input_folder, folder, f)
genome_seq = merge_multiple_fasta_to_single_record(gz_fna, delim)
# write out fasta record for current folder
if genome_seq:
oh.write(">"+genome_name+"\n")
oh.write(str(genome_seq)+"\n")
else:
print "No genome fna was found for %s"%genome_name
def extractGenbankGenomicFna(args):
""" given the genbank genome records downloaded using getGenomesFromGenbank,
extract the *_genomic.fna records from each genome subfolder, copy them to the
specified folder.
"""
input_folder = args.input_GenbankGenomes
output_folder = args.outdir
subprocess.check_call(['mkdir', '-p', output_folder])
# input folder has a lot of subfolders, each subfolder contains its genome information
for folder in os.listdir(input_folder):
# take the fodler name as its genome name
genome_name = folder
print "[extractGenbankGenomicFna]: Extracting genomic fna: "+genome_name
output_file = os.path.join(output_folder, genome_name+".fna")
for f in os.listdir(os.path.join(input_folder, folder)):
if f.endswith("_genomic.fna.gz") and ("_from_" not in f):
gz_fna = os.path.join(input_folder, folder, f)
#print(gz_fna)
try:
with gzip.open(gz_fna, "r") as ih, open(output_file, "w") as oh:
for line in ih:
oh.write(line)
except Exception as e:
err = "[extractGenbankGenomicFna]: Error raised: %s"%e
raise SubCommandError(err)
break
else:
print "[extractGenbankGenomicFna]: No protein file was found for %s"%genome_name
continue
def extractGenbankProteins(args):
""" given the genbank genome records downloaded using getGenomesFromGenbank,
extract the protein records from each genome subfolder, copy them to the
specified folder.
"""
input_folder = args.input_GenbankGenomes
output_folder = args.outdir
subprocess.check_call(['mkdir', '-p', output_folder])
# input folder has a lot of subfolders, each subfolder contains its genome information
for folder in os.listdir(input_folder):
# take the fodler name as its genome name
genome_name = folder
print "[extractGenbankProteins]: Extracting genome proteins: "+genome_name
output_file = os.path.join(output_folder, genome_name+".faa")
for f in os.listdir(os.path.join(input_folder, folder)):
if f.endswith("_protein.faa.gz"):
gz_faa = os.path.join(input_folder, folder, f)
try:
with gzip.open(gz_faa, "r") as ih, open(output_file, "w") as oh:
for line in ih:
oh.write(line)
except Exception as e:
err = "[extractGenbankProteins]: Error raised: %s"%e
raise SubCommandError(err)
else:
print "[extractGenbankProteins]: No protein file was found for %s"%genome_name
continue
def renameFastaFile(args):
"""given an input folder with fasta files inside, rename each fasta file and
it's header name, accroding to the lookup table
"""
input_folder = args.inputFastaFolder
lookupTable = args.lookupTable
output_folder = args.outdir
length = args.length
subprocess.check_call(['mkdir', '-p', output_folder])
# read in lookup table
lookup_dict = {}
newNames = set()
with open(lookupTable, "r") as ih:
for line in ih:
if line.startswith('#'):
continue
line = line.strip().split("\t")
k, v = line[0], line[1]
assert v not in newNames, "repeat new name was not allowed: %s"%v
assert len(v) <= length, "the length of %s was more than %d!"%(v, length)
newNames.add(v)
lookup_dict[k] = v
# rename fasta records
for fa in os.listdir(input_folder):
if fa in lookup_dict:
newName = lookup_dict[fa]
print "[renameFastaFile]: renaming %s to %s"%(fa, newName)
with open(os.path.join(output_folder, newName+".fa"), 'w') as oh:
fa_rec = SeqIO.read(os.path.join(input_folder, fa), 'fasta')
#fa_rec.name = newName
#SeqIO.write(fa_rec, oh, 'fasta')
oh.write(">"+newName+"\n")
oh.write(str(fa_rec.seq)+"\n")
else:
print "[renameFastaFile]: Warining, no new name was found for %s, no change"%fa
subprocess.check_call(['mkdir', '-p', os.path.join(output_folder, "nochange")])
subprocess.check_call(['cp', os.path.join(input_folder, fa),
os.path.join(output_folder, "nochange", fa)
])
def get_fasta_from_hmmsearch_result(hmmsearch_domtblout, query_fasta, output_folder):
"""given a query fasta, and hmmsearch_domtblout, find the best hit for each
marker, and write it out to output_folder with marker name as prefix
"""
marker2queries = {} # {marker:[query1, query2]}
query_marker_evalue = {} # {query: {marker: evalue}}
query_match = {} # {marker:{query1:query_match_len}}
marker_match = {} # {marker:{query1:marker_match_len}}
len_dict = {} # {query: query_length, marker:marker_length}
with open(hmmsearch_domtblout, "r") as ih:
for line in ih:
#print line
if line.startswith("#"):
continue
line = line.strip().split()
query = line[0]
query_len = int(line[2])
marker = line[3]
marker_len = int(line[5])
full_evalue = float(line[6])
#dom_evalue = float(line[11])
marker_match_len = int(line[16]) - int(line[15])
query_match_len = int(line[18]) - int(line[17])
len_dict[query] = query_len
len_dict[marker] = marker_len
if marker not in query_match:
query_match[marker] = {}
query_match[marker][query] = query_match_len
else:
if not query in query_match[marker]:
query_match[marker][query] = query_match_len
else:
query_match[marker][query] += query_match_len
if marker not in marker_match:
marker_match[marker] = {}
marker_match[marker][query] = marker_match_len
else:
if not query in marker_match[marker]:
marker_match[marker][query] = marker_match_len
else:
marker_match[marker][query] += marker_match_len
if not query in query_marker_evalue:
query_marker_evalue[query] = {}
query_marker_evalue[query][marker] = full_evalue
else:
if marker in query_marker_evalue[query]:
if full_evalue < query_marker_evalue[query]:
query_marker_evalue[query][marker] = full_evalue
if marker not in marker2queries:
marker2queries[marker] = [query]
else:
marker2queries[marker].append(query)
marker2bestquery = {} # {marker: best_query}
for marker, queries in marker2queries.items():
#print marker, queries
for query in queries:
query_len_pct = query_match[marker][query] / float(len_dict[query])
marker_len_pct = marker_match[marker][query] / float(len_dict[marker])
if (query_len_pct < 0.7) and (marker_len_pct < 0.7):
continue
else:
if marker not in marker2bestquery:
marker2bestquery[marker] = query
else:
old_query = marker2bestquery[marker]
if query_marker_evalue[query] < query_marker_evalue[old_query]:
marker2bestquery[marker] = query
# query2fasta
query2fasta = {}
for query_rec in SeqIO.parse(query_fasta, 'fasta'):
query = query_rec.name
query2fasta[query] = query_rec.seq
for marker, query in marker2bestquery.items():
#print marker, query
with open(os.path.join(output_folder, marker+".faa"), "w") as oh:
name = query.rsplit('_', 1)[0]
oh.write(">"+name+"\n")
oh.write(str(query2fasta[query])+"\n")
def AMPHORE2HmmMarkerScanner(args):
""" a simplified python implementation of MarkerScanner.pl of AMPHORE2, given
a folder with fasta sequences inside, a hmm file with hmm models inside,
search markers within input fasta sequences
"""
fasta_folder = args.inputFastaFolder
input_marker = args.inputMarkerFile
force = args.force
evalue = args.evalue
output_folder = args.outdir
subprocess.check_call(['mkdir', '-p', output_folder])
for fa in os.listdir(fasta_folder):
if not (fa.endswith('.fa') or fa.endswith('.fasta') or fa.endswith('.fas')):
print "[AMPHORE2HmmMarkerScanner]: %s is not a fasta file, pass"%fa
continue
filestem = os.path.splitext(fa)[0]
curr_outfolder = os.path.join(output_folder, filestem)
subprocess.check_call(['mkdir', '-p', curr_outfolder])
input_fa = os.path.join(fasta_folder, fa)
out_orf = os.path.join(curr_outfolder, filestem+".orf")
out_hmm = os.path.join(curr_outfolder, filestem+".hmmsearch")
# get orf
print "[AMPHORE2HmmMarkerScanner]: getting orfs of %s ..."%fa
if os.path.exists(out_orf):
if not force:
print "[AMPHORE2HmmMarkerScanner]: orf exists, nothing to do, \
use --force to override"
continue
subprocess.check_call(['getorf',
'-sequence', input_fa,
'-outseq', out_orf,
'-table', '1',
'-minsize', '100'
])
# run hmmsearch
print "[AMPHORE2HmmMarkerScanner]: hmmsearch of markers for %s ..."%fa
subprocess.check_call(['hmmsearch',
'-Z', '5000',
'--cpu', '40',
'-E', str(evalue),
'--domE', str(evalue),
'--domtblout', out_hmm,
'-o', '/dev/null',
input_marker,
out_orf
])
# get protein fasta belonging to each marker
print "[AMPHORE2HmmMarkerScanner]: write marker proteins for %s ..."%fa
get_fasta_from_hmmsearch_result(out_hmm, out_orf, curr_outfolder)
def rnammerMarkerScanner(args):
""" a wrapper function to call rnammer for rRNA scanning, for each genome
sequence in the input folder, it will scan the genome and found all rRNAs
"""
rnammer = args.rnammer
fasta_folder = args.inputFastaFolder
kingdom = args.superkingdom
type = args.molecule
force = args.force
output_folder = args.outdir
prefix = args.prefix
subprocess.check_call(['mkdir', '-p', output_folder])
for fa in os.listdir(fasta_folder):
if not (fa.endswith('.fa') or fa.endswith('.fasta') or fa.endswith('.fas')):
print "[rnammerMarkerScanner]: %s is not a fasta file, pass"%fa
continue
if prefix:
filestem = prefix
else:
filestem = os.path.splitext(fa)[0]
curr_outfolder = os.path.join(output_folder, filestem)
subprocess.check_call(['mkdir', '-p', curr_outfolder])
input_fa = os.path.join(fasta_folder, fa)
out_rRNA = os.path.join(curr_outfolder, filestem+".fa")
out_gff = os.path.join(curr_outfolder, filestem+".gff")
out_hmm = os.path.join(curr_outfolder, filestem+".hmmsearch")
# get orf
print "[rnammerMarkerScanner]: getting rRNAs for %s ..."%fa
if os.path.exists(out_rRNA):
if not force:
print "[rnammerMarkerScanner]: rRNA file exists, nothing to do, \
use --force to override"
continue
else:
print "[rnammerMarkerScanner]: rRNA file for %s was override"%fa
subprocess.check_call([rnammer,
'-S', kingdom,
'-m', type,
'-f', out_rRNA,
'-gff', out_gff,
'-h', out_hmm,
input_fa
])
def mergeEachMarkerProteinSet(args):
""" given a list of marker names, concate all found marker proteins in each
genome.
"""
input_folder = args.inputUmbrellaFolder
marker_list = args.inputMarkerList
output_folder = args.outdir
subprocess.check_call(['mkdir', '-p', output_folder])
# read in marker names
markers = []
with open(marker_list, "r") as ih:
for line in ih:
if line.startswith('#'):
continue
line = line.strip()
markers.append(line)
# concate marker proteins
for marker in markers:
with open(os.path.join(output_folder, marker+".mfa"), "w") as oh:
for genome in os.listdir(input_folder):
for f in os.listdir(os.path.join(input_folder, genome)):
if marker in f and (f.endswith(".faa") or f.endswith(".fa") or f.endswith(".fasta")):
fasta = SeqIO.read(os.path.join(input_folder, genome, f), 'fasta')
oh.write(">"+genome+"\n")
oh.write(str(fasta.seq)+"\n")
def mergerRNASet(args):
""" given the name of marker rRNA, concate this marker rRNA for all genomes.
"""
input_folder = args.inputUmbrellaFolder
marker = args.molecule+"_rRNA"
strategy = args.strategy # all or best
output_folder = args.outdir
length_cutoff = args.length_cutoff
numberOfNs = args.numberOfNs
subprocess.check_call(['mkdir', '-p', output_folder])
# merge marker rRNAs
with open(os.path.join(output_folder, marker+"_"+strategy+".fa"), "w") as oh:
for genome_folder in os.listdir(input_folder):
if not os.path.isdir(os.path.join(input_folder, genome_folder)):
continue
for f in os.listdir(os.path.join(input_folder, genome_folder)):
if f.endswith(".fa") or f.endswith(".fasta") or f.endswith(".fas"):
genome_name = genome_folder
marker_seq = os.path.join(input_folder, genome_folder, f)
if strategy == 'best':
best = None
best_score = 0
for rec in SeqIO.parse(marker_seq, "fasta"):
if marker in rec.description:
if len(rec.seq) < length_cutoff:
continue
if rec.seq.count('N') + rec.seq.count('n') >= numberOfNs:
continue
if not best:
best = rec
best_score = float(rec.description.split("=")[-1])
else:
current = rec
current_score = float(rec.description.split("=")[-1])
if current_score > best_score:
best_score = current_score
best = current
if best:
if best.seq.count('N') >=50:
continue
oh.write(">"+genome_name+"\n")
oh.write(str(best.seq)+"\n")
else:
count = 0
for rec in SeqIO.parse(marker_seq, 'fasta'):
seq_set = set()
if marker in rec.description:
if rec.seq.count('N') + rec.seq.count('n') >= numberOfNs:
continue
if len(rec.seq) < length_cutoff:
continue
if str(rec.seq) in seq_set:
continue
count += 1
seq_set.add(str(rec.seq))
oh.write(">"+genome_name+"_"+str(count)+"\n")
oh.write(str(rec.seq)+"\n")
def runSSUAlignAndMask(args):
""" run ssu-align and ssu-mask for input sequences
"""
ssu_align = args.ssu_align
ssu_mask = args.ssu_mask
input_fasta = args.inputFastaFile
force = args.force
prefix = args.prefix
kingdom = args.kingdom
no_align = args.no_align
no_search = args.no_search
output_RNA = args.output_RNA
output_stk = args.output_stk
if prefix:
output_folder = prefix+"_ssu-align"
else:
basename = os.path.basename(input_fasta)
filestem = os.path.splitext(basename)[0]
output_folder = filestem + "_ssu-align"
# do ssu alignment
cmd = [ssu_align, '-n', kingdom]
if force:
cmd.append('-f')
if not output_RNA:
cmd.append('--dna')
if no_align:
cmd.append('--no-align')
if no_search:
cmd.append('--no-search')
cmd.append(input_fasta)
cmd.append(output_folder)
try:
p = subprocess.Popen(cmd)
p.wait()
except Exception as e:
print "Error raised at running runSSUAlignAndMask:%s"%e
sys.exit(1)
# do mask
cmd=[ssu_mask]
if not output_RNA:
cmd.append('--dna')
if not output_stk:
cmd.append('--afa')
cmd.append(output_folder)
try:
p = subprocess.Popen(cmd)
p.wait()
except Exception as e:
err = "[runSSUAlignAndMask]: Error raised at running runSSUAlignAndMask:%s"%e
raise SubCommandError(err)
def runMSAonMarkers(args):
""" given an input folder contains marker protein sequences, run multiple sequence alignment
on each marker
"""
input_folder = args.inputMarkerProteinsFolder
suffix = args.suffix
if args.outdir:
output_folder = args.outdir
else:
output_folder = input_folder
for f in os.listdir(input_folder):
if f.endswith(suffix):
filestem = os.path.splitext(f)[0]
# run msa
try:
align = subprocess.Popen([ 'muscle',
'-in', os.path.join(input_folder, f),
'-out', os.path.join(output_folder, filestem+".msf")
])
align.wait()
trim = subprocess.Popen(['trimal',
'-in', os.path.join(output_folder, filestem+".msf"),
'-out', os.path.join(output_folder, filestem+"_trimAl.msf"),
'-gt', '0.9',
'-cons', '60',
'-w', '3'
])
trim.wait()
except Exception as e:
err = "[runMSAonMarkers]: Error raised at runMSAonMarkers: %s"%e
raise SubCommandError(err)
def concatMarkersforEachGenome(args):
"""given a folder with marker protein alignments inside, concatenate markers
for each genome, gaps will be placed if no marker found for one genome.
"""
input_folder = args.inputMarkerAlignmentFolder
suffix = args.suffix
if not args.outdir:
output_folder = input_folder
else:
output_folder = args.outdir
# summarize all markers, marker length and all genome names
genome2markers = {} # {genome:{marker1:seq1, marker2:seq2, ...}}
marker2length = {} # {marker1:length1, marker2:length2}
for f in os.listdir(input_folder):
if f.endswith(suffix):
marker = f.split("_")[0]
aln_records = SeqIO.parse(os.path.join(input_folder, f), 'fasta')
for aln in aln_records:
genome = aln.name
seq = str(aln.seq).strip()
marker_len = len(seq)
if marker not in marker2length:
marker2length[marker] = marker_len
else:
assert marker_len == marker2length[marker], "Error: Are you\
sure the input file is an alignment? The length is different!"
if genome not in genome2markers:
genome2markers[genome] = {marker:seq}
else:
genome2markers[genome].update({marker:seq})
# write out concated markers for each genome
all_genomes = sorted(genome2markers.keys())
all_markers = sorted(marker2length.keys())
with open(os.path.join(output_folder, "concatMarkersforEachGenome_supermatrix.msf"), 'w') as oh:
for genome in all_genomes:
header = ">"+genome+"\n"
seq = ""
for marker in all_markers:
marker_seq = genome2markers[genome].get(marker, None)
if not marker_seq:
marker_seq = "-"*marker2length[marker]
seq += marker_seq.strip()
oh.write(header)
oh.write(seq+"\n")
# write out marker partitions
with open(os.path.join(output_folder, "concatMarkersforEachGenome_supermatrix.coord"), 'w') as coord:
accum_start = 0
for marker in all_markers:
start = accum_start + 1
stop = start + marker2length[marker]-1
coord.write("%s = %d-%d\n"%(marker, start, stop))
accum_start = stop
def convertFasta2Phylip(input_fasta, output_prefix, format='phylip-sequential'):
""" convert input_fasta to phylip-relaxed format using Biopython
"""
AlignIO.convert(input_fasta, 'fasta', output_prefix+".phy", format)
def convertMSFtoPhylip(args):
"""convert multiple sequence alignments in multiple fasta format to phylip-sequential format
"""
input_folder = args.inputSuperAlignmentFolder
suffix = args.suffix
phylip_relaxed = args.phylip_relaxed
if args.outdir:
output_folder = args.outdir
else:
output_folder = input_folder
outfmt = "phylip-sequential"
if phylip_relaxed:
outfmt = 'phylip-relaxed'
# convert to phylip-relaxed format
for f in os.listdir(input_folder):
if f.endswith(suffix):
input_fasta = os.path.join(input_folder, f)
basename = os.path.splitext(f)[0]
output_prefix = os.path.join(output_folder, basename)
print "[convertMSFtoPhylip]: converting %s to %s format ..."%(f, outfmt)
# it's compulsory to use phylip-sequential format for PartitionFinderProtein
convertFasta2Phylip(input_fasta, output_prefix, format=outfmt)
def runPartitionFinderProtein(args):
""" run PartitionFinderProtein on concatenated super alignment, find the best
protein models and partitions.
"""
alignment = args.inputPhylipSuperAlign
coordinate = args.inputCoordFile
workingdir = args.newWorkingFolder
program = args.PartitionFinderProtein
partition_finder_template = """## ALIGNMENT FILE ##
alignment = %s;
## BRANCHLENGTHS: linked | unlinked ##
branchlengths = linked;
## MODELS OF EVOLUTION for PartitionFinder: all | raxml | mrbayes | <list> ##
## for PartitionFinderProtein: all_protein | <list> ##
models = all_protein;
# MODEL SELECCTION: AIC | AICc | BIC #
model_selection = BIC;
## DATA BLOCKS: see manual for how to define ##
[data_blocks]
%s
## SCHEMES, search: all | user | greedy ##
[schemes]
search = greedy;
#user schemes go here if search=user. See manual for how to define.#
"""
# copy alignment to workingdir
subprocess.check_call(['mkdir', '-p', workingdir])