Bug fixes in seqpos of dna timeseries

biobb_dna/dna/dna_timeseries.py

@@ -3,6 +3,7 @@
 """Module containing the HelParTimeSeries class and the command line interface."""
 import argparse
 import zipfile
+import re
 from pathlib import Path
 
 import pandas as pd
@@ -27,9 +28,9 @@ class HelParTimeSeries(BiobbObject):
         properties (dict):
             * **sequence** (*str*) - (None) Nucleic acid sequence corresponding to the input .ser file. Length of sequence is expected to be the same as the total number of columns in the .ser file, minus the index column (even if later on a subset of columns is selected with the *usecols* option).
             * **bins** (*int*) - (None) Bins for histogram. Parameter has same options as matplotlib.pyplot.hist.
-            * **helpar_name** (*str*) - (Optional) helical parameter name.
+            * **helpar_name** (*str*) - (None) Helical parameter name. It must match the name of the helical parameter in the .ser input file. Values: majd, majw, mind, minw, inclin, tip, xdisp, ydisp, shear, stretch, stagger, buckle, propel, opening, rise, roll, twist, shift, slide, tilt, alphaC, alphaW, betaC, betaW, gammaC, gammaW, deltaC, deltaW, epsilC, epsilW, zetaC, zetaW, chiC, chiW, phaseC, phaseW.
             * **stride** (*int*) - (1000) granularity of the number of snapshots for plotting time series.
-            * **seqpos** (*list*) - (None) list of sequence positions (columns indices starting by 0) to analyze.  If not specified it will analyse the complete sequence.
+            * **seqpos** (*list*) - (None) list of sequence positions (columns indices starting by 1) to analyze.  If not specified it will analyse the complete sequence.
             * **remove_tmp** (*bool*) - (True) [WF property] Remove temporal files.
             * **restart** (*bool*) - (False) [WF property] Do not execute if output files exist.
             * **sandbox_path** (*str*) - ("./") [WF property] Parent path to the sandbox directory.
@@ -101,11 +102,12 @@ def __init__(self, input_ser_path, output_zip_path,
                     "Helical parameter name is invalid! "
                     f"Options: {constants.helical_parameters}")
 
-        # get base length and unit from helical parameter name
+        # get base length from helical parameter name
         if self.helpar_name.lower() in constants.hp_singlebases:
             self.baselen = 0
         else:
             self.baselen = 1
+        # get unit from helical parameter name
         if self.helpar_name in constants.hp_angular:
             self.hp_unit = "Degrees"
         else:
@@ -128,31 +130,71 @@ def launch(self) -> int:
         if self.sequence is None or len(self.sequence) < 2:
             raise ValueError("sequence is null or too short!")
 
-        # check seqpos
+        # calculate cols with 0 index
         if self.seqpos is not None:
-            if (max(self.seqpos) > len(self.sequence) - 2) or (min(self.seqpos) < 1):
+            cols = [i-1 for i in self.seqpos]
+        else:
+            cols = list(range(len(self.sequence)))
+
+        # sort cols in ascending order
+        cols.sort()
+
+        # check seqpos for base pairs
+        if self.seqpos is not None and self.helpar_name in constants.hp_basepairs:
+            if (max(cols) > len(self.sequence) - 2) or (min(cols) < 0):
+                raise ValueError(
+                    f"seqpos values must be between 1 and {len(self.sequence) - 1}")
+            if not (isinstance(self.seqpos, list) and len(self.seqpos) > 1):
+                raise ValueError(
+                    "seqpos must be a list of at least two integers")
+        # check seqpos for non base pairs
+        elif self.seqpos is not None and self.helpar_name not in constants.hp_basepairs:
+            if (max(cols) > len(self.sequence) - 1) or (min(cols) < 0):
                 raise ValueError(
-                    f"seqpos values must be between 1 and {len(self.sequence) - 2}")
+                    f"seqpos values must be between 1 and {len(self.sequence)}")
             if not (isinstance(self.seqpos, list) and len(self.seqpos) > 1):
                 raise ValueError(
                     "seqpos must be a list of at least two integers")
 
-        # read input .ser file
-        ser_data = read_series(
-            self.stage_io_dict['in']['input_ser_path'],
-            usecols=self.seqpos)
-        if self.seqpos is None:
-            ser_data = ser_data[ser_data.columns[1:-1]]
+        if self.helpar_name in constants.hp_basepairs:
+            # remove first and last base pairs from cols if they match 0 and len(sequence)
+            if min(cols) == 0:
+                cols.pop(0)
+            if max(cols) == len(self.sequence) - 1:
+                cols.pop(-1)
+
             # discard first and last base(pairs) from sequence
-            sequence = self.sequence[1:]
-            subunits = [
-                f"{i+1}_{sequence[i:i+1+self.baselen]}"
-                for i in range(len(ser_data.columns))]
+            sequence = self.sequence[1:-1]
+            # create indices list
+            indices = cols.copy()
+            # create subunits list from cols
+            subunits = [f"{i+1}_{sequence[i-1:i+self.baselen]}" for i in cols]
+            # clean subunits (leave only basepairs)
+            pattern = re.compile(r'\d+_[A-Za-z]{2}')
+            # get removed items
+            removed_items = [s for s in subunits if not pattern.fullmatch(s)]
+            # get indices of removed items (in integer format and starting from 0)
+            removed_numbers = [re.match(r'\d+', item).group() for item in removed_items]
+            removed_numbers = list(map(int, removed_numbers))
+            removed_numbers = [i-1 for i in removed_numbers]
+            # remove non basepairs from subunits and indices
+            subunits = [s for s in subunits if pattern.fullmatch(s)]
+            indices = [i for i in indices if i not in removed_numbers]
         else:
             sequence = self.sequence
-            subunits = [
-                f"{i+1}_{sequence[i:i+1+self.baselen]}"
-                for i in self.seqpos]
+            # create indices list
+            indices = cols.copy()
+            # trick for getting the index column from the .ser file
+            indices.insert(0, 0)
+            # create subunits list from cols
+            subunits = [f"{i+1}_{sequence[i:i+1+self.baselen]}" for i in cols]
+
+        # read input .ser file
+        ser_data = read_series(
+            self.stage_io_dict['in']['input_ser_path'],
+            usecols=indices)
+
+        # get columns for selected bases
         ser_data.columns = subunits
 
         # write output files for all selected bases (one per column)

biobb_dna/dna/dna_timeseries_unzip.py

@@ -6,6 +6,7 @@
 import shutil
 import argparse
 
+from biobb_dna.utils import constants
 from biobb_common.generic.biobb_object import BiobbObject
 from biobb_common.configuration import settings
 from biobb_common.tools import file_utils as fu
@@ -27,7 +28,7 @@ class DnaTimeseriesUnzip(BiobbObject):
             * **type** (*str*) - (None) Type of analysis, series or histogram. Values: series, hist.
             * **parameter** (*str*) - (None) Type of parameter. Values: majd, majw, mind, minw, inclin, tip, xdisp, ydisp, shear, stretch, stagger, buckle, propel, opening, rise, roll, twist, shift, slide, tilt, alphaC, alphaW, betaC, betaW, gammaC, gammaW, deltaC, deltaW, epsilC, epsilW, zetaC, zetaW, chiC, chiW, phaseC, phaseW.
             * **sequence** (*str*) - (None) Nucleic acid sequence used for generating dna_timeseries output file.
-            * **index** (*int*) - (0) Base pair index in the parameter 'sequence', starting from 0.
+            * **index** (*int*) - (1) Base pair index in the parameter 'sequence', starting from 1.
             * **remove_tmp** (*bool*) - (True) [WF property] Remove temporal files.
             * **restart** (*bool*) - (False) [WF property] Do not execute if output files exist.
             * **sandbox_path** (*str*) - ("./") [WF property] Parent path to the sandbox directory.
@@ -79,7 +80,7 @@ def __init__(self, input_zip_file,
         self.type = properties.get('type', None)
         self.parameter = properties.get('parameter', None)
         self.sequence = properties.get('sequence', None)
-        self.index = properties.get('index', 0)
+        self.index = properties.get('index', 1)
         self.properties = properties
 
         # Check the properties
@@ -108,11 +109,8 @@ def launch(self) -> int:
             exit(1)
 
         # Check that the parameter is valid
-        if self.parameter not in ["majd", "majw", "mind", "minw", "inclin", "tip", "xdisp", "ydisp", "shear", "stretch",
-                                  "stagger", "buckle", "propel", "opening", "rise", "roll", "twist", "shift", "slide",
-                                  "tilt", "alphaC", "alphaW", "betaC", "betaW", "gammaC", "gammaW", "deltaC", "deltaW",
-                                  "epsilC", "epsilW", "zetaC", "zetaW", "chiC", "chiW", "phaseC", "phaseW"]:
-            fu.log(f"Parameter {self.parameter} not valid. Valid parameters are: majd, majw, mind, minw, inclin, tip, xdisp, ydisp, shear, stretch, stagger, buckle, propel, opening, rise, roll, twist, shift, slide, tilt, alphaC, alphaW, betaC, betaW, gammaC, gammaW, deltaC, deltaW, epsilC, epsilW, zetaC, zetaW, chiC, chiW, phaseC, phaseW.",
+        if self.parameter not in constants.helical_parameters:
+            fu.log(f"Parameter {self.parameter} not valid. Valid parameters are: {constants.helical_parameters}.",
                    self.out_log, self.global_log)
             exit(1)
 
@@ -124,16 +122,16 @@ def launch(self) -> int:
             exit(1)
 
         # Check that the index is valid
-        if self.index < 0 or self.index >= len(self.sequence) - 1:
+        if self.index < 1 or self.index >= len(self.sequence) - 1:
             fu.log(f"Index {self.index} not valid. It should be between 0 and {len(self.sequence) - 2}.",
                    self.out_log, self.global_log)
             exit(1)
 
         # Get index sequence base and next base
-        bp = self.sequence[self.index] + self.sequence[self.index + 1]
+        bp = self.sequence[self.index-1] + self.sequence[self.index]
 
         # Get the filename
-        filename = f"{self.type}_{self.parameter}_{self.index + 1}_{bp}"
+        filename = f"{self.type}_{self.parameter}_{self.index}_{bp}"
         csv_file = f"{filename}.csv"
         jpg_file = f"{filename}.jpg"
 

biobb_dna/test/conf.yml

@@ -56,7 +56,7 @@ dna_timeseries:
     output_zip_path: timeseries_output.zip
     ref_output_zip_path: file:test_reference_dir/dna/timeseries_output.zip
   properties:
-    sequence: "CGCGAATTCGCG"
+    sequence: CGCGAATTCGCG
     seqpos: [4,5]
 
 dna_timeseries_unzip:

biobb_dna/test/reference/dna/dna_timeseries_unzip.csv

@@ -1,52 +1,48 @@
 shift,density
--1.87,2.0
--1.789607843137255,1.0
--1.70921568627451,2.0
--1.6288235294117648,3.0
--1.5484313725490197,5.0
--1.4680392156862747,7.0
--1.3876470588235295,4.0
--1.3072549019607844,15.0
--1.2268627450980394,18.0
--1.1464705882352944,26.0
--1.0660784313725493,36.0
--0.9856862745098041,64.0
--0.905294117647059,58.0
--0.824901960784314,91.0
--0.7445098039215687,123.0
--0.6641176470588237,164.0
--0.5837254901960787,194.0
--0.5033333333333336,217.0
--0.4229411764705886,241.0
--0.34254901960784334,291.0
--0.2621568627450983,264.0
--0.18176470588235327,296.0
--0.10137254901960802,311.0
--0.02098039215686298,330.0
-0.05941176470588205,278.0
-0.1398039215686273,317.0
-0.22019607843137212,279.0
-0.3005882352941174,227.0
-0.38098039215686264,223.0
-0.46137254901960745,188.0
-0.5417647058823527,164.0
-0.6221568627450975,144.0
-0.7025490196078428,101.0
-0.782941176470588,84.0
-0.8633333333333328,54.0
-0.9437254901960781,55.0
-1.024117647058823,31.0
-1.1045098039215682,29.0
-1.1849019607843134,14.0
-1.2652941176470582,11.0
-1.3456862745098035,8.0
-1.4260784313725487,11.0
-1.5064705882352936,8.0
-1.5868627450980388,2.0
-1.667254901960784,1.0
-1.7476470588235289,3.0
-1.8280392156862741,2.0
-1.908431372549019,1.0
-1.9888235294117642,0.0
-2.0692156862745095,1.0
-2.1496078431372547,1.0
+-2.52,1.0
+-2.42,1.0
+-2.32,3.0
+-2.2199999999999998,6.0
+-2.12,8.0
+-2.02,12.0
+-1.92,19.0
+-1.8199999999999998,28.0
+-1.72,40.0
+-1.62,39.0
+-1.52,60.0
+-1.42,97.0
+-1.3199999999999998,81.0
+-1.22,99.0
+-1.1199999999999999,108.0
+-1.02,163.0
+-0.9199999999999999,143.0
+-0.8199999999999998,187.0
+-0.72,245.0
+-0.6199999999999999,196.0
+-0.52,333.0
+-0.41999999999999993,258.0
+-0.31999999999999984,308.0
+-0.21999999999999975,339.0
+-0.11999999999999966,269.0
+-0.020000000000000018,362.0
+0.08000000000000007,298.0
+0.18000000000000016,278.0
+0.28000000000000025,264.0
+0.38000000000000034,177.0
+0.48,161.0
+0.5800000000000001,124.0
+0.6800000000000002,101.0
+0.7800000000000002,59.0
+0.8800000000000003,40.0
+0.98,32.0
+1.08,17.0
+1.1800000000000002,14.0
+1.2800000000000002,8.0
+1.3800000000000003,11.0
+1.48,6.0
+1.5800000000000005,1.0
+1.6800000000000002,1.0
+1.7799999999999998,2.0
+1.8800000000000003,0.0
+1.98,0.0
+2.0800000000000005,1.0

biobb_dna/test/reference/dna/dna_timeseries_unzip.txt

@@ -1,8 +1,8 @@
+series_shift_4_GA.csv
+series_shift_4_GA.jpg
+hist_shift_4_GA.csv
+hist_shift_4_GA.jpg
 series_shift_5_AA.csv
 series_shift_5_AA.jpg
 hist_shift_5_AA.csv
 hist_shift_5_AA.jpg
-series_shift_6_AT.csv
-series_shift_6_AT.jpg
-hist_shift_6_AT.csv
-hist_shift_6_AT.jpg