187 incorrect numpy rng endpoint off by one low high

README.md

@@ -76,6 +76,12 @@ $ poetry install
 ```
 
 Notes: If any packages are struggling to resolve, check the channels and try to manually pip install the package to see if that helps (but note that NEAT is not tested on the pip versions.)
+If poetry hangs for you, try the following fix (from https://github.com/python-poetry/poetry/issues/8623):
+
+```
+export PYTHON_KEYRING_BACKEND=keyring.backends.null.Keyring
+```
+then re-run `poetry install`
 
 Test your install by running:
 ```

config_template/simple_template.yml

@@ -28,6 +28,6 @@ overwrite_output: .
 
 mode: .
 size: .
-jobs: .
+threads: .
 cleanup_splits: .
 reuse_splits: .

neat/models/mutation_model.py

@@ -49,14 +49,14 @@ class MutationModel(SnvModel, InsertionModel, DeletionModel):
     def __init__(self,
                  avg_mut_rate: float = default_avg_mut_rate,
                  homozygous_freq: float = default_homozygous_freq,
-                 variant_probs: dict[variants: float, ...] = default_variant_probs,
+                 variant_probs: dict[VariantTypes, float] = default_variant_probs,
                  transition_matrix: np.ndarray = default_mutation_sub_matrix,
                  is_cancer: bool = False,
                  # Any new parameters needed for new models should go below
                  trinuc_trans_matrices: np.ndarray = default_trinuc_trans_matrices,
                  trinuc_mut_bias: np.ndarray = default_trinuc_mut_bias,
-                 insert_len_model: dict[int: float] = default_insertion_len_model,
-                 deletion_len_model: dict[int: float] = default_deletion_len_model):
+                 insert_len_model: dict[int, float] = default_insertion_len_model,
+                 deletion_len_model: dict[int, float] = default_deletion_len_model):
 
         # Any new mutation types will need to be instantiated in the mutation model here
         SnvModel.__init__(self,
@@ -78,7 +78,7 @@ def __init__(self,
         self.all_dels = []
         self.all_ins = []
 
-    def get_mutation_type(self, rng: Generator) -> variants:
+    def get_mutation_type(self, rng: Generator) -> VariantTypes:
         """
         Picks one of the mutation types at random using a weighted list from the model.
         Note that the order of mutation types is Insertion, Deletion, SNV. To update the model selection if any
@@ -122,9 +122,9 @@ def generate_insertion(self, location: int, ref: Seq, rng: Generator) -> Inserti
         """
         This method generates an insertion object, based on the insertion model
 
-        :param location: The location of the variant, relative to the reference
-        :param ref: The reference for which to generate the variant
-        :param rng: The random number generator for the run
+        :param location: The location of the variant, relative to the reference.
+        :param ref: The reference for which to generate the variant.
+        :param rng: The random number generator for the run.
         :return:
         """
         # Note that insertion length model is based on the number of bases inserted. We add 1 to the length

neat/models/variant_models.py

@@ -46,7 +46,7 @@ class InsertionModel(VariantModel):
     _type = Insertion
     _description = "An insertion of N nucleotides into a chromosome."
 
-    def __init__(self, insert_len_model: dict[int: float, ...]):
+    def __init__(self, insert_len_model: dict[int, float]):
         # Creating probabilities from the weights
         tot = sum(insert_len_model.values())
         self.insertion_len_model = {key: val / tot for key, val in insert_len_model.items()}
@@ -78,12 +78,12 @@ class DeletionModel(VariantModel):
     _type = Deletion
     _description = "A deletion of a random number of bases"
 
-    def __init__(self, deletion_len_model: dict[int, float, ...]):
+    def __init__(self, deletion_len_model: dict[int, float]):
         # Creating probabilities from the weights
         tot = sum(deletion_len_model.values())
         self.deletion_len_model = {key: val/tot for key, val in deletion_len_model.items()}
 
-    def get_deletion_length(self, rng: Generator, size: int = None) -> int | list[int, ...]:
+    def get_deletion_length(self, rng: Generator, size: int = None) -> int | list[int]:
         """
         Get size number of inserts lengths. Size == 1 results in an int return, else a list of ints.
 

neat/read_simulator/runner.py

@@ -2,6 +2,7 @@
 Runner for generate_reads task
 """
 import logging
+import os
 import time
 import multiprocessing as mp
 
@@ -130,7 +131,8 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str):
     output_opts: list = []
     output_files: list = []
     # a dict with contig keys, then the thread index, and finally the applicable contig variants as the value
-    all_variants: dict[str, dict[int, ContigVariants]] = {chrom: {} for chrom in reference_index.keys()}
+    # TODO Remove if not needed
+    # all_variants: dict[str, dict[int, ContigVariants]] = {chrom: {} for chrom in reference_index.keys()}
     thread_idx = 1
     contig_list = list(reference_keys_with_lens.keys())
     contig_dict = {contig: contig_list.index(contig) for contig in reference_keys_with_lens.keys()}
@@ -173,7 +175,8 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str):
                     mutation_rate_dict[contig],
                 )
                 _LOG.info(f"Completed simulating contig {contig}.")
-                all_variants[contig][idx] = local_variants
+                # TODO Remove if not needed
+                # all_variants[contig][idx] = local_variants
                 output_files.append((thread_idx, files_written))
             else:
                 thread_input_variants = filter_thread_variants(input_variants_dict[contig], (start, start+length))
@@ -205,7 +208,8 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str):
 
         # Need to organize the results, as above
         for idx, contig, local_variants, files_written in results.get():
-            all_variants[contig][idx] = local_variants
+            # TODO Remove if not needed
+            # all_variants[contig][idx] = local_variants
             output_files.append((idx, files_written))
 
     _LOG.info("Processing complete, writing output")
@@ -229,7 +233,6 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str):
         elif "fastq" in file.name:
             continue
         else:
-            # Must be a vcf
             pysam.tabix_index(str(file), preset="vcf", force=force)
 
     _LOG.info(f"Read simulator complete in {time.time() - analysis_start} s")

neat/read_simulator/single_runner.py

@@ -80,14 +80,17 @@ def read_simulator_single(
     local_ref_name = list(local_ref_index.keys())[0]
     local_seq_record = local_ref_index[local_ref_name]
 
+    if len(local_seq_record) < local_options.read_len:
+        _LOG.debug("Record too small for processing")
+
+
     coords = (block_start, block_start+len(local_seq_record))
     mutation_rate_regions = recalibrate_mutation_regions(mutation_regions, coords, mut_model.avg_mut_rate)
 
     # Creates files and sets up objects for files that can be written to as needed.
     # Also creates headers for bam and vcf.
     # We'll also keep track here of what files we are producing.
     # We don't really need to write out the VCF. We should be able to store it in memory
-    local_options.produce_vcf = False
     local_output_file_writer = OutputFileWriter(options=local_options, header=bam_header)
     """
     Begin Analysis
@@ -98,7 +101,7 @@ def read_simulator_single(
         reference=local_seq_record,
         ref_start=block_start,
         mutation_rate_regions=mutation_rate_regions,
-        existing_variants=input_variants_local,
+        input_variants=input_variants_local,
         mutation_model=mut_model,
         max_qual_score=max_qual_score,
         options=local_options,
@@ -146,7 +149,8 @@ def read_simulator_single(
             os.rename(str(sorted_bam), str(local_output_file_writer.bam))
             _LOG.info(f"bam for thread {thread_idx} written")
 
-    write_block_vcf(local_variants, contig_name, local_ref_index, local_output_file_writer)
+    if local_options.produce_vcf:
+        write_block_vcf(local_variants, contig_name, block_start, local_ref_index, local_output_file_writer)
 
     local_output_file_writer.flush_and_close_files(False)
     file_dict = {
@@ -165,14 +169,22 @@ def read_simulator_single(
 def write_block_vcf(
         local_variants: ContigVariants,
         contig: str,
+        ref_start: int,
         ref_index: dict,
         ofw: OutputFileWriter,
 ):
+    """
+    :param local_variants: The ContigVariants object for this block
+    :param contig: The name of the contig this block comes from
+    :param ref_start: The reference start position, relative to the contig
+    :param ref_index: The local index for this reference block
+    :param ofw: The OutputFileWriter object for this block
+    """
     # take contig name from ref index because we know it is in the proper order
     locations = sorted(local_variants.variant_locations)
     for location in locations:
         for variant in local_variants[location]:
-            ref, alt = local_variants.get_ref_alt(variant, ref_index[contig])
+            ref, alt = local_variants.get_ref_alt(variant, ref_index[contig], ref_start)
             sample = local_variants.get_sample_info(variant)
             # +1 to position because the VCF uses 1-based coordinates
             #          .id should give the more complete name

neat/read_simulator/utils/generate_variants.py

@@ -52,12 +52,11 @@ def map_non_n_regions(sequence) -> np.ndarray:
 
     return base_check
 
-
 def generate_variants(
         reference: SeqRecord,
         ref_start: int,
         mutation_rate_regions: list[tuple[int, int, float]],
-        existing_variants: ContigVariants,
+        input_variants: ContigVariants,
         mutation_model: MutationModel,
         options: Options,
         max_qual_score: int,
@@ -67,13 +66,16 @@ def generate_variants(
     :param reference: The reference to generate variants off of. This should be a SeqRecord object
     :param ref_start: where on the reference this SeqRecord above starts.
     :param mutation_rate_regions: Genome segments with associated mutation rates
-    :param existing_variants: Any input variants or overlaps to pick up from the previous contig
+    :param input_variants: Any input variants or overlaps to pick up from the previous contig
     :param mutation_model: The mutation model for the dataset
     :param options: The options for the run
     :param max_qual_score: the maximum quality score for the run.
     """
-    return_variants = existing_variants
-    # Step 1: Create a VCF of mutations
+    return_variants = ContigVariants()
+    for variant_location in input_variants.variant_locations:
+        if ref_start <= variant_location < ref_start + len(reference):
+            for variant in input_variants.contig_variants[variant_location]:
+                return_variants.add_variant(variant)
 
     # pase out the mutation rates
     mutation_rates = np.array([x[2] for x in mutation_rate_regions])
@@ -121,50 +123,67 @@ def generate_variants(
         # For no input mutation regions bed, this will return the entire sequence.
         mut_region = options.rng.choice(a=local_mut_regions, p=probability_rates)
         mut_region_offset = (int(mut_region[0]-ref_start), int(mut_region[1]-ref_start), mut_region[2])
+
         # Pick a random starting place. Randint is inclusive of endpoints, so we subtract 1
-        window_start = options.rng.integers(mut_region_offset[0], mut_region_offset[1] - 1, dtype=int)
+        if mut_region_offset[1] <= mut_region_offset[0]:
+            _LOG.error("Invalid Mutation Region")
+            raise ValueError
 
+        window_start = options.rng.integers(mut_region_offset[0], mut_region_offset[1] - 1, dtype=int)
         found = False
+        if window_start > len(reference):
+            _LOG.error("Invalid Mutation Region")
+            raise ValueError
         if reference[window_start] not in ALLOWED_NUCL:
-            # pick a random location to the right
-            plus = options.rng.integers(window_start + 1, mut_region_offset[1] - 1, dtype=int)
-            if reference[plus] in ALLOWED_NUCL:
-                found = True
-                window_start = plus
-            else:
-                # If that didn't work pick a random location to the left
-                if window_start - 1 > mut_region_offset[0]:
-                    minus = options.rng.integers(mut_region_offset[0], window_start - 1, dtype=int)
-                    if reference[minus] in ALLOWED_NUCL:
-                        found = True
-                        window_start = minus
+            # Ensure we are not too close to the edge
+            if window_start + 1 < mut_region_offset[1] - 1 and mut_region_offset[1] - 1 > 0:
+                # pick a random location to the right
+                plus = options.rng.integers(window_start + 1, mut_region_offset[1] - 1, dtype=int)
+                if reference[plus] in ALLOWED_NUCL:
+                    found = True
+                    window_start = plus
+                else:
+                    # If that didn't work pick a random location to the left
+                    if mut_region_offset[0] < window_start - 1 and window_start - 1 > 0:
+                        minus = options.rng.integers(mut_region_offset[0], window_start - 1, dtype=int)
+                        if reference[minus] in ALLOWED_NUCL:
+                            found = True
+                            window_start = minus
         else:
             found = True
 
         # If we couldn't find a spot, try again
         if not found:
             continue
 
-        end_point = min(
-            options.rng.integers(
-                window_start,
-                min(mut_region_offset[1], window_start+max_window_size)-1,
-                dtype=int
-            ),
-            # Don't go past the barrier
-            len(reference)-1
-        )
-        if reference[end_point] not in ALLOWED_NUCL:
-            # Didn't find it to the right, look to the left
-            end_point = options.rng.integers(
-                max(window_start - max_window_size, mut_region_offset[0]),
-                window_start,
-                dtype=int
+        if window_start < min(mut_region_offset[1], window_start+max_window_size)-1:
+            end_point = min(
+                options.rng.integers(
+                    window_start,
+                    min(mut_region_offset[1], window_start+max_window_size)-1,
+                    dtype=int
+                ),
+                # Don't go past the barrier
+                len(reference)-1
             )
             if reference[end_point] not in ALLOWED_NUCL:
-                # No suitable end_point, so we try again
-                continue
+                # Didn't find it to the right, look to the left
+                if max(window_start - max_window_size, mut_region_offset[0]) > window_start:
+                    # ensure we are still in bounds
+                    continue
+                end_point = options.rng.integers(
+                    max(window_start - max_window_size, mut_region_offset[0]),
+                    window_start,
+                    dtype=int
+                )
+                if reference[end_point] not in ALLOWED_NUCL:
+                    # No suitable end_point, so we try again
+                    continue
+        else:
+            end_point = None
 
+        if not window_start or not end_point:
+            continue
         # Sorting assures that wherever we found the end point, the coordinates will be in the correct order for slicing
         mutation_slice = sorted([window_start, end_point])
         slice_distance = mutation_slice[1] - mutation_slice[0]
@@ -203,21 +222,24 @@ def generate_variants(
                 # Because the slice is a substring, we add the start point onto the relative location to get the
                 # location relative to the reference.
                 position = find_random_non_n(options.rng, n_gaps)  # position in slice
-                location = position + mutation_slice[0]  # location relative to reference
+                location = position + mutation_slice[0] + ref_start  # location relative to overall contig
                 if variant_type == Insertion:
+                    # Want to insert the variant in the overall context
                     temp_variant = mutation_model.generate_insertion(location, subsequence[position], options.rng)
                 else:
+                    # Want to insert the variant in the overall context
                     temp_variant = mutation_model.generate_deletion(location, options.rng)
 
             # Case 2: SNV
             elif variant_type == SingleNucleotideVariant:
                 # We'll sample for the location within this slice
                 # It's a relative location, so we add the start point of the subsequence to that.
                 position = mutation_model.sample_trinucs(options.rng)  # position in slice
-                location = position + mutation_slice[0]  # location relative to reference
-                if location == 0:
+                local_location = position + mutation_slice[0]  # location relative to slice
+                location = local_location + ref_start  # relative to overall contig
+                if local_location == 0:
                     continue
-                trinuc = reference[location: location+3].seq.upper()
+                trinuc = reference[local_location: local_location+3].seq.upper()
                 disallowed_chars = False
                 for letter in trinuc:
                     if letter not in ALLOWED_NUCL:
@@ -284,7 +306,7 @@ def generate_variants(
             # The count will tell us how many we actually added v how many we were trying to add
             n_added += 1
 
-    # _LOG.debug(f'Finished generating chunk random mutations in {(time.time() - start_time)/60:.2f} minutes')
+    _LOG.debug(f'Finished generating chunk random mutations in {(time.time() - start_time)/60:.2f} minutes')
 
     return return_variants
 

neat/read_simulator/utils/output_file_writer.py

@@ -14,7 +14,6 @@
 import logging
 from typing import Any
 
-import pysam
 from Bio import bgzf
 from pathlib import Path