-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathreprophylo.py
6263 lines (5319 loc) · 275 KB
/
reprophylo.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
reprophyloversion=1.0
############################################################################################
if False:
"""
ReproPhylo version 1
General purpose phylogenetics package for reproducible and experimental analysis
Amir Szitenebrg
David H Lunt
EvoHull.org
University of Hull
Developed with:
CPython 2.7.6
IPython 1.2.1
ete2 2.2rev1056
biopython 1.64
dendropy 3.12.0
cloud 2.8.5
numpy 1.8.2
matplotlib 1.3.1
pandas
RAxML 8
Phylobayes 3
Trimal 1
Muscle
Mafft 7
Pal2nal 14
"""
##############################################################################################
from Bio import SeqIO
import os, csv, sys, dendropy, re, time, random, glob, platform, warnings, rpgit, ast, gb_syn,css
import HTML, inspect, shutil
import subprocess as sub
from Bio.Seq import Seq
import numpy as np
import matplotlib.pyplot as plt
from Bio.Alphabet import IUPAC
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation, CompoundLocation
from Bio.Align.Applications import MafftCommandline, MuscleCommandline
from StringIO import StringIO
from Bio import AlignIO
from Bio.Phylo.Applications import RaxmlCommandline
from Bio.Align import MultipleSeqAlignment
from Bio.SeqUtils import GC
from ete2 import *
from collections import Counter
#import pandas as pd
import math
import __builtin__
##############################################################################################
class Locus:
##############################################################################################
""" Configure the loci stored in the ReproPhylo Project.
>>> locus = Locus('dna', 'CDS', 'coi', ['cox1','COX1','coi','COI','CoI'])
>>> print(locus)
Locus(char_type=dna, feature_type=CDS, name=coi, aliases=cox1; COX1; coi; COI; CoI)
"""
char_type = 'NotSet'
feature_type = 'NotSet'
name = 'NotSet'
aliases = []
def __init__(self, char_type=char_type, feature_type=feature_type,
name=name, aliases=aliases):
self.char_type = char_type
self.feature_type = feature_type
self.name = name
self.aliases = aliases
valid = ['dna','prot']
if not self.char_type in valid:
raise ValueError('self.char_type should be \'dna\' or \'prot\'')
if not type(self.feature_type) is str:
raise ValueError('self.feature_type should be a string')
if not type(self.name) is str:
raise ValueError('self.name should be a string')
if not type(self.aliases) is list:
raise ValueError('self.aliases should be a list')
else:
for a in self.aliases:
if not type(a) is str:
raise ValueError('aliases in self.aliases have to be strings')
def __str__(self):
aliases_str = ('; ').join(self.aliases)
return ('Locus(char_type='+self.char_type+', feature_type='+self.feature_type+
', name='+self.name+', aliases='+aliases_str+')')
##############################################################################################
class Concatenation:
##############################################################################################
"""This class is used to configure concatenations given loci and rules.
>>> coi = Locus('dna', 'CDS', 'coi', ['cox1','COX1','coi','COI','CoI'])
>>> ssu = Locus('dna', 'rRNA', '18S', ['18S rRNA','SSU rRNA'])
>>> bssu = Locus('dna', 'rRNA', '16S', ['16S rRNA'])
>>> lsu = Locus('dna', 'rRNA', '28S', ['28S rRNA', 'LSU rRNA'])
>>> alg11 = Locus('dna', 'CDS', 'ALG11', ['ALG11'])
>>> loci = [coi, ssu, bssu, lsu, alg11]
>>> concatenation = Concatenation(name='combined', loci=loci,
... otu_meta='OTU_name',
... otu_must_have_all_of=['coi'],
... otu_must_have_one_of =[['16S','28S'],['ALG11','18S']],
... define_trimmed_alns=["MuscleDefaults@dummyTrimMethod"])
>>> print(str(concatenation))
Concatenation named combined, with loci coi,18S,16S,28S,ALG11,
of which coi must exist for all species
and at least one of each group of [ 16S 28S ][ ALG11 18S ] is represented.
Alignments with the following names: MuscleDefaults@dummyTrimMethod are prefered
"""
otu_must_have_all_of = []
otu_must_have_one_of = 'any'
define_trimmed_alns = [] #should be Locus_name@Alignment_method_name@Trimming_mathod_name
feature_id_dict = {}
def __init__(self,
name,
loci,
otu_meta,
otu_must_have_all_of = otu_must_have_all_of,
otu_must_have_one_of = otu_must_have_one_of,
define_trimmed_alns = define_trimmed_alns):
self.name = name
self.loci = loci
self.otu_meta = otu_meta
self.otu_must_have_all_of = otu_must_have_all_of
self.otu_must_have_one_of = otu_must_have_one_of
if isinstance(otu_must_have_all_of,str):
raise IOError('The keyword \'otu_must_have_all_of\' has to be a list')
if isinstance(otu_must_have_one_of[0],str) and not otu_must_have_one_of == 'any':
raise IOError('The keyword \'otu_must_have_one_of\' has to be a list of lists')
if self.otu_must_have_one_of == 'any':
self.otu_must_have_one_of = [[l.name for l in self.loci]]
self.feature_id_dict = {} # Will hold the feature_id list for each otu
self.define_trimmed_alns = define_trimmed_alns # To choose between alternative
# alignments of the same locus
self.used_trimmed_alns = {} #To hold the alignment chosen for each locus
# Validate loci list
seen = []
for locus in loci:
if not isinstance(locus, Locus):
raise TypeError("Expecting Locus object in loci list")
if locus.name in seen:
raise NameError('Locus ' + locus.name + ' appears more than once in self.loci')
else:
seen.append(locus.name)
def __str__(self):
loci_names = [i.name for i in self.loci]
loci_string = ''
for l in loci_names:
loci_string += l+','
loci_string = loci_string[:-1]
must_have = ''
for i in self.otu_must_have_all_of:
must_have += i+','
must_have = must_have[:-1]
trimmed_alignmnets_spec = ''
one_of = ''
for i in self.otu_must_have_one_of:
one_of += '[ '
for j in i:
one_of += j+' '
one_of += ']'
if (self.define_trimmed_alns) > 0:
for i in self.define_trimmed_alns:
trimmed_alignmnets_spec += i
return ("Concatenation named %s, with loci %s,\n"
"of which %s must exist for all species\n"
"and at least one of each group of %s is represented.\n"
"Alignments with the following names: %s are prefered"
% (self.name, loci_string, must_have, one_of, trimmed_alignmnets_spec))
##############################################################################################
if False:
"""
Reprophylo Project Utilities
Used in the Project class but are not in the classe's methods
"""
##############################################################################################
## Git management
__builtin__.git = False
# git log template
gitline = "<<<<\n%s\nSTDOUT:\n%s\nSTDERR:%s\n>>>>\n"
def undate_git_log(pj, out, err):
if not err:
err = 'None'
if not out:
out = 'None'
pj.git_log += gitline%(str(time.asctime()),str(out), str(err))
def start_git(pj):
__builtin__.git = True # flag it on
cwd = os.getcwd()
if os.path.isdir(cwd + '/.git'):
# a repo exists, check it belongs to this project by checking the description
try:
assert open(cwd + '/.git/description').read().strip().rstrip() == pj.pickle_name.strip().rstrip()
warnings.warn('Git repository exists for this Project')
except:
raise RuntimeError('The Git repository in the CWD does not belong to this project. Either the pickle'+
' moved, or this is a preexsisting repo. Try one of the following: Delete the local '+
' .Git dir if you don\'t need it, move the pickle and the notebook to a new work dir,'+
' or if possible, move them back to their original location. You may also disable Git'+
' by with stop_git().')
else:
# start a rep
out, err = rpgit.gitInit()
undate_git_log(pj, out, err)
# write the pickle name as the repo description
hndl = open(cwd + '/.git/description', 'wt')
hndl.write(pj.pickle_name.strip().rstrip())
hndl.close()
warnings.warn('The new repository is called %s.'%open(cwd + '/.git/description', 'r').read().rstrip())
# list scripts and notebooks
import fnmatch
matches = []
for root, dirnames, filenames in os.walk(cwd):
for filename in fnmatch.filter(filenames, '*.py'):
matches.append(os.path.join(root, filename))
for filename in fnmatch.filter(filenames, '*.ipynb'):
matches.append(os.path.join(root, filename))
# git add scripts and notebooks
for match in matches:
out, err = rpgit.gitAdd(match)
undate_git_log(pj, out, err)
# commit scripts and notebooks
comment = "%i script file(s) from %s" % (len(matches), time.asctime())
out, err = rpgit.gitCommit(comment)
undate_git_log(pj, out, err)
def stop_git():
__builtin__.git = False # flad it off
cwd = os.getcwd()
# list, git add and git commit scripts and notebooks
# list
import fnmatch
matches = []
for root, dirnames, filenames in os.walk(cwd):
for filename in fnmatch.filter(filenames, '*.py'):
matches.append(os.path.join(root, filename))
for filename in fnmatch.filter(filenames, '*.ipynb'):
matches.append(os.path.join(root, filename))
# add
for match in matches:
out, err = rpgit.gitAdd(match)
undate_git_log(pj, out, err)
# commit
comment = "%i script file(s) from %s" % (len(matches), time.asctime())
out, err = rpgit.gitCommit(comment)
undate_git_log(pj, out, err)
## end git management
def platform_report():
"""
Prints machine specs, os specs and dependencies at time of execution
>>> isinstance(platform_report(), list)
True
"""
import pkg_resources
modules = [] # and their versions
for i in ('ete2','biopython','dendropy','cloud'):
try:
modules.append(i+' version: '+
pkg_resources.get_distribution(i).version)
except:
pass
modules.append('reprophylo version %s'%str(reprophyloversion))
return(['Platform: '+platform.platform(aliased=0, terse=0),
'Processor: '+platform.processor(),
'Python build: '+platform.python_build()[0] + platform.python_build()[1],
'Python compiler: '+platform.python_compiler(),
'Python implementation: ' +platform.python_implementation(),
'Python version: ' + platform.python_version()]+
modules+
['User: ' +platform.uname()[1]])
def write_alns(pj, format = 'fasta'):
"""
Writes untrimmed sequence alignment files that are in pj in a biopython format
"""
if len(pj.alignments.keys()) == 0:
raise IOError('Align the records first')
else:
for key in pj.alignments:
AlignIO.write(pj.alignments[key], key+'_aln.'+format, format)
def keep_feature(feature, loci):
""" Returns true if a feature's type is in one of the loci and if the gene
or product qualifiers is in the aliases of one of the loci, for data collection
from a genbank or embl file
# making a dummy feature
>>> coi = Locus('dna','CDS','coi', ['cox1','COX1','coi','COI','CoI'])
>>> location = FeatureLocation(1,100)
>>> feature = SeqFeature()
>>> feature.location = location
>>> feature.type = 'CDS'
>>> feature.qualifiers['gene'] = ['CoI']
# testing if fits any of the Project Locus objects
>>> a = keep_feature(feature, [coi])
>>> print(a)
True"""
keep = 0
for g in loci:
if not g.name in g.aliases:
g.aliases.append(g.name)
if feature.type == 'source':
keep = 1
elif feature.type == g.feature_type:
qual = None
if 'gene' in feature.qualifiers.keys():
qual = 'gene'
elif 'product' in feature.qualifiers.keys():
qual = 'product'
if qual and feature.qualifiers[qual][0] in g.aliases:
keep = 1
if keep == 1:
return True
else:
return False
return
def dwindle_record(record, loci):
"""
Retains only features that are called by Locus objects and records with features that are
called by Locus objects
# Making a dummy locus
>>> coi = Locus('dna','CDS','coi', ['cox1','COX1','coi','COI','CoI'])
# Making a dummy record with a feature that fits a Locus object (kept_feature)
# and a feature that does not (dwindled_feature)
>>> location = FeatureLocation(1,100)
>>> kept_feature = SeqFeature()
>>> kept_feature.location = location
>>> kept_feature.type = 'CDS'
>>> kept_feature.qualifiers['gene'] = ['CoI']
>>> dwindled_feature = SeqFeature()
>>> dwindled_feature.location = location
>>> dwindled_feature.type = 'rRNA'
>>> dwindled_feature.qualifiers['gene'] = ['LSU']
>>> s = 'atgc'*1000
>>> record = SeqRecord(seq=Seq(s, IUPAC.ambiguous_dna), id='1', description='spam')
>>> record.features.append(kept_feature)
>>> record.features.append(dwindled_feature)
>>> print(len(record.features))
2
# Dwindling the record
>>> a = dwindle_record(record, [coi])
>>> print(len(record.features))
1
"""
dwindled_features = []
feature_count = 0
for feature in record.features:
if keep_feature(feature, loci)== True:
# determine feature id
if feature.type == 'source' and not 'feature_id' in feature.qualifiers.keys():
feature.qualifiers['feature_id'] = [record.id + '_source']
elif not 'feature_id' in feature.qualifiers.keys():
feature.qualifiers['feature_id'] = [record.id + '_f' + str(feature_count)]
feature_count += 1
# determine prop ambiguity and GC content
if not feature.type == 'source':
feature_seq = feature.extract(record.seq)
degen = len(feature_seq)
for i in ['A','T','G','C','U','a','t','g','c','u']:
degen -= feature_seq.count(i)
feature.qualifiers['GC_content'] = [str(GC(feature_seq))]
feature.qualifiers['nuc_degen_prop'] = [str(float(degen)/len(feature_seq))]
if 'translation' in feature.qualifiers.keys():
transl = feature.qualifiers['translation'][0]
degen = 0
for i in ['B', 'X', 'Z', 'b', 'x', 'z']:
degen += transl.count(i)
feature.qualifiers['prot_degen_prop'] = [str(float(degen)/len(transl))]
dwindled_features.append(feature)
record.features = dwindled_features
return record
def is_embl_or_gb(input_filename):
suffixes = ['.gb','.embl']
gb = False
for s in suffixes:
if s in input_filename:
gb = True
return gb
def parse_input(input_filename, fmt):
return SeqIO.parse(input_filename, fmt)
def list_to_string(List):
"""
Handles list printing as a nice string in the pj.write(format="csv") method
>>> L = ['a','b','b']
>>> print(list_to_string(L))
a;b;b
"""
string = ''
for i in List:
if type(i) is str and '\n' in i:
string += lines_to_line(i).rstrip()+';'
else:
string += str(i)+';'
return string[:-1]
def lines_to_line(lines):
"""
Replaces newline with space in the pj.write(format="csv") method
"""
lines = lines.split('\n')
return (' ').join(lines)
def type_to_single_line_str(var):
"""
Returns any type as a one line string for the pj.write(format="csv") method
"""
if type(var) is str and '\n' in var:
return lines_to_line(var)
elif type(var) is str or type(var) is int or type(var) is float:
return str(var)
elif type(var) is list and len(var) == 1:
return str(var[0])
elif type(var) is list and len(var) > 0:
return list_to_string(var)
else:
return var
def get_qualifiers_dictionary(project, feature_id):
"""
Takes sequence record annotation, source qualifiers and feature qualifiers and puts them
in a flat dictionary
This is being replaced by __get_qualifiers_dictionary__ which deals with the records as a dict
and is much faster. Eventually, records will be handled as a dict throughout, instead of as
a list.
# Making a dummy locus
>>> coi = Locus('dna','CDS','coi', ['cox1','COX1','coi','COI','CoI'])
# Making a dummy Project
>>> pj = Project([coi], git=False)
# making a dummy record
>>> s = 'atgc'*1000
>>> location = FeatureLocation(1,100)
>>> feature = SeqFeature()
>>> feature.location = location
>>> feature.type = 'CDS'
>>> feature.qualifiers['gene'] = ['CoI']
>>> feature.qualifiers['feature_id'] = ['1_f0']
>>> source = SeqFeature()
>>> source.location = FeatureLocation(0,3999)
>>> source.type = 'source'
>>> source.qualifiers['organism'] = ['Tetillda radiata']
>>> record = SeqRecord(seq=Seq(s, IUPAC.ambiguous_dna), id='1', description='spam')
>>> record.features.append(feature)
>>> record.features.append(source)
>>> record.annotations["evidence"] = 'made up'
>>> pj.records = [record]
# executing get_qualifiers_dictionary()
>>> qual_dict = get_qualifiers_dictionary(pj, '1_f0')
>>> qual_items = qual_dict.items()
>>> qual_items.sort(key = lambda i: i[0])
>>> for key, val in qual_items: print(key.ljust(20,' ') + val.ljust(20,' '))
annotation_evidence made up
feature_id 1_f0
gene CoI
record_id 1
source_organism Tetillda radiata
"""
if type(feature_id) is list and len(feature_id) > 1:
raise IOError('get_qualifiers_dictionary takes one feature_id at a time')
if type(feature_id) is list:
feature_id = feature_id[0]
record_id = feature_id.rpartition('_')[0]
qualifiers_dictionary={'record_id': record_id}
for record in project.records:
if record.id in feature_id:
for annotation in record.annotations.keys():
qualifiers_dictionary['annotation_'+annotation]=record.annotations[annotation]
for feature in record.features:
if feature.type == 'source':
for qualifier in feature.qualifiers.keys():
qualifiers_dictionary['source_'+qualifier]=feature.qualifiers[qualifier][0]
elif feature.qualifiers['feature_id'][0] == feature_id:
for qualifier in feature.qualifiers.keys():
qualifiers_dictionary[qualifier]=feature.qualifiers[qualifier][0]
return qualifiers_dictionary
def __get_qualifiers_dictionary__(project, feature_id):
"""
This will replace the public version. It uses the Project._records_dict to pull
the record using the record id, instead of iterating Project.records, which is very slow.
It requires Project.__records_list_to_dict__() to execute beforehand.
"""
if type(feature_id) is list and len(feature_id) > 1:
raise IOError('get_qualifiers_dictionary takes one feature_id at a time')
if type(feature_id) is list:
feature_id = feature_id[0]
record_id = feature_id.rpartition('_')[0]
record = project._records_dict[record_id]
qualifiers_dictionary={'record_id': record_id}
for annotation in record.annotations.keys():
qualifiers_dictionary['annotation_'+annotation]=record.annotations[annotation]
for feature in record.features:
if feature.type == 'source':
for qualifier in feature.qualifiers.keys():
qualifiers_dictionary['source_'+qualifier]=feature.qualifiers[qualifier][0]
elif feature.qualifiers['feature_id'][0] == feature_id:
for qualifier in feature.qualifiers.keys():
qualifiers_dictionary[qualifier]=feature.qualifiers[qualifier][0]
return qualifiers_dictionary
def seq_format_from_suffix(suffix):
"""
Guesses input format from suffix
>>> print(seq_format_from_suffix('gb'))
genbank
"""
suffixes = {'fasta': ['fas','fasta','fa','fna'],
'genbank': ['gb','genbank'],
'embl': ['embl']}
found = False
for key in suffixes.keys():
if suffix in suffixes[key]:
found = True
return key
if not found:
raise RuntimeError(suffix+' is not a recognised suffix of an unaligned sequence file')
def read_feature_quals_from_tab_csv(csv_filename):
"""
This is used to update feature qualifiers from a tab delimited file
"""
import re
header = open(csv_filename, 'r').readlines()[0].rstrip().split('\t')
feature_id_col = header.index('feature_id')
taxonomy_col = header.index('taxonomy')
seq_col = header.index('seq')
translation_col = None
if 'translation' in header:
translation_col = header.index('translation')
csv_info = {}
for line in [l.rstrip().split('\t') for l in open(csv_filename, 'r').readlines()[1:]]:
if not line[0] in csv_info.keys():
csv_info[line[0]] = {'source':{},
'taxonomy':[],
'features':{}
}
if csv_info[line[0]]['taxonomy'] == []:
csv_info[line[0]]['taxonomy'] = line[taxonomy_col].split(';')
csv_info[line[0]]['features'][line[feature_id_col]] = {}
get_source = False
if csv_info[line[0]]['source'] == {}:
get_source = True
for i in range(len(header)):
if get_source and 'source:_' in header[i]:
qual_name = re.sub('source:_','',header[i])
if not line[i] == 'null' and not line[i] == '' and line[i]:
csv_info[line[0]]['source'][qual_name] = line[i].split(';')
elif (not 'source:_' in header[i] and not line[i] == 'null' and not line[i] == '' and line[i] and
not i in [seq_col, translation_col, taxonomy_col, feature_id_col]):
csv_info[line[0]]['features'][line[feature_id_col]][header[i]] = line[i].split(';')
return csv_info
## Alignment statistics
def count_positions(aln_column):
counts = {}
for i in aln_column:
if i in counts.keys():
counts[i] += 1
else:
counts[i] = 1
return counts
def global_aln_stats(aln_obj):
total_gaps = 0
prop_list = []
non_uniform_count = aln_obj.get_alignment_length()
parsimony_informative = 0
for i in range(aln_obj.get_alignment_length()):
total_gaps += aln_obj[:, i].count('-')
prop_list.append(aln_obj[:, i].count('-')/float(len(aln_obj)))
if len(count_positions(aln_obj[:, i]).keys()) == 1:
non_uniform_count -= 1
elif (len(count_positions(aln_obj[:, i]).keys()) == 2 and
'-' in count_positions(aln_obj[:, i]).keys()):
non_uniform_count -= 1
if len([p for p in count_positions(aln_obj[:, i]).keys() if (p != '-' and count_positions(aln_obj[:, i])[p] > 1)]) > 1:
parsimony_informative += 1
mean_gap_prop = sum(prop_list)/aln_obj.get_alignment_length()
return (mean_gap_prop, non_uniform_count, parsimony_informative)
def count_undetermined_lines(aln_obj, cutoff=0):
count = 0
ids = []
if aln_obj.get_alignment_length() < cutoff*2:
warnings.warn('The cutoff to exclude a sequence is more than half of the alignmnet length')
elif aln_obj.get_alignment_length() <= cutoff:
raise RuntimeWarning('The cutoff to exclude a sequence is as long or longer than the alignment')
for seq in aln_obj:
if str(seq.seq).count('-') >= aln_obj.get_alignment_length()-cutoff:
count += 1
ids.append(seq.id)
return count, ids
def count_collapsed_aln_seqs(aln_obj):
count = 1
seen_seqs = [str(aln_obj[0].seq)]
for seq in aln_obj[1:]:
str_seq = str(seq.seq)
if len([s for s in seen_seqs if (str_seq in s or s in str_seq or s == str_seq)]) == 0:
count += 1
seen_seqs.append(str_seq)
return count
def aln_summary(aln_obj, cutoff=0):
lines = ["Alignment length: %i" % aln_obj.get_alignment_length(),
"Number of rows: %i" % len(aln_obj),
"Unique sequences: %i"%count_collapsed_aln_seqs(aln_obj),
"Average gap prop.: %f\nVariable columns: %i\nParsimony informative: %i"
%global_aln_stats(aln_obj),
"Undetermined sequences: %i"%(count_undetermined_lines(aln_obj, cutoff=cutoff)[0]),
"Undetermined sequence cutoff: %i"%cutoff
]
return [lines, len(aln_obj), count_undetermined_lines(aln_obj, cutoff=cutoff), count_collapsed_aln_seqs(aln_obj)]
##
def loci_list_from_csv(loci):
"""
Parse the loci csv file given to Project
"""
# verify format
if any(len(line.split(',')) >= 4 for line in open(loci, 'r').readlines()):
pass
else:
raise IOError("File %s has no valid loci of format char_type,feature_type,name,aliases"%loci)
loci_dict = {}
loci_list = []
for line in [line.rstrip() for line in open(loci, 'r').readlines() if len(line.rstrip()) > 0]:
# verify format
if len(line.split(',')) < 4:
raise IOError("The line %s in file %s is missing arguments. Needs at least char_type,feature_type,name,aliases"%
(line.rstrip(), loci))
# look for synonyms
else:
group = None
try:
group = int(line.rstrip().split(',')[-1])
except:
pass
if group:
locus_exists = False
for name in loci_dict:
if 'group' in loci_dict[name].keys() and loci_dict[name]['group'] == group:
loci_dict[name]['aliases'] += line.split(',')[3:-1]
locus_exists = True
if not locus_exists:
loci_dict[line.split(',')[2]] = {'group': int(line.rstrip().split(',')[-1]),
'char_type': line.split(',')[0],
'feature_type': line.split(',')[1],
'aliases': line.split(',')[3:-1]
}
else:
loci_dict[line.split(',')[2]] = {'group': None,
'char_type': line.split(',')[0],
'feature_type': line.split(',')[1],
'aliases': line.split(',')[3:]
}
for name in loci_dict:
loci_list.append(Locus(loci_dict[name]['char_type'],
loci_dict[name]['feature_type'],
name,
loci_dict[name]['aliases']))
return loci_list
def parse_paup_charset(nexus_filename):
"""
Takes a nexus file with PAUP style charset commands.
Returns a dictionary with partition names as keys and a list of
integers representing the start and end of the partition as a value.
Position count starts from 0.
Handles paup commands of the following format:
CHARSET locus_name=1-129;
or
charset locus_name = 1 - 129 ;
"""
try:
AlignIO.read(nexus_filename, 'nexus')
except:
n = len(list(AlignIO.parse(nexus_filename, 'nexus')))
raise IOError('Cannot handle more then one matrix in %s. Got %i matrices'%
(nexus_filename, n))
charsets = {}
charset_lines = [l for l in open(nexus_filename,'r').readlines() if
(l.startswith('CHARSET') or l.startswith('charset'))]
if len(charset_lines) == 0:
raise IOError("There are no CHARSET commands in %s"%nexus_filename)
for line in charset_lines:
try:
info = line.split()[1].split(';')[0]
locus_name, range = info.split('=')
locus_name = locus_name.strip().rstrip()
start = int(range.split('-')[0].strip().rstrip())-1
end = int(range.split('-')[1].strip().rstrip())-1
charsets[locus_name] = [start,end]
except:
raise IOError('Expects "charset set_name = start_int - end_int;"'+
' (case insensitive, spaces around the "=" or "-" not mandatory). Got %s'%line)
return charsets
def pj_from_nexus_w_charset(nexus_filename, output_dir, char_type,
feature_type, project=False, pickle=False, git=False):
"""
Takes a nexus file with PAUP style charset commands as input.
Creates a separate fasta file for each alignment partition
Returns a list of fasta filenames and a list of Locus objects
If project==True, returns a Project instance with the loci, alignments and records instead
"""
from reprophylo import Locus
from Bio import AlignIO
charsets = parse_paup_charset(nexus_filename)
alignment = AlignIO.read(nexus_filename, 'nexus')
filenames = []
loci_list = []
for locus_name in charsets:
s = charsets[locus_name][0]
e = charsets[locus_name][1]
outname = "%s/%s.fasta"%(output_dir,locus_name)
AlignIO.write(alignment[:, s:e], outname, 'fasta')
filenames.append(outname)
loci_list.append(Locus(char_type, feature_type, locus_name, [locus_name]))
if project:
from reprophylo import Project
pj = Project(loci_list, pickle=pickle, git=git)
i=1
for f in filenames:
locus_name = f.split('/')[-1].split('.')[0]
print '%i/%i reading %s'%(i,len(filenames), locus_name)
i += 1
pj.read_alignment(f, char_type, feature_type, locus_name)
return pj
else:
return filenames, loci_list
##############################################################################################
class Project:
##############################################################################################
"""
The Project class contians all the data and has methods to analyze it. It allows for
experimental analysis by running alternative analyses and formally comparing the
outputs. The pickle_pj() function allows to pickle the project, including the data,
intermediates and results, as well as a description of the methods.It allows for a rerun
of the whole analysis as is, as well as for a reconfiguration of the analysis or addition
of data. If git is installed, it can be called by 'import rpgit'. As a result, git can be
initiated using start_git(). A git repository will be created in the CWD, if it doesn't already exist.
Input datasets, .py, .ipynb and .pkl files in the CWD will be version controlled.
Version control can be paused in the middle of the script
by calling stop_git() and restarted by calling start_git() again.
"""
def __init__(self, loci, pickle=None, git=True):
"""
# making dummy loci
>>> coi = Locus('dna','CDS','coi',['COX1','cox1'])
>>> ssu = Locus('dna','rRNA','18S',['18S','SSU'])
# Making a Project object
>>> pj = Project([coi,ssu], git=False)
>>> print(str(pj))
Project object with the loci coi,18S,
"""
self.records = []
self._records_dict = {}
self.starttime = str(time.asctime())
self.user = None
if os.path.isfile('USER'):
self.user = []
for line in open('USER','r').readlines():
key, arg = line.rstrip().split('=')
self.user.append([key, arg])
self.loci = loci
self.records_by_locus = {}
self.concatenations = []
self.alignments = {}
self.trimmed_alignments = {}
self.trees = {}
self.used_methods = {}
self.sets = {}
self.git_log = ''
self.pickle_name=pickle
if self.pickle_name and os.path.exists(self.pickle_name):
raise IOError('Pickle %s exists. If you want to keep using it do pj=unpickle_pj(\'%s\') instead.'%(self.pickle_name,self.pickle_name))
if git and not self.pickle_name:
raise IOError('Must have pickle to run Git. Either specify pickle or git=False')
elif git:
start_git(self)
self.defaults = {'raxmlHPC': programspath+'raxmlHPC-PTHREADS-SSE3',
'mafft': 'mafft',
'muscle': programspath+'muscle',
'trimal': programspath+'trimal',
'pb': programspath+'pb',
'bpcomp': programspath+'bpcomp',
'tracecomp': programspath+'tracecomp',
'fasttree': programspath+'FastTreeMP',
'pal2nal': programspath+'pal2nal.pl',
# PROGRAM PLUG
# 'program name': programpath+'the basic command'
}
seen = []
if isinstance(loci,list):
for locus in loci:
if not isinstance(locus, Locus):
raise TypeError("Expecting Locus object in loci list. "+locus+
" not a Locus object")
if locus.name in seen:
raise NameError('Locus ' + locus.name + ' apears more than once in self.loci')
else:
seen.append(locus.name)
elif isinstance(loci,str):
self.loci = loci_list_from_csv(loci)
#print 'Read the following loci from file %s:'%loci
#for l in self.loci:
#print str(l)
self.aln_summaries = []
if self.pickle_name:
pickle_pj(self, self.pickle_name, track=False)
if __builtin__.git and self.pickle_name:
import rpgit
comment = "%s from %s" % (str(self), time.asctime())
out, err = rpgit.gitAdd(self.pickle_name)
undate_git_log(self, out, err)
cwd = os.getcwd()
import fnmatch
matches = []
for root, dirnames, filenames in os.walk(cwd):
for filename in fnmatch.filter(filenames, '*.py'):
matches.append(os.path.join(root, filename))
for filename in fnmatch.filter(filenames, '*.ipynb'):
matches.append(os.path.join(root, filename))
for match in matches:
out, err = rpgit.gitAdd(match)
undate_git_log(self, out, err)
out, err = rpgit.gitCommit(comment)
undate_git_log(self, out, err)
def __str__(self):
loci_string = ''
for i in self.loci:
loci_string += i.name+','
return 'Project object with the loci '+loci_string
def __records_list_to_dict__(self):
self._records_dict = SeqIO.to_dict(self.records)
def last_git_log(self):
print self.git_log.split('<<<<')[-1]
def show_commits(self):
print rpgit.gitLog()[0]
###################################
# Project methods for reading data
###################################
def read_embl_genbank(self, input_filenames_list):
"""
Read a file from Genbank of EMBL