updated the saturation function to now force the y-intercept to be ze…

…ro when calculating slope. scikit-learn is a new dependency.

Pipfile

@@ -8,6 +8,7 @@ verify_ssl = true
 [packages]
 biopython = "==1.76"
 scipy = "==1.4.1"
+scikit-learn = "==1.4.2"
 
 [requires]
 python_version = "3.9"

change_log.txt

@@ -1,5 +1,8 @@
 Major changes to PhyKIT are summarized here.
 
+1.19.4
+  - Saturation function forces y-intercept to be zero when calculating slope
+
 1.19.3
   - Saturation function now uses uncorrected distances instead of pairwise identities
 

docs/change_log/index.rst

@@ -8,6 +8,9 @@ Change log
 
 Major changes to PhyKIT are summarized here.
 
+**1.19.4**:
+Saturation function forces y-intercept to be zero when calculating slope
+
 **1.19.3**:
 Saturation function now uses uncorrected distances instead of pairwise identities
 

phykit/services/tree/saturation.py

@@ -2,7 +2,8 @@
 import itertools
 from typing import Tuple
 
-import scipy
+import numpy as np
+from sklearn.linear_model import LinearRegression
 
 from .base import Tree
 from ...helpers.files import get_alignment_and_format as get_alignment_and_format_helper
@@ -43,14 +44,19 @@ def run(self):
             uncorrected_distances,
         ) = self.loop_through_combos_and_calculate_pds_and_pis(combos, alignment, tree)
 
-        # calculate linear regression
-        _, _, r_value, _, _ = scipy.stats.linregress(
-            uncorrected_distances, patristic_distances
-        )
 
+
+        # calculate slope and fit the y-intercept to zero
+        # Fitting the y-intercept to zero follows Jeffroy et al.
+        # See fig 2 https://www.cell.com/trends/genetics/fulltext/S0168-9525(06)00051-5
+        model = LinearRegression(fit_intercept=False)
+        model.fit(
+            np.array(patristic_distances).reshape(-1, 1),
+            np.array(uncorrected_distances)
+        )
         # report res
         self.print_res(
-            self.verbose, combos, uncorrected_distances, patristic_distances, r_value
+            self.verbose, combos, uncorrected_distances, patristic_distances, model.coef_[0]
         )
 
     def process_args(self, args):
@@ -107,7 +113,7 @@ def print_res(
         combos: list,
         uncorrected_distances: list,
         patristic_distances: list,
-        r_value: float,
+        slope: float,
     ) -> None:
         """
         print results to stdout
@@ -121,6 +127,6 @@ def print_res(
                         f"{combo[0]}-{combo[1]}\t{round(dist,4)}\t{round(patristic_distance, 4)}"
                     )
             else:
-                print(round(r_value**2, 4))
+                print(round(slope, 4))
         except BrokenPipeError:
             pass

phykit/version.py

@@ -1 +1 @@
-__version__ = "1.19.3"
+__version__ = "1.19.4"

requirements.txt

@@ -1,4 +1,5 @@
 biopython>=1.81
 numpy>=1.24.0
 scipy>=1.11.3
+scikit-learn>=1.4.2
 cython

setup.py

@@ -21,6 +21,7 @@
     "biopython>=1.81",
     "numpy>=1.24.0",
     "scipy>=1.11.3",
+    "scikit-learn>=1.4.2",
     "cython"
 ]
 

tests/integration/tree/test_saturation_integration.py

@@ -14,7 +14,7 @@
 class TestSaturation(object):
     @patch("builtins.print")
     def test_saturation(self, mocked_print):
-        expected_result = 0.8451
+        expected_result = 0.4919
         testargs = [
             "phykit",
             "saturation",
@@ -29,7 +29,7 @@ def test_saturation(self, mocked_print):
 
     @patch("builtins.print")
     def test_saturation_alias(self, mocked_print):
-        expected_result = 0.8451
+        expected_result = 0.4919
         testargs = [
             "phykit",
             "sat",