Vectorize percentile scoring

src/kinase_library/enrichment/binary_enrichment.py

@@ -12,6 +12,7 @@
 
 from ..utils import _global_vars, exceptions, utils
 from ..modules import data, enrichment
+from ..objects import core
 from ..objects import phosphoproteomics as pps
 
 #%%
@@ -277,21 +278,22 @@ def kinase_enrichment(self, kin_type, kl_method, kl_thresh,
                 print('Not all kinase scores were provided. Re-calculating scores for background data')
                 self.bg_pps.score(kin_type=kin_type,kinases=kinases)
         elif kl_method in ['percentile','percentile_rank']:
+            scored_phosprot = core.ScoredPhosphoProteome(phosprot_name=_global_vars.phosprot_name)
             if not hasattr(self.fg_pps, kin_type + '_' + data_att) or rescore:
                 print('\nCalculating percentiles for foreground data')
-                self.fg_pps.percentile(kin_type=kin_type,kinases=kinases)
+                self.fg_pps.percentile(kin_type=kin_type,kinases=kinases,customized_scored_phosprot=scored_phosprot)
                 self.phosprot_name = _global_vars.phosprot_name
             elif not set(kinases)<=set(getattr(self.fg_pps, kin_type + '_' + data_att)):
                 print('Not all kinase percentiles were provided. Re-calculating percentiles for foreground data')
-                self.fg_pps.percentile(kin_type=kin_type,kinases=kinases)
+                self.fg_pps.percentile(kin_type=kin_type,kinases=kinases,customized_scored_phosprot=scored_phosprot)
                 self.phosprot_name = _global_vars.phosprot_name
             if not hasattr(self.bg_pps, kin_type + '_' + data_att) or rescore:
                 print('\nCalculating percentiles for background data')
-                self.bg_pps.percentile(kin_type=kin_type,kinases=kinases)
+                self.bg_pps.percentile(kin_type=kin_type,kinases=kinases,customized_scored_phosprot=scored_phosprot)
                 self.phosprot_name = _global_vars.phosprot_name
             elif not set(kinases)<=set(getattr(self.bg_pps, kin_type + '_' + data_att)):
                 print('Not all kinase percentiles were provided. Re-calculating percentiles for background data')
-                self.bg_pps.percentile(kin_type=kin_type,kinases=kinases)
+                self.bg_pps.percentile(kin_type=kin_type,kinases=kinases,customized_scored_phosprot=scored_phosprot)
                 self.phosprot_name = _global_vars.phosprot_name
 
         fg_score_data = getattr(self.fg_pps, kin_type + '_' + data_att)

src/kinase_library/enrichment/differential_phosphorylation.py

@@ -12,6 +12,7 @@
 
 from ..utils import _global_vars, exceptions, utils
 from ..modules import data, enrichment
+from ..objects import core
 from ..objects import phosphoproteomics as pps
 from ..enrichment import binary_enrichment as be
 from ..logger import logger
@@ -296,15 +297,16 @@ def kinase_enrichment(self, kin_type, kl_method, kl_thresh,
                     hasattr(self.upreg_sites_pps, kin_type+'_percentile_ranks') and
                     hasattr(self.downreg_sites_pps, kin_type+'_percentile_ranks') and
                     hasattr(self.unreg_sites_pps, kin_type+'_percentile_ranks')) or rescore:
+                scored_phosprot = core.ScoredPhosphoProteome(phosprot_name=_global_vars.phosprot_name)
                 print('\nCalculating percentiles for upregulated sites ({} substrates)'.format(len(self.upreg_sites_data)))
                 logger.info('Calculating percentiles for upregulated sites ({} substrates)'.format(len(self.upreg_sites_data)))
-                upreg_sites_percentile = self.upreg_sites_pps.percentile(kin_type=kin_type, kinases=kinases, non_canonical=non_canonical, values_only=True)
+                upreg_sites_percentile = self.upreg_sites_pps.percentile(kin_type=kin_type, kinases=kinases, non_canonical=non_canonical, values_only=True, customized_scored_phosprot=scored_phosprot)
                 print('\nCalculating percentiles for downregulated sites ({} substrates)'.format(len(self.downreg_sites_data)))
                 logger.info('Calculating percentiles for downregulated sites ({} substrates)'.format(len(self.downreg_sites_data)))
-                downreg_sites_percentile = self.downreg_sites_pps.percentile(kin_type=kin_type, kinases=kinases, non_canonical=non_canonical, values_only=True)
+                downreg_sites_percentile = self.downreg_sites_pps.percentile(kin_type=kin_type, kinases=kinases, non_canonical=non_canonical, values_only=True, customized_scored_phosprot=scored_phosprot)
                 print('\nCalculating percentiles for background (unregulated) sites ({} substrates)'.format(len(self.unreg_sites_data)))
                 logger.info('Calculating percentiles for background (unregulated) sites ({} substrates)'.format(len(self.unreg_sites_data)))
-                unreg_sites_percentile = self.unreg_sites_pps.percentile(kin_type=kin_type, kinases=kinases, non_canonical=non_canonical, values_only=True)
+                unreg_sites_percentile = self.unreg_sites_pps.percentile(kin_type=kin_type, kinases=kinases, non_canonical=non_canonical, values_only=True, customized_scored_phosprot=scored_phosprot)
                 self.phosprot_name = _global_vars.phosprot_name
             else:
                 upreg_sites_percentile = getattr(self.upreg_sites_pps, kin_type+'_percentiles')

src/kinase_library/objects/core.py

@@ -1286,6 +1286,18 @@ def __init__(self, phosprot_name, kin_type=None, phosprot_file=None, phosprot_pa
 
         self.log2_values = log2_values
 
+    def get_sorted_scores(self, kin_type):
+        """Return column-sorted score array, computing and caching on first call."""
+        cache_attr = '_' + kin_type + '_sorted_scores'
+        if not hasattr(self, cache_attr):
+            scores = getattr(self, kin_type + '_scores')
+            if not self.log2_values:
+                setattr(self, cache_attr, np.sort(np.log2(scores.values), axis=0))
+            else:
+                setattr(self, cache_attr, np.sort(scores.values, axis=0))
+            setattr(self, '_' + kin_type + '_sorted_columns', list(scores.columns))
+        return getattr(self, cache_attr), getattr(self, '_' + kin_type + '_sorted_columns')
+
     def merge_data_scores(self, kin_type, data_seq_col=None):
         """
         Merging phosphoproteome data and scores

src/kinase_library/objects/phosphoproteomics.py

@@ -6,6 +6,7 @@
 import numpy as np
 import pandas as pd
 import pyarrow.parquet as pq
+from tqdm import tqdm
 
 from ..utils import _global_vars, exceptions, utils
 from ..modules import data
@@ -441,22 +442,22 @@ def percentile(self, kin_type=None, kinases=None,
             exceptions.check_kin_type(kin_type)
         exceptions.check_kin_list_type(kinases, kin_type=kin_type)
 
-        if kin_type == 'ser_thr':
-            scored_phosprot = all_scored_phosprot.ser_thr_scores
-        elif kin_type == 'tyrosine':
-            scored_phosprot = all_scored_phosprot.tyrosine_scores
-        else:
+        if kin_type not in ('ser_thr', 'tyrosine'):
             raise ValueError('Wrong kinase type.')
-        scored_phosprot = scored_phosprot.loc[:,kinases] # only for requested kinases if subset
 
-        # If scored phopshoproteome is linear values - converting it to log2 values
-        if not all_scored_phosprot.log2_values:
-            scored_phosprot = np.log2(scored_phosprot)
+        sorted_scores, sorted_columns = all_scored_phosprot.get_sorted_scores(kin_type)
+        col_indices = [sorted_columns.index(k) for k in kinases]
+        sorted_subset = sorted_scores[:, col_indices]
+        n_phosprot = sorted_subset.shape[0]
+        score_vals = score[kinases].values
 
         print('Calculating percentile for '+str(len(getattr(self,kin_type+'_substrates')))+' '+kin_type+' substrates')
         logger.info('Calculating percentile for '+str(len(getattr(self,kin_type+'_substrates')))+' '+kin_type+' substrates')
-        percent_output = scored_phosprot.progress_apply(lambda x: x.sort_values().searchsorted(score[x.name], side='right'))/len(scored_phosprot)*100
-        percent_output.index = score.index
+        result = np.empty_like(score_vals, dtype=float)
+        for j in tqdm(range(sorted_subset.shape[1]), desc='Percentile scoring'):
+            result[:, j] = np.searchsorted(sorted_subset[:, j], score_vals[:, j], side='right')
+        percent_output = pd.DataFrame(result / n_phosprot * 100,
+                                      index=score.index, columns=kinases)
 
         percent_output = percent_output.round(round_digits)
         percent_rank_output = percent_output.rank(method='min', ascending=False, axis=1).astype(int)