06 Apr 2026 15 min read computational biology

從 DMS 到機器學習：預測藥物是否能拯救突變蛋白的完整 Pipeline

這篇文章在講什麼？

想像你有一個蛋白質，因為基因突變而「摺疊失敗」——它無法正確摺疊，所以被細胞當作垃圾丟掉，無法到達它應該工作的位置（例如細胞表面），這就是很多罕見疾病的成因。

現在有一種叫做 藥理伴護蛋白（Pharmacological Chaperone, PC） 的小分子藥物，它可以跟蛋白質結合、幫助它正確摺疊、讓它順利被運送到該去的地方，但問題是：不是每個突變都能被救回來。

那麼，能不能用機器學習來預測「哪些突變可以被藥物拯救」？

這就是這篇論文在做的事，在 V2R 內部交叉驗證中，最佳模型 AUROC 為 0.751；將模型凍結後直接遷移到另一個蛋白質 Rhodopsin 時，AUROC 達到 0.843。

Part 1：背景知識-給剛入門的你

什麼是 GPCR？

GPCR（G protein-coupled receptor，G 蛋白偶聯受體）是人體中最大的受體蛋白家族，大約有 800 多個成員，它們橫跨在細胞膜上，負責接收外界訊號（荷爾蒙、神經傳導物質、光等等），然後把訊號傳進細胞內部。

你可能聽過的 GPCR 包括：多巴胺受體（精神科藥物的標靶）、腎上腺素受體（氣喘藥物的標靶）、以及本研究的主角——V2R（血管加壓素 2 型受體） 和 Rhodopsin（視紫紅質）。

V2R 負責腎臟的水分回收，V2R 突變會導致腎性尿崩症（nephrogenic diabetes insipidus）——患者的腎臟無法正常濃縮尿液。

Rhodopsin 是視網膜中負責感光的受體，Rhodopsin 突變會導致色素性視網膜炎（retinitis pigmentosa），最終可能失明。

什麼是 Deep Mutational Scanning（DMS）？

傳統上，研究一個突變的效果需要一個一個做實驗，但 DMS 技術可以一次測量一個蛋白質上幾乎所有可能的單點突變的效果。

以 V2R 為例：V2R 有 371 個胺基酸位置，每個位置可以突變成 19 種其他胺基酸，所以理論上有 371 × 19 = 7,049 個可能的錯義突變，DMS 實驗可以一次測量所有這些突變的表面表現量，不過經過原始資料的品質控制（排除 nonsense 突變、同義突變、以及缺少基線或藥物處理數值的變異），最終納入分析的是 6,649 個錯義突變。

2025 年的兩篇重要研究提供了本論文的資料基礎：

Mighell & Lehner (2025)：測量了 V2R 所有突變在有/無 tolvaptan（一種 PC）處理下的表面表現量
Manian et al. (2025)：測量了 Rhodopsin 所有突變在有/無 YC-001（另一種 PC）處理下的運輸效率

什麼是 Pharmacological Chaperone（藥理伴護蛋白）？

想像蛋白質摺疊就像摺紙：正常的蛋白質可以順利摺成正確的形狀，但突變可能讓某些步驟出錯，導致蛋白質被困在錯誤的中間狀態。

PC 就像一個「模具」——它跟蛋白質結合後，穩定了正確的構象，幫助蛋白質完成摺疊並通過細胞的品質管制系統。

在已上市的藥物中，Migalastat（治療法布瑞氏症，GLA 突變）是最典型的 pharmacological chaperone 範例，在更廣義的蛋白質穩定化與 corrector 脈絡中，也常提到 Elexacaftor/Tezacaftor（CFTR correctors，治療囊性纖維化）與 Tafamidis（transthyretin stabilizer，治療類澱粉心肌病），這些藥物的作用機制不完全相同，但都涉及透過小分子穩定蛋白質構象來恢復功能。

Part 2：把生物問題轉成 ML 問題

定義標籤：什麼是「有缺陷」？什麼是「被拯救」？

這是整個 pipeline 最關鍵的設計決策之一，我們需要把連續的實驗數值轉換成二元標籤。

Step 1：定義「有缺陷」（defective）的突變

V2R 的 DMS 資料中，每個突變都有一個 baseline surface expression score，野生型（正常）的分數被標準化為 1.0，完全失去功能的 nonsense 突變約為 0。

我選擇 0.7 作為閾值：baseline < 0.7 的突變被定義為「有缺陷」。

這個閾值的來源：V2R 同義突變（synonymous variants，不改變蛋白質序列的突變）的中位數是 1.001，0.7 大約是中位數減去 2 個標準差，這代表表現量顯著低於正常的突變。

Step 2：定義「被拯救」（rescued）

在有缺陷的突變中，如果加了藥物後表現量回到 0.7 以上，就標記為「rescued」（label = 1），否則就是「not rescued」（label = 0）。

V2R：6,649 個錯義突變
  → 2,133 個有缺陷（baseline < 0.7）
    → 1,477 個被拯救（69.2%）
    → 656 個未被拯救（30.8%）

閾值選擇會影響結果嗎？我在論文中做了 threshold sensitivity analysis，把閾值從 0.5 測到 0.8，AUROC 範圍是 0.734–0.753，結論不變，這種 robustness check 在投稿時幾乎是必做的。

選擇 Features：我們知道哪些跟「能不能被救」有關的資訊？

Feature engineering 是這個專案的核心，我把 features 分成兩層：

Layer 1（蛋白質無關的，可跨蛋白質遷移）：

Feature	它在衡量什麼	來源
AlphaMissense score	這個突變有多「致病」	DeepMind 的預訓練模型
RaSP ΔΔG	預測的熱力學穩定性變化	結構預測工具
ThermoMPNN ΔΔG	另一個穩定性預測	深度學習工具
ESM1b score	蛋白質語言模型的突變效果評分	Meta 的蛋白質語言模型
delta_hydro	疏水性變化	Kyte-Doolittle 量表
Baseline expression	未加藥時的表現量	DMS 實驗數據

Layer 2（蛋白質特異的，只用於機制分析）：

Feature	它在衡量什麼
到 AVP 結合位點的距離	突變離藥物結合口袋有多遠
到 Gs 偶聯界面的距離	突變離訊號傳導界面有多遠

為什麼分兩層？ 因為 Layer 2 的 features 是 V2R 特有的——Rhodopsin 用的是完全不同的藥物和結合位點，如果我們想做跨蛋白質預測，就只能用 Layer 1 的 features，這個「可遷移性」的設計從一開始就是刻意的。

Part 3：模型訓練與評估

Progressive Baseline：一步一步加 feature

很多論文只報告最終模型的效能，但這樣你無法理解「每個 feature 到底貢獻了多少」，我採用 progressive feature addition 的方式：

B0: 隨機猜測                    → AUROC 0.500
B1: AlphaMissense               → AUROC 0.575  (+0.075)
B2: + RaSP, ThermoMPNN          → AUROC 0.648  (+0.073)
B3: + ESM1b                     → AUROC 0.654  (+0.006)
B4: + delta_hydro               → AUROC 0.676  (+0.022)
B5: + baseline expression       → AUROC 0.751  (+0.075)

關鍵觀察：

AlphaMissense 單獨用效果很差（0.575），這告訴我們一件重要的事：「致病性」≠「可拯救性」，一個突變可以是致病的，但我們不知道它能不能被藥物救回來。
穩定性預測器（RaSP, ThermoMPNN）貢獻很大（+0.073）。 蛋白質的熱力學穩定性確實跟能不能被 PC 拯救有關。
ESM1b 幾乎沒有額外貢獻（+0.006），這可能是因為 ESM1b 的訊號已經被 AlphaMissense 和穩定性預測器捕捉了。
Baseline expression 是最強的單一 feature（+0.075），越嚴重的缺陷越難被拯救，這很直覺。

AUROC 是什麼？

AUROC（Area Under the Receiver Operating Characteristic curve）衡量的是分類器區分兩個類別的能力，0.5 代表跟隨機猜一樣，1.0 代表完美區分。0.75 大約代表「還不錯，有實用價值」。

你可以這樣想像：如果隨機抽一個被拯救的突變和一個未被拯救的突變，模型正確判斷誰是誰的機率就是 AUROC，0.751 代表大約 75% 的機率能判斷正確。

為什麼用 Random Forest？

我同時訓練了 Random Forest（RF）和 Logistic Regression（LR）做比較，RF 在 4/5 個 feature set 上都優於 LR，特別是在純計算 features（B4）時差距最大（+0.100）。

這代表 features 之間存在非線性交互作用，例如：一個中度不穩定的突變如果在藥物結合位點附近，可能就救不回來；但同樣中度不穩定的突變如果在其他位置，就可能被救回來，RF 可以捕捉這種條件式的規則，LR 不行。

# 核心訓練程式碼（簡化版）
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_predict

rf = RandomForestClassifier(
    n_estimators=200,
    class_weight='balanced',  # 處理類別不平衡
    random_state=42
)

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
y_pred = cross_val_predict(rf, X, y, cv=cv, method='predict_proba')[:, 1]  # 取正類（rescued）的機率

class_weight='balanced' 很重要。 我們的資料是 69.2% rescued vs 30.8% not-rescued，如果不處理這個不平衡，模型會傾向於把所有東西都預測成 rescued，balanced 會自動根據類別比例調整權重。

Positional Leakage Test：你的模型真的學到了什麼？

這是一個很多人會忽略但非常重要的驗證步驟。

問題是這樣的：蛋白質上不同位置有不同的結構角色，如果 cross-validation 的 fold 裡，訓練集和測試集包含了同一個位置的不同突變，模型可能只是學到了「哪些位置容易被救」，而不是真正學到了 variant-level 的特徵。

Leave-positions-out cross-validation 把整個 residue position 作為一個單位來分 fold——如果某個位置的突變在測試集中，那這個位置的所有突變都不會出現在訓練集中。

結果：

Standard 5-fold CV:          AUROC = 0.751
Leave-positions-out 10-fold: AUROC = 0.745
差距：0.006

差距幾乎為零，代表模型確實在學習 variant-level 的特徵，而不是在背位置資訊。

Part 4：機制分析——為什麼有些突變救不回來？

這是我最喜歡的部分，因為這不只是數字，而是可以用生物學解讀的分析。

根據特徵分布的差異，656 個未被拯救的突變可以區分成兩組，它們的特徵模式支持兩類主要 failure mode 的解釋：

第一類：藥物結合位點被破壞（NR-binding, n=150）

這些突變位於 tolvaptan 的結合位點（由原始論文通過統計分析鑑定的 32 個位置）。

特徵：

中度的 baseline expression（0.184）
小的疏水性變化（-0.46）
低的不穩定性分數

解釋： 蛋白質的整體結構沒有被嚴重破壞，但藥物結合的口袋被改變了，藥物進不去，自然就救不了。

第二類：蛋白質太嚴重地被破壞（NR-non-binding, n=506）

這些突變不在結合位點，但就是救不回來。

特徵：

極低的 baseline expression（0.064）
大的疏水性變化（-3.09）
高的不穩定性分數

解釋： 蛋白質的摺疊被嚴重破壞了，破壞程度超過了藥物所能提供的穩定化能力。就像一棟已經塌了一半的房子，光靠加一根柱子是撐不起來的。

這個觀察為什麼值得注意？ 如果這個解釋成立，它預測了一件事：不同的藥物-蛋白質組合，這兩種失敗模式的比例會不同，如果藥物的結合位點很深、很保守，第一類失敗就會比較少；如果結合位點很淺、很靈活，第一類失敗就會比較多，這是可以在未來實驗中驗證的預測。

Part 5：跨蛋白質遷移——外部驗證

為什麼這一步很重要？

如果模型只能預測 V2R 的突變，那它的價值有限——因為 V2R 的 DMS 資料已經有實驗結果了，模型的真正價值在於：能不能預測我們還沒有實驗資料的蛋白質？

我把 V2R 訓練好的模型完全凍結（frozen），不做任何調整，直接拿去預測 Rhodopsin 的突變能不能被 YC-001 拯救。

注意這兩個系統有多不同：

不同的蛋白質（V2R vs. Rhodopsin）
不同的藥物（tolvaptan vs. YC-001）
不同的實驗方法（flow cytometry vs. transcriptional reporter）
不同的疾病（腎性尿崩症 vs. 色素性視網膜炎）

只用四個 feature 的原因

六個 Layer 1 features 中，有兩個（RaSP 和 ThermoMPNN）因為技術問題無法為 Rhodopsin 計算：

RaSP 的 GitHub repository 已經不公開了
ThermoMPNN 有無法解決的 dependency 衝突

所以跨蛋白質預測只用了四個 features：AlphaMissense、baseline expression、delta_hydro、ESM1b。

提醒： 工具不可用是真實研究中常見的問題，誠實報告這個限制比假裝它不存在更重要，這也是為什麼我在論文的 Limitations 部分詳細說明了這一點。

結果：

Model A (AlphaMissense only):                    AUROC = 0.677
Model B (AlphaMissense + baseline):              AUROC = 0.788
Model D (+ hydrophobicity + ESM1b):              AUROC = 0.843

0.843 的 AUROC 甚至高於 V2R 內部的 0.751，不過要注意：兩個數字不能直接比較，因為實驗系統不同、標籤的雜訊特性不同，Rhodopsin 使用的 composite score 是一個 meta-analytic 指標，可能比 V2R 的單一實驗量測更乾淨。

Permutation Test：確認結果不是運氣

我做了 1,000 次 label permutation test：把 Rhodopsin 的標籤隨機打亂 1,000 次，每次都重新計算 AUROC。

null distribution 的平均值是 0.500（跟預期一致），而我們觀察到的 0.843 遠遠超過所有 permuted AUROC（p < 0.001）。

# Permutation test（簡化版）
import numpy as np
from sklearn.metrics import roc_auc_score

# y_pred_rho 為 frozen V2R model 對 Rhodopsin 預測的正類（rescued）機率
observed_auroc = roc_auc_score(y_rho, y_pred_rho)
n_permutations = 1000
permuted_aurocs = []

for _ in range(n_permutations):
    y_perm = np.random.permutation(y_rho)
    perm_auroc = roc_auc_score(y_perm, y_pred_rho)
    permuted_aurocs.append(perm_auroc)

p_value = np.mean(np.array(permuted_aurocs) >= observed_auroc)

附錄 A：Feature Engineering 實戰細節（進階）

以下內容適合有 Python 和 ML 經驗、想深入了解或重現此研究的讀者，如果你只想了解研究的整體思路和結果，可以跳過。

ESM1b Score 的計算

ESM1b 是 Meta 開發的蛋白質語言模型，類似 NLP 中的 BERT，但訓練在蛋白質序列上。它可以給每個位置的每個胺基酸一個「合理性」分數。

# ESM1b log-likelihood ratio 計算
import torch
import esm

# 載入模型
model, alphabet = esm.pretrained.esm1b_t33_650M_UR50S()
batch_converter = alphabet.get_batch_converter()
model.eval()

# 準備序列
data = [("rhodopsin", wild_type_sequence)]
batch_labels, batch_strs, batch_tokens = batch_converter(data)

# 取得每個位置每個胺基酸的 log probability
with torch.no_grad():
    results = model(batch_tokens, repr_layers=[33])
    logits = results["logits"]  # shape: [1, seq_len, vocab_size]

# 計算 log-likelihood ratio
# LLR = log P(mutant | context) - log P(wildtype | context)
log_probs = torch.log_softmax(logits, dim=-1)

解釋： 你可以把 ESM1b 想成一個「讀過」幾百萬條蛋白質序列的模型，當你問它「在這個位置放一個 Glycine 合理嗎？」，它會根據周圍的胺基酸序列給你一個分數，如果這個位置通常是疏水性胺基酸，你放一個親水性的胺基酸，分數就會很低。

Hydrophobicity Change 的計算

這是最簡單的 feature，但貢獻不小。

# Kyte-Doolittle hydrophobicity scale
kd_scale = {
    'A': 1.8, 'R': -4.5, 'N': -3.5, 'D': -3.5, 'C': 2.5,
    'Q': -3.5, 'E': -3.5, 'G': -0.4, 'H': -3.2, 'I': 4.5,
    'L': 3.8, 'K': -3.9, 'M': 1.9, 'F': 2.8, 'P': -1.6,
    'S': -0.8, 'T': -0.7, 'W': -0.9, 'Y': -1.3, 'V': 4.2
}

delta_hydro = kd_scale[mutant_aa] - kd_scale[wildtype_aa]

一個大的負 delta_hydro 通常代表疏水性胺基酸被替換成親水性胺基酸——如果這發生在蛋白質的疏水核心，就可能嚴重破壞摺疊。

附錄 B：重現這個研究（進階）

環境設置：

# 建議使用 conda 環境
conda create -n rescue python=3.10
conda activate rescue

# 安裝核心套件
pip install pandas numpy scikit-learn openpyxl matplotlib seaborn

# 如果要計算 ESM1b features
pip install fair-esm torch

資料取得：

V2R 資料： Mighell & Lehner (2025) 的 Supplementary Table S1（DOI: 10.1038/s41594-025-01659-6）
Rhodopsin 資料： Manian et al. (2025) 的 GitHub（https://github.com/octantbio/rho-dms）
分析程式碼： https://github.com/pcc402-art/project-rescue

重現步驟：

git clone https://github.com/pcc402-art/project-rescue.git
cd project-rescue

# Phase 1: 基礎模型訓練
python phase1_baseline.py

# Phase 3: 完整重跑（含所有 robustness checks）
python phase3_v2_full_rerun.py

# Phase 5: 最終投稿版分析
python phase5_presubmission.py

# 生成圖表
python figures_v2.py

關於作者

我是Po-Chun Chen，獨立研究者，研究方向為計算生物學與資訊安全，這是我以獨立研究者身份發表的第一篇論文。

如果你有問題、想討論，歡迎透過 [email protected] 聯繫我。

論文連結： Chen, P.-C. (2026). Computational prediction of pharmacological chaperone rescuability in GPCR missense variants with cross-protein transfer from V2R to Rhodopsin. Zenodo. https://doi.org/10.5281/zenodo.19434570

程式碼： https://github.com/pcc402-art/project-rescue