|
1 | 1 | import json |
2 | 2 | import logging |
3 | 3 | from pathlib import Path |
4 | | -from typing import Any, List, Tuple, Union, cast |
| 4 | +from typing import Any, Dict, List, Optional, Tuple, Union, cast |
5 | 5 |
|
6 | 6 | import numpy as np |
7 | 7 | import pandas as pd |
| 8 | +from sklearn.neighbors import KNeighborsRegressor |
8 | 9 | import torch |
9 | 10 | from numpy.typing import NDArray |
10 | 11 | from pandas import DataFrame |
|
14 | 15 | from scipy.special import softmax |
15 | 16 | from sklearn.metrics import recall_score |
16 | 17 |
|
| 18 | +from app.helpers.sequence_util import is_homolog_seq_id |
17 | 19 | from folde.types import FolDEModelConfig, ModelDiff, ModelEvaluation |
18 | 20 |
|
19 | 21 | DMS_SHORTNAMES = { |
20 | | - # 'A0A140D2T1_ZIKV_Sourisseau_2019': "A0A140D2T1_ZIKV", |
21 | | - # 'A0A2Z5U3Z0_9INFA_Doud_2016': "A0A2Z5U3Z0_9INFA", |
22 | | - # 'ADRB2_HUMAN_Jones_2020': "ADRB2", |
23 | | - # 'BLAT_ECOLX_Stiffler_2015': "BLAT_ECOLX", |
24 | | - # 'C6KNH7_9INFA_Lee_2018': "C6KNH7", |
25 | | - # 'IF1_ECOLI_Kelsic_2016': "IF1_ECOLI", |
26 | | - # 'MK01_HUMAN_Brenan_2016': "MK01_HUMAN", |
27 | | - # 'P53_HUMAN_Giacomelli_2018_Null_Etoposide': "P53", |
28 | | - # 'PHOT_CHLRE_Chen_2023': "PHOT_CHLRE", |
29 | 22 | "BLAT_ECOLX_Firnberg_2014": "BLAT_ECOLX", |
30 | 23 | "ANCSZ_Hobbs_2022": "ANCSZ", |
31 | 24 | "HXK4_HUMAN_Gersing_2022_activity": "HXK4_HUMAN", |
|
36 | 29 | "P53_HUMAN_Giacomelli_2018_Null_Nutlin": "P53_Null", |
37 | 30 | "HSP82_YEAST_Flynn_2019": "HSP82_YEAST", |
38 | 31 | "P53_HUMAN_Giacomelli_2018_WT_Nutlin": "P53_WT", |
39 | | - # "MK01_HUMAN_Brenan_2016": "MK01_HUMAN", |
40 | 32 | "HEM3_HUMAN_Loggerenberg_2023": "HEM3_HUMAN", |
41 | 33 | "PPM1D_HUMAN_Miller_2022": "PPM1D_HUMAN", |
42 | 34 | "SPG1_STRSG_Olson_2014": "SPG1", |
43 | | - # VALIDATION TARGETS |
44 | 35 | "ADRB2_HUMAN_Jones_2020": "ADRB2_HUMAN", |
45 | 36 | "P53_HUMAN_Giacomelli_2018_Null_Nutlin": "P53_HUMAN_Null", |
46 | 37 | "P53_HUMAN_Giacomelli_2018_WT_Nutlin": "P53_HUMAN_WT", |
47 | 38 | "MK01_HUMAN_Brenan_2016": "MK01_HUMAN", |
48 | 39 | "KCNJ2_MOUSE_Coyote-Maestas_2022_function": "KCNJ2_MOUSE", |
49 | 40 | "CAS9_STRP1_Spencer_2017_positive": "CAS9_STRP1", |
50 | 41 | "SC6A4_HUMAN_Young_2021": "SC6A4_HUMAN", |
51 | | - # 'OXDA_RHOTO_Vanella_2023_expression', |
52 | | - # 'HSP82_YEAST_Mishra_2016', |
53 | 42 | "PTEN_HUMAN_Mighell_2018": "PTEN_HUMAN", |
54 | 43 | "S22A1_HUMAN_Yee_2023_activity": "S22A1_HUMAN", |
55 | 44 | "KKA2_KLEPN_Melnikov_2014": "KKA2_KLEPN", |
56 | | - # Include some "easy to engineer" targets. |
57 | 45 | "PPARG_HUMAN_Majithia_2016": "PPARG_HUMAN", |
58 | | - # 'P53_HUMAN_Giacomelli_2018_Null_Etoposide', |
59 | 46 | "MET_HUMAN_Estevam_2023": "MET_HUMAN", |
60 | 47 | "MTHR_HUMAN_Weile_2021": "MTHR_HUMAN", |
61 | 48 | "LGK_LIPST_Klesmith_2015": "LGK_LIPST", |
@@ -595,3 +582,110 @@ def load_checkpointed_model_eval( |
595 | 582 | raise ValueError(f"Failed to find any matching files.") |
596 | 583 |
|
597 | 584 | return composite_result |
| 585 | + |
| 586 | + |
| 587 | +class NaturalnessImputer(object): |
| 588 | + def __init__(self): |
| 589 | + self.is_pretrained = False |
| 590 | + self.single_mutant_naturalness_df: Optional[pd.DataFrame] = None |
| 591 | + self.knns: Dict[str, KNeighborsRegressor] = {} |
| 592 | + |
| 593 | + def pretrain(self, naturalness_df: pd.DataFrame, embedding_series: pd.Series): |
| 594 | + if self.is_pretrained: |
| 595 | + raise ValueError("Model is already pretrained.") |
| 596 | + |
| 597 | + self.knns = {} |
| 598 | + for naturalness_column in naturalness_df.columns: |
| 599 | + naturalness_series = naturalness_df[naturalness_column] |
| 600 | + assert naturalness_series.index.equals(embedding_series.index) |
| 601 | + assert embedding_series.index.is_unique, "embedding_series contains duplicate indices" |
| 602 | + |
| 603 | + assert not naturalness_series.isna().any(), "naturalness_series contains NANs" |
| 604 | + |
| 605 | + X = np.array([np.array(emb) for emb in embedding_series.values]) |
| 606 | + y = naturalness_series.values |
| 607 | + |
| 608 | + # Fit KNN regressor |
| 609 | + knn = KNeighborsRegressor(n_neighbors=5) |
| 610 | + knn.fit(X, y) |
| 611 | + self.knns[naturalness_column] = knn |
| 612 | + |
| 613 | + # Also store the naturalness series for prediction later. |
| 614 | + self.single_mutant_naturalness_df = naturalness_df |
| 615 | + self.is_pretrained = True |
| 616 | + |
| 617 | + def impute(self, naturalness_df: pd.DataFrame, embedding_series: pd.Series) -> List[pd.Series]: |
| 618 | + assert self.single_mutant_naturalness_df is not None |
| 619 | + assert len(naturalness_df.columns) == len(self.knns) |
| 620 | + assert set(naturalness_df.columns) == set(self.single_mutant_naturalness_df.columns) |
| 621 | + |
| 622 | + # Get the ensemble of naturalness scores with missing values imputed. |
| 623 | + |
| 624 | + ensemble_of_computed_naturalness: List[pd.Series] = [] |
| 625 | + |
| 626 | + for naturalness_column in naturalness_df.columns: |
| 627 | + |
| 628 | + naturalness_series = naturalness_df[naturalness_column] |
| 629 | + single_mutant_naturalness_series = self.single_mutant_naturalness_df[naturalness_column] |
| 630 | + knn = self.knns[naturalness_column] |
| 631 | + |
| 632 | + def get_naturalness(seq_id, direct_naturalness) -> float: |
| 633 | + """Try computing naturalness for mutants even if none was provided by extrapolating for multimutants.""" |
| 634 | + if direct_naturalness is not None and not pd.isna(direct_naturalness): |
| 635 | + return direct_naturalness |
| 636 | + |
| 637 | + if is_homolog_seq_id(seq_id): |
| 638 | + return np.nan |
| 639 | + |
| 640 | + # Break it down into single mutants. |
| 641 | + seq_id_parts = seq_id.split("_") |
| 642 | + |
| 643 | + # For multimutants, we compute naturalness as the product of the naturalness of the single mutants. |
| 644 | + computed_naturalness = single_mutant_naturalness_series.loc[seq_id_parts].sum() |
| 645 | + if pd.isna(computed_naturalness): |
| 646 | + raise ValueError( |
| 647 | + f"Computed naturalness is NAN for {seq_id} with parts {seq_id_parts}" |
| 648 | + ) |
| 649 | + return computed_naturalness |
| 650 | + |
| 651 | + naturalness_series.index.name = "seq_id" |
| 652 | + computed_naturalness_series = naturalness_series.reset_index(name="wt_marginal").apply( |
| 653 | + lambda r: get_naturalness(r.seq_id, r.wt_marginal), axis=1 |
| 654 | + ) |
| 655 | + computed_naturalness_series.index = naturalness_series.index |
| 656 | + |
| 657 | + # Do KNN imputation to fill in NANs from homologs. |
| 658 | + if computed_naturalness_series.isna().any(): |
| 659 | + if not self.is_pretrained: |
| 660 | + raise ValueError( |
| 661 | + "Model is not pretrained, so cannot fill in NANs from homologs." |
| 662 | + ) |
| 663 | + |
| 664 | + logging.info( |
| 665 | + f"Filling in NANs from homologs for {computed_naturalness_series.isna().sum()}/{len(computed_naturalness_series)} naturalness values." |
| 666 | + ) |
| 667 | + assert embedding_series is not None |
| 668 | + embedding_array = np.array([np.array(emb) for emb in embedding_series.values]) |
| 669 | + naturalness_array = computed_naturalness_series.values |
| 670 | + |
| 671 | + # Find indices for known and missing |
| 672 | + missing_mask = computed_naturalness_series.isna().to_numpy() |
| 673 | + X_missing = embedding_array[missing_mask] |
| 674 | + |
| 675 | + imputed_values = knn.predict(X_missing) |
| 676 | + |
| 677 | + # Fill in the missing values |
| 678 | + imputed_naturalness = naturalness_array.copy() |
| 679 | + imputed_naturalness[missing_mask] = imputed_values |
| 680 | + |
| 681 | + # Convert back to Series |
| 682 | + computed_naturalness_series = pd.Series( |
| 683 | + imputed_naturalness, index=naturalness_series.index |
| 684 | + ) |
| 685 | + |
| 686 | + if computed_naturalness_series.isna().any(): |
| 687 | + raise ValueError( |
| 688 | + f"Computed naturalness series still has NANs: {computed_naturalness_series.isna().sum()}/{len(computed_naturalness_series)}" |
| 689 | + ) |
| 690 | + ensemble_of_computed_naturalness.append(computed_naturalness_series) |
| 691 | + return ensemble_of_computed_naturalness |
0 commit comments