minor fixes to 0.1.8

Author: Yun Deng

releases/singer-0.1.8-beta-linux-x86_64.tar.gz

@@ -5,21 +5,26 @@ import os
 import subprocess
 import argparse
 import random
-import uuid
+# import uuid
+import numpy as np
+# import pysam
 
 def random_string():
     unique_id = uuid.uuid4().hex
     return unique_id
 
+'''
 def write_rate_map(output_prefix, rate, x, y):
     rate_table = np.array([[0, y, rate]])
     rand_suffix = random_string()
     map_filename = f"{output_prefix}_{rand_suffix}.txt"
+    print(rate_table)
     np.savetxt(map_filename, rate_table, delimiter=' ')
     return map_filename
+'''
 
 def compute_average_rate(map_file, x, y):
-    data = np.loadtxt(map_file)
+    data = np.loadtxt(map_file, ndmin=2)
     mask = (data[:, 1] > x) & (data[:, 0] < y)
     selected_data = data[mask]
 
@@ -38,28 +43,21 @@ def compute_average_rate(map_file, x, y):
     average_rate = weighted_sum / total_length
     return average_rate
 
-def tabix_vcf(vcf_file):
-    index_file = vcf_file + '.tbi'
-    if not os.path.exists(index_file):
-        print(f"Indexing VCF file: {vcf_file}")
-        subprocess.run(['tabix', '-p', 'vcf', vcf_file], check=True)
-    else:
-        print(f"Index file already exists: {index_file}")
-
+'''
 def calculate_diversity_with_filter(vcf_file, start, end):
-    tabix_vcf(vcf_file)
+    # Open the VCF file using pysam
     vcf = pysam.VariantFile(vcf_file)
-    try:
-        first_variant = next(vcf.fetch())
-        chrom = first_variant.chrom
-    except StopIteration:
-        raise ValueError("The VCF file is empty or does not contain any variants.")
-
+    
     num_variants = 0
     total_diversity = 0.0
     snp_positions = []
-
-    for record in vcf.fetch(chrom, start, end):
+    
+    # Iterate through all records and filter based on position
+    for record in vcf:
+        if record.pos < start:
+            continue
+        if record.pos > end:
+            break
         snp_positions.append(record.pos)
         allele_counts = [0] * len(record.alleles)
         num_samples = 0
@@ -80,34 +78,97 @@ def calculate_diversity_with_filter(vcf_file, start, end):
         total_diversity += variant_diversity
         num_variants += 1
 
+    # Calculate SNP distances
     snp_distances = np.diff(snp_positions)
     
     if len(snp_distances) == 0:
         return 0.0
-    snp_distances = np.append(snp_positions[0] - start, snp_distances, )
+    
+    # Append the distances to the start and end positions
+    snp_distances = np.append(snp_positions[0] - start, snp_distances)
     snp_distances = np.append(snp_distances, end - snp_positions[-1])
+    
+    # Sort the distances and filter out the longest 5 distances
     sorted_snp_distances = np.sort(snp_distances)
-    filtered_length = np.sum(sorted_snp_distances[0:-5])
+    filtered_length = np.sum(sorted_snp_distances[:-5])
+    
+    # Calculate nucleotide diversity
     nucleotide_diversity = total_diversity / filtered_length
+    
+    return nucleotide_diversity
+'''
+
+def calculate_diversity(vcf_file, start, end):
+    total_diversity = 0.0
+    segment_length = end - start
+    num_variants = 0
+
+    # Open the VCF file and read it line by line
+    with open(vcf_file, 'r') as vcf:
+        for line in vcf:
+            # Skip header lines
+            if line.startswith("#"):
+                continue
+
+            # Split the VCF record into fields
+            fields = line.strip().split('\t')
+            pos = int(fields[1])  # Get the position of the variant
+            
+            # Check if the variant is within the target region
+            if pos < start:
+                continue
+            if pos > end:
+                break
+            
+            # Extract genotype information for each sample (after the 9th column)
+            alleles = [fields[3]] + fields[4].split(",")  # REF and ALT alleles
+            allele_counts = [0] * len(alleles)
+            num_samples = 0
+
+            # Iterate over the genotype columns to count allele occurrences
+            for sample_field in fields[9:]:
+                gt = sample_field.split(':')[0]  # Extract the genotype (e.g., "0/1")
+                if '.' in gt:  # Missing data
+                    continue
+                
+                # Genotype fields are usually formatted like "0/1", "1/1", etc.
+                gt_alleles = gt.replace('|', '/').split('/')  # Handle phased/unphased
+                num_samples += 1
+                for allele in gt_alleles:
+                    allele_counts[int(allele)] += 1
+
+            # Calculate allele frequencies
+            if num_samples == 0:
+                continue  # Skip if there are no valid samples
+            allele_freqs = [ac / (2 * num_samples) for ac in allele_counts]
+            
+            # Calculate diversity for the current variant
+            variant_diversity = sum(f * (1 - f) for f in allele_freqs if f > 0 and f < 1)
+            
+            # Add to total diversity
+            total_diversity += variant_diversity
+            num_variants += 1
+
+    nucleotide_diversity = total_diversity / segment_length
     return nucleotide_diversity
 
 def compute_Ne(vcf_prefix, start, end, m):
-    vcf_filename = vcf + ".vcf"
-    effective_diversity = calculate_diversity_with_filter(vcf_filename, start, end)
-    Ne = effective_diversity/4/m
+    vcf_filename = vcf_prefix + ".vcf"
+    diversity = calculate_diversity(vcf_filename, start, end)
+    Ne = diversity/4/m
     return Ne 
 
-def run_singer(Ne, mut_map_filename, recomb_map_filename, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
+def run_rate_singer(Ne, m, r, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
     attempts = 0
     max_attempts = 100
     random.seed(seed)
-    random_seeds = [random.randint(0, 2**30 - 1) for _ in range(max_attempts)]   
-   
+    random_seeds = [random.randint(0, 2**30 - 1) for _ in range(max_attempts)]
+
     script_dir = os.path.dirname(os.path.realpath(__file__))
-    singer_executable = os.path.join(script_dir, "singer") 
-    
-    start_cmd = f"{singer_executable} -Ne {Ne} -mut_map {mut_map_filename} -recomb_map {recomb_map_filename} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin}"    
-    debug_cmd = f"{singer_executable} -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"    
+    singer_executable = os.path.join(script_dir, "singer")
+
+    start_cmd = f"{singer_executable} -Ne {Ne} -m {m} -r {r} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin}"
+    debug_cmd = f"{singer_executable} -Ne {Ne} -m {m} -r {r} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"
 
     seeded_start_cmd = start_cmd + f" -seed {seed}"
     print(seeded_start_cmd)
@@ -123,7 +184,7 @@ def run_singer(Ne, mut_map_filename, recomb_map_filename, start, end, vcf_prefix
         print("Auto-debug failed. Contact the author for help: [email protected]")
         sys.exit(1)
 
-def run_haps_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start, end, haps_prefix, output_prefix, num_iters, thin, polar, seed):
+def resume_rate_singer(Ne, m, r, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
     attempts = 0
     max_attempts = 100
     random.seed(seed)
@@ -132,12 +193,12 @@ def run_haps_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start, e
     script_dir = os.path.dirname(os.path.realpath(__file__))
     singer_executable = os.path.join(script_dir, "singer")
 
-    start_cmd = f"{singer_executable} -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -haps {haps_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin}"
-    debug_cmd = f"{singer_executable} -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -haps {haps_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"
+    resume_cmd = f"{singer_executable} -Ne {Ne} -m {m} -r {r} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -resume"
+    debug_cmd = f"{singer_executable} -Ne {Ne} -m {m} -r {r} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"
 
-    seeded_start_cmd = start_cmd + f" -seed {seed}"
-    print(seeded_start_cmd)
-    process = subprocess.run(seeded_start_cmd.split(), check=False)
+    seeded_resume_cmd = resume_cmd + f" -seed {seed}"
+    print(seeded_resume_cmd)
+    process = subprocess.run(seeded_resume_cmd.split(), check=False)
     while process.returncode != 0 and attempts < max_attempts:
         print(f"Auto-debug iteration: {attempts}")
         seeded_debug_cmd = debug_cmd + f" -seed {random_seeds[attempts]}"
@@ -149,21 +210,21 @@ def run_haps_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start, e
         print("Auto-debug failed. Contact the author for help: [email protected]")
         sys.exit(1)
 
-def resume_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
+def run_map_singer(Ne, mut_map_filename, recomb_map_filename, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
     attempts = 0
     max_attempts = 100
     random.seed(seed)
-    random_seeds = [random.randint(0, 2**30 - 1) for _ in range(max_attempts)]
-
+    random_seeds = [random.randint(0, 2**30 - 1) for _ in range(max_attempts)]   
+   
     script_dir = os.path.dirname(os.path.realpath(__file__))
-    singer_executable = os.path.join(script_dir, "singer")
+    singer_executable = os.path.join(script_dir, "singer") 
+    
+    start_cmd = f"{singer_executable} -Ne {Ne} -mut_map {mut_map_filename} -recomb_map {recomb_map_filename} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin}"    
+    debug_cmd = f"{singer_executable} -Ne {Ne} -mut_map {mut_map_filename} -recomb_map {recomb_map_filename} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"    
 
-    resume_cmd = f"{singer_executable} -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -resume" 
-    debug_cmd = f"{singer_executable} -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"
-      
-    seeded_resume_cmd = resume_cmd + f" -seed {seed}"
-    print(seeded_resume_cmd)
-    process = subprocess.run(seeded_resume_cmd.split(), check=False)
+    seeded_start_cmd = start_cmd + f" -seed {seed}"
+    print(seeded_start_cmd)
+    process = subprocess.run(seeded_start_cmd.split(), check=False)
     while process.returncode != 0 and attempts < max_attempts:
         print(f"Auto-debug iteration: {attempts}")
         seeded_debug_cmd = debug_cmd + f" -seed {random_seeds[attempts]}"
@@ -175,8 +236,8 @@ def resume_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start, end
         print("Auto-debug failed. Contact the author for help: [email protected]")
         sys.exit(1)
 
-      
-def resume_fast_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
+
+def resume_map_singer(Ne, mut_map_filename, recomb_map_filename, start, end, vcf_prefix, output_prefix, num_iters, thin, polar, seed):
     attempts = 0
     max_attempts = 100
     random.seed(seed)
@@ -185,9 +246,9 @@ def resume_fast_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start
     script_dir = os.path.dirname(os.path.realpath(__file__))
     singer_executable = os.path.join(script_dir, "singer")
 
-    resume_cmd = f"{singer_executable} -fast -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -resume"
-    debug_cmd = f"{singer_executable} -fast -Ne {Ne} -m {mutation_rate} -r {mutation_rate * recombination_to_mutation_ratio} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"
-
+    resume_cmd = f"{singer_executable} -Ne {Ne} -mut_map {mut_map_filename} -recomb_map {recomb_map_filename} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -resume" 
+    debug_cmd = f"{singer_executable} -Ne {Ne} -mut_map {mut_map_filename} -recomb_filename {recomb_map_filename} -input {vcf_prefix} -output {output_prefix} -start {start} -end {end} -polar {polar} -n {num_iters} -thin {thin} -debug"
+      
     seeded_resume_cmd = resume_cmd + f" -seed {seed}"
     print(seeded_resume_cmd)
     process = subprocess.run(seeded_resume_cmd.split(), check=False)
@@ -202,19 +263,19 @@ def resume_fast_singer(Ne, mutation_rate, recombination_to_mutation_ratio, start
         print("Auto-debug failed. Contact the author for help: [email protected]")
         sys.exit(1)
 
-
+      
 def main():
     parser = argparse.ArgumentParser(description='Sample and infer ARG from genetic variation data')
 
     parser.add_argument('-Ne', type=float, default=-1, help='Effective population size.')
-    parser.add_argument('-m', type=float, required=True, help='Mutation rate.')
+    parser.add_argument('-m', type=float, help='Mutation rate.')
     parser.add_argument('-ratio', type=float, default=1, help='Recombination to mutation ratio. Default: 1.')
     parser.add_argument('-recomb_map', type=str, help='Filename for the recombination rate map.')
     parser.add_argument('-mut_map', type=str, help='Filename for the mutation rate map.')
     parser.add_argument('-vcf', type=str, help='VCF file prefix (without .vcf extension).')
     parser.add_argument('-output', type=str, required=True, help='Output file prefix.')
-    parser.add_argument('-start', type=str, required=True, help='Start position.')
-    parser.add_argument('-end', type=str, required=True, help='End position.')
+    parser.add_argument('-start', type=float, required=True, help='Start position.')
+    parser.add_argument('-end', type=float, required=True, help='End position.')
     parser.add_argument('-n', type=int, default=100, help='Number of MCMC samples. Default: 100.')
     parser.add_argument('-thin', type=int, default=20, help='Thinning interval length. Default: 20.')
     parser.add_argument('-polar', type=float, default=0.5, help='Site flip probability. Default: 0.5.')
@@ -230,28 +291,31 @@ def main():
     if not args.m and not args.mut_map:
         parser.error("You must provide either -m or -mut_map")
 
-    if not args.ratio and not args.recomb_map:
-        parser.error("You must provide either -ratio or -recomb_map") 
-
 
     if args.m:
-        args.mut_map = write_rate_map(args.m, args.output, args.start, args.end)
-
-    if args.ratio:
-        if not args.m:
-            parser.error("You must provide -m when using the -ratio flag")
+        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:
-            args.recomb_map = write_rate_map(args.ratio*args.m, args.vcf, args.start, args.end, args.output)
-    
-    if not args.Ne:
-        args.m = compute_average_rate(args.mut_map, args.start, args.end)
-        args.Ne = compute_Ne(args.vcf, args.start, args.end, args.m)
-
-    if (args.resume):
-        resume_singer(args.Ne, args.mut_map, args.recomb_map, args.start, args.end, args.vcf, args.output, args.n, args.thin, args.polar, args.seed)
-    else:
-        run_singer(args.Ne, args.mut_map, args.recomb_map, args.start, args.end, args.vcf, args.output, args.n, args.thin, args.polar, args.seed)
+            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 
+
+    else if args.mut_map and not args.m:
+        args.m = compute_average_rate(args.mut_map, args.start, args.end) 
+        if not args.recomb_map:
+            parser.error("You must provide -recomb_map when using the -mut_map flag")
 
+        if args.Ne < 0:
+            args.m = compute_average_rate(args.mut_map, args.start, args.end)
+            args.Ne = compute_Ne(args.vcf, args.start, args.end, args.m)
+
+        if (args.resume):
+            resume_map_singer(args.Ne, args.mut_map, args.recomb_map, args.start, args.end, args.vcf, args.output, args.n, args.thin, args.polar, args.seed)
+        else: 
+            run_map_singer(args.Ne, args.mut_map, args.recomb_map, args.start, args.end, args.vcf, args.output, args.n, args.thin, args.polar, args.seed)
+    return
+
 
 if __name__ == "__main__":
     main()