Merge pull request #124 from ncsa/develop

Develop

environment.yml

@@ -17,4 +17,4 @@ dependencies:
   - pip:
     - pysam
     - frozendict
-    - poetry>=1.3.*
+    - poetry==1.3.*

neat/cli/commands/model_sequencing_error.py

@@ -48,12 +48,9 @@ def add_arguments(self, parser: argparse.ArgumentParser):
                             dest="quality_scores",
                             nargs='+',
                             required=False,
-                            default=[2, 11, 25, 37],
-                            help="Quality score bins. Enter as a space separeted list for binned scores "
-                                 "(-Q 2 11 24 37), or enter a single maximum score for a full range (-Q 42 gives all"
-                                 "scores from 1-42). The default is binned quality scores: [2, 11, 24, 37]. Note that"
-                                 "using quality score bins on an unbinned fastq will result in a binned model, at the"
-                                 "cost of some inaccuracy.")
+                            default=42,
+                            help="Quality score max. The default 42. The lowest possible score is 1. To used binned"
+                                 "scoring, enter a space separated list of scores, e.g., -Q 2 15 23 37")
 
         parser.add_argument('-m',
                             type=int,

neat/gen_mut_model/runner.py

@@ -85,18 +85,23 @@ def runner(reference_index,
 
     # Pre-parsing to find all the matching chromosomes between ref and vcf
     _LOG.info('Processing VCF file...')
-    matching_variants, matching_chromosomes = read_and_filter_variants(vcf_to_process,
-                                                                       reference_index,
-                                                                       ignore)
+    matching_variants, matching_chromosomes = read_and_filter_variants(
+        vcf_to_process,
+        reference_index,
+        ignore
+    )
 
     if not matching_variants or not matching_chromosomes:
         _LOG.error("No valid variants detected. Check names in vcf versus reference and/or bed.")
         sys.exit(1)
 
-    trinuc_ref_count, bed_track_length = count_trinucleotides(reference_index,
-                                                              bed,
-                                                              outcounts_file,
-                                                              matching_chromosomes)
+    trinuc_ref_count, bed_track_length = count_trinucleotides(
+        reference_index,
+        bed,
+        outcounts_file,
+        matching_chromosomes,
+        save_trinuc
+    )
 
     if not trinuc_ref_count:
         _LOG.error("No valid trinucleotides detected in reference.")
@@ -296,47 +301,6 @@ def runner(reference_index,
         for k in sorted(snp_trans_freq):
             _LOG.info(f'p({k[0]} --> {k[1]} | SNP occurs) = {snp_trans_freq[k]}')
 
-    # Save counts, if requested
-    if save_trinuc:
-        trinuc_output_data = {
-            'p(snp)': average_snp_freq,
-            'p(insertion)': average_insertion_frequency,
-            'p(deletion)': average_deletion_frequency,
-            'overall average mut rate': avg_mut_rate,
-            'total variants processed': total_var,
-            'homozygous probability': homozygous_frequency
-        }
-
-        for k in sorted(trinuc_mut_prob):
-            trinuc_output_data[f'p({k} mutates)'] = trinuc_mut_prob[k]
-
-        for k in sorted(trinuc_trans_probs):
-            trinuc_output_data[f'p({k[0]} --> {k[1]} | {k[0]} mutates)'] = trinuc_trans_probs[k]
-
-        for k in sorted(deletion_counts):
-            trinuc_output_data[f'p(del length = {abs(k)} | deletion occurs)'] = deletion_frequency[k]
-
-        for k in sorted(insertion_counts):
-            trinuc_output_data[f'p(insert length = {abs(k)} | insertion occurs)'] = insertion_freqency[k]
-
-        for k in sorted(snp_trans_freq):
-            trinuc_output_data[f'p({k[0]} --> {k[1]} | SNP occurs)'] = snp_trans_freq[k]
-
-        trinuc_output_path = Path(output)
-        trinuc_output_path = trinuc_output_path.with_suffix('')
-        trinuc_output_path = trinuc_output_path.with_suffix('.trinuc.pickle.gz')
-        _LOG.info(f'Saving trinucleotide counts to: {trinuc_output_path}')
-
-        with open_output(trinuc_output_path, 'w+') as trinuc_outfile:
-
-            pickle.dump(trinuc_output_data, trinuc_outfile)
-
-            # Human-readable content of file below
-            trinuc_outfile.write('\n')
-
-            for key, value in trinuc_output_data.items():
-                trinuc_outfile.write(f'{key}: {value}\n')
-
     _LOG.info(f'p(snp)   = {average_snp_freq}')
     _LOG.info(f'p(insertion) = {average_insertion_frequency}')
     _LOG.info(f'p(deletion) = {average_deletion_frequency}')
@@ -360,8 +324,8 @@ def runner(reference_index,
         transition_matrix=snp_transition_bias,
         trinuc_trans_matrices=trinuc_transition_bias,
         trinuc_mut_bias=trinuc_mutation_probs,
-        insert_len_model=insertion_counts,
-        deletion_len_model=deletion_counts
+        insert_len_model=insertion_freqency,
+        deletion_len_model=deletion_frequency
     )
 
     print('\nSaving model...')
@@ -405,6 +369,8 @@ def compute_mut_runner(reference,
     if outcounts:
         validate_input_path(outcounts)
         outcounts = Path(outcounts)
+    elif save_trinuc:
+        outcounts = Path(output + '.trinuc.pickle.gz')
 
     print('Processing reference...')
     reference_index = SeqIO.index(reference, 'fasta')
@@ -477,8 +443,18 @@ def compute_mut_runner(reference,
     output = Path(output + '.pickle.gz')
     validate_output_path(output, overwrite=overwrite_output)
 
-    runner(reference_index, vcf_to_process, vcf_columns, outcounts, show_trinuc, save_trinuc,
-           output, bed, human_sample, skip_common)
+    runner(
+        reference_index,
+        vcf_to_process,
+        vcf_columns,
+        outcounts,
+        show_trinuc,
+        save_trinuc,
+        output,
+        bed,
+        human_sample,
+        skip_common
+    )
 
     if os.path.exists('temp.vcf'):
         os.remove('temp.vcf')

neat/gen_mut_model/utils.py

@@ -4,6 +4,8 @@
 
 import json
 import sys
+import pickle
+import gzip
 
 import numpy as np
 
@@ -165,14 +167,16 @@ def cluster_list(list_to_cluster, delta):
 def count_trinucleotides(reference_index,
                          bed,
                          trinuc_counts,
-                         matching_chroms):
+                         matching_chroms,
+                         save_trinuc_file):
     """
     Counts the frequency of the various trinucleotide combinations in the dataset
 
     :param dict reference_index: The index of the fasta, from SeqIO
     :param Path bed: Full path to bed file, if using
     :param Path trinuc_counts: Full path to existing trinculeotide counts file, if input, or the output path
     :param list matching_chroms: List of matching chromosomes
+    :param save_trinuc_file: Boolean determines whether to save trinucleotide counts file.
     :return dict, int: A dictionary of trinculeotides and counts, the number of bases spanned
     """
     # Data format: TRINUC: COUNT (e.g., "AAA": 12), where COUNT is the number of times trinuc was observed
@@ -183,11 +187,11 @@ def count_trinucleotides(reference_index,
     track_len = 0
 
     trinuc_counts_exists = False
-    save_trinuc_file = False
     if trinuc_counts:
         if trinuc_counts.is_file():
+            _LOG.info("Trinucleotide counts file exists, skipping save to avoid overwriting data.")
             trinuc_counts_exists = True
-        save_trinuc_file = True
+            save_trinuc_file = False
 
     if bed:
         _LOG.info("Counting trinucleotide combinations in bed regions")
@@ -207,7 +211,6 @@ def count_trinucleotides(reference_index,
                         trinuc_ref_count = count_trinuc(trinuc, trinuc_ref_count)
 
     # Solution to attribute error (needs to be checked)
-    # TODO remove ref_name from this dict
     elif not trinuc_counts_exists:
         _LOG.info("Counting trinucleotide combinations in reference.")
         for ref_name in matching_chroms:
@@ -220,13 +223,12 @@ def count_trinucleotides(reference_index,
             with open_output(trinuc_counts, 'w') as countfile:
                 _LOG.info('Saving trinuc counts to file...')
                 # Convert all values to writable formats
-                trinuc_ref_count = {str(x): int(y) for x, y in trinuc_ref_count.items()}
-                countfile.write(json.dumps(trinuc_ref_count))
+                pickle.dump(trinuc_ref_count, countfile)
 
     else:
         _LOG.info(f'Loading file: {trinuc_counts}.')
-        with open_input(trinuc_counts) as counts:
-            trinuc_ref_count = json.load(counts)
+        with gzip.open(trinuc_counts, 'rb') as counts:
+            trinuc_ref_count = pickle.load(counts)
         if save_trinuc_file:
             _LOG.warning('Existing trinucelotide file will not be changed.')
 

neat/model_sequencing_error/runner.py

@@ -62,9 +62,11 @@ def model_seq_err_runner(
     _LOG.debug(f"Quality offset: {offset}")
 
     final_quality_scores: list
+    binned_scores = False
     if len(qual_scores) == 1:
         final_quality_scores = list(range(1, qual_scores[0] + 1))
     else:
+        binned_scores = True
         final_quality_scores = sorted(qual_scores)
 
     _LOG.debug(f'Quality scores: {final_quality_scores}')
@@ -109,11 +111,14 @@ def model_seq_err_runner(
     for file in files:
         file_num += 1
         _LOG.info(f'Reading file {file_num} of {len(files)}')
-        parameters_by_position, file_avg_error, read_length = parse_file(file,
-                                                                         final_quality_scores,
-                                                                         num_records_to_process,
-                                                                         offset,
-                                                                         read_length)
+        parameters_by_position, file_avg_error, read_length = parse_file(
+            file,
+            final_quality_scores,
+            num_records_to_process,
+            offset,
+            read_length,
+            binned_scores
+        )
 
         read_parameters.append(parameters_by_position)
         average_errors.append(file_avg_error)

neat/model_sequencing_error/utils.py

@@ -10,7 +10,6 @@
 
 from scipy.stats import mode
 from ..common import open_input
-from ..models import take_closest
 
 __all__ = [
     "parse_file"
@@ -139,10 +138,7 @@ def parse_file(input_file: str, quality_scores: list, max_reads: int, qual_offse
 
             for j in range(read_length):
                 # The qualities of each read_position_scores
-                quality_bin = take_closest(quality_scores, qualities_to_check[j])
-                bin_index = quality_scores.index(quality_bin)
-                temp_q_count[j][bin_index] += 1
-                qual_score_counter[quality_bin] += 1
+                qual_score_counter[qualities_to_check[j]] += 1
 
             if records_read % quarters == 0:
                 _LOG.info(f'reading data: {(records_read / max_reads) * 100:.0f}%')

neat/models/__init__.py

@@ -1,2 +1 @@
 from .models import *
-from .utils import *

neat/models/default_sequencing_error_model.py

@@ -17,8 +17,8 @@
      [0.2505, 0.2552, 0.4943, 0.0]]
 )
 
-# This list may not be the final list
-default_quality_scores = np.array([2, 11, 25, 37])
+# All numbers from 1-42
+default_quality_scores = np.arange(1, 43)
 
 # This puts a high probability toward getting a maximum quality score. The current values
 # should be considered temporary. We're working on final values.

neat/models/models.py

@@ -18,7 +18,6 @@
 from ..common import TRINUC_IND, ALLOWED_NUCL, NUC_IND, DINUC_IND
 from .default_mutation_model import *
 from .default_sequencing_error_model import *
-from .utils import bin_scores, take_closest
 
 __all__ = [
     "MutationModel",
@@ -60,7 +59,9 @@ class InsertionModel(VariantModel):
     def __init__(self,
                  insert_len_model: dict[int: float, ...],
                  rng: Generator = None):
-        self.insert_len_model = insert_len_model
+        # 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()}
         self.rng = rng
 
     def get_insertion_length(self, size: int = None) -> int | list[int, ...]:
@@ -71,8 +72,8 @@ def get_insertion_length(self, size: int = None) -> int | list[int, ...]:
                      Greater than 1 returns a list of ints.
         :return: int or list of ints.
         """
-        return self.rng.choice(a=list(self.insert_len_model),
-                               p=[*self.insert_len_model.values()],
+        return self.rng.choice(a=list(self.insertion_len_model),
+                               p=[*self.insertion_len_model.values()],
                                size=size, shuffle=False)
 
 
@@ -91,7 +92,9 @@ class DeletionModel(VariantModel):
     def __init__(self,
                  deletion_len_model: dict[int: float, ...],
                  rng: Generator = None):
-        self.deletion_len_model = deletion_len_model
+        # 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()}
         self.rng = rng
 
     def get_deletion_length(self, size: int = None) -> int | list[int, ...]:
@@ -281,6 +284,9 @@ def generate_snv(self, trinucleotide: Seq, reference_location: int) -> SingleNuc
         transition_matrix = self.trinuc_trans_matrices[DINUC_IND[trinucleotide[0] + "_" + trinucleotide[2]]]
         # then determine the trans probs based on the middle nucleotide
         transition_probs = transition_matrix[NUC_IND[trinucleotide[1]]]
+        # Creating probabilities from the weights
+        transition_sum = sum(transition_probs)
+        transition_probs = [x/transition_sum for x in transition_probs]
         # Now pick a random alternate, weighted by the probabilities
         alt = self.rng.choice(ALLOWED_NUCL, p=transition_probs)
         temp_snv = SingleNucleotideVariant(reference_location, alt=alt)
@@ -376,11 +382,8 @@ def __init__(self,
         self.insertion_model = insertion_model
         self.uniform_quality_score = None
         if self.is_uniform:
-            # bin scores returns a list, so we need the first (only) element of the list
-            converted_avg_err = bin_scores(self.quality_scores,
-                                           [int(-10. * np.log10(self.average_error))])[0]
-            # Set score to the lowest of the max of the quality scores and the bin closest to the input avg error.
-            self.uniform_quality_score = min([max(self.quality_scores), converted_avg_err])
+            # Set score to the lowest of the max of the quality scores and the input avg error.
+            self.uniform_quality_score = min([max(self.quality_scores), int(-10. * np.log10(self.average_error) + 0.5)])
         self.rng = rng
 
     def get_sequencing_errors(self,
@@ -498,7 +501,13 @@ def get_quality_scores(self,
             for i in quality_index_map:
                 score = self.rng.normal(self.quality_score_probabilities[i][0],
                                         scale=self.quality_score_probabilities[i][1])
-                score = take_closest(self.quality_scores, score)
+                # make sure score is in range and an int
+                score = round(score)
+                if score > 42:
+                    score = 42
+                if score < 1:
+                    score = 1
+
                 temp_qual_array.append(score)
 
         if self.rescale_qualities:
@@ -509,9 +518,9 @@ def get_quality_scores(self,
                                                           self.quality_score_error_rate[n]) + 0.5)])
                               for n in temp_qual_array]
             # Now rebin the quality scores.
-            temp_qual_array = np.array(bin_scores(self.quality_scores, rescaled_quals))
+            temp_qual_array = np.array(rescaled_quals)
         else:
-            temp_qual_array = np.array(bin_scores(self.quality_scores, temp_qual_array))
+            temp_qual_array = np.array(temp_qual_array)
 
         return temp_qual_array[:input_read_length]