Skip to content

Generalize recombination map support to all scripts #24

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dcdehaas wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Generalize recombination map support to all scripts #24

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
31a0f2a
Add -recomb_map option for parallel_singer
dcdatcu Sep 23, 2024
832f759
Rate_map has to take the -start option into account
dcdehaas Nov 7, 2024
8f55d02
"Support" regions of 0 recombination; Hg38 has a few of these
dcdehaas Nov 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions SINGER/SINGER/Sampler.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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) {
Copy link
Author

@dcdehaas dcdehaas Oct 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this was always true, because of lines 20 and 21 above.

if (!has_map) {
arg.compute_rhos_thetas(recomb_rate, mut_rate);
} else {
arg.compute_rhos_thetas(recomb_map, mut_map);
Expand Down Expand Up @@ -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);
Expand Down
1 change: 1 addition & 0 deletions SINGER/SINGER/Sampler.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down
17 changes: 6 additions & 11 deletions SINGER/SINGER/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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();
Expand Down
32 changes: 23 additions & 9 deletions SINGER/SINGER/parallel_singer
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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}")
Expand All @@ -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__":
Expand Down
7 changes: 3 additions & 4 deletions SINGER/SINGER/singer_master
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand Down