Fix bugs

docs/source/usage.rst

@@ -34,6 +34,9 @@ Input Arguments
                           The last marker to consider. Default: all markers considered.
       -main_metafounder MAIN_METAFOUNDER
                           The metafounder to use where parents are unknown with input "0". Default: MF_1.
+      -map_file MAP_FILE  
+                          A map file for all loci.
+
 
 |Software| requires a pedigree file (``-ped_file``) and one or more genomic data files to run the analysis.
 
@@ -158,8 +161,7 @@ Hybrid peeling arguments
     Single locus arguments:
       -seg_file SEG_FILE    A segregation probabilities file for hybrid peeling.
       -seg_map_file SEG_MAP_FILE
-                            A map file for loci in the segregation probabilities file.
-      -map_file MAP_FILE    A map file for all loci in hybrid peeling.
+                            A map file for loci in the segregation probabilities file in hybrid peeling.
 
 In order to run hybrid peeling the user needs to supply a ``-map_file`` which gives the genetic positions for the SNPs in the sequence allele read counts data supplied, a ``-seg_map_file`` which gives the genetic position for the SNPs in the segregation file, and a ``-seg_file`` which gives the segregation values generated via multi-locus iterative peeling. These arguments are not required for running in multi-locus mode.
 

example/runScript.sh

@@ -6,32 +6,33 @@ mkdir -p outputs
 AlphaPeel
 
 # Example 1: Performing multi-locus peeling with genotype data:
-AlphaPeel -genotypes data/genotypes.txt \
-         -pedigree data/pedigree.txt \
-         -out outputs/multilocus \
-         -nCycles 5 \
-         -runType multi \
-         -maxthreads 6
+AlphaPeel -geno_file data/genotypes.txt \
+         -ped_file data/pedigree.txt \
+         -out_file outputs/multilocus \
+         -n_cycle 5 \
+         -method multi \
+         -n_thread 6 \
+         -seg_prob
 
 # Example 1b: Performing multi-locus peeling with genotype data and calling the values with a threshold of 0.98
-AlphaPeel -genotypes data/genotypes.txt \
-         -pedigree data/pedigree.txt \
-         -out outputs/multilocus_with_phase \
-         -nCycles 5 \
-         -runType multi \
-         -maxthreads 6 \
+AlphaPeel -geno_file data/genotypes.txt \
+         -ped_file data/pedigree.txt \
+         -out_file outputs/multilocus_with_phase \
+         -n_cycle 5 \
+         -method multi \
+         -n_thread 6 \
          -geno_threshold 0.98 \
          -hap_threshold 0.98 \
          -geno \
          -hap
 
 # Example 2: Performing single-locus "hybrid" peeling with sequence data and pre-computed segregation estimates (generated from Example 1).
-AlphaPeel -seqfile data/sequence.txt \
-         -pedigree data/pedigree.txt \
-         -mapfile data/genotypes-map.txt\
-         -out outputs/hybrid \
-         -runType single \
-         -segmapfile data/segregation-map.txt \
-         -segfile outputs/multilocus.seg \
-         -nCycles 5 \
-         -maxthreads 6
+AlphaPeel -seq_file data/sequence.txt \
+         -ped_file data/pedigree.txt \
+         -map_file data/genotypes-map.txt\
+         -out_file outputs/hybrid \
+         -method single \
+         -seg_map_file data/segregation-map.txt \
+         -seg_file outputs/multilocus.seg_prob.txt \
+         -n_cycle 5 \
+         -n_thread 6

pyproject.toml

@@ -29,7 +29,9 @@ AlphaPeel = "tinypeel.tinypeel:main"
 
 [project.urls]
 "Homepage" = "https://github.com/AlphaGenes/AlphaPeel"
+"Documentation" = "https://alphapeel.readthedocs.io/en/latest/"
 "Bug Tracker" = "https://github.com/AlphaGenes/AlphaPeel/issues"
+"Changelog" = "https://github.com/AlphaGenes/AlphaPeel/blob/devel/Changelog.md"
 
 [tool.setuptools]
 package-dir = {"" = "src"}

src/tinypeel/Peeling/PeelingInfo.py

@@ -5,6 +5,7 @@
 
 from ..tinyhouse import ProbMath
 from ..tinyhouse import HaplotypeOperations
+from ..tinyhouse import InputOutput
 
 
 #####################################################################
@@ -36,9 +37,15 @@ def createPeelingInfo(pedigree, args, createSeg=True, phaseFounder=False):
     peelingInfo = jit_peelingInformation(
         nInd=pedigree.maxIdn, nFam=pedigree.maxFam, nLoci=nLoci, createSeg=createSeg
     )
+
     peelingInfo.isXChr = args.x_chr
     # Information about the peeling positions are handled elsewhere.
     peelingInfo.positions = None
+    if args.map_file:
+        peelingInfo.positions = np.array(
+            InputOutput.readMapFile(args.map_file, args.startsnp, args.stopsnp)[2],
+            dtype=np.int64,
+        )
 
     mutation_rate = args.mutation_rate
     # Generate the segregation tensors.
@@ -149,9 +156,9 @@ def setupTransmission(length, peelingInfo):
     if peelingInfo.positions is None:
         localMap = np.linspace(0, 1, num=peelingInfo.nLoci, dtype=np.float32)
     else:
-        localMap = (
-            peelingInfo.positions / peelingInfo.positions[-1]
-        )  # This should be sorted. Need to add in code to check.
+        localMap = (peelingInfo.positions - peelingInfo.positions[0]) / (
+            peelingInfo.positions[-1] - peelingInfo.positions[0]
+        )
     for i in range(peelingInfo.nLoci - 1):
         distance = localMap[i + 1] - localMap[i]
         distance = distance * length
@@ -300,7 +307,7 @@ def getHetMidpoint(geno):
 spec["transmissionRate"] = optional(float32[:])
 spec["maf"] = optional(float32[:])
 
-spec["positions"] = optional(float32[:])  # Not sure we use this.
+spec["positions"] = optional(int64[:])
 spec["iteration"] = int64
 
 

src/tinypeel/tinypeel.py

@@ -331,9 +331,7 @@ def generateSingleLocusSegregation(peelingInfo, pedigree, args):
     """
     if args.segfile is not None:
         # This just gets the locations in the map files.
-        snpMap = np.array(
-            InputOutput.readMapFile(args.map_file, args.startsnp, args.stopsnp)[2]
-        )
+        snpMap = peelingInfo.positions
         segMap = np.array(InputOutput.readMapFile(args.seg_map_file)[2])
 
         loci, distance = getLociAndDistance(snpMap, segMap)
@@ -614,6 +612,14 @@ def getArgs():
             "main_metafounder",
         ],
     )
+    input_parser.add_argument(
+        "-map_file",
+        default=None,
+        required=False,
+        type=str,
+        help="A map file for all loci.",
+    )
+
     # Output options
     output_parser = parser.add_argument_group("Output Options")
 
@@ -790,19 +796,12 @@ def getArgs():
     )
 
     singleLocus_parser = parser.add_argument_group("Hybrid peeling arguments")
-    singleLocus_parser.add_argument(
-        "-map_file",
-        default=None,
-        required=False,
-        type=str,
-        help="a map file for all loci in hybrid peeling.",
-    )
     singleLocus_parser.add_argument(
         "-seg_map_file",
         default=None,
         required=False,
         type=str,
-        help="a map file for loci in the segregation probabilities file.",
+        help="A map file for loci in the segregation probabilities file in hybrid peeling.",
     )
     singleLocus_parser.add_argument(
         "-seg_file",

tests/accuracy_tests/run_accu_test.py

@@ -192,7 +192,7 @@ def assess_peeling(sim_path, get_params, output_path, name, method, metafounder,
         else:
             print(f"File: {file}")
 
-        Marker_accu = [str(get_marker_accu(new_file[:, 1:], true_file[:, 1:]))]
+        Marker_accu = [str(get_marker_accu(new_file[:, :], true_file[:, :]))]
         for gen in range(nGen):
             Marker_accu.append(
                 str(
@@ -214,7 +214,7 @@ def assess_peeling(sim_path, get_params, output_path, name, method, metafounder,
         print("Marker_accuracies", " ".join(Marker_accu))
 
         Ind_accu = [
-            str(get_ind_accu(new_file[:, 1:], true_file[:, 1:], nIndPerGen, None))
+            str(get_ind_accu(new_file[:, :], true_file[:, :], nIndPerGen, None))
         ]
         for gen in range(nGen):
             Ind_accu.append(

tests/accuracy_tests/sim_for_alphapeel_accu_test/sim_for_accu.R

@@ -258,7 +258,7 @@ for (ind in (nIndPerGen + 1):nInd) {
 
 # ----- Realised recombination rate -----
 
-rec_count <- matrix(data = 0, nrow = nLociAll, ncol = 1)
+rec_count <- matrix(data = 0, nrow = nLociAll - 1, ncol = 1)
 for (ind in ((nIndPerGen + 1):nInd)) {
   for (chr in (1:nChr)) {
     indRecHist <- recHist[[ind]][[chr]]

tests/accuracy_tests/sim_for_alphapeel_accu_test/sim_for_sex_chr_accu.R

@@ -325,7 +325,7 @@ for (ind in (nIndPerGen + 1):nInd) {
 
 # ----- Realised recombination rate -----
 
-rec_count <- matrix(data = 0, nrow = nLociAll, ncol = 1)
+rec_count <- matrix(data = 0, nrow = nLociAll - 1, ncol = 1)
 for (ind in ((nIndPerGen + 1):nInd)) {
   for (chr in (1:nChr)) {
     indRecHist <- recHist[[ind]][[chr]]