22import polars as pl
33from pathlib import Path
44from Bio import SeqIO
5- from ._rs import rs_preprocess , rs_match , rs_discontinuous
5+ from ._rs import rs_preprocess , rs_match , rs_discontinuous , rs_metadata
66
77VALID_OUTPUT_FORMATS = ['dataframe' , 'csv' , 'tsv' , 'xlsx' , 'json' ]
88FASTA_EXTENSIONS = {
@@ -160,7 +160,17 @@ def _search(self, k, max_mismatches, peptides=None):
160160
161161 def best_match_search (self ):
162162 peptides_remaining = self .query .copy ()
163- all_results = []
163+ acc = tuple ([] for _ in range (8 )) # 8 columnar accumulators
164+
165+ def collect_matched (cols ):
166+ matched_ids = set ()
167+ matched_col = cols [2 ]
168+ for i in range (len (cols [0 ])):
169+ if matched_col [i ] is not None :
170+ for j in range (8 ):
171+ acc [j ].append (cols [j ][i ])
172+ matched_ids .add (cols [0 ][i ])
173+ return matched_ids
164174
165175 initial_k = min (len (seq ) for _ , seq in peptides_remaining )
166176 ks = [initial_k ]
@@ -172,52 +182,34 @@ def best_match_search(self):
172182 for k in ks :
173183 if not peptides_remaining :
174184 break
175-
176185 min_len = min (len (seq ) for _ , seq in peptides_remaining )
177186 max_mm = (min_len // k ) - 1
178187 if max_mm < 0 :
179188 max_mm = 0
180-
181- results = self ._search (k , max_mm , peptides_remaining )
182-
183- matched_ids = set ()
184- for row in results :
185- if row [2 ] != '' :
186- all_results .append (row )
187- matched_ids .add (row [0 ])
188-
189+ matched_ids = collect_matched (self ._search (k , max_mm , peptides_remaining ))
189190 peptides_remaining = [
190191 (qid , seq ) for qid , seq in peptides_remaining if qid not in matched_ids
191192 ]
192193 print (f" -> k={ k } , mismatches<={ max_mm } : { len (peptides_remaining )} remaining" )
193194
194195 if peptides_remaining :
195196 max_mm = min (len (seq ) for _ , seq in peptides_remaining ) // 2
196-
197197 while peptides_remaining :
198198 shortest_len = min (len (seq ) for _ , seq in peptides_remaining )
199199 if max_mm >= shortest_len :
200200 break
201-
202201 max_mm += 1
203- results = self ._search (2 , max_mm , peptides_remaining )
204-
205- matched_ids = set ()
206- for row in results :
207- if row [2 ] != '' :
208- all_results .append (row )
209- matched_ids .add (row [0 ])
210-
202+ matched_ids = collect_matched (self ._search (2 , max_mm , peptides_remaining ))
211203 peptides_remaining = [
212204 (qid , seq ) for qid , seq in peptides_remaining if qid not in matched_ids
213205 ]
214206 print (f" -> k=2, mismatches<={ max_mm } : { len (peptides_remaining )} remaining" )
215207
216- if peptides_remaining :
217- for qid , seq in peptides_remaining :
218- all_results .append ([ qid , seq ] + [ '' ] * 14 )
208+ for qid , seq in peptides_remaining : # leftover unmatched
209+ for j , val in enumerate (( qid , seq , None , None , None , "[]" , None , None )) :
210+ acc [ j ] .append (val )
219211
220- df = self ._to_dataframe (all_results )
212+ df = self ._to_dataframe (acc )
221213 return self ._best_match_filter (df )
222214
223215 def _best_match_filter (self , df ):
@@ -274,46 +266,63 @@ def _best_match_filter(self, df):
274266
275267 return pl .concat ([matched_df , unmatched_df ], how = "vertical" )
276268
277- def _to_dataframe (self , results ):
278- schema = {
279- 'Query ID' : pl .Utf8 ,
280- 'Query Sequence' : pl .Utf8 ,
281- 'Matched Sequence' : pl .Utf8 ,
282- 'Protein ID' : pl .Utf8 ,
283- 'Protein Name' : pl .Utf8 ,
284- 'Species' : pl .Utf8 ,
285- 'Taxon ID' : pl .Utf8 ,
286- 'Gene' : pl .Utf8 ,
287- 'Mismatches' : pl .Utf8 ,
288- 'Mutated Positions' : pl .Utf8 ,
289- 'Index start' : pl .Utf8 ,
290- 'Index end' : pl .Utf8 ,
291- 'Protein Existence Level' : pl .Utf8 ,
292- 'Sequence Version' : pl .Utf8 ,
293- 'Gene Priority' : pl .Utf8 ,
294- 'SwissProt Reviewed' : pl .Utf8 ,
295- }
296-
297- if not results :
298- return pl .DataFrame (schema = schema )
299-
300- df = pl .DataFrame (results , schema = schema , orient = 'row' )
301-
302- df = df .with_columns (
303- pl .when (pl .col (col ) == '' ).then (None ).otherwise (pl .col (col )).alias (col )
304- for col in df .columns
269+ def _metadata_table (self ) -> pl .DataFrame :
270+ """Per-protein metadata (built once from this proteome's index) for the edge join."""
271+ import glob
272+ pattern = os .path .join (self .preprocessed_files_path , f'{ self .proteome_name } _*mers.pepidx' )
273+ idx_files = glob .glob (pattern )
274+ if not idx_files :
275+ raise FileNotFoundError (f'No .pepidx for { self .proteome_name } in { self .preprocessed_files_path } ' )
276+ pnum , pid , name , species , taxon , gene , exist , seqver , geneprio , swiss = rs_metadata (idx_files [0 ])
277+ m = pl .DataFrame ({
278+ 'protein_num' : pl .Series (pnum , dtype = pl .UInt32 ),
279+ 'Protein ID' : pid , 'Protein Name' : name , 'Species' : species , 'Taxon ID' : taxon ,
280+ 'Gene' : gene , 'Protein Existence Level' : exist , 'Sequence Version' : seqver ,
281+ 'Gene Priority' : geneprio , 'SwissProt Reviewed' : swiss ,
282+ })
283+ str_cols = ['Protein ID' ,'Protein Name' ,'Species' ,'Taxon ID' ,'Gene' ,
284+ 'Protein Existence Level' ,'Sequence Version' ,'Gene Priority' ,'SwissProt Reviewed' ]
285+ m = m .with_columns (
286+ pl .when (pl .col (c ) == '' ).then (None ).otherwise (pl .col (c )).alias (c ) for c in str_cols
305287 )
306-
307- df = df .with_columns ([
308- pl .col ('Mismatches' ).cast (pl .Int64 ),
309- pl .col ('Index start' ).cast (pl .Int64 ),
310- pl .col ('Index end' ).cast (pl .Int64 ),
288+ return m .with_columns ([
311289 pl .col ('Protein Existence Level' ).cast (pl .Int64 ),
312290 pl .col ('Sequence Version' ).cast (pl .Int64 ),
313291 pl .col ('Gene Priority' ).cast (pl .Int64 ),
314292 pl .col ('SwissProt Reviewed' ).cast (pl .Int64 ).cast (pl .Boolean ),
315293 ])
316294
295+ FINAL_COLUMNS = [
296+ 'Query ID' ,'Query Sequence' ,'Matched Sequence' ,'Protein ID' ,'Protein Name' ,'Species' ,
297+ 'Taxon ID' ,'Gene' ,'Mismatches' ,'Mutated Positions' ,'Index start' ,'Index end' ,
298+ 'Protein Existence Level' ,'Gene Priority' ,'SwissProt Reviewed' ,
299+ ]
300+
301+ def _to_dataframe (self , cols ):
302+ """Build the result DataFrame from columnar Rust output, reconstructing protein
303+ metadata via a single join (instead of cloning it into every hit row)."""
304+ qid , qseq , matched , pnum , mm , mutated , istart , iend = cols
305+
306+ if not qid :
307+ schema = {c : pl .Utf8 for c in self .FINAL_COLUMNS }
308+ for c in ('Mismatches' ,'Index start' ,'Index end' ,'Protein Existence Level' ,'Gene Priority' ):
309+ schema [c ] = pl .Int64
310+ schema ['SwissProt Reviewed' ] = pl .Boolean
311+ return pl .DataFrame (schema = schema )
312+
313+ base = pl .DataFrame ({
314+ 'Query ID' : qid ,
315+ 'Query Sequence' : qseq ,
316+ 'Matched Sequence' : matched ,
317+ 'protein_num' : pl .Series (pnum , dtype = pl .UInt32 ),
318+ 'Mismatches' : pl .Series (mm , dtype = pl .Int64 ),
319+ 'Mutated Positions' : mutated ,
320+ 'Index start' : pl .Series (istart , dtype = pl .Int64 ),
321+ 'Index end' : pl .Series (iend , dtype = pl .Int64 ),
322+ })
323+
324+ df = base .join (self ._metadata_table (), on = 'protein_num' , how = 'left' ).drop ('protein_num' )
325+
317326 if self .sequence_version :
318327 df = df .with_columns (
319328 pl .when (pl .col ('Sequence Version' ).is_not_null ())
@@ -322,5 +331,4 @@ def _to_dataframe(self, results):
322331 .alias ('Protein ID' )
323332 )
324333
325- df = df .drop ('Sequence Version' )
326- return df
334+ return df .drop ('Sequence Version' ).select (self .FINAL_COLUMNS )
0 commit comments