簡介本資源為2021年華為杯數(shù)學(xué)建模競(jìng)賽D題「抗乳腺癌候選藥物的優(yōu)化建?！沟耐暾獯鸢嫦騾⒓訑?shù)學(xué)建模競(jìng)賽的高校學(xué)生及對(duì)生物信息學(xué)、醫(yī)療數(shù)據(jù)分析感興趣的進(jìn)階學(xué)習(xí)者。內(nèi)容圍繞特征選擇、回歸預(yù)測(cè)、二分類建模、最優(yōu)化求解、模型訓(xùn)練與驗(yàn)證、數(shù)據(jù)預(yù)處理、結(jié)果可視化等核心環(huán)節(jié)展開可幫助讀者理解如何將機(jī)器學(xué)習(xí)方法落地到藥物研發(fā)場(chǎng)景。壓縮包共83個(gè)文件以32個(gè)Python腳本、22個(gè)CSV數(shù)據(jù)文件、14張PNG圖表為主另含XLSX、XML、Markdown與DOCX等輔助材料整體約31.82MB目錄按代碼、數(shù)據(jù)與文檔分層組織便于按模塊查閱。目前已有3477人學(xué)習(xí)下載適合需要復(fù)盤賽題思路、參考建模流程與代碼實(shí)現(xiàn)、查漏補(bǔ)缺的參賽者與研究者。1. 華為杯D題抗乳腺癌候選藥物優(yōu)化建模這道題到底在考什么2021年華為杯研究生數(shù)學(xué)建模競(jìng)賽D題題目全稱是「抗乳腺癌候選藥物的優(yōu)化建模」屬于典型的「數(shù)據(jù)驅(qū)動(dòng) 多目標(biāo)優(yōu)化」復(fù)合型賽題。它給出的是一批化合物分子描述符數(shù)據(jù)要求參賽隊(duì)圍繞生物活性、ADMET性質(zhì)吸收、分布、代謝、排泄、毒性建立預(yù)測(cè)模型再在此基礎(chǔ)上做分子結(jié)構(gòu)的優(yōu)化篩選。很多隊(duì)伍拿到題的第一反應(yīng)是直接上機(jī)器學(xué)習(xí)調(diào)包結(jié)果發(fā)現(xiàn)數(shù)據(jù)維度高、樣本量小、指標(biāo)之間還互相打架最后論文寫得像實(shí)驗(yàn)報(bào)告模型精度也上不去。這道題真正考的不是某個(gè)算法會(huì)不會(huì)用而是你能不能把「預(yù)測(cè)」和「優(yōu)化」串成一條閉環(huán)先用回歸模型把pIC50等活性指標(biāo)預(yù)測(cè)準(zhǔn)再用分類模型判斷化合物是否滿足ADMET約束最后在滿足約束的前提下搜索活性更高的分子描述符組合。適合已經(jīng)做過一兩次數(shù)學(xué)建模、想沖擊華為杯優(yōu)秀論文的研究生也適合想系統(tǒng)練一遍「特征工程 多目標(biāo)優(yōu)化」落地流程的工程師。下面我按自己復(fù)現(xiàn)這道題時(shí)的實(shí)際路徑把選型理由、參數(shù)設(shè)置和踩過的坑講清楚。2. 數(shù)據(jù)預(yù)處理與特征篩選從729個(gè)描述符里挑出能用的2.1 先搞清楚數(shù)據(jù)長什么樣再動(dòng)手這道題官方給的數(shù)據(jù)通常包含訓(xùn)練集和測(cè)試集兩部分訓(xùn)練集里每個(gè)化合物有分子描述符MD、pIC50值、以及若干ADMET標(biāo)簽。分子描述符數(shù)量在729個(gè)左右涵蓋拓?fù)渲笖?shù)、分子指紋、物化性質(zhì)等。樣本量通常只有幾百到一千出頭屬于典型的「寬而淺」數(shù)據(jù)。直接把這729列全丟進(jìn)模型結(jié)果一定是過擬合交叉驗(yàn)證分?jǐn)?shù)好看但測(cè)試集崩掉。我一般會(huì)先做三件事檢查缺失值分布、看每列方差、算一下描述符與pIC50的皮爾遜相關(guān)系數(shù)。方差接近0的列直接刪因?yàn)閷?duì)區(qū)分樣本沒貢獻(xiàn)相關(guān)系數(shù)絕對(duì)值低于某個(gè)閾值的列也先放一邊但不要急著永久刪除后面做特征重要性時(shí)還要回頭看。import pandas as pd import numpy as np from sklearn.feature_selection import VarianceThreshold # 讀取訓(xùn)練數(shù)據(jù)假設(shè)描述符列從第3列開始 train pd.read_csv(train.csv) desc_cols train.columns[2:] # 根據(jù)實(shí)際列名調(diào)整 # 1. 刪除方差極低的描述符 selector VarianceThreshold(threshold0.01) X_var selector.fit_transform(train[desc_cols]) kept_cols desc_cols[selector.get_support()] print(f方差篩選后剩余特征數(shù): {len(kept_cols)}) # 2. 計(jì)算與pIC50的相關(guān)系數(shù) corr train[kept_cols].corrwith(train[pIC50]).abs() corr_selected corr[corr 0.1].index.tolist() print(f相關(guān)性篩選后剩余特征數(shù): {len(corr_selected)})這段代碼的邏輯是先用方差閾值砍掉近似常數(shù)列再用相關(guān)系數(shù)做粗篩。threshold0.01這個(gè)值不是固定的如果數(shù)據(jù)做過標(biāo)準(zhǔn)化可以適當(dāng)調(diào)高到0.05如果沒標(biāo)準(zhǔn)化方差本身量綱差異大建議先做StandardScaler再篩。相關(guān)系數(shù)閾值0.1也是經(jīng)驗(yàn)值目的是保留弱相關(guān)但可能非線性有用的特征真正精細(xì)的篩選交給后面的模型。2.2 用隨機(jī)森林做嵌入法特征重要性排序粗篩之后特征數(shù)可能還有兩三百個(gè)這時(shí)候用嵌入法讓模型自己挑。隨機(jī)森林的特征重要性排序穩(wěn)定、對(duì)共線性不敏感適合這種高維小樣本場(chǎng)景。我一般會(huì)跑一次RandomForestRegressor取importance累計(jì)貢獻(xiàn)達(dá)到95%的前N個(gè)特征作為最終輸入。from sklearn.ensemble import RandomForestRegressor rf RandomForestRegressor(n_estimators500, max_depth8, random_state42, n_jobs-1) rf.fit(train[corr_selected], train[pIC50]) imp pd.Series(rf.feature_importances_, indexcorr_selected).sort_values(ascendingFalse) cum_imp imp.cumsum() / imp.sum() final_features cum_imp[cum_imp 0.95].index.tolist() print(f最終入模特征數(shù): {len(final_features)})n_estimators500是為了讓重要性估計(jì)更穩(wěn)定max_depth8是防止單棵樹過深導(dǎo)致重要性被噪聲主導(dǎo)。累計(jì)貢獻(xiàn)閾值0.95意味著保留能解釋95%重要性的特征剩下的5%大概率是噪聲。這一步做完特征數(shù)通常能壓到50以內(nèi)后面建模會(huì)輕很多。注意特征篩選必須在交叉驗(yàn)證的每一折內(nèi)部獨(dú)立做不能在全量數(shù)據(jù)上篩完再交叉驗(yàn)證否則會(huì)信息泄漏交叉驗(yàn)證分?jǐn)?shù)虛高。我見過不少論文在這里翻車。3. 活性預(yù)測(cè)模型pIC50回歸怎么選、怎么調(diào)3.1 為什么我首選梯度提升樹而不是深度學(xué)習(xí)這道題的樣本量決定了深度學(xué)習(xí)不是好選擇。幾百個(gè)樣本、幾十個(gè)特征神經(jīng)網(wǎng)絡(luò)參數(shù)量隨便就超過樣本量訓(xùn)練集loss降得再低也是背答案。梯度提升樹XGBoost、LightGBM、CatBoost在這個(gè)量級(jí)上表現(xiàn)穩(wěn)定對(duì)特征縮放不敏感還能直接輸出特征重要性做解釋。我一般先用LightGBM跑基線再用XGBoost對(duì)比最后用貝葉斯優(yōu)化調(diào)參。import lightgbm as lgb from sklearn.model_selection import KFold, cross_val_score model lgb.LGBMRegressor( n_estimators800, learning_rate0.03, num_leaves15, max_depth6, min_child_samples10, subsample0.8, colsample_bytree0.7, reg_alpha0.1, reg_lambda0.1, random_state42 ) kf KFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(model, train[final_features], train[pIC50], cvkf, scoringneg_mean_squared_error) rmse np.sqrt(-scores.mean()) print(f5折交叉驗(yàn)證RMSE: {rmse:.4f})num_leaves15和max_depth6是控制模型復(fù)雜度的關(guān)鍵小樣本下葉子數(shù)太多必然過擬合。learning_rate0.03配合n_estimators800是慢學(xué)習(xí)率多棵樹的經(jīng)典組合比0.1配200棵更穩(wěn)。min_child_samples10保證每個(gè)葉子至少有10個(gè)樣本防止模型記住個(gè)別離群點(diǎn)。這些參數(shù)不是拍腦袋是我在類似規(guī)模數(shù)據(jù)集上反復(fù)試出來的區(qū)間你可以在這個(gè)基礎(chǔ)上用Optuna做精細(xì)搜索。3.2 評(píng)估指標(biāo)不能只看RMSE回歸任務(wù)里RMSE是基礎(chǔ)但這道題最終要服務(wù)于「篩選高活性化合物」所以排序能力比絕對(duì)誤差更重要。我建議同時(shí)看Spearman相關(guān)系數(shù)和Top-K命中率把預(yù)測(cè)pIC50最高的前10%化合物挑出來看真實(shí)pIC50也排在前10%的比例有多少。這個(gè)指標(biāo)直接對(duì)應(yīng)后續(xù)優(yōu)化環(huán)節(jié)的實(shí)用性。from scipy.stats import spearmanr # 用交叉驗(yàn)證的預(yù)測(cè)結(jié)果計(jì)算排序指標(biāo) from sklearn.model_selection import cross_val_predict pred cross_val_predict(model, train[final_features], train[pIC50], cvkf) spearman spearmanr(pred, train[pIC50]).correlation print(fSpearman相關(guān)系數(shù): {spearman:.4f}) # Top-K命中率 k int(len(pred) * 0.1) top_pred_idx np.argsort(pred)[-k:] top_true_idx np.argsort(train[pIC50].values)[-k:] hit_rate len(set(top_pred_idx) set(top_true_idx)) / k print(fTop-10%命中率: {hit_rate:.4f})Spearman能到0.7以上、Top-10%命中率能到0.5以上這個(gè)模型就算可用了。如果Spearman低于0.6優(yōu)先檢查特征篩選是不是漏掉了關(guān)鍵描述符而不是急著換模型。4. ADMET分類與多目標(biāo)約束把「能用」和「好用」分開4.1 ADMET標(biāo)簽不平衡怎么處理ADMET性質(zhì)通常是二分類標(biāo)簽比如是否具有肝毒性、是否高血漿蛋白結(jié)合率。這類數(shù)據(jù)往往正負(fù)樣本比例懸殊直接訓(xùn)練分類器會(huì)偏向多數(shù)類。我的做法是用SMOTE做少數(shù)類過采樣但只在訓(xùn)練折內(nèi)部做驗(yàn)證折保持原始分布。分類器選LightGBM的class_weightbalanced也能緩解但SMOTE在小樣本下更直接。from imblearn.over_sampling import SMOTE from sklearn.model_selection import StratifiedKFold from sklearn.metrics import f1_score skf StratifiedKFold(n_splits5, shuffleTrue, random_state42) f1_scores [] for train_idx, val_idx in skf.split(train[final_features], train[ADMET_label]): X_tr, X_val train[final_features].iloc[train_idx], train[final_features].iloc[val_idx] y_tr, y_val train[ADMET_label].iloc[train_idx], train[ADMET_label].iloc[val_idx] smote SMOTE(random_state42, k_neighbors3) X_res, y_res smote.fit_resample(X_tr, y_tr) clf lgb.LGBMClassifier(n_estimators300, learning_rate0.05, num_leaves15, random_state42) clf.fit(X_res, y_res) pred_val clf.predict(X_val) f1_scores.append(f1_score(y_val, pred_val)) print(fADMET分類平均F1: {np.mean(f1_scores):.4f})k_neighbors3是因?yàn)樾颖鞠锣従犹鄷?huì)生成不真實(shí)的合成樣本。F1比準(zhǔn)確率更適合評(píng)估不平衡分類如果F1低于0.6說明ADMET標(biāo)簽本身噪聲大或者特征區(qū)分度不夠這時(shí)候不要硬調(diào)模型回去檢查標(biāo)簽定義和特征覆蓋。4.2 多目標(biāo)優(yōu)化活性最大化和ADMET約束怎么同時(shí)滿足這道題的核心輸出是一組推薦化合物要求活性高且ADMET性質(zhì)達(dá)標(biāo)。我把它建模成帶約束的單目標(biāo)優(yōu)化目標(biāo)函數(shù)是預(yù)測(cè)pIC50最大化約束是各ADMET分類器預(yù)測(cè)概率超過閾值。搜索空間是分子描述符的可行域但描述符之間不獨(dú)立不能簡單隨機(jī)采樣。常見做法是在現(xiàn)有化合物庫中做篩選排序而不是從頭生成分子。# 對(duì)測(cè)試集化合物打分活性預(yù)測(cè) ADMET約束過濾 test_pred_activity model.predict(test[final_features]) admet_probs clf.predict_proba(test[final_features])[:, 1] # 約束ADMET達(dá)標(biāo)概率 0.7 mask admet_probs 0.7 candidates test[mask].copy() candidates[pred_pIC50] test_pred_activity[mask] candidates candidates.sort_values(pred_pIC50, ascendingFalse) print(f滿足ADMET約束的候選化合物數(shù): {len(candidates)}) print(candidates[[compound_id, pred_pIC50]].head(10))閾值0.7是平衡召回和精度的經(jīng)驗(yàn)值如果候選太少可以降到0.6候選太多就提到0.8。這一步的關(guān)鍵是先把約束卡死再在可行域內(nèi)排序而不是把兩個(gè)目標(biāo)加權(quán)求和——加權(quán)系數(shù)很難解釋評(píng)審也不買賬。5. 避坑與排查這道題最容易翻車的五個(gè)地方5.1 特征篩選在交叉驗(yàn)證外面做導(dǎo)致分?jǐn)?shù)虛高現(xiàn)象交叉驗(yàn)證RMSE只有0.3測(cè)試集RMSE飆到0.8以上。原因在全量數(shù)據(jù)上做特征篩選驗(yàn)證折的信息泄漏到了訓(xùn)練過程。解決把方差篩選、相關(guān)性篩選、隨機(jī)森林重要性全部封裝進(jìn)Pipeline在每折訓(xùn)練集上fit驗(yàn)證集上transform。5.2 把ADMET標(biāo)簽當(dāng)回歸做現(xiàn)象用回歸模型預(yù)測(cè)ADMET數(shù)值結(jié)果全是0.5左右的無效輸出。原因ADMET標(biāo)簽是二分類回歸模型優(yōu)化MSE會(huì)傾向于預(yù)測(cè)均值。解決改用分類模型評(píng)估指標(biāo)換成F1或AUC不要看RMSE。5.3 優(yōu)化環(huán)節(jié)直接對(duì)描述符做梯度上升現(xiàn)象生成的描述符組合在化學(xué)上不可能存在評(píng)審一眼看出是編的。原因描述符之間有物理約束獨(dú)立擾動(dòng)會(huì)破壞分子合法性。解決在現(xiàn)有化合物庫內(nèi)做篩選排序或者用遺傳算法在合法分子空間搜索不要對(duì)描述符向量直接做梯度優(yōu)化。5.4 忽略描述符的量綱差異現(xiàn)象基于距離的模型KNN、SVM效果遠(yuǎn)差于樹模型。原因描述符量綱從0.001到10000都有距離計(jì)算被大量綱特征主導(dǎo)。解決用樹模型可以跳過標(biāo)準(zhǔn)化但如果要用SVM或KNN必須先做StandardScaler或MinMaxScaler。5.5 論文里只寫最終模型不寫選型過程現(xiàn)象評(píng)審質(zhì)疑為什么用LightGBM不用隨機(jī)森林答不上來。原因只跑了最終模型沒做基線對(duì)比。解決至少跑三個(gè)基線線性回歸、隨機(jī)森林、XGBoost用表格列出交叉驗(yàn)證RMSE和Spearman選型理由自然就有了。6. 從復(fù)現(xiàn)到拿獎(jiǎng)一個(gè)被低估的提分技巧這道題拿高分的關(guān)鍵往往不在模型本身而在「可解釋性」和「閉環(huán)驗(yàn)證」這兩塊。我復(fù)盤過幾篇華為杯優(yōu)秀論文發(fā)現(xiàn)它們的共同點(diǎn)是不僅給出了預(yù)測(cè)模型還把特征重要性映射回化學(xué)意義比如哪些拓?fù)渲笖?shù)對(duì)應(yīng)分子柔性、哪些指紋對(duì)應(yīng)疏水性然后解釋為什么這些性質(zhì)影響抗乳腺癌活性。評(píng)審看的是你有沒有把數(shù)學(xué)建模和領(lǐng)域知識(shí)接上。具體操作上我習(xí)慣在最終論文里加一張「關(guān)鍵描述符-生物活性」對(duì)照表用SHAP值量化每個(gè)特征對(duì)pIC50的貢獻(xiàn)方向。SHAP比特征重要性多了一層方向信息能直接說「這個(gè)描述符增大時(shí)活性上升還是下降」。import shap explainer shap.TreeExplainer(model) shap_values explainer.shap_values(train[final_features]) # 輸出每個(gè)特征的平均絕對(duì)SHAP值按降序排列 shap_importance pd.Series( np.abs(shap_values).mean(axis0), indexfinal_features ).sort_values(ascendingFalse) print(shap_importance.head(15))這張表出來之后挑前5個(gè)特征去查化學(xué)數(shù)據(jù)庫確認(rèn)它們對(duì)應(yīng)的分子性質(zhì)寫進(jìn)論文的「結(jié)果分析」部分。這一步花不了多少時(shí)間但能讓論文從「調(diào)包報(bào)告」變成「有洞察的建模工作」。另一個(gè)提分點(diǎn)是做敏感性分析把ADMET約束閾值從0.6到0.9逐檔變化看推薦化合物數(shù)量和平均預(yù)測(cè)活性的變化曲線。這條曲線能說明你的優(yōu)化方案在不同嚴(yán)格程度下都穩(wěn)定而不是卡在一個(gè)特定閾值上才有效。我一般會(huì)跑5檔畫一張雙軸圖左邊是候選數(shù)量、右邊是平均pIC50評(píng)審看到這種分析基本會(huì)給加分。最后說個(gè)血淚經(jīng)驗(yàn)這道題的數(shù)據(jù)預(yù)處理和特征篩選至少占整個(gè)工作量的60%建模和優(yōu)化各占20%。很多隊(duì)伍反過來80%時(shí)間調(diào)模型結(jié)果特征沒選好怎么調(diào)都上不去。我現(xiàn)在的習(xí)慣是先把特征篩選的Pipeline搭穩(wěn)交叉驗(yàn)證分?jǐn)?shù)穩(wěn)定了再動(dòng)模型參數(shù)。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取