Add -recomb_map option for parallel_singer

singer_master already supports it.

Author: dcdatcu

SINGER/SINGER/Sampler.cpp

@@ -21,6 +21,7 @@ Sampler::Sampler(double pop_size, Rate_map &rm, Rate_map &mm) {
     recomb_rate = rm.mean_rate()*Ne;
     recomb_map = rm;
     mut_map = mm;
+    has_map = true;
 }
 
 void Sampler::set_precision(double c, double q) {
@@ -306,7 +307,7 @@ void Sampler::build_singleton_arg() {
     arg = ARG(Ne, sequence_length);
     arg.discretize(bin_size);
     arg.build_singleton_arg(n);
-    if (mut_rate > 0 and recomb_rate > 0) {
+    if (!has_map) {
         arg.compute_rhos_thetas(recomb_rate, mut_rate);
     } else {
         arg.compute_rhos_thetas(recomb_map, mut_map);
@@ -787,7 +788,7 @@ void Sampler::load_resume_arg() {
     coord_file = output_prefix + "_coordinates.txt";
     arg.read(node_file, branch_file, recomb_file, mut_file);
     arg.read_coordinates(coord_file);
-    if (mut_rate > 0 and recomb_rate > 0) {
+    if (!has_map) {
         arg.compute_rhos_thetas(recomb_rate, mut_rate);
     } else {
         arg.compute_rhos_thetas(recomb_map, mut_map);

SINGER/SINGER/Sampler.hpp

@@ -48,6 +48,7 @@ class Sampler {
     vector<Node_ptr> ordered_sample_nodes = {};
     unordered_map<double, set<Node_ptr>> carriers = {};
     unordered_map<Node_ptr, set<double>> mutation_sets = {};
+    bool has_map = false;
 
     Sampler();
 

SINGER/SINGER/main.cpp

@@ -224,16 +224,11 @@ int main(int argc, const char * argv[]) {
             exit(1);
         }
     }
-    if (r < 0) {
-        cerr << "-r flag missing or invalid value. " << endl;
-        exit(1);
-    }
-    if (m < 0) {
-        cerr << "-m flag missing or invalid value. " << endl;
-        exit(1);
-    }
-    if (Ne < 0) {
-        cerr << "-Ne flag missing or invalid value. " << endl;
+    bool bothRates = (r > 0 && m > 0);
+    bool bothMaps = !(recomb_map_filename.empty() || mut_map_filename.empty());
+    if (!(bothRates || bothMaps)) {
+        cerr << "Must provide either valid values for -r and -m, OR provide a valid "
+             << "mutation map and recombination map. " << endl;
         exit(1);
     }
     if (input_filename.size() == 0) {
@@ -253,7 +248,7 @@ int main(int argc, const char * argv[]) {
         exit(1);
     }
     Sampler sampler;
-    if (r > 0 and m > 0) {
+    if (bothRates) {
         sampler = Sampler(Ne, r, m);
     } else {
         Rate_map recomb_map = Rate_map();

SINGER/SINGER/parallel_singer

@@ -14,7 +14,7 @@ def index_vcf(vcf_prefix, block_length):
     subprocess.run(["python", indexer, vcf_prefix, str(block_length)])
 
 
-def run_singer_in_parallel(vcf_prefix, output_prefix, mutation_rate, ratio, block_length, num_iters, thinning_interval, Ne, polar, num_cores):
+def run_singer_in_parallel(vcf_prefix, output_prefix, mutation_rate, ratio, block_length, num_iters, thinning_interval, Ne, polar, num_cores, recomb_map, mut_map):
     """Run singer in parallel using the specified parameters."""
 
     # Read the breakpoints from the index file
@@ -30,8 +30,18 @@ def run_singer_in_parallel(vcf_prefix, output_prefix, mutation_rate, ratio, bloc
         start = breakpoints[i]
 
         # Base command
-        cmd = f"{singer_master_executable} -Ne {Ne} -m {mutation_rate} -ratio {ratio} -vcf {vcf_prefix} -output {output_prefix}_{i}_{i+1} -start {start} -end {start + block_length} -n {num_iters} -thin {thinning_interval} -polar {polar}"
-        
+        cmd = f"{singer_master_executable} -vcf {vcf_prefix} -output {output_prefix}_{i}_{i+1} -start {start} -end {start + block_length} -n {num_iters} -thin {thinning_interval} -polar {polar}"
+        if Ne is not None:
+            cmd += f" -Ne {Ne}"
+        if recomb_map:
+            cmd += f" -recomb_map \"{recomb_map}\""
+        else:
+            cmd += f" -ratio {ratio}"
+        if mut_map is not None:
+            cmd += f" -mut_map {mut_map}"
+        else:
+            cmd += f" -m {mutation_rate}"
+
         cmd_list.append(cmd)
 
     # Execute the commands in parallel
@@ -47,28 +57,32 @@ def convert_long_ARG(vcf_prefix, output_prefix, num_iters, freq):
 
 def main():
     parser = argparse.ArgumentParser(description="Parallelize singer runs by cutting the chromosome")
-    parser.add_argument('-Ne', type=float, required=True, help='Effective population size. Default: 1e4.')
-    parser.add_argument("-m", type=float, required=True, help="Mutation rate.")
+    parser.add_argument('-Ne', type=float, required=False, help='Effective population size. Default: 1e4.')
+    parser.add_argument("-m", type=float, required=False, help="Mutation rate.")
+    parser.add_argument("-mut_map", type=str, help="Filename for the mutation rate map.")
     parser.add_argument("-ratio", type=float, default=1, help="Recombination to mutation ratio. Default: 1.")
+    parser.add_argument("-recomb_map", default=None, required=False, help="Filename for the recombination rate map.")
     parser.add_argument("-L", type=int, default=int(1e6), help="Block length. Default: 1e6.")
     parser.add_argument("-vcf", type=str, required=True, help="VCF file prefix (without .vcf or .vcf.gz extension).")
     parser.add_argument("-output", type=str, required=True, help="Output file prefix.")
     parser.add_argument("-n", type=int, required=True, help="Number of MCMC samples.")
     parser.add_argument("-thin", type=int, required=True, help="Thinning interval length.")
     parser.add_argument("-polar", type=float, default=0.5, required=False, help="Site flip probability. Default: 0.5.")
     parser.add_argument("-freq", type=float, default=1, required=False, help="Convert to tskit every {freq} samples. Default: 1.")
-    parser.add_argument("-num_cores", type=int, default=20, required=False, help="Number of cores. Default: 20.")    
-        
+    parser.add_argument("-num_cores", type=int, default=20, required=False, help="Number of cores. Default: 20.")
+
     if len(sys.argv) == 1:
         parser.print_help(sys.stderr)
         sys.exit(1)
 
     args = parser.parse_args()
 
     print("Parameters:")
-    print(f"Effective population size: {args.Ne}")
+    print(f"Effective population size: {args.Ne if args.Ne is not None else 'estimated'}")
     print(f"Mutation rate: {args.m}")
+    print(f"Mutation map: {args.mut_map}")
     print(f"Recombination to mutation ratio: {args.ratio}")
+    print(f"Recombination map: {args.recomb_map}")
     print(f"Block length: {args.L}")
     print(f"VCF file prefix: {args.vcf}")
     print(f"Output file prefix: {args.output}")
@@ -79,7 +93,7 @@ def main():
     print(f"Number of cores: {args.num_cores}") 
 
     index_vcf(args.vcf, args.L)
-    run_singer_in_parallel(args.vcf, args.output, args.m, args.ratio, args.L, args.n, args.thin, args.Ne, args.polar, args.num_cores)
+    run_singer_in_parallel(args.vcf, args.output, args.m, args.ratio, args.L, args.n, args.thin, args.Ne, args.polar, args.num_cores, args.recomb_map, args.mut_map)
     convert_long_ARG(args.vcf, args.output, args.n, args.freq)
 
 if __name__ == "__main__":

SINGER/SINGER/singer_master

@@ -293,16 +293,15 @@ def main():
 
 
     if args.m:
+        if args.recomb_map:
+            parser.error("-recomb_map is not compatible with -m; you must use -mut_map")
         if args.Ne < 0:
             args.Ne = compute_Ne(args.vcf, args.start, args.end, args.m) 
         if (args.resume):
             resume_rate_singer(args.Ne, args.m, args.m*args.ratio, args.start, args.end, args.vcf, args.output, args.n, args.thin, args.polar, args.seed)
         else:
             run_rate_singer(args.Ne, args.m, args.m*args.ratio, args.start, args.end, args.vcf, args.output, args.n, args.thin, args.polar, args.seed) 
-    return 
-
-    if args.mut_map and not args.m:
-        args.m = compute_average_rate(args.mut_map, args.start, args.end) 
+    elif args.mut_map:
         if not args.recomb_map:
             parser.error("You must provide -recomb_map when using the -mut_map flag")