1. 從數(shù)學(xué)建模到真實(shí)世界一次抗乳腺癌藥物研究的實(shí)戰(zhàn)復(fù)盤去年我參與指導(dǎo)了一個(gè)由學(xué)生團(tuán)隊(duì)完成的數(shù)學(xué)建模項(xiàng)目題目恰好是“基于大數(shù)據(jù)與機(jī)器學(xué)習(xí)的抗乳腺癌藥物研究”。這聽(tīng)起來(lái)像是一個(gè)典型的、充滿學(xué)術(shù)氣息的競(jìng)賽題目但當(dāng)我們真正扎進(jìn)去把那些抽象的“模型”、“算法”和“數(shù)據(jù)”變成可運(yùn)行的代碼和可解釋的結(jié)果時(shí)整個(gè)過(guò)程遠(yuǎn)比想象中復(fù)雜和有趣。它不再是一道簡(jiǎn)單的“題”而是一個(gè)微縮版的真實(shí)世界科研項(xiàng)目。今天我想拋開(kāi)那些華麗的獲獎(jiǎng)證書和論文摘要以一個(gè)親歷者的視角復(fù)盤我們是如何一步步將“數(shù)學(xué)建模C題大數(shù)據(jù)機(jī)器學(xué)習(xí)”這些宏大概念落地成一個(gè)有血有肉、有坑有收獲的具體項(xiàng)目。如果你也對(duì)如何將機(jī)器學(xué)習(xí)應(yīng)用于生物醫(yī)學(xué)數(shù)據(jù)分析或者對(duì)數(shù)學(xué)建模競(jìng)賽的實(shí)戰(zhàn)過(guò)程感到好奇那么這篇復(fù)盤或許能給你一些不一樣的啟發(fā)。這個(gè)項(xiàng)目的核心目標(biāo)很明確利用公開(kāi)的乳腺癌相關(guān)生物醫(yī)學(xué)大數(shù)據(jù)如基因表達(dá)數(shù)據(jù)、藥物敏感性數(shù)據(jù)、臨床病理數(shù)據(jù)構(gòu)建機(jī)器學(xué)習(xí)模型來(lái)預(yù)測(cè)藥物對(duì)特定乳腺癌亞型的療效并嘗試挖掘潛在的生物標(biāo)志物或藥物作用機(jī)制。這本質(zhì)上是一個(gè)高維、小樣本、強(qiáng)噪聲的預(yù)測(cè)與關(guān)聯(lián)分析問(wèn)題在生物信息學(xué)和計(jì)算生物學(xué)領(lǐng)域非常典型。我們的工作流大致可以拆解為數(shù)據(jù)獲取與理解、特征工程與降維、模型構(gòu)建與優(yōu)化、結(jié)果解釋與生物意義挖掘這四個(gè)環(huán)環(huán)相扣的階段。下面我就按這個(gè)邏輯把每個(gè)環(huán)節(jié)的思考、操作和踩過(guò)的坑詳細(xì)道來(lái)。2. 數(shù)據(jù)戰(zhàn)場(chǎng)尋找、清洗與理解你的“彈藥”任何數(shù)據(jù)科學(xué)項(xiàng)目的起點(diǎn)和基石都是數(shù)據(jù)。對(duì)于“抗乳腺癌藥物”這個(gè)主題理想的數(shù)據(jù)應(yīng)該包含至少兩個(gè)維度患者/細(xì)胞系的分子特征如基因表達(dá)、突變和對(duì)應(yīng)的藥物反應(yīng)數(shù)據(jù)如IC50、敏感性標(biāo)簽。我們當(dāng)時(shí)主要使用了兩個(gè)公開(kāi)數(shù)據(jù)庫(kù)癌癥基因組圖譜TCGA和癌癥細(xì)胞系百科全書CCLE。TCGA提供了大量乳腺癌患者的臨床信息和多組學(xué)數(shù)據(jù)而CCLE則包含了眾多癌細(xì)胞系對(duì)多種藥物的敏感性數(shù)據(jù)。2.1 數(shù)據(jù)獲取與初步探索的實(shí)戰(zhàn)細(xì)節(jié)直接從官網(wǎng)下載原始數(shù)據(jù)只是第一步更關(guān)鍵的是理解數(shù)據(jù)的結(jié)構(gòu)和含義。例如TCGA的基因表達(dá)數(shù)據(jù)通常是RNA-Seq的FPKM或TPM值這是一個(gè)巨大的矩陣行是樣本列是基因。我們的第一個(gè)操作就是用Python的pandas和numpy加載數(shù)據(jù)并立即進(jìn)行探索性數(shù)據(jù)分析EDA。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns # 假設(shè)已下載并解壓了TCGA-BRCA的基因表達(dá)數(shù)據(jù)文件 # 這里用一個(gè)簡(jiǎn)化的示例說(shuō)明流程 expression_data pd.read_csv(tcga_brca_expression.csv, index_col0) # 行索引為樣本ID列索引為基因名 clinical_data pd.read_csv(tcga_brca_clinical.csv, index_col0) # 包含樣本的亞型、分期等信息 print(f表達(dá)數(shù)據(jù)形狀: {expression_data.shape}) # 例如 (1000, 20000) 意味著1000個(gè)樣本2萬(wàn)個(gè)基因 print(f臨床數(shù)據(jù)形狀: {clinical_data.shape}) # 查看數(shù)據(jù)前幾行和基本信息 print(expression_data.head()) print(clinical_data[PAM50_subtype].value_counts()) # 查看乳腺癌亞型分布第一個(gè)坑數(shù)據(jù)對(duì)齊與樣本匹配。TCGA和CCLE的數(shù)據(jù)樣本ID體系完全不同且一個(gè)基于患者組織一個(gè)基于實(shí)驗(yàn)室細(xì)胞系。直接合并建模是行不通的。我們的策略是以CCLE的細(xì)胞系數(shù)據(jù)作為訓(xùn)練集因?yàn)橛兴幬锓磻?yīng)標(biāo)簽利用CCLE中細(xì)胞系的基因表達(dá)數(shù)據(jù)構(gòu)建預(yù)測(cè)模型然后嘗試將訓(xùn)練好的模型應(yīng)用于TCGA的患者數(shù)據(jù)進(jìn)行“跨域”預(yù)測(cè)和生物標(biāo)志物發(fā)現(xiàn)。這要求我們?cè)谔卣骰蛏媳仨殞?duì)齊因此需要取兩個(gè)數(shù)據(jù)集的基因交集。第二個(gè)坑缺失值與異常值。生物數(shù)據(jù)缺失嚴(yán)重特別是臨床數(shù)據(jù)中的某些字段。對(duì)于基因表達(dá)數(shù)據(jù)我們采用了相對(duì)保守的策略刪除在超過(guò)50%樣本中表達(dá)量為0或缺失的基因。對(duì)于數(shù)值型缺失根據(jù)情況用中位數(shù)或KNN插補(bǔ)。異常值檢測(cè)使用了箱線圖和Z-score方法但對(duì)于基因表達(dá)這種通常符合對(duì)數(shù)正態(tài)分布的數(shù)據(jù)我們更關(guān)注極端高表達(dá)值是否具有生物學(xué)意義如致癌基因擴(kuò)增而非簡(jiǎn)單剔除。注意在生物醫(yī)學(xué)領(lǐng)域粗暴地刪除或填充數(shù)據(jù)可能會(huì)抹掉重要的生物學(xué)信號(hào)如某個(gè)基因在特定亞型中普遍不表達(dá)。每一次數(shù)據(jù)清洗決策最好都能結(jié)合一些基本的生物學(xué)知識(shí)進(jìn)行判斷。2.2 特征工程的生物醫(yī)學(xué)視角拿到清洗后的基因表達(dá)矩陣比如1000個(gè)細(xì)胞系 x 15000個(gè)基因直接扔進(jìn)模型是災(zāi)難性的維度災(zāi)難、過(guò)擬合。特征工程的目標(biāo)是降維和提煉信息。方差過(guò)濾這是最直接的一步刪除在所有樣本中表達(dá)量幾乎沒(méi)有變化的基因。這些基因攜帶的信息量極少。from sklearn.feature_selection import VarianceThreshold selector VarianceThreshold(threshold0.1) # 設(shè)定一個(gè)方差閾值 expression_data_reduced selector.fit_transform(expression_data) print(f方差過(guò)濾后特征數(shù): {expression_data_reduced.shape[1]})基于生物學(xué)知識(shí)的過(guò)濾我們查閱了文獻(xiàn)和KEGG、GO等數(shù)據(jù)庫(kù)重點(diǎn)關(guān)注與乳腺癌通路相關(guān)的基因集如PI3K-Akt信號(hào)通路、細(xì)胞周期、DNA損傷修復(fù)等相關(guān)基因。這相當(dāng)于引入“領(lǐng)域先驗(yàn)知識(shí)”大幅縮小特征范圍提高模型的可解釋性。統(tǒng)計(jì)方法篩選對(duì)于回歸問(wèn)題預(yù)測(cè)IC50值我們使用與藥物反應(yīng)相關(guān)性分析如皮爾遜相關(guān)系數(shù)對(duì)于分類問(wèn)題敏感/耐藥使用方差分析ANOVA或卡方檢驗(yàn)篩選出在不同反應(yīng)組間差異顯著的基因。高級(jí)降維經(jīng)過(guò)上述篩選特征維度可能仍在數(shù)百到數(shù)千。我們進(jìn)一步使用了主成分分析PCA和t-SNE。PCA主要用于后續(xù)的模型輸入降維而t-SNE用于可視化觀察樣本在低維空間是否能夠按藥物敏感性或乳腺癌亞型自然聚類。from sklearn.decomposition import PCA pca PCA(n_components50) # 保留前50個(gè)主成分 expression_pca pca.fit_transform(expression_data_scaled) # 注意必須先標(biāo)準(zhǔn)化 print(fPCA累計(jì)方差貢獻(xiàn)率: {np.sum(pca.explained_variance_ratio_)})核心心得特征工程不是純粹的數(shù)學(xué)游戲。在生物醫(yī)學(xué)項(xiàng)目中每一步篩選最好都能對(duì)應(yīng)一個(gè)潛在的生物學(xué)假設(shè)。例如我們最終保留的特征基因列表應(yīng)該能回答“為什么是這些基因”這個(gè)問(wèn)題。這為后續(xù)的結(jié)果解釋打下了堅(jiān)實(shí)基礎(chǔ)。3. 模型構(gòu)建選擇、訓(xùn)練與驗(yàn)證的博弈特征準(zhǔn)備好后就進(jìn)入了核心的建模環(huán)節(jié)。我們的預(yù)測(cè)任務(wù)主要有兩類回歸預(yù)測(cè)連續(xù)的IC50值和分類預(yù)測(cè)敏感/耐藥二分類或多分類如對(duì)不同藥物的響應(yīng)。3.1 模型選型與集成策略我們沒(méi)有押寶單一模型而是構(gòu)建了一個(gè)模型池進(jìn)行對(duì)比和集成這在實(shí)際研究和競(jìng)賽中都非常常見(jiàn)。基礎(chǔ)模型線性回歸 / 邏輯回歸作為基線模型。雖然簡(jiǎn)單但如果特征經(jīng)過(guò)良好篩選線性模型往往能有不錯(cuò)的表現(xiàn)且解釋性極強(qiáng)。我們可以直接得到特征的系數(shù)將其視為對(duì)藥物敏感性的“貢獻(xiàn)度”。支持向量機(jī)SVM對(duì)于高維數(shù)據(jù)特別是當(dāng)特征數(shù)大于樣本數(shù)時(shí)線性SVC或使用RBF核的SVM有時(shí)能表現(xiàn)出優(yōu)勢(shì)。我們使用GridSearchCV來(lái)優(yōu)化懲罰系數(shù)C和核函數(shù)參數(shù)gamma。隨機(jī)森林Random Forest這是我們的主力模型之一。它能處理高維數(shù)據(jù)對(duì)缺失值和異常值不敏感并能提供特征重要性排序。我們用它來(lái)做分類和回歸。梯度提升樹(shù)如XGBoost, LightGBM另一個(gè)主力模型。通常在表格數(shù)據(jù)上表現(xiàn)優(yōu)于隨機(jī)森林但需要更仔細(xì)的參數(shù)調(diào)優(yōu)。我們使用LightGBM因?yàn)樗?xùn)練速度快且對(duì)類別不平衡問(wèn)題有較好的處理能力。集成方法Stacking我們嘗試了將隨機(jī)森林、LightGBM和SVM的預(yù)測(cè)結(jié)果作為新特征輸入到一個(gè)元學(xué)習(xí)器如邏輯回歸或線性回歸中進(jìn)行二次訓(xùn)練。這種方法有時(shí)能提升1-2%的精度但增加了復(fù)雜度。Voting / Averaging對(duì)于分類問(wèn)題采用硬投票或軟投票對(duì)于回歸問(wèn)題直接對(duì)多個(gè)模型的預(yù)測(cè)結(jié)果取平均。這是一個(gè)簡(jiǎn)單有效的提升魯棒性的方法。關(guān)鍵一步交叉驗(yàn)證與數(shù)據(jù)劃分。生物醫(yī)學(xué)數(shù)據(jù)常有批次效應(yīng)或樣本依賴性問(wèn)題。我們嚴(yán)格使用分層K折交叉驗(yàn)證確保每一折中各類別如不同亞型、不同敏感性的比例與整體數(shù)據(jù)集一致。絕對(duì)避免在訓(xùn)練集中出現(xiàn)測(cè)試集的任何信息。from sklearn.model_selection import StratifiedKFold, cross_val_score from lightgbm import LGBMClassifier model LGBMClassifier(random_state42, n_jobs-1) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(model, X_train, y_train, cvcv, scoringroc_auc) print(f5折交叉驗(yàn)證AUC均值: {scores.mean():.3f} (/- {scores.std():.3f}))3.2 超參數(shù)調(diào)優(yōu)與過(guò)擬合對(duì)抗機(jī)器學(xué)習(xí)競(jìng)賽和實(shí)際應(yīng)用的一大區(qū)別在于對(duì)過(guò)擬合的警惕程度。在有限的數(shù)據(jù)集上一個(gè)在交叉驗(yàn)證中表現(xiàn)優(yōu)異的模型可能在獨(dú)立測(cè)試集或真實(shí)世界中泛化能力很差。調(diào)優(yōu)工具我們使用了GridSearchCV和RandomizedSearchCV。對(duì)于像LightGBM這種參數(shù)較多的模型RandomizedSearchCV隨機(jī)搜索效率更高。貝葉斯優(yōu)化如optuna庫(kù)也是很好的選擇但當(dāng)時(shí)由于時(shí)間關(guān)系沒(méi)有深入。早停法Early Stopping對(duì)于梯度提升樹(shù)這類迭代模型早停法是防止過(guò)擬合的利器。我們?cè)诿恳徽劢徊骝?yàn)證的訓(xùn)練集中再劃分一個(gè)驗(yàn)證集用于監(jiān)控驗(yàn)證集性能當(dāng)性能不再提升時(shí)停止訓(xùn)練。from sklearn.model_selection import train_test_split X_train_part, X_val, y_train_part, y_val train_test_split(X_train, y_train, test_size0.2, stratifyy_train, random_state42) lgb_train lgb.Dataset(X_train_part, y_train_part) lgb_eval lgb.Dataset(X_val, y_val, referencelgb_train) params {...} # 模型參數(shù) gbm lgb.train(params, lgb_train, num_boost_round1000, valid_sets[lgb_eval], callbacks[lgb.early_stopping(stopping_rounds50)]) # 早停正則化在模型參數(shù)中顯式地加入L1或L2正則化項(xiàng)如邏輯回歸的penaltyl1 LightGBM的lambda_l1,lambda_l2限制模型復(fù)雜度。簡(jiǎn)化模型有時(shí)候最好的防過(guò)擬合策略是使用更簡(jiǎn)單的模型如線性模型或減少特征數(shù)量。我們對(duì)比了使用全基因集、通路基因集和統(tǒng)計(jì)篩選基因集構(gòu)建的模型發(fā)現(xiàn)后者雖然特征少但測(cè)試集性能往往更穩(wěn)定。踩坑實(shí)錄我們?cè)欢仍谟?xùn)練集上達(dá)到了接近0.95的AUC欣喜若狂。但當(dāng)用從未參與過(guò)任何訓(xùn)練過(guò)程的完全獨(dú)立的測(cè)試集來(lái)自另一個(gè)數(shù)據(jù)源評(píng)估時(shí)AUC驟降到0.65左右。這就是典型的過(guò)擬合。復(fù)盤發(fā)現(xiàn)問(wèn)題出在特征工程環(huán)節(jié)我們?cè)诤Y選差異基因時(shí)使用了整個(gè)數(shù)據(jù)集包含未來(lái)的測(cè)試集的統(tǒng)計(jì)量導(dǎo)致信息泄露。正確的做法是在交叉驗(yàn)證的每一折中僅使用該折的訓(xùn)練集數(shù)據(jù)來(lái)進(jìn)行特征篩選然后用篩選出的特征來(lái)轉(zhuǎn)換該折的訓(xùn)練集和驗(yàn)證集。這是一個(gè)極其重要且容易忽略的細(xì)節(jié)。4. 從預(yù)測(cè)到洞見(jiàn)模型解釋與生物意義挖掘得到一個(gè)高精度的“黑箱”模型并不是終點(diǎn)尤其是在生物醫(yī)學(xué)領(lǐng)域。醫(yī)生和生物學(xué)家更關(guān)心的是模型依據(jù)什么做出判斷哪些基因或通路是關(guān)鍵這背后暗示了怎樣的生物學(xué)機(jī)制4.1 模型可解釋性技術(shù)應(yīng)用特征重要性樹(shù)模型隨機(jī)森林、LightGBM天然提供特征重要性如基尼重要性、分裂增益。我們將其排序列出Top N的基因。importances gbm.feature_importances_ indices np.argsort(importances)[::-1] top_genes [gene_names[i] for i in indices[:20]] print(Top 20重要基因:, top_genes)SHAP值分析這是當(dāng)前最流行的模型解釋工具之一。SHAP值可以量化每個(gè)特征對(duì)單個(gè)樣本預(yù)測(cè)結(jié)果的貢獻(xiàn)度既能看全局重要性也能看局部單個(gè)樣本解釋。import shap explainer shap.TreeExplainer(gbm) shap_values explainer.shap_values(X_test) # 全局摘要圖 shap.summary_plot(shap_values, X_test, feature_namesgene_names) # 對(duì)某個(gè)特定樣本例如一個(gè)對(duì)藥物敏感的患者的解釋 shap.force_plot(explainer.expected_value, shap_values[0,:], X_test.iloc[0,:], feature_namesgene_names)通過(guò)SHAP圖我們可以清晰地看到對(duì)于預(yù)測(cè)為“敏感”的樣本是哪些基因的高表達(dá)或低表達(dá)推動(dòng)了這一預(yù)測(cè)。部分依賴圖PDP與個(gè)體條件期望圖ICE用于可視化單個(gè)或兩個(gè)特征對(duì)模型預(yù)測(cè)結(jié)果的邊際效應(yīng)。例如我們可以觀察某個(gè)關(guān)鍵基因的表達(dá)量從低到高變化時(shí)模型預(yù)測(cè)的敏感性概率如何變化。4.2 生物信息學(xué)富集分析拿到一列重要的基因名單后下一步就是回答它們的生物學(xué)功能。我們使用在線工具如DAVID、Metascape或R/Python包如clusterProfiler進(jìn)行基因本體GO富集分析和京都基因與基因組百科全書KEGG通路富集分析。這個(gè)過(guò)程通常是將模型篩選出的Top 200個(gè)重要基因作為“輸入基因列表”以人類所有基因?yàn)楸尘坝?jì)算這些基因在哪些生物學(xué)過(guò)程、分子功能、細(xì)胞組分或信號(hào)通路上顯著富集。結(jié)果會(huì)以p值或錯(cuò)誤發(fā)現(xiàn)率FDR排序。我們的一次關(guān)鍵發(fā)現(xiàn)模型識(shí)別出的重要基因在“雌激素反應(yīng)”、“G2/M細(xì)胞周期檢查點(diǎn)”等通路上顯著富集。這與已知的乳腺癌生物學(xué)尤其是激素受體陽(yáng)性型乳腺癌高度吻合并且提示我們模型的預(yù)測(cè)可能部分依賴于細(xì)胞增殖相關(guān)的通路。這極大地增強(qiáng)了我們模型結(jié)果的可信度和生物學(xué)意義。5. 項(xiàng)目復(fù)盤超越競(jìng)賽的思考與實(shí)用建議回顧整個(gè)項(xiàng)目從一道競(jìng)賽題目出發(fā)我們實(shí)際上走完了一個(gè)小型生物信息學(xué)研究的完整流程。在這個(gè)過(guò)程中技術(shù)上的挑戰(zhàn)固然很多但更深的體會(huì)來(lái)自于對(duì)數(shù)據(jù)科學(xué)項(xiàng)目本質(zhì)的思考。首先問(wèn)題定義比模型選擇更重要。最初我們糾結(jié)于用哪個(gè)高級(jí)模型。后來(lái)明白清晰定義要解決的具體問(wèn)題是預(yù)測(cè)IC50值還是分類敏感/耐藥是預(yù)測(cè)單一藥物還是多藥聯(lián)合以及如何評(píng)估用什么指標(biāo)AUC RMSE 還是特異性/敏感性才是首要的。這直接決定了數(shù)據(jù)如何準(zhǔn)備、特征如何構(gòu)造、模型如何設(shè)計(jì)。其次數(shù)據(jù)質(zhì)量決定天花板特征工程決定接近天花板的速度。我們花了超過(guò)60%的時(shí)間在數(shù)據(jù)獲取、清洗、理解和特征工程上。一個(gè)干凈、信息量大的特征集即使用簡(jiǎn)單的線性模型也能得到不錯(cuò)的結(jié)果。反之再?gòu)?fù)雜的模型在垃圾數(shù)據(jù)上也只會(huì)產(chǎn)出垃圾結(jié)果。再者可解釋性是生物醫(yī)學(xué)AI的“生命線”。在醫(yī)療領(lǐng)域一個(gè)無(wú)法解釋的“黑箱”模型無(wú)論其交叉驗(yàn)證分?jǐn)?shù)多高都很難獲得臨床信任。SHAP、LIME等可解釋性技術(shù)結(jié)合傳統(tǒng)的生物信息學(xué)富集分析是連接機(jī)器學(xué)習(xí)預(yù)測(cè)與生物學(xué)理解的橋梁。在論文或報(bào)告中這部分內(nèi)容的深度往往比模型精度本身更能體現(xiàn)工作價(jià)值。最后工程化思維與協(xié)作。真實(shí)的項(xiàng)目涉及多環(huán)節(jié)、多工具。我們使用Git進(jìn)行版本控制用Jupyter Notebook進(jìn)行探索性分析用Python腳本封裝數(shù)據(jù)處理和模型訓(xùn)練流水線用Docker或Conda管理復(fù)現(xiàn)環(huán)境。良好的代碼結(jié)構(gòu)和文檔習(xí)慣不僅是為了比賽更是為了未來(lái)自己或他人能夠復(fù)現(xiàn)和拓展這項(xiàng)工作。對(duì)于想要參加類似數(shù)學(xué)建模競(jìng)賽或從事相關(guān)交叉學(xué)科研究的朋友我的建議是不要只盯著最新的神經(jīng)網(wǎng)絡(luò)模型。從扎實(shí)的數(shù)據(jù)處理、統(tǒng)計(jì)基礎(chǔ)、和經(jīng)典的機(jī)器學(xué)習(xí)模型如樹(shù)模型開(kāi)始理解每一個(gè)步驟背后的“為什么”并始終將你的工作與領(lǐng)域知識(shí)如乳腺癌的分子分型、常見(jiàn)治療靶點(diǎn)緊密結(jié)合。當(dāng)你能夠清晰地向一位生物學(xué)家解釋你的模型為什么認(rèn)為某個(gè)患者可能對(duì)某藥敏感并指出可能與哪些已知通路相關(guān)時(shí)你的項(xiàng)目就真正產(chǎn)生了價(jià)值。這次“華為杯”的經(jīng)歷與其說(shuō)是一次競(jìng)賽不如說(shuō)是一次將數(shù)學(xué)、計(jì)算機(jī)與生命科學(xué)知識(shí)融會(huì)貫通的寶貴實(shí)踐。