forked from mjhubisz/argweaver
-
Notifications
You must be signed in to change notification settings - Fork 0
/
arg-sample.cpp
2092 lines (1856 loc) · 77.9 KB
/
arg-sample.cpp
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
// C/C++ includes
#ifdef ARGWEAVER_MPI
#include "mpi.h"
#endif
#include <time.h>
#include <memory>
#include <sys/stat.h>
#include <unistd.h>
// arghmm includes
#include "argweaver/compress.h"
#include "argweaver/ConfigParam.h"
#include "argweaver/emit.h"
#include "argweaver/fs.h"
#include "argweaver/logging.h"
#include "argweaver/mem.h"
#include "argweaver/parsing.h"
#include "argweaver/sample_arg.h"
#include "argweaver/sequences.h"
#include "argweaver/total_prob.h"
#include "argweaver/track.h"
#include "argweaver/est_popsize.h"
#include "argweaver/mcmcmc.h"
#include "argweaver/coal_records.h"
#include "argweaver/recomb.h"
using namespace argweaver;
// version info
#define VERSION_TEXT "1.0"
#define VERSION_INFO "\
ARGweaver " VERSION_TEXT " \n\
Matt Rasmussen\n\
Sampler for large ancestral recombination graphs\n\
"
// file extensions
const char *SMC_SUFFIX = ".smc";
const char *SITES_SUFFIX = ".sites";
const char *STATS_SUFFIX = ".stats";
const char *LOG_SUFFIX = ".log";
const char *COAL_RECORDS_SUFFIX = ".cr";
// help categories
const int ADVANCED_OPT = 1;
const int POPMODEL_OPT = 2;
const int EXPERIMENTAL_OPT = 3;
const int EXIT_ERROR = 1;
// parsing command-line options
class Config
{
public:
Config()
{
make_parser();
resample_region[0] = -1;
resample_region[1] = -1;
}
void make_parser()
{
config.clear();
config.add(new ConfigParamComment("Input/output options"));
// input/output
config.add(new ConfigParam<string>
("-o", "--output", "<output prefix>", &out_prefix,
"arg-sample",
"prefix for all output filenames (default='arg-sample')"));
config.add(new ConfigParam<string>
("-s", "--sites", "<sites alignment>", &sites_file,
"sequence alignment in sites format"));
config.add(new ConfigParam<string>
("-f", "--fasta", "<fasta alignment>", &fasta_file,
"sequence alignment in FASTA format"));
config.add(new ConfigParam<string>
("", "--vcf", "<.vcf.gz file>", &vcf_file,
"sequence alignment in gzipped vcf format. Must also supply"
" --region in format chr:start-end, and tabix index file"
" (.vcf.gz.tbi) must also be present. tabix program needs to"
" be installed and in PATH or path must be provided with"
" --tabix-dir. Assumes samples are diploid and unphased;"
" each individual will have _1 and _2 appended to its name for its"
" two haploid lineages. Indel-type variants are skipped."
" Any positions not specified in VCF are assumed to be"
" invariant"));
config.add(new ConfigParam<string>
("", "--vcf-files", "<vcf_file_list.txt>", &vcf_list_file,
"Same as --vcf but loads all individuals from multiple .vcf.gz files."
" The argument should be a file containing list of files to read"
" (one per line). All files should be aligned to same reference"
" genome"));
config.add(new ConfigParam<string>
("", "--rename-seqs", "<name_map_file.txt>", &rename_file,
"Used to rename sequences (usually from cryptic names in VCF"
" files to more meaningful names). The file should have two"
" columns, with the old sequence name, followed by the new"
" sequence name. Names can be diploid or haploid, and any old"
" names not matching the current sequence set will be ignored."));
config.add(new ConfigParam<string>
("", "--tabix-dir", "<directory>", &tabix_dir,
" path to tabix executable. May be required if using --vcf"
" and tabix is not in PATH"));
config.add(new ConfigParam<string>
("", "--age-file", "<age file>", &age_file,
" file giving age for any ancient samples (two-columns, "
" first column is sample name, second age in generations)."
" Age will be rounded to nearest discrete time point in"
" the model, with a maximum of the second-to-oldest point."));
config.add(new ConfigSwitch
("", "--use-genotype-probs", &use_genotype_probs,
"(for VCF input) Use genotype probabilities given by the"
" PL, GL, or PP scores in the VCF file. Default: ignore these scores and"
" treat assigned genotypes (after any filters) as"
" certain.", ADVANCED_OPT));
config.add(new ConfigParam<string>
("", "--region", "<start>-<end>",
&subregion_str, "",
"sample ARG for only a region of the sites (optional). Note"
"the [chr:] prefix should be added ONLY if alignment is in"
"vcf format. If this option is given, the vcf must also"
"be indexed with tabix"));
config.add(new ConfigParam<string>
("", "--subsites", "<subsites file>", &subsites_file,
"file listing NAMES from sites file (or sequences from"
" fasta) to keep; others will not be used. May be diploid"
" or haploid names (i.e., ind will have same effect as"
" ind_1 and ind_2)"));
config.add(new ConfigParam<string>
("-a", "--arg", "<SMC file>", &arg_file, "",
"initial ARG file (*.smc) for resampling (optional)"));
/* config.add(new ConfigParam<string>
("", "--cr", "<CR file>", &cr_file, "",
"initial ARGfile (*.cf) for resampling (optional)"));*/
config.add(new ConfigSwitch
("", "--unphased", &unphased,
"data is unphased (will integrate over phasings)."));
config.add(new ConfigSwitch
("", "--phased-vcf", &phased_vcf,
"Treat data as phased even with VCF input (warning: there is no"
" support to read phasing information from VCF. This option is"
" for special case where VCF input is for inbred individuals"
" and the redundant haplotype for each individiual is removed"
" with --subsites"));
config.add(new ConfigParam<string>
("", "--unphased-file", "<filename>", &unphased_file, "",
"use this file to identify haplotype pairs (file should"
" have two sequence names per line). By default, will"
" first check for the naming convention <ind>_1 <ind>_2 to"
" determine pairs. If the convention is not used, will assume"
" haplotype pairs are adjacent to each other in sequence"
" file. This option should not be used with VCF files."));
config.add(new ConfigParam<double>
("", "--randomize-phase", "<frac_random>", &randomize_phase,
0.0, "randomize phasings at start", ADVANCED_OPT));
config.add(new ConfigParam<double>
("", "--add-switch-errors", "<switch_error_rate>", &switch_error_rate,
0.0, "add switch errors at this rate", ADVANCED_OPT));
config.add(new ConfigSwitch
("", "--no-sample-phase", &no_sample_phase,
"Do not sample phase. Otherwise, if data is unphased, phase"
" will be sampled at same frequency as ARGs", ADVANCED_OPT));
config.add(new ConfigParam<int>
("", "--sample-phase", "<niters>", &sample_phase_step, 0,
"output phasings every <niters> samples. The default"
" is to sample at same frequency as ARGs (see --sample-step),"
" unless --no-sample-phase is given",
ADVANCED_OPT));
config.add(new ConfigParamComment("Masking options"));
config.add(new ConfigParam<string>
("", "--maskmap", "<sites mask>",
&maskmap, "",
"mask map file (optional)"));
config.add(new ConfigParam<string>
("", "--ind-maskmap", "<ind mask file>",
&ind_maskmap, "",
"File with two columns giving sample name and file with"
" mask for that sample. Sample name can refer to haploid or"
" diploid. If they are diploid, will append _1 and _2 to"
" names to search for relevant lineages to mask"));
config.add(new ConfigSwitch
("", "--expand-ind-mask",
&expand_ind_mask,
"This option is used to deal with Ns in a SITES file in the"
" case when the sites file was written by arg-sample run"
" with compression. In this case the N should apply to the"
" entire compressed site, and not a single base. Hint: to"
" run arg-sample and initialize the sites from a previous"
" run, use --sites <prev_outroot>[.rep].sites.gz"
" --maskmap <prev_outroot>.masked_regions.bed"
" --expand-ind-mask. The compression rate needs to be the"
" same as in the previous run. No other masking options are"
" necessary",
ADVANCED_OPT));
config.add(new ConfigParam<int>
("", "--mask-Ns", "<num>",
&maskN, -1,
"mask out any site in which at least num lineages have "
"missing data (this may increase efficiency)"));
config.add(new ConfigParam<string>
("", "--mask-cluster", "<numsnp,window>", &mask_cluster, "",
"mask any windows of size <window> which have at least"
" <numsnp> variants (at least two alleles other than N"
" must be observed at each site to count as a variant;"
" windows are determined after applying --subsites but"
" before any masking or compression is performed)"));
config.add(new ConfigParam<string>
("", "--vcf-genotype-filter", "<filter string>", &vcf_filter,
"String describing filtering for individual genotypes in VCF file."
" For example: \"GQ<5;DP<4;DP>30\" will mask genotypes with"
" GQ < 5 or DP < 4 or DP > 30. Currently the comparison operator"
" can only be < or >. All filtered genotypes will be masked."));
config.add(new ConfigParam<double>
("", "--vcf-min-qual", "<minQualScore>", &vcf_min_qual, 0.0,
"Minimum QUAL score for variants read from VCF. Others will be"
" masked. Default=0"));
config.add(new ConfigParam<double>
("", "--mask-uncertain", "<cutoff>", &mask_uncertain, 0.0,
"(for use --use-genotype-probs)"
" Mask any genotype if probability of most likely call < cutoff"
" according to PL score in VCF file", ADVANCED_OPT));
#ifdef ARGWEAVER_MPI
config.add(new ConfigSwitch
("", "--mpi", &mpi, "this is an mpi run, add <rank>.sites"
" to sites file name and <rank>. to out root, and"
" <rank>.smc.gz to --arg option (if given)"));
#endif
// model parameters
config.add(new ConfigParamComment("Model parameters"));
config.add(new ConfigParam<double>
("-m", "--mutrate", "<mutation rate>", &mu, 2.5e-8,
"mutations per site per generation (default=2.5e-8)"));
config.add(new ConfigParam<string>
("-M", "--mutmap", "<mutation rate map file>", &mutmap, "",
"mutation map file (optional)"));
config.add(new ConfigParam<double>
("-r", "--recombrate", "<recombination rate>", &rho, 1.5e-8,
"recombination per site per generation (default=1.5e-8)"));
config.add(new ConfigParam<string>
("-R", "--recombmap", "<recombination rate map file>",
&recombmap, "",
"recombination map file (optional)"));
config.add(new ConfigParam<int>
("-t", "--ntimes", "<ntimes>", &ntimes, 20,
"number of time points (default=20)"));
config.add(new ConfigParam<double>
("", "--maxtime", "<maxtime>", &maxtime, 200e3,
"maximum time point in generations (default=200e3)"));
config.add(new ConfigParam<double>
("", "--time-step", "<time>", &time_step, 0,
"linear time step in generations (optional)"));
config.add(new ConfigParam<double>
("", "--delta", "<delta>", &delta, 0.01,
"delta value for choosing log times (bigger value-> more"
" dense time points at leaves (default=0.01)"));
config.add(new ConfigParam<string>
("", "--popsize-file", "<popsize filename>", &popsize_file, "",
"file containing population sizes for each time span (optional)"));
config.add(new ConfigParam<string>
("", "--times-file", "<times filename>", ×_file, "",
"file containing time points (optional)"));
config.add(new ConfigParam<string>
("-N", "--popsize", "<population size>", &popsize_str,
"10000",
"effective population size (default=1e4)"));
config.add(new ConfigParam<string>
("", "--popsize-file", "<popsize filename>", &popsize_file, "",
"Two column file with columns time, popsize."
" The time indicates the time (in generations from present) that"
" a population changes to the given diploid size."
" The rows should be sorted in order of increasing time, with"
" the first row having t=0.\n"
" If using the multiple population model, a third colum indicates"
" the population number, with the first population numbered zero"));
config.add(new ConfigSwitch
("", "--smc-prime", &model.smc_prime,
"Use SMC' model (this option is generally slower"
" but represents a more accurate model). It is"
" recommended with the --pop-tree-file option, or for"
" ARGs with small numbers of highly diverged samples"));
// note that --invisible-recombs has value=true,
// so that providing the switch makes the variable false
config.add(new ConfigSwitch
("", "--invisible-recombs", &invisible_recombs,
"Output invisible recombinations which do not affect tree"
" topology (this will provide more accurate count of"
" recombination events, but may increase the runtime",
ADVANCED_OPT, false));
config.add(new ConfigSwitch
("", "--smc-orig", &do_nothing,
"This option is deprecated and does nothing; it is kept"
" for backwards compatibility. The original SMC model is the"
" default; see --smc-prime", ADVANCED_OPT));
// Population model options
config.add(new ConfigParam<string>
("", "--pop-file", "<population assigment file>", &pop_file,
"file assigning each haplotype to a population index",
POPMODEL_OPT));
config.add(new ConfigParam<string>
("-P", "--pop-tree-file", "<population file>", &pop_tree_file,
"", "File describing population tree (for multiple populations)"
" (Note --smc-prime is recommended with this option)",
POPMODEL_OPT));
config.add(new ConfigParam<int>
("", "--max-migs", "<max_migs>", &max_migrations,
1, "For use with --pop-tree-file, do not thread lineages with more"
" than this many migrations. Setting this to any value other"
" than one is experimental, and may significantly increase"
" run-time and cause issues with MCMC mixing",
POPMODEL_OPT));
config.add(new ConfigSwitch
("", "--no-resample-mig", &no_resample_mig,
"Do not perform migration-specific resampling",
POPMODEL_OPT));
config.add(new ConfigParam<int>
("", "--start-mig", "<start_mig>", &start_mig_iter,
0, "Start allowing migration events at this MCMC iteration"
" <start_mig>", POPMODEL_OPT));
// Population size sampling options
config.add(new ConfigParamComment("Population size sampling", EXPERIMENTAL_OPT));
config.add(new ConfigParam<int>
("", "--sample-popsize", "<num>", &sample_popsize_num, 0,
"sample population size using Hamiltonian Monte Carlo every"
"<num> threading operations (default=0 means do not sample)",
EXPERIMENTAL_OPT));
config.add(new ConfigParam<int>
("", "--popsize-em", "<n>", &popsize_em, 0,
"Do EM update of popsizes after every n threading operations",
EXPERIMENTAL_OPT));
config.add(new ConfigParam<double>
("", "--popsize-em-min-event", "<num>", &popsize_em_min_event,
2000.0,
"Minimum number of events per time interval; time intervals with"
" fewer events will be combined with previous time interval for"
" EM computations", EXPERIMENTAL_OPT));
config.add(new ConfigParam<int>
("", "--popsize-config", "<num>", &popsize_config, 0,
"Choose configuration for population sizes:\n"
" 0: constant population size across time (default)\n"
" 1: different population size in each population\n"
" 2: most fine-grained; different population size in each time\n"
" interval and population",
EXPERIMENTAL_OPT));
config.add(new ConfigSwitch
("", "--sample-popsize-const", &sample_popsize_const,
"update popsize but keep constant across times/populations",
EXPERIMENTAL_OPT));
config.add(new ConfigParam<string>
("", "--sample-popsize-config", "<popsize config file>",
&popsize_config_file, "",
"optional, for use with --sample-popsize.\n"
" Overrides popsize-config option, which provides some common configs\n"
" This file should contain up to 2-5 entries per line, in the format:"
" <param_name> <time_idx> <pop_idx=0> <optimize=1> <init_val=N>\n"
" The last three values are optional. All entries with the same\n"
" param_name will be constrained to have the same population size\n"
" Population sizes are initialized to the default value of N unless\n"
" specified here.", EXPERIMENTAL_OPT));
config.add(new ConfigParam<double>
("", "--epsilon", "<val>", &epsilon,
0.01, "(for use with --sample-popsize) epsilon value for"
"Hamiltonian population size updates", EXPERIMENTAL_OPT));
config.add(new ConfigParam<double>
("", "--pseudocount", "<val>", &pseudocount,
1.0, "(for use with --sample-popsize) gives weight to prior",
EXPERIMENTAL_OPT));
#ifdef ARGWEAVER_MPI
config.add(new ConfigParam<int>
("", "--mcmcmc", "<int>", &mcmcmc_numgroup,
1, "number of mcmcmc threads",
EXPERIMENTAL_OPT));
config.add(new ConfigParam<double>
("", "--mcmcmc-heat", "<val>", &mcmcmc_heat,
0.05, "heat interval for each thread in (MC)^3 group",
EXPERIMENTAL_OPT));
#endif
config.add(new ConfigSwitch
("", "--init-popsize-random", &init_popsize_random,
"(for use with --sample-popsize). Initialize each"
" population size parameter to a random number sampled"
" uniformly in [5000, 50000]", EXPERIMENTAL_OPT));
// sampling
config.add(new ConfigParamComment("Sampling"));
config.add(new ConfigParam<int>
("-n", "--iters", "<# of iterations>", &niters, 1000,
"(default=1000)"));
config.add(new ConfigParam<string>
("", "--resample-region", "<start>-<end>",
&resample_region_str, "",
"region to resample of input ARG (optional)"));
config.add(new ConfigSwitch
("", "--resume", &resume, "resume a previous run"));
config.add(new ConfigSwitch
("", "--overwrite", &overwrite,
"force an overwrite of a previous run"));
config.add(new ConfigSwitch
("", "--no-sample-arg", &no_sample_arg,
"Do not sample the ARG; hold at initial value "
"(usually for use with --sample-popsize)",
EXPERIMENTAL_OPT));
// misc
config.add(new ConfigParamComment("Miscellaneous"));
config.add(new ConfigParam<int>
("-c", "--compress-seq", "<compression factor>",
&compress_seq, 1,
"alignment compression factor (default=1)"));
config.add(new ConfigParam<int>
("", "--climb", "<# of climb iterations>", &nclimb, 0,
"(default=0)", ADVANCED_OPT));
config.add(new ConfigParam<int>
("", "--num-buildup", "<# of buildup iterations>", &num_buildup,
1, "(default=0)", ADVANCED_OPT));
config.add(new ConfigParam<int>
("", "--sample-step", "<sample step size>", &sample_step,
10, "number of iterations between steps (default=10)"));
config.add(new ConfigSwitch
("", "--no-compress-output", &no_compress_output,
"do not gzip output files"));
config.add(new ConfigParam<int>
("-x", "--randseed", "<random seed>", &randseed, 0,
"seed for random number generator (default=current time)"));
config.add(new ConfigSwitch
("", "--write-sites", &write_sites,
"Write sites (after all masking options) to <outroot>.sites.gz"));
config.add(new ConfigSwitch
("", "--write-sites-only", &write_sites_only,
"Same as --write-sites, but exit after reading and writing sites"));
config.add(new ConfigSwitch
("", "--include-masked-sites", &write_masked_sites,
"Relevant for --write-sites, --write-sites-only, or when sampling phase:"
" include sites that are masked in all samples in the output. By default"
" these are not included, but the set of masked sites is written to a"
" bed file named <outroot>.masked_sites.bed"));
// advanced options
config.add(new ConfigParamComment("Advanced Options", ADVANCED_OPT));
config.add(new ConfigSwitch
("", "--gibbs", &gibbs,
"use Gibbs sampling", ADVANCED_OPT));
config.add(new ConfigParam<double>
("", "--prob-path-switch", "<probability>",
&prob_path_switch, .1,
"removal path switch (default=.1)", ADVANCED_OPT));
config.add(new ConfigSwitch
("", "--infsites", &infsites,
"assume infinite sites model (at most one mutation per site)",
ADVANCED_OPT));
config.add(new ConfigParam<int>
("", "--resample-window", "<window size>",
&resample_window, 100000,
"sliding window for resampling (default=100000)",
ADVANCED_OPT));
config.add(new ConfigParam<int>
("", "--resample-window-iters", "<iterations>",
&resample_window_iters, 10,
"number of iterations per sliding window for resampling"
" (default=10)", ADVANCED_OPT));
// help information
config.add(new ConfigParamComment("Information"));
config.add(new ConfigParam<int>
("-V", "--verbose", "<verbosity level>",
&verbose, LOG_LOW,
"verbosity level 0=quiet, 1=low, 2=medium, 3=high"));
config.add(new ConfigSwitch
("-q", "--quiet", &quiet, "suppress logging to stderr"));
config.add(new ConfigSwitch
("-v", "--version", &version, "display version information"));
config.add(new ConfigSwitch
("-h", "--help", &help,
"display help information"));
config.add(new ConfigSwitch
("", "--help-advanced", &help_advanced,
"display help information about advanced options"));
config.add(new ConfigSwitch
("", "--help-experimental", &help_experimental,
"display help information about experimental options",
ADVANCED_OPT));
config.add(new ConfigSwitch
("-d", "--help-popmodel", &help_popmodel,
"display help information for population-based models (ARGweaver-D)"));
}
int parse_args(int argc, char **argv)
{
// parse arguments
if (!config.parse(argc, (const char**) argv)) {
if (argc < 2)
config.printHelp();
return EXIT_ERROR;
}
// display help
if (help) {
config.printHelp();
return EXIT_ERROR;
}
// display debug help
if (help_advanced) {
config.printHelp(stderr, ADVANCED_OPT);
return EXIT_ERROR;
}
if (help_popmodel) {
config.printHelp(stderr, POPMODEL_OPT);
return EXIT_ERROR;
}
if (help_experimental) {
config.printHelp(stderr, EXPERIMENTAL_OPT);
return EXIT_ERROR;
}
// display version info
if (version) {
printf(VERSION_INFO);
return EXIT_ERROR;
}
#ifdef ARGWEAVER_MPI
mcmcmc_group = 0;
int groupsize = MPI::COMM_WORLD.Get_size() / mcmcmc_numgroup;
mcmcmc_group = MPI::COMM_WORLD.Get_rank() / groupsize;
if (mcmcmc_group != 0) {
char tmp[1000];
sprintf(tmp, ".%i", mcmcmc_group);
mcmcmc_prefix = string(tmp);
}
printLog(LOG_LOW, "mcmcmc_prefix = %s\n", mcmcmc_prefix.c_str());
printLog(LOG_LOW, "mcmcmc_group=%i\n", mcmcmc_group);
#endif
return 0;
}
ConfigParser config;
// input/output
string fasta_file;
string sites_file;
string vcf_file;
string vcf_list_file;
string rename_file;
string vcf_filter;
double vcf_min_qual;
string subsites_file;
string out_prefix;
string arg_file;
// string cr_file;
string subregion_str;
string age_file;
string maskmap;
string ind_maskmap;
bool expand_ind_mask;
int maskN;
string mask_cluster;
string tabix_dir;
double mask_uncertain;
bool use_genotype_probs;
// model parameters
string popsize_str;
string pop_tree_file;
int max_migrations;
string pop_file;
double mu;
double rho;
int ntimes;
double maxtime;
double time_step;
double delta;
string popsize_file;
string times_file;
string mutmap;
string recombmap;
ArgModel model;
int popsize_em;
double popsize_em_min_event;
bool popsize_prior_neighbor;
bool init_popsize_random;
int popsize_config;
string popsize_config_file;
int sample_popsize_num;
bool sample_popsize_const;
bool invisible_recombs;
double epsilon;
double pseudocount;
#ifdef ARGWEAVER_MPI
double mcmcmc_heat;
int mcmcmc_group;
int mcmcmc_numgroup;
bool mpi;
#endif
string mcmcmc_prefix;
bool no_sample_arg;
bool no_resample_mig;
int start_mig_iter;
// search
int nclimb;
int num_buildup;
int niters;
string resample_region_str;
int resample_region[2];
bool resume;
bool overwrite;
string resume_stage;
int resume_iter;
int resample_window;
int resample_window_iters;
bool gibbs;
// misc
int compress_seq;
int sample_step;
bool no_compress_output;
int randseed;
double prob_path_switch;
bool infsites;
bool unphased;
string unphased_file;
bool phased_vcf;
double randomize_phase;
double switch_error_rate;
bool no_sample_phase;
int sample_phase_step;
bool all_masked;
bool write_sites;
bool write_sites_only;
bool write_masked_sites;
// help/information
bool quiet;
int verbose;
bool version;
bool help;
bool help_advanced;
bool help_experimental;
bool help_popmodel;
bool do_nothing;
// logging
FILE *stats_file;
};
bool parse_region(const char *region, int *start, int *end)
{
int start_idx=0;
for (int i=0; i < (int)strlen(region); i++)
if (region[i]==':') {
start_idx = i+1;
break;
}
return sscanf(®ion[start_idx], "%d-%d", start, end) == 2;
}
//=============================================================================
// logging
void set_up_logging(const Config &c, int level, const char *log_mode) {
Logger *logger;
/* printf("set_up_logging %s %s %s %i %i\n",
c.out_prefix.c_str(), c.mcmcmc_prefix.c_str(), LOG_SUFFIX, c.quiet,
level);*/
setLogLevel(level);
if (c.quiet) {
// log only to file
logger = &g_logger;
} else {
// log to both stdout and file
logger = new Logger(NULL, level);
g_logger.setChain(logger);
}
string log_filename = c.out_prefix + c.mcmcmc_prefix + LOG_SUFFIX;
if (!logger->openLogFile(log_filename.c_str(), log_mode)) {
printError("Could not open log file '%s'", log_filename.c_str());
abort();
}
}
// log the program version and start time
void log_intro(int level)
{
time_t t = time(NULL);
printLog(level, "arg-sample " VERSION_TEXT "\n");
printLog(level, "start time: %s", ctime(&t)); // newline is in ctime
}
// log the command used
void log_prog_commands(int level, int argc, char **argv)
{
printLog(level, "command:");
for (int i=0; i<argc; i++) {
printLog(level, " %s", argv[i]);
}
printLog(level, "\n");
}
//=============================================================================
// statistics output
void print_stats_header(Config *config) {
fprintf(config->stats_file, "stage\titer\tprior\tprior2\tlikelihood\tjoint\t"
"recombs\tnoncompats\targlen");
if (config->invisible_recombs)
fprintf(config->stats_file, "\tinvis_recombs");
if (config->model.popsize_config.sample) {
list<PopsizeConfigParam> l = config->model.popsize_config.params;
for (list<PopsizeConfigParam>::iterator it=l.begin();
it != l.end(); ++it) {
fprintf(config->stats_file, "\t%s", it->name.c_str());
}
} else if (config->popsize_em) {
for (int pop=0; pop < config->model.num_pops(); pop++) {
for (int i=0; i < config->model.ntimes-1; i++) {
char str[100];
if (config->model.num_pops() == 1)
sprintf(str, "N%i", i);
else sprintf(str, "N%i.%i", pop, i);
fprintf(config->stats_file, "\t%s", str);
}
}
}
if (config->model.pop_tree != NULL) {
for (unsigned int i=0; i < config->model.pop_tree->mig_params.size(); i++)
fprintf(config->stats_file, "\t%s", config->model.pop_tree->mig_params[i].name.c_str());
}
fprintf(config->stats_file, "\n");
}
void print_stats(FILE *stats_file, const char *stage, int iter,
ArgModel *model,
const Sequences *sequences, LocalTrees *trees,
const SitesMapping* sites_mapping, const Config *config,
const TrackNullValue *maskmap_uncompressed,
const vector<int> &invisible_recomb_pos0=vector<int>(),
const vector<Spr> &invisible_recombs=vector<Spr>())
{
// calculate number of recombinations
int nrecombs = trees->get_num_trees() - 1;
// calculate number of non-compatiable sites
int nseqs = trees->get_num_leaves();
char *seqs[nseqs];
for (int i=0; i<nseqs; i++)
seqs[i] = sequences->seqs[trees->seqids[i]];
int noncompats = count_noncompat(trees, seqs, nseqs, sequences->length());
// get memory usage in MB
double maxrss = get_max_memory_usage() / 1000.0;
// calculate likelihood, prior, and joint probabilities
// uncompressed local trees
vector<int> invisible_recomb_pos;
if (sites_mapping) {
uncompress_local_trees(trees, sites_mapping);
uncompress_model(model, sites_mapping, config->compress_seq);
sites_mapping->uncompress(invisible_recomb_pos0, invisible_recomb_pos);
}
double prior = calc_arg_prior(model, trees, NULL, NULL, -1, -1,
invisible_recomb_pos, invisible_recombs);
double prior2 = calc_arg_prior_recomb_integrate(model, trees, NULL, NULL, NULL);
double likelihood = config->all_masked ? 0.0 :
calc_arg_likelihood(model, sequences, trees,
sites_mapping,
maskmap_uncompressed);
double joint = prior + likelihood;
double arglen = get_arglen(trees, model->times);
// recompress local trees
if (sites_mapping) {
compress_local_trees(trees, sites_mapping);
compress_model(model, sites_mapping, config->compress_seq);
}
// output stats
fprintf(stats_file, "%s\t%d\t%f\t%f\t%f\t%f\t%d\t%d\t%f",
stage, iter,
prior, prior2, likelihood, joint, nrecombs, noncompats, arglen);
if (config->invisible_recombs)
fprintf(stats_file, "\t%i", (int)invisible_recomb_pos.size());
if (model->popsize_config.sample) {
list<PopsizeConfigParam> l=model->popsize_config.params;
for (list<PopsizeConfigParam>::iterator it=l.begin();
it != l.end(); ++it) {
set<PopTime>::iterator it2 = it->intervals.begin();
fprintf(stats_file, "\t%f", model->popsizes[it2->pop][it2->time]);
it2++;
// just checking here
set<PopTime>::iterator it3 = it2;
while (it3 != it->intervals.end()) {
assert(model->popsizes[it2->pop][it2->time] ==
model->popsizes[it3->pop][it3->time]);
it3++;
}
}
} else if (config->popsize_em) {
for (int pop=0; pop < model->num_pops(); pop++)
for (int i=0; i < model->ntimes-1; i++)
fprintf(stats_file, "\t%f", model->popsizes[pop][2*i]);
}
if (model->pop_tree != NULL) {
for (unsigned int i=0; i < model->pop_tree->mig_params.size(); i++) {
MigParam mp = model->pop_tree->mig_params[i];
fprintf(stats_file, "\t%e",
model->pop_tree->mig_matrix[mp.time_idx].get(mp.from_pop, mp.to_pop));
}
}
fprintf(stats_file, "\n");
fflush(stats_file);
printLog(LOG_LOW, "\n"
"prior: %f\n"
"likelihood: %f\n"
"joint: %f\n"
"nrecombs: %d\n"
"noncompats: %d\n"
"arglen: %f\n"
"max memory: %.1f MB\n\n",
prior, likelihood, joint, nrecombs, noncompats, arglen, maxrss);
}
//=============================================================================
// sample output
// Returns the iteration-specific ARG filename
string get_out_arg_file(const Config &config, int iter)
{
char iterstr[10];
snprintf(iterstr, 10, ".%d", iter);
return config.out_prefix + config.mcmcmc_prefix + iterstr + SMC_SUFFIX;
}
/*string get_out_cr_file(const Config &config, int iter)
{
char iterstr[10];
snprintf(iterstr, 10, ".%d", iter);
return config.out_prefix + config.mcmcmc_prefix + iterstr + COAL_RECORDS_SUFFIX;
}*/
string get_out_sites_file(const Config &config, int iter)
{
char iterstr[12];
if (iter < 0) {
iterstr[0]='\0';
} else {
snprintf(iterstr, 12, ".%d", iter);
}
string sitesfile = config.out_prefix + config.mcmcmc_prefix + iterstr + SITES_SUFFIX;
if (!config.no_compress_output)
sitesfile = sitesfile + ".gz";
return sitesfile;
}
bool log_sequences(string chrom, const Sequences *sequences,
const Config *config,
const SitesMapping *sites_mapping, int iter) {
Sites sites(chrom);
string out_sites_file = get_out_sites_file(*config, iter);
make_sites_from_sequences(sequences, &sites);
if (sites_mapping)
uncompress_sites(&sites, sites_mapping);
CompressStream stream(out_sites_file.c_str(), "w");
if (!stream.stream) {
printError("cannot write '%s'", out_sites_file.c_str());
return false;
}
write_sites(stream.stream, &sites, config->write_masked_sites);
return true;
}
bool log_local_trees(const ArgModel *model, const Sequences *sequences,
LocalTrees *trees, const SitesMapping* sites_mapping,
const Config *config, int iter,
const vector<int> &self_recomb_pos0=vector<int>(),
const vector<Spr> &self_recombs=vector<Spr>())
{
string out_arg_file = get_out_arg_file(*config, iter);
const vector<int> *self_recomb_ptr;
if (!config->no_compress_output)
out_arg_file += ".gz";
// write local trees uncompressed
vector<int> self_recomb_pos1;
if (sites_mapping) {
uncompress_local_trees(trees, sites_mapping);
sites_mapping->uncompress(self_recomb_pos0, self_recomb_pos1);
self_recomb_ptr = &self_recomb_pos1;
} else self_recomb_ptr = &self_recomb_pos0;
// setup output stream
CompressStream stream(out_arg_file.c_str(), "w");
if (!stream.stream) {
printError("cannot write '%s'", out_arg_file.c_str());
return false;
}
write_local_trees(stream.stream, trees, sequences, model->times,
model->pop_tree != NULL,
*self_recomb_ptr, self_recombs);
// testing for now; output coal records version
/* string out_cr_file = get_out_cr_file(*config, iter);
if (!config->no_compress_output)
out_cr_file += ".gz";
CompressStream stream2(out_cr_file.c_str(), "w");
if (!stream2.stream) {
printError("cannot write '%s'", out_cr_file.c_str());
return false;
}
write_coal_records(stream2.stream, model, trees, sequences); */
if (sites_mapping)
compress_local_trees(trees, sites_mapping);
return true;
}
//=============================================================================
// Read a list of doubles from a file.
bool read_doubles(const char *filename, vector<double> &values) {
FILE *infile = fopen(filename, "r");
if (infile == NULL) {
return false;
}
// Read all doubles and sort them.
double value;
while (fscanf(infile, "%lf", &value) == 1)
values.push_back(value);
fclose(infile);
return true;
}
/*bool read_init_arg_cr(const char *cr_file, const ArgModel *model,
LocalTrees *trees, vector<string> &seqnames)
{
CompressStream stream(cr_file, "r");
if (!stream.stream) {
printError("cannot read '%s'", cr_file);
return false;
}
return read_coal_records(stream.stream, model, trees, seqnames);
}*/
bool read_init_arg(const char *arg_file, const ArgModel *model,
LocalTrees *trees, vector<string> &seqnames)
{
CompressStream stream(arg_file, "r");
if (!stream.stream) {
printError("cannot read '%s'", arg_file);
return false;
}
return read_local_trees(stream.stream, model->times, model->ntimes,
trees, seqnames);
}