Merge branch 'main' into help_and_docs

Author: prihoda

tracetrack/entities/record.py

@@ -106,12 +106,14 @@ class TraceSeqRecord(Record):
     """
     Class for storing information loaded from trace files (.ab1 files).
     """
-    def __init__(self, seq, quality, id=None, traces=None, base_locations=None, reference=None, reverse: bool = None):
+    def __init__(self, seq, quality, f, id=None, traces=None, base_locations=None, reference=None,
+                 reverse: bool = None):
         """
         This constructor is only called by read function below or from other functions starting with an instance, so
         the parameters need not be understood in too much details - reflects structure of a trace file.
         :param seq: Seq object with DNA sequence
         :param quality: list of int representing per base quality
+        :param f: mixed peaks detection threshold
         :param traces: dict containing a list of trace values for each base
         :param base_locations: list of positions (int) on the x axis, where each base is called
         :param reference: reference sequence in the DB, only added after instantiating
@@ -122,7 +124,8 @@ def __init__(self, seq, quality, id=None, traces=None, base_locations=None, refe
         self.traces = traces
         self.base_locations = base_locations
         self.quality = quality
-        self.mixed_peaks = self.find_mixed_peaks()
+        self.f = f
+        self.mixed_peaks = self.find_mixed_peaks(f)
         # should the original sequence be stored because of finding mixed peaks? Or just ignore al N's...
         self.reference = reference
         self.reverse = reverse
@@ -148,8 +151,8 @@ def read(cls, path):
         # get locations in trace array for all bases
         base_locations = list(record.annotations['abif_raw']["PLOC1"])
 
-        return cls(record.seq, record.letter_annotations['phred_quality'], id=name, traces=traces,
-                   base_locations=base_locations)
+        return cls(record.seq, record.letter_annotations['phred_quality'], f=0.15, id=name,
+                   traces=traces, base_locations=base_locations)
 
     def filter_sequence_by_quality(self, threshold, end_threshold):
         """
@@ -167,7 +170,8 @@ def filter_sequence_by_quality(self, threshold, end_threshold):
             traces=self.traces,
             base_locations=self.base_locations,
             reference=self.reference,
-            reverse=self.reverse
+            reverse=self.reverse,
+            f=self.f
         )
 
     def reverse_complement(self, **kwargs):
@@ -183,7 +187,8 @@ def reverse_complement(self, **kwargs):
             traces={str(Seq(base).reverse_complement()): values[::-1] for base, values in self.traces.items()},
             base_locations=[num_locations - i - 1 for i in self.base_locations[::-1]],
             reverse=False,
-            reference=self.reference
+            reference=self.reference,
+            f=self.f
         )
 
     def has_base_above_threshold(self):
@@ -242,7 +247,7 @@ def peak_borders(self, i: int):
             end = min(pos + width - 2, len(self.traces['A']))
         return start, end
 
-    def find_mixed_peaks(self, fraction: float = 0.15):
+    def find_mixed_peaks(self, fraction: float):
         """
         For each position of the sequence, determine if the peak in the chromatogram is "mixed".
         Disregard mixed signals in regions with low signal to noise ratio (generally bad quality region)
@@ -253,12 +258,12 @@ def find_mixed_peaks(self, fraction: float = 0.15):
         stn = [self.signal_to_noise(i) for i in range(len(self.base_locations))]
         avg_stn = sum(stn) / len(stn)
         mixed_peaks = []
+        threshold = max(25, avg_stn * 1.35)
 
         for i, pos in enumerate(self.base_locations):
             stn_local = stn[i-10:i] + stn[i+1:i+10]
             signal_to_noise = sum(stn_local) / 20
 
-            threshold = max(25, avg_stn * 1.35)
             if signal_to_noise < threshold:
                 continue
                 # bad StN ratio -> disregard potential mixed positions
@@ -269,13 +274,22 @@ def find_mixed_peaks(self, fraction: float = 0.15):
             peaks = {base: values[pos] for base, values in self.traces.items()}
             if base != "N":
                 main_peak = peaks[base.upper()]
-                for letter, area in areas.items():
-                    # check for both area and height of peak
-                    if base != letter and area > (areas[base.upper()] * fraction) and peaks[letter] > (main_peak * fraction) \
-                            and self.is_concave(pos, letter):
-                        mixed_peaks.append(i)
+            else:
+                # If the called base is "N", we use the highest peak present as the main peak
+                main_peak = max(peaks.values())
+                base = max(peaks, key=peaks.get)
+            for letter, area in areas.items():
+                # check for both area and height of peak
+                if base != letter and area > (areas[base.upper()] * fraction) and peaks[letter] > (main_peak * fraction) \
+                        and self.is_concave(pos, letter):
+                    mixed_peaks.append(i)
         return mixed_peaks
 
+    def re_find_mixed_peaks(self, f):
+        """Function to re-calculate mixed peaks after object has been constructed"""
+        self.mixed_peaks = self.find_mixed_peaks(f)
+        self.f = f
+
     def signal_to_noise(self, i: int):
         """
         Calculate the signal to noise ratio of a position as the ratio of the primary peak area to the sum of all other

tracetrack/flask_server.py

@@ -238,6 +238,16 @@ def alignment_settings_post(task_id, settings_task=None):
     threshold = int(request.form['threshold'])
     end_threshold = int(request.form['end_threshold'])
     separate = 'separate' in request.form
+    fraction = request.form['fraction']
+    messages = []
+    try:
+        fraction_f = float(fraction)
+        if fraction_f <= 0 or fraction_f > 1:
+            messages.append(f"The number provided for 'f' ({fraction}) is outside of the [0, 1] interval. Using default threshold 0.15.")
+            fraction_f = 0.15
+    except ValueError:
+        messages.append(f"The input for 'f' ({fraction}) could not be converted to a number. Using default threshold 0.15.")
+        fraction_f = 0.15
 
     result = scheduler.get_result(task_id)
     if isinstance(result, Exception):
@@ -256,8 +266,9 @@ def alignment_settings_post(task_id, settings_task=None):
             direction = str(request.form[f'dir_for_{seq_id}'])
             dir_flag = True if direction == "Rev" else False
             seq.flag_as_reverse(dir_flag)
+            seq.re_find_mixed_peaks(fraction_f)
 
-    new_id = schedule_tasks(sequences, population_names, db, separate, threshold, end_threshold)
+    new_id = schedule_tasks(sequences, population_names, db, separate, threshold, end_threshold, messages)
     return redirect("/results/" + task_id + "/" + new_id)
 
 

tracetrack/server_utils.py

@@ -558,13 +558,18 @@ def get_color_bin(value, thresholds):
 
 
 def trace_record_dict(trace: AlignedTrace) -> dict:
-    return {
+    d = {
         'sequence': str(trace.record.seq),
         'traces': [{'base': base, 'values': trace.record.traces[base]} for base in 'GATC'],
         'locations': trace.record.base_locations,
         'alignedPositions': trace.get_aligned_positions(),
         'id': trace.record.id
     }
+    last = max(0, d["locations"][d["alignedPositions"][0]] - 112)  # 112 points to the left of start of reference are still shown
+    for i in range(last):
+        for tr in d["traces"]:
+            tr["values"][i] = 0
+    return d
 
 
 def position_dict(alignment: Alignment, position: int) -> dict:
@@ -586,13 +591,14 @@ def hash_trace_name(pop_number, seq_number):
     return f"{pop_number}_{seq_number}"
 
 
-def schedule_tasks(sequences, population_names, db, separate, threshold, end_threshold):
+def schedule_tasks(sequences, population_names, db, separate, threshold, end_threshold, messages):
     settings = Settings(threshold, end_threshold, separate)
     inputs = []
     inputs_sorted = []
     for population_seqs, population_name in zip(sequences, population_names):
         new_seqlist = []
-        warnings = []
+        warnings = [] + messages
+        messages = []
         for record in population_seqs:
             new_record = record.filter_sequence_by_quality(threshold, end_threshold)
             if new_record.has_base_above_threshold():

tracetrack/templates/help.html

@@ -37,6 +37,13 @@ <h3>End Trimming Score Threshold</h3>
         As the quality towards the ends of trace files tends to decrease, both ends of sequences are discarded until three bases in a row pass the
         trimming threshold. This can also be set by the user and should be higher than the score threshold, otherwise it has no effect.
     </p>
+    <h3>Threshold for calling mixed peaks (<i>f</i>)</h3>
+    <p>
+        The application detects positions, where more than one trace peak is present, resulting in one primary and one or more secondary base calls.
+        This threshold specifies the minimal fraction of the area of the primary peak that another peak has to attain in order to be considered a secondary peak.
+    </p>
+
+
     <h3>Matching to Reference</h3>
     <p>
         Trace files are matched to reference sequences automatically. TraceTrack first checks if the reference ID is contained in the trace file name,

tracetrack/templates/index.html

@@ -50,6 +50,9 @@ <h4 class="badged-header"><span class="badge badge-primary">2</span>References</
                 <div class="example-used" style="display: none;">
                     <h4>Using example files</h4>
                 </div>
+                <div >
+
+                </div>
 
                 <br>
                 <div class="form-group">

tracetrack/templates/settings.html

@@ -38,6 +38,13 @@ <h4 class="badged-header"><span class="badge badge-primary">3</span>Alignment Se
                 </div>
             </div>
             <br>
+            <div class="form-row">
+                <div class="form-group">
+                    <input type="text" class="form-control" value="0.15" id="fraction" name="fraction" placeholder="0.15">
+                    <label class="formGroupExampleInput" for="fraction">Threshold for calling mixed peaks (<i>f</i>)</label>
+                </div>
+            </div>
+            <br>
             <div class="form-row">
                 <h4 class="badged-header"><span class="badge badge-primary">4</span>Reference Assignment</h4>
                 {% for trace_key, trace in ref_assignment.items() %}