-
Notifications
You must be signed in to change notification settings - Fork 2
/
design.py
1728 lines (1595 loc) · 86.9 KB
/
design.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 python3
"""Design guides for diagnostics."""
import argparse
from collections import defaultdict
import logging
import os
import random
import re
import shutil
import sys
import tempfile
import numpy as np
from adapt import alignment
from adapt import guide_search
from adapt.prepare import align
from adapt.prepare import ncbi_neighbors
from adapt.prepare import prepare_alignment
from adapt import primer_search
from adapt.specificity import alignment_query
from adapt import target_search
from adapt.utils import log
from adapt.utils import predict_activity
from adapt.utils import seq_io
from adapt.utils import year_cover
from adapt.utils import mutate
from adapt.utils import weight
from adapt.utils import thermo
from adapt.utils.version import get_project_path, get_latest_model_version
try:
import boto3
from botocore.exceptions import ClientError
except ImportError:
cloud = False
else:
cloud = True
__author__ = 'Hayden Metsky <hmetsky@broadinstitute.org>, Priya P. Pillai <ppillai@broadinstitute.org>'
logger = logging.getLogger(__name__)
# Define defaults for parameters associated with different objectives
# The reason to define them here, rather than in argparse, is that
# defining them outside of argparse lets us check whether the
# argument has been set (versus just set to its default)
OBJ_PARAM_DEFAULTS = {
'minimize-guides': {
'guide_mismatches': 0,
'guide_cover_frac': 1.0
},
'maximize-activity': {
'guide_mismatches': 0,
'soft_guide_constraint': 1,
'hard_guide_constraint': 5,
'penalty_strength': 0.25,
'maximization_algorithm': 'random-greedy'
}
}
def check_obj_args(args):
"""Check argument depending on the objective function and set defaults.
Raise critical messages if arguments are set for the wrong objective
function than what is specified, and set defaults for these arguments if
they are not specified
Args:
args: namespace of arguments provided to this executable
"""
if args.obj == 'minimize-guides':
if (args.soft_guide_constraint or args.hard_guide_constraint or
args.penalty_strength or args.maximization_algorithm):
logger.critical(("When --obj is 'minimize-guides', the "
"following arguments are not used: --soft-guide-constraint, "
"--hard-guide-constraint, --penalty-strength, "
"--maximization-algorithm"))
if args.use_simple_binary_activity_prediction:
raise Exception(("Cannot use --use-simple-binary-activity-prediction "
"when --obj is 'minimize-guides'"))
if ((args.predict_activity_model_path or
args.predict_cas13a_activity_model is not None) and
args.guide_mismatches is None):
logger.critical(("When --obj is 'minimize-guides', the argument "
"-gm/--guide-mismatches is used even when a predictive "
"model is provided. The argument improves runtime by lowering "
"the number of calls to the model: to detect a sequence, a "
"guide must satisfy --guide-mismatches *and* the predicted "
"activity threshold. This argument has not been provided and "
"the default value (%d) might be too restrictive; to ignore "
"a requirement on mismatches, consider setting "
"--guide-mismatches to a sufficiently high value."),
OBJ_PARAM_DEFAULTS['minimize-guides']['guide_mismatches'])
for arg in OBJ_PARAM_DEFAULTS['minimize-guides'].keys():
if vars(args)[arg] is None:
vars(args)[arg] = OBJ_PARAM_DEFAULTS['minimize-guides'][arg]
if args.obj == 'maximize-activity':
if args.use_simple_binary_activity_prediction:
if (args.guide_cover_frac or args.cover_by_year_decay):
logger.critical(("When --obj is 'maximize-activity', the "
"following arguments are not used: "
"--guide-cover-frac, --cover-by-year-decay"))
if args.predict_activity_model_path:
logger.critical(("When --use-simple-binary-activity-prediction "
"is set, --predict-activity-model-path is not used"))
if args.predict_cas13a_activity_model:
logger.critical(("When --use-simple-binary-activity-prediction "
"is set, --predict-cas13a-activity-model is not used"))
else:
if (args.guide_mismatches or args.guide_cover_frac or
args.cover_by_year_decay):
logger.critical(("When --obj is 'maximize-activity', the "
"following arguments are not used: --guide-mismatches, "
"--guide-cover-frac, --cover-by-year-decay"))
for arg in OBJ_PARAM_DEFAULTS['maximize-activity'].keys():
if vars(args)[arg] is None:
vars(args)[arg] = OBJ_PARAM_DEFAULTS['maximize-activity'][arg]
def seqs_grouped_by_year(seqs, args, sequence_weights=None):
"""Group sequences according to their year and assigned partial covers.
Args:
seqs: dict mapping sequence name to sequence, as read from a FASTA
args: namespace of arguments provided to this executable
Returns:
tuple (aln, years_idx, cover_frac) where aln is an
alignment.Alignment object from seqs; years_idx is a dict
mapping each year to the set of indices in aln representing
sequences for that year; and cover_frac is a dict mapping each
year to the desired partial cover of sequences from that year,
as determined by args.cover_by_year_decay (if args.search_cmd
is 'complete-targets', then cover_frac is a tuple (g, p) where
g is the previously described dict calculated from
args.guide_cover_frac and p is the same calculated from
args.primer_cover_frac)
"""
years_fn, year_highest_cover, year_cover_decay = args.cover_by_year_decay
norm_weights = weight.normalize(sequence_weights, seqs.keys())
# Map sequence names to index in alignment, and construct alignment
seq_list = []
seq_norm_weights = []
seq_idx = {}
for i, (name, seq) in enumerate(seqs.items()):
seq_idx[name] = i
seq_list.append(seq)
seq_norm_weights.append(norm_weights[name])
aln = alignment.Alignment.from_list_of_seqs(
seq_list, seq_norm_weights=seq_norm_weights)
# Read sequences for each year, and check that every sequence has
# a year
years = year_cover.read_years(years_fn)
all_seqs_with_year = set.union(*years.values())
for seq in seq_idx.keys():
if seq not in all_seqs_with_year:
raise Exception("Unknown year for sequence '%s'" % seq)
# Convert years dict to map to indices rather than sequence names
years_idx = {}
for year in years.keys():
# Skip names not in seq_idx because the years file may contain
# sequences that are not in seqs
years_idx[year] = set(seq_idx[name] for name in years[year]
if name in seq_idx)
# Construct desired partial cover for each year
guide_cover_frac = year_cover.construct_partial_covers(
years.keys(), year_highest_cover, args.guide_cover_frac, year_cover_decay)
if args.search_cmd == 'complete-targets':
primer_cover_frac = year_cover.construct_partial_covers(
years.keys(), year_highest_cover, args.primer_cover_frac, year_cover_decay)
cover_frac = (guide_cover_frac, primer_cover_frac)
else:
cover_frac = guide_cover_frac
return aln, years_idx, cover_frac
def parse_required_guides_and_ignore(args):
"""Parse files giving required guides and ignored sequences.
Args:
args: namespace of arguments provided to this executable
Returns:
tuple (required_guides, ignored_ranges, disallowed_kmers) where
required_guides is a representation of data in the args.required_guides
file; ignored_ranges is a representation of data in the
args.ignored_ranges file; and disallowed_kmers is a
representation of data in the args.disallowed_kmers file
"""
num_aln = len(args.in_fasta)
# Read required guides, if provided
if args.required_guides:
required_guides = seq_io.read_required_guides(
args.required_guides, args.guide_length, num_aln)
else:
required_guides = [{} for _ in range(num_aln)]
# Read ignored ranges, if provided
if args.ignored_ranges:
ignored_ranges = seq_io.read_ignored_ranges(
args.ignored_ranges, num_aln)
else:
ignored_ranges = [set() for _ in range(num_aln)]
# Read disallowed kmers, if provided
if args.disallowed_kmers:
disallowed_kmers = seq_io.read_disallowed_kmers(
args.disallowed_kmers,
min_len_warning=5,
max_len_warning=args.guide_length)
else:
disallowed_kmers = set()
return required_guides, ignored_ranges, disallowed_kmers
def prepare_alignments(args):
"""Download, curate, and align sequences for input.
Args:
args: namespace of arguments provided to this executable
Returns:
tuple (in_fasta, taxid_for_fasta, years_tsv, aln_tmp_dirs, out_tsv, design_for)
in which in_fasta is a list of paths to fasta files each containing an
alignment, taxid_for_fasta[i] gives a taxon id for in_fasta[i],
years_tsv gives a tempfile storing a tsv file containing a year for
each sequence across all the fasta files (only if
args.cover_by_year_decay is set), aln_tmp_dirs is a list of temp
directories that need to be cleaned up, out_tsv[i] is a path to the
file at which to write the output for in_fasta[i], and design_for[i]
indicates whether to actually design for in_fasta[i] (True/False)
"""
logger.info(("Setting up to prepare alignments"))
# Set the path to mafft
align.set_mafft_exec(args.mafft_path)
# Initialize variables to None
am = None
asm = None
sequences_to_use = None
accessions_to_use = None
taxs_to_design_for = None
# Read list of taxonomies
if args.input_type == 'auto-from-args':
s = None if args.segment == 'None' else args.segment
args.tax_id = ncbi_neighbors.determine_current_taxid(args.tax_id)
ref_accs = ncbi_neighbors.construct_references(
args.tax_id, args.segment) if not args.ref_accs else args.ref_accs
meta_filt = None
meta_filt_against = None
if args.metadata_filter:
meta_filt = seq_io.read_metadata_filters(args.metadata_filter)
if args.specific_against_metadata_filter:
meta_filt_against = seq_io.read_metadata_filters(args.specific_against_metadata_filter)
taxs = [(None, args.tax_id, s, ref_accs, meta_filt, meta_filt_against)]
elif args.input_type == 'auto-from-file':
taxs = seq_io.read_taxonomies(args.in_tsv)
elif args.input_type == 'fasta':
sequences_to_use = {}
taxs = []
args.sample_seqs = None
args.write_annotation = False
args.write_input_seqs = False
args.write_weights = False
args.weight_by_log_size_of_subtaxa = False
for i, unaligned_fasta in enumerate(args.in_fasta):
sequences_to_use[(0, unaligned_fasta)] = seq_io.read_fasta(unaligned_fasta)
taxs.append((None, 0, unaligned_fasta, [], None, None))
else:
raise Exception(("Unknown input type '%s'") % args.input_type)
if args.input_type in ['auto-from-args', 'auto-from-file']:
# Setup alignment and alignment stat memoizers
if args.prep_memoize_dir:
if args.prep_memoize_dir[:5] == "s3://":
bucket = args.prep_memoize_dir.split("/")[2]
try:
if args.aws_access_key_id is not None and args.aws_secret_access_key is not None:
S3 = boto3.client("s3", aws_access_key_id=args.aws_access_key_id,
aws_secret_access_key=args.aws_secret_access_key)
else:
S3 = boto3.client("s3")
S3.head_bucket(Bucket=bucket)
except ClientError as e:
if e.response['Error']['Code'] == "404":
S3.create_bucket(Bucket=bucket)
elif e.response['Error']['Code'] == "403":
raise Exception(("Incorrect AWS Access Key ID or Secret Access Key")) from e
else:
raise
except ConnectionError as e:
raise Exception(("Cannot connect to Amazon S3")) from e
if args.prep_memoize_dir[-1] == "/":
align_memoize_dir = '%saln' % args.prep_memoize_dir
align_stat_memoize_dir = '%sstats' % args.prep_memoize_dir
else:
align_memoize_dir = '%s/aln' % args.prep_memoize_dir
align_stat_memoize_dir = '%s/stats' % args.prep_memoize_dir
else:
if not os.path.isdir(args.prep_memoize_dir):
raise Exception(("Path '%s' does not exist") %
args.prep_memoize_dir)
align_memoize_dir = os.path.join(args.prep_memoize_dir, 'aln')
if not os.path.exists(align_memoize_dir):
os.makedirs(align_memoize_dir)
align_stat_memoize_dir = os.path.join(args.prep_memoize_dir, 'stats')
if not os.path.exists(align_stat_memoize_dir):
os.makedirs(align_stat_memoize_dir)
am = align.AlignmentMemoizer(align_memoize_dir,
aws_access_key_id=args.aws_access_key_id,
aws_secret_access_key=args.aws_secret_access_key)
asm = align.AlignmentStatMemoizer(align_stat_memoize_dir,
aws_access_key_id=args.aws_access_key_id,
aws_secret_access_key=args.aws_secret_access_key)
# Read specified accessions, if provided
if args.use_accessions:
accessions_to_use = seq_io.read_accessions_for_taxonomies(
args.use_accessions)
# Read specified sequences, if provided
if args.use_fasta:
sequences_to_use = seq_io.read_sequences_for_taxonomies(
args.use_fasta)
# Only design for certain taxonomies, if provided
if args.only_design_for:
taxs_to_design_for = seq_io.read_taxonomies_to_design_for(
args.only_design_for)
# Construct alignments for each taxonomy
in_fasta = []
taxid_for_fasta = []
years_tsv_per_aln = []
aln_tmp_dirs = []
out_tsv = []
design_for = []
specific_against_metadata_accs = []
annotations = []
sequence_weights = []
for label, tax_id, segment, ref_accs, meta_filt, meta_filt_against in taxs:
aln_file_dir = tempfile.TemporaryDirectory()
if args.cover_by_year_decay:
years_tsv_tmp = tempfile.NamedTemporaryFile()
years_tsv_tmp_name = years_tsv_tmp.name
else:
years_tsv_tmp = None
years_tsv_tmp_name = None
if accessions_to_use is not None:
if (tax_id, segment) in accessions_to_use:
accessions_to_use_for_tax = accessions_to_use[(tax_id, segment)]
else:
accessions_to_use_for_tax = None
else:
accessions_to_use_for_tax = None
if sequences_to_use is not None:
if (tax_id, segment) in sequences_to_use:
sequences_to_use_for_tax = sequences_to_use[(tax_id, segment)]
else:
sequences_to_use_for_tax = None
else:
sequences_to_use_for_tax = None
if (accessions_to_use_for_tax is not None and
sequences_to_use_for_tax is not None):
raise Exception(("Cannot use both --use-accessions and "
"--use-fasta for the same taxonomy"))
annotation_tsv = os.path.join(args.out_tsv_dir, label) if label \
else args.write_annotation
nc, specific_against_metadata_acc, annotation, sequence_weight = \
prepare_alignment.prepare_for(
tax_id, segment, ref_accs,
aln_file_dir.name, aln_memoizer=am, aln_stat_memoizer=asm,
sample_seqs=args.sample_seqs,
years_tsv=years_tsv_tmp_name,
annotation_tsv=annotation_tsv,
cluster_threshold=args.cluster_threshold,
accessions_to_use=accessions_to_use_for_tax,
sequences_to_use=sequences_to_use_for_tax,
meta_filt=meta_filt,
meta_filt_against=meta_filt_against,
subtaxa_weight=args.weight_by_log_size_of_subtaxa)
for i in range(nc):
in_fasta += [os.path.join(aln_file_dir.name, str(i) + '.fasta')]
taxid_for_fasta += [tax_id]
specific_against_metadata_accs.append(specific_against_metadata_acc)
annotations.append(annotation[i])
sequence_weights.append(sequence_weight[i])
if taxs_to_design_for is None:
design_for += [True]
else:
design_for += [(tax_id, segment) in taxs_to_design_for]
years_tsv_per_aln += [years_tsv_tmp]
aln_tmp_dirs += [aln_file_dir]
if label is None:
if args.input_type == 'fasta':
# If the input is a FASTA, the 'tax_id' is actually an index number
out_tsv.extend([args.out_tsv[tax_id] + '.' + str(i) + '.tsv' for i in range(nc)])
else:
out_tsv.extend([args.out_tsv + '.' + str(i) + '.tsv' for i in range(nc)])
else:
for i in range(nc):
out_name = label + '.' + str(i) + '.tsv'
out_tsv.append(os.path.join(args.out_tsv_dir, out_name))
if args.write_input_seqs:
# Write the sequences that are in the alignment being used
# as input
all_seq_names = []
for i in range(nc):
fn = os.path.join(aln_file_dir.name, str(i) + '.fasta')
seqs = seq_io.read_fasta(fn)
all_seq_names += list(seqs.keys())
all_seq_names = sorted(all_seq_names)
if label is None:
# args.write_input_seqs gives the path to where to write
# the list
out_file = args.write_input_seqs + '.txt'
else:
# Determine where to write the sequence names based on
# the label and args.out_tsv_dir
out_name = label + '.input-sequences.txt'
out_file = os.path.join(args.out_tsv_dir, out_name)
with open(out_file, 'w') as fw:
for name in all_seq_names:
fw.write(name + '\n')
if args.write_input_aln:
# Write the alignments being used as input
for i in range(nc):
fn = os.path.join(aln_file_dir.name, str(i) + '.fasta')
if label is None:
if args.input_type == 'fasta':
# If the input is a FASTA, the 'tax_id' is actually an index number
copy_path = args.write_input_aln[tax_id] + '.' + str(i) + '.fasta'
else:
# args.write_input_aln gives the prefix of the path to
# which to write the alignment
copy_path = args.write_input_aln + '.' + str(i) + '.fasta'
else:
# Determine where to write the alignment based on the
# label and args.out_tsv_dir
out_name = label + '.' + str(i) + '.fasta'
copy_path = os.path.join(args.out_tsv_dir, out_name)
shutil.copyfile(fn, copy_path)
if args.write_weights:
# Write the weights being used
for i in range(nc):
if label is None:
# args.write_input_seqs gives the path to where to write
# the list
out_file = args.write_weights + '.' + str(i) + '.weights.tsv'
else:
# Determine where to write the sequence names based on
# the label and args.out_tsv_dir
out_name = label + '.' + str(i) + '.weights.tsv'
out_file = os.path.join(args.out_tsv_dir, out_name)
with open(out_file, 'w') as fw:
for seq_name, seq_weight in sequence_weight[i].items():
fw.write("%s\t%f\n" %(seq_name, seq_weight))
# Combine all years tsv (there is one per fasta file)
if any(f is not None for f in years_tsv_per_aln):
years_tsv = tempfile.NamedTemporaryFile()
with open(years_tsv.name, 'w') as fw:
for tf in years_tsv_per_aln:
if tf is not None:
with open(tf.name) as fin:
for line in fin:
fw.write(line)
tf.close()
else:
years_tsv = None
return (in_fasta, taxid_for_fasta, years_tsv, aln_tmp_dirs, out_tsv,
design_for, specific_against_metadata_accs, annotations,
sequence_weights)
def design_for_id(args):
"""Design guides for differential identification across targets.
Args:
args: namespace of arguments provided to this executable
"""
# Create an alignment object for each input
# If obj is not 'minimize-guides', guide_cover_frac and seq_groups
# will be set to None
alns = []
seq_groups_per_input = []
guide_cover_frac_per_input = []
primer_cover_frac_per_input = []
for i, in_fasta in enumerate(args.in_fasta):
seqs = seq_io.read_fasta(in_fasta)
# Warn if some, but not all, sequences have weights set
if len(args.sequence_weights[i]) > 0 and len(args.sequence_weights[i]) < len(seqs):
logger.warning(("Only %i sequences of %i have weights set for "
"alignment %i; the rest will be given a default weight of 1 "
"pre-normalization."
%(len(args.sequence_weights[i]), len(seqs), i)))
if args.cover_by_year_decay:
aln, seq_groups, cover_frac = seqs_grouped_by_year(seqs, args,
sequence_weights=args.sequence_weights[i])
if args.search_cmd == 'complete-targets':
guide_cover_frac, primer_cover_frac = cover_frac
else:
guide_cover_frac = cover_frac
primer_cover_frac = None
else:
seq_names = list(seqs.keys())
seq_list = [seqs[seq_name] for seq_name in seq_names]
norm_weights = weight.normalize(args.sequence_weights[i], seq_names)
seq_norm_weights = [norm_weights[seq_name]
for seq_name in seq_names]
aln = alignment.Alignment.from_list_of_seqs(
seq_list, seq_norm_weights=seq_norm_weights)
seq_groups = None
guide_cover_frac = args.guide_cover_frac
if args.search_cmd == 'complete-targets':
primer_cover_frac = args.primer_cover_frac
else:
primer_cover_frac = None
alns += [aln]
seq_groups_per_input += [seq_groups]
guide_cover_frac_per_input += [guide_cover_frac]
primer_cover_frac_per_input += [primer_cover_frac]
# Keep track of how many items in alns to use for design (the first N of them)
num_aln_for_design = len(args.in_fasta)
# Add sequence lists to alns to be specific against
# Note that although these are being place in alns, they may not
# be actual alignments (they can be unaligned sequences)
# If specific-against-metadata-filter is specified, add metadata filtered accessions to alns
# and store what index it is located at in alns. Also, store start and end indices of
# specific_against_metadata sequence lists
specific_against_metadata_indices = {}
specific_against_metadata_start = len(alns)
for i in range(num_aln_for_design):
if len(args.specific_against_metadata_accs[i]) > 0:
logger.info(("Fetching %d sequences within taxon %d to be specific against"),
len(args.specific_against_metadata_accs[i]), args.taxid_for_fasta[i])
seqs = prepare_alignment.fetch_sequences_for_acc_list(list(
args.specific_against_metadata_accs[i]))
seqs_list = alignment.SequenceList(list(seqs.values()))
alns += [seqs_list]
specific_against_metadata_indices[i] = len(alns) - 1
specific_against_metadata_end = len(alns)
# Also add the specific_against_fastas sequences into alns
for specific_against_fasta in args.specific_against_fastas:
seqs = seq_io.read_fasta(specific_against_fasta)
seqs_list = alignment.SequenceList(list(seqs.values()))
alns += [seqs_list]
# Finally, add specific_against_taxa sequences into alns after
# downloading them
if args.specific_against_taxa:
taxs_to_be_specific_against = seq_io.read_taxonomies_to_design_for(
args.specific_against_taxa)
for given_taxid, segment in taxs_to_be_specific_against:
logger.info(("Fetching sequences to be specific against tax "
"%d (segment: %s)"), given_taxid, segment)
taxid = ncbi_neighbors.determine_current_taxid(given_taxid)
seqs = prepare_alignment.fetch_sequences_for_taxonomy(
taxid, segment)
seqs_list = alignment.SequenceList(list(seqs.values()))
alns += [seqs_list]
specific_against_exists = ((len(args.specific_against_fastas) > 0 or
args.specific_against_taxa is not None) or
(specific_against_metadata_end - specific_against_metadata_start) > 0)
required_guides, ignored_ranges, disallowed_kmers = \
parse_required_guides_and_ignore(args)
required_flanking_seqs = (args.require_flanking5, args.require_flanking3)
# Allow G-U base pairing, unless it is explicitly disallowed
allow_gu_pairs = not args.do_not_allow_gu_pairing
# Assign an id in [0, 1, 2, ...] for each taxid to design for
# Find all alignments with each taxid
aln_with_taxid = defaultdict(set)
for i, taxid in enumerate(args.taxid_for_fasta):
aln_with_taxid[taxid].add(i)
num_taxa = len(aln_with_taxid)
logger.info(("Designing for %d taxa"), num_taxa)
# Read taxonomies to ignore for specificity, if specified
tax_ignore = {}
if (args.input_type in ['auto-from-file', 'auto-from-args'] and
args.taxa_to_ignore_for_specificity):
tax_ignore = seq_io.read_taxonomy_specificity_ignore(
args.taxa_to_ignore_for_specificity)
if args.specific_against_fastas or args.specific_against_taxa:
logger.warning(("Taxa to ignore for specificity cannot from "
"--specific-against-*"))
# Construct the data structure for guide queries to perform
# differential identification
if num_taxa > 1 or specific_against_exists:
logger.info(("Constructing data structure to permit guide queries for "
"differential identification"))
if args.diff_id_method == "lshnn":
aq = alignment_query.AlignmentQuerierWithLSHNearNeighbor(alns,
args.guide_length, args.diff_id_mismatches, allow_gu_pairs)
elif args.diff_id_method == "shard":
aq = alignment_query.AlignmentQuerierWithKmerSharding(alns,
args.guide_length, args.diff_id_mismatches, allow_gu_pairs)
else:
raise Exception(("Unknown method for querying specificity: '%s'" %
args.diff_id_method))
aq.setup()
else:
logger.info(("Only one taxon was provided, so not constructing "
"data structure to permit queries for differential "
"identification"))
aq = None
for i in range(num_aln_for_design):
taxid = args.taxid_for_fasta[i]
logger.info(("Finding guides for alignment %d (of %d), which is in "
"taxon %d"), i + 1, num_aln_for_design, taxid)
if args.design_for is not None and args.design_for[i] is False:
logger.info("Skipping design for this alignment")
continue
aln = alns[i]
seq_groups = seq_groups_per_input[i]
guide_cover_frac = guide_cover_frac_per_input[i]
primer_cover_frac = primer_cover_frac_per_input[i]
required_guides_for_aln = required_guides[i]
ignored_ranges_for_aln = ignored_ranges[i]
alns_in_same_taxon = aln_with_taxid[taxid]
# For metadata filtering, we only want to be specific against the
# accessions in alns[specific_against_metadata_index]
specific_against_metadata_index = specific_against_metadata_indices[i] \
if i in specific_against_metadata_indices else None
if aq is not None:
guide_is_specific = aq.guide_is_specific_to_alns_fn(
alns_in_same_taxon, args.diff_id_frac,
do_not_memoize=args.do_not_memoize_oligo_computations)
else:
# No specificity to check
guide_is_specific = lambda guide: True
def guide_is_suitable(guide):
# Return True iff the guide does not contain a disallowed
# k-mer and is specific to aln
# Return False if the guide contains a disallowed k-mer
for kmer in disallowed_kmers:
if kmer in guide:
return False
# Return True if guide does not hit too many sequences in
# alignments other than aln
return guide_is_specific(guide)
# Mask alignments from this taxon from being reported in queries
# because we will likely get many guide sequences that hit its
# alignments, but we do not care about these for checking specificity
if aq is not None:
for j in range(num_aln_for_design):
if j in alns_in_same_taxon:
aq.mask_aln(j)
logger.info(("Masking alignment %d as it is taxon %d"),
j + 1, taxid)
# Also mask any specific against metadata sequence lists that do not
# match the one for this alignment
if specific_against_metadata_index is not None:
logger.info(("Masking all sequence lists for specificity except the one "
"for alignment %d (taxon %d)"), i + 1, taxid)
for j in range(specific_against_metadata_start, specific_against_metadata_end):
if specific_against_metadata_index != j:
aq.mask_aln(specific_against_metadata_index)
# Also mask taxonomies to ignore when determining specificity
# of taxid
if taxid in tax_ignore:
for ignore_taxid in tax_ignore[taxid]:
for j in aln_with_taxid[ignore_taxid]:
logger.info(("Masking alignment %d (from taxon %d) "
"from specificity queries"), j + 1, ignore_taxid)
aq.mask_aln(j)
# Construct activity predictor
if args.use_simple_binary_activity_prediction:
predictor = predict_activity.SimpleBinaryPredictor(
args.guide_mismatches,
allow_gu_pairs,
required_flanking_seqs=required_flanking_seqs)
elif (args.predict_activity_model_path or
args.predict_cas13a_activity_model is not None):
if args.predict_activity_model_path:
cla_path, reg_path = args.predict_activity_model_path
else:
dir_path = get_project_path()
cla_path_all = os.path.join(dir_path, 'models', 'classify',
'cas13a')
reg_path_all = os.path.join(dir_path, 'models', 'regress',
'cas13a')
if len(args.predict_cas13a_activity_model) not in (0,2):
raise Exception(("If setting versions for "
"--predict-cas13a-activity-model, both a version for "
"the classifier and the regressor must be set."))
if (len(args.predict_cas13a_activity_model) == 0 or
args.predict_cas13a_activity_model[0] == 'latest'):
cla_version = get_latest_model_version(cla_path_all)
else:
cla_version = args.predict_cas13a_activity_model[0]
if (len(args.predict_cas13a_activity_model) == 0 or
args.predict_cas13a_activity_model[1] == 'latest'):
reg_version = get_latest_model_version(reg_path_all)
else:
reg_version = args.predict_cas13a_activity_model[1]
cla_path = os.path.join(cla_path_all, cla_version)
reg_path = os.path.join(reg_path_all, reg_version)
if args.predict_activity_thres:
# Use specified thresholds on classification and regression
cla_thres, reg_thres = args.predict_activity_thres
else:
# Use default thresholds specified with the model
cla_thres, reg_thres = None, None
predictor = predict_activity.Predictor(cla_path, reg_path,
classification_threshold=cla_thres,
regression_threshold=reg_thres)
else:
if args.predict_activity_thres:
raise Exception(("Cannot set --predict-activity-thres without "
"setting --predict-activity-model-path or "
"--predict-cas13a-activity-model"))
if args.obj == 'maximize-activity':
raise Exception(("Either --predict-activity-model-path or "
"--predict-cas13a-activity-model must be specified if "
"--obj is 'maximize-activity'. Alternatively, if you "
"cannot use our pre-trained Cas13a model nor another "
"model, you can set "
"--use-simple-binary-activity-prediction"))
if args.predict_activity_degradation:
raise Exception(("--predict-activity-model-path must be "
"specified if --predict-activity-degradation is set "
"(unless --use-simple-binary-activity-prediction is set)"))
# Do not predict activity
predictor = None
# Construct mutator, if necessary
if args.predict_activity_degradation:
mutator = mutate.GTRSubstitutionMutator(
aln, *args.predict_activity_degradation,
args.predict_activity_degradation_mu,
args.predict_activity_degradation_t,
args.predict_activity_degradation_n)
else:
mutator = None
# Find an optimal set of guides for each window in the genome,
# and write them to a file; ensure that the selected guides are
# specific to this alignment
if args.obj == 'minimize-guides':
gs = guide_search.GuideSearcherMinimizeGuides(
aln,
args.guide_length,
args.guide_mismatches,
guide_cover_frac,
args.missing_thres,
seq_groups=seq_groups,
pre_filter_fns=[guide_is_suitable],
required_oligos=required_guides_for_aln,
ignored_ranges=ignored_ranges_for_aln,
allow_gu_pairs=allow_gu_pairs,
required_flanking_seqs=required_flanking_seqs,
predictor=predictor,
do_not_memoize=args.do_not_memoize_oligo_computations)
elif args.obj == 'maximize-activity':
gs = guide_search.GuideSearcherMaximizeActivity(
aln,
args.guide_length,
args.soft_guide_constraint,
args.hard_guide_constraint,
args.penalty_strength,
args.missing_thres,
algorithm=args.maximization_algorithm,
pre_filter_fns=[guide_is_suitable],
required_oligos=required_guides_for_aln,
ignored_ranges=ignored_ranges_for_aln,
allow_gu_pairs=allow_gu_pairs,
required_flanking_seqs=required_flanking_seqs,
predictor=predictor,
do_not_memoize=args.do_not_memoize_oligo_computations)
if args.search_cmd == 'sliding-window':
# Find an optimal set of guides for each window in the genome,
# and write them to a file
gs.find_guides_with_sliding_window(args.window_size,
args.out_tsv[i],
window_step=args.window_step,
sort=args.sort_out,
print_analysis=(args.log_level<=logging.INFO))
elif args.search_cmd == 'complete-targets':
# Find optimal targets (primer and guide set combinations),
# and write them to a file
if args.primer_gc_content_bounds is None:
primer_gc_content_bounds = None
else:
primer_gc_content_bounds = tuple(args.primer_gc_content_bounds)
if args.primer_thermo:
shared_memo = {}
conditions = thermo.Conditions(sodium=args.pcr_sodium_conc,
magnesium=args.pcr_magnesium_conc, dNTP=args.pcr_dntp_conc,
oligo_concentration=args.pcr_oligo_conc,
target_concentration=args.pcr_target_conc)
left_primer_predictor = predict_activity.TmPredictor(
args.ideal_primer_melting_temperature +
thermo.CELSIUS_TO_KELVIN,
conditions, False, shared_memo=shared_memo)
right_primer_predictor = predict_activity.TmPredictor(
args.ideal_primer_melting_temperature +
thermo.CELSIUS_TO_KELVIN,
conditions, True, shared_memo=shared_memo)
post_filter_primers = [lambda oligo:
thermo.has_no_secondary_structure(oligo, conditions)]
# Since TmPredictor returns negative values, the algorithm must
# be greedy.
lps = primer_search.PrimerSearcherMaximizeActivity(
aln, args.primer_length, args.primer_length,
args.soft_primer_constraint, args.hard_primer_constraint,
args.primer_penalty_strength, args.missing_thres,
algorithm='greedy',
primer_gc_content_bounds=primer_gc_content_bounds,
post_filter_fns=post_filter_primers,
predictor=left_primer_predictor)
rps = primer_search.PrimerSearcherMaximizeActivity(
aln, args.primer_length, args.primer_length,
args.soft_primer_constraint, args.hard_primer_constraint,
args.primer_penalty_strength, args.missing_thres,
algorithm='greedy',
primer_gc_content_bounds=primer_gc_content_bounds,
post_filter_fns=post_filter_primers,
predictor=right_primer_predictor)
primer_set_filter_fns = [lambda oligo_set:
thermo.has_no_heterodimers(oligo_set, conditions)]
else:
lps = primer_search.PrimerSearcherMinimizePrimers(
aln, args.primer_length, args.primer_mismatches,
primer_cover_frac, args.missing_thres, seq_groups=seq_groups,
primer_gc_content_bounds=primer_gc_content_bounds)
rps = lps
ts = target_search.TargetSearcher(lps, rps, gs,
max_primers_at_site=args.max_primers_at_site,
max_target_length=args.max_target_length,
obj_weights=args.obj_fn_weights,
only_account_for_amplified_seqs=args.only_account_for_amplified_seqs,
halt_early=args.halt_search_early, mutator=mutator,
primer_set_filter_fns=[])
ts.find_and_write_targets(args.out_tsv[i],
best_n=args.best_n_targets, no_overlap=args.do_not_overlap,
annotations=args.annotations[i])
else:
raise Exception("Unknown search subcommand '%s'" % args.search_cmd)
# i should no longer be masked from queries
if aq is not None:
aq.unmask_all_aln()
def run(args):
logger = logging.getLogger(__name__)
# Set random seed for entire program
random.seed(args.seed)
np.random.seed(args.seed)
check_obj_args(args)
logger.info("Running design.py with arguments: %s", args)
if args.input_type in ['auto-from-file', 'auto-from-args']:
if args.ncbi_api_key:
ncbi_neighbors.set_ncbi_api_key(args.ncbi_api_key)
if args.input_type == 'auto-from-file':
if not os.path.isdir(args.out_tsv_dir):
raise Exception(("Output directory '%s' does not exist") %
args.out_tsv_dir)
# Prepare input alignments, stored in temp fasta files
in_fasta, taxid_for_fasta, years_tsv, aln_tmp_dirs, out_tsv, \
design_for, specific_against_metadata_accs, annotations, \
sequence_weights = prepare_alignments(args)
args.in_fasta = in_fasta
args.taxid_for_fasta = taxid_for_fasta
args.out_tsv = out_tsv
args.design_for = design_for
args.specific_against_metadata_accs = specific_against_metadata_accs
args.annotations = annotations
args.sequence_weights = sequence_weights
if args.cover_by_year_decay:
# args.cover_by_year_decay contains two parameters: the year
# with the highest cover and the decay; add in (to the beginning)
# the file listing the years
year_highest_cover, year_cover_decay = args.cover_by_year_decay
args.cover_by_year_decay = (years_tsv.name, year_highest_cover,
year_cover_decay)
elif args.input_type == 'fasta':
if len(args.in_fasta) != len(args.out_tsv):
raise Exception(("Number output TSVs must match number of input "
"FASTAs"))
if args.write_input_aln:
if len(args.in_fasta) != len(args.write_input_aln):
raise Exception(("Number input alignment filenames must match "
"number of input FASTAs"))
if args.weight_sequences:
if len(args.in_fasta) != len(args.weight_sequences):
raise Exception(("Number of sequence weight TSVs must match "
"number of input FASTAs"))
args.sequence_weights = [seq_io.read_sequence_weights(fn) for
fn in args.weight_sequences]
else:
args.sequence_weights = [defaultdict(lambda: 1) for _ in
args.in_fasta]
args.design_for = None
args.taxid_for_fasta = list(range(len(args.in_fasta)))
args.specific_against_metadata_accs = [[] for _ in range(len(args.in_fasta))]
args.annotations = [[] for _ in range(len(args.in_fasta))]
aln_tmp_dirs = []
if args.unaligned:
in_fasta, taxid_for_fasta, years_tsv, aln_tmp_dirs, out_tsv, \
design_for, specific_against_metadata_accs, annotations, \
sequence_weights = prepare_alignments(args)
args.in_fasta = in_fasta
args.taxid_for_fasta = taxid_for_fasta
args.out_tsv = out_tsv
args.design_for = design_for
args.specific_against_metadata_accs = specific_against_metadata_accs
args.annotations = annotations
else:
args.out_tsv = [out_tsv + '.tsv' for out_tsv in args.out_tsv]
else:
raise Exception("Unknown input type subcommand '%s'" % args.input_type)
design_for_id(args)
# Close temporary files storing alignments
for td in aln_tmp_dirs:
td.cleanup()
if args.input_type in ['auto-from-file', 'auto-from-args']:
if years_tsv is not None:
years_tsv.close()
def argv_to_args(argv):
parser = argparse.ArgumentParser()
###########################################################################
# OPTIONS AVAILABLE ACROSS ALL SUBCOMMANDS
###########################################################################
base_subparser = argparse.ArgumentParser(add_help=False)
# Guide length
base_subparser.add_argument('-gl', '--guide-length', type=int, default=28,
help="Length of guide to construct")
# Objective function
base_subparser.add_argument('--obj',
choices=['maximize-activity', 'minimize-guides'],
default='maximize-activity',
help=(("Objective function to solve. 'maximize-activity' maximizes "
"the expected activity of the guide set of the target genomes "
"subject to soft and hard constraints on the size of the guide "
"set. 'minimize-guides' minimizes the number of guides in the "
"guide set subject to coverage constraints across the target "