☰
混合STR图谱解析:人数判别、比例估计与基因型推断可复现
2026/10/27 20:11:29 网站建设 项目流程

简介:面向2025年深圳杯数学建模竞赛D题参赛者与数学建模研究者的完整分析论文,围绕赛题重述、问题分析、模型假设与符号定义展开,并给出问题一至问题三的建模与求解思路,可用于赛前研读、赛后复盘及同类分类与组合优化问题的方案参考。压缩包共1个PDF文件,约2.71MB,正文按章节组织,包含特征工程设计、分类模型结构与数学表达、模型总体性能对比、不同贡献者数量下的表现分析,以及贡献者兼容性筛选、混合比例优化估计和组合匹配评分函数等内容,并附可运行代码与相关数据,便于对照理解模型实现与评估指标。目录层级清晰,涵盖数据分布概览、模型性能比较和重要特征分析,适合需要完整赛题论文框架与可复现材料的读者查阅。已有330人学习。

1. 从一份赛题论文说起:四个问题、一套可运行代码、一份数据

第一次打开这份资源是在一个赶进度的晚上,队友把一份 PDF 甩过来,说模型都写完了,代码跑不起来。这是某地区数学建模竞赛 D 题的完整分析论文,主题是混合 STR 图谱的解析——给定一堆峰高、等位基因、片段长度,判断样本里混了几个人、各自占多少比例、每个人的基因型长什么样,最后再做一轮去噪提精度。它的价值不在于论文写得多漂亮,而在于四问配了可运行代码和配套数据,是一条能从头跑到尾的分析链路。适合两类人:正在啃这类赛题、需要一份能照着复现的参考答案的建模选手;以及手上有一堆多源信号、想做分解和归属判断的开发者。

2. 项目骨架与技术选型:四问如何拆成四个可运行的模块

2.1 论文目录与代码目录的映射关系

论文的章节顺序和代码的执行顺序是错位的,这点一定要先理清楚,否则你会在「问题三的求解」里找「问题一的特征工程」找半天。论文是按「问题重述 → 问题分析 → 模型假设 → 符号定义 → 问题一到四 → 模型推广」组织的,而代码是按「数据 → 特征 → 分类 → 推断 → 去噪 → 评估」组织的。我把它对齐成下面这张表:

论文章节代码模块输入输出
三、模型假设 / 四、符号定义config.py无全局常量、峰高阈值、基因座列表
五、问题一(贡献者人数判别)features.py+train_num.py峰表peaks.csv每样本特征向量、人数分类模型
六、问题二(DNA 比例估计)train_ratio.py带比例标签的样本按人数分桶的比例分类器
七、问题三(基因型推断)infer_genotype.py混合样本 + 参考库候选组合排序、推断基因型
八、问题四(去噪)denoise.py原始峰高信号去噪后峰表、筛峰结果
九、模型推广无独立代码—论文里的文字讨论

这份资源的目录结构大致是这样,我按自己复现时的习惯整理过,模块名对齐上表:

model_d_paper/ ├── data/ │ ├── peaks.csv # 峰表:SampleID, Marker, Allele, Height, Size, IsOL │ ├── samples.csv # 样本表:SampleID, ContributorNum, RatioLabel │ └── reference.csv # 候选个体参考基因型 ├── config.py # 阈值、基因座、路径常量 ├── features.py # 基因座级 + 样本级特征 ├── train_num.py # 问题一:人数分类 ├── train_ratio.py # 问题二:比例分类 ├── infer_genotype.py # 问题三:兼容性评分 + 组合搜索 ├── denoise.py # 问题四:小波去噪 + 伪峰筛除 ├── evaluate.py # 混淆矩阵、准确率/召回率 └── requirements.txt

提示:如果你拿到的资源里文件名和这里不一致,先按「哪一章的求解结果需要哪个函数的输出」倒着找一遍调用关系,比按文件名猜快得多。

2.2 环境与依赖装起来

技术栈很轻,纯 Python 数据科学栈,没有 GPU、没有深度学习框架的必要。论文里提到的四类模型(随机森林、GBDT、MLP、SVM)全在 scikit-learn 里,小波去噪用 PyWavelets,加权最小二乘的求解用 scipy 的nnls就够。

# 建议单独建虚拟环境,避免和你机器上的科学计算包打架 python -m venv venv source venv/bin/activate # Windows 用 venv\Scripts\activate pip install numpy pandas scikit-learn scipy PyWavelets \ matplotlib seaborn openpyxl joblib

版本上我一般会锁一下,因为HistGradientBoostingClassifier的参数在旧版本里名字不一样,pywt.threshold的mode参数在某些小版本里对空数组处理有差异。稳妥写法是把下面这行也放进requirements.txt:

numpy>=1.24 pandas>=2.0 scikit-learn>=1.3 scipy>=1.10 PyWavelets>=1.4 joblib>=1.3

2.3 数据读取与字段规范化

这一步是整个链路的地基,字段类型错了,后面所有特征都是错的。峰表里最容易翻车的是Height被读成字符串或者带千分位逗号,IsOL被读成True/False之外的字符串。我习惯在读取之后立刻做一次类型收敛:

import pandas as pd import numpy as np def load_peaks(path="data/peaks.csv") -> pd.DataFrame: df = pd.read_csv(path) # 峰高必须是数值,读成字符串会在后面的 mean/std 上静默出错 df["Height"] = pd.to_numeric(df["Height"], errors="coerce") df["Size"] = pd.to_numeric(df["Size"], errors="coerce") # OL 是"无法与标准梯度比对"的标记,统一成布尔 df["IsOL"] = df["IsOL"].astype(str).str.lower().isin(["1", "true", "yes"]) # 丢掉峰高缺失的行,保留其余 df = df.dropna(subset=["Height"]).reset_index(drop=True) # 关键:等位基因标签统一成字符串,避免 12 和 "12" 被当成两个不同等位基因 df["Allele"] = df["Allele"].astype(str) return df

这里的逻辑说明几句。errors="coerce"把非法值变成 NaN 再统一丢弃,比让 pandas 抛异常更适合脏数据场景;Allele强制转成字符串是因为后续做集合交集时,12和"12"在 Python 里不相等,会出现「明明峰在那儿却匹配不上」的玄学问题。IsOL的三态判断覆盖了常见的导出格式,多写这一行能省掉后面一堆调试。

3. 特征工程与人数判别:GBDT 在问题一里凭什么跑到 0.96

3.1 从峰表到样本级特征矩阵

问题一的目标是判断混合样本里有几个贡献者,是个标准的多类分类。论文的特征设计分四组:等位基因结构特征、峰高统计特征、片段长度分布特征、逻辑与结构标记。落到代码上,要分两层做——先算基因座级特征,再聚合成样本级特征。

import numpy as np import pandas as pd def locus_features(peaks: pd.DataFrame, height_thr: float = 50.0) -> dict: """单个样本、单个基因座下的峰 -> 一组特征""" h = peaks["Height"].to_numpy(dtype=float) h = h[h >= height_thr] # 低于阈值的当噪声丢掉,阈值可配 if h.size == 0: return {} hs = np.sort(h)[::-1] # 降序,方便取 top-k mean, std = h.mean(), h.std(ddof=0) top1 = hs[0] top2 = hs[1] if hs.size > 1 else 0.0 top4 = hs[:4] return { "allele_count": h.size, # 有效等位基因数 "non_ol_count": int((~peaks.loc[peaks["Height"] >= height_thr, "IsOL"]).sum()), "height_mean": mean, "height_std": std, "height_cv": std / (mean + 1e-9), # 变异系数,衡量峰高离散度 "height_min": h.min(), "height_max": h.max(), "top1_top2_ratio": top1 / (top2 + 1e-9), # 主次峰比,判断是否有主导贡献者 "top4_ratio": top4.sum() / (h.sum() + 1e-9), # 前四高峰占全部峰高的比例 "size_range": peaks["Size"].max() - peaks["Size"].min(), "size_std": peaks["Size"].std(ddof=0), }

这段函数的参数就两个:peaks是某个样本某个基因座的全部峰,height_thr是峰高阈值。为什么要有阈值?原始图谱里低于 50 RFU 的峰大概率是基线噪声,不滤掉会让allele_count虚高,直接把「2 人」判成「5 人」。这个阈值在不同批次的试剂盒下会浮动,所以一定要参数化,别硬编码。

样本级聚合就是把 16 个基因座的特征按样本汇总。求和类的特征(等位基因数)取 sum 和 max,比值类的取 mean,这样既有总量信息也有单点极值信息:

def sample_features(peaks_df: pd.DataFrame) -> pd.DataFrame: rows = [] for (sid, marker), g in peaks_df.groupby(["SampleID", "Marker"]): f = locus_features(g) if f: f.update(SampleID=sid, Marker=marker) rows.append(f) lf = pd.DataFrame(rows) agg = lf.groupby("SampleID").agg( allele_count_sum=("allele_count", "sum"), allele_count_max=("allele_count", "max"), non_ol_sum=("non_ol_count", "sum"), height_cv_mean=("height_cv", "mean"), top1_top2_mean=("top1_top2_ratio", "mean"), top4_ratio_mean=("top4_ratio", "mean"), height_max_max=("height_max", "max"), size_range_max=("size_range", "max"), ).reset_index() return agg

allele_count_sum和allele_count_max是两个不同维度的信息:前者反映整张图谱有多"拥挤",后者反映最拥挤的那个基因座有多拥挤。论文里expected_allele_ratio(等位基因数与理论人数比值)拿到了最高特征重要性,本质就是这两类信息被标签信息放大后的结果,这一点在第 3.3 节展开。

3.2 四个分类器的对照与训练参数

论文对比了随机森林、GBDT、MLP、SVM 四类模型,最终 GBDT 在测试集上到 0.9608,交叉验证 0.9584±0.0136,稳得不像话。对照组里 MLP 训练集 0.75、交叉验证只有 0.36,是典型的过拟合;SVM 全线低于 0.6,说明这类特征空间下线性/近似线性的边界不够用。

复现时我建议直接用HistGradientBoostingClassifier,因为它在几万行特征表上比传统GradientBoostingClassifier快一个数量级,效果基本一致:

from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier from sklearn.model_selection import StratifiedGroupKFold, cross_val_score from sklearn.metrics import classification_report, confusion_matrix from sklearn.preprocessing import StandardScaler from sklearn.pipeline import make_pipeline from sklearn.neural_network import MLPClassifier from sklearn.svm import SVC def build_models(): return { "gbdt": HistGradientBoostingClassifier( max_iter=300, # 迭代轮数,论文里用的是 300 量级 learning_rate=0.08, # 学习率压小一点,配合较多轮数更稳 max_leaf_nodes=31, # 叶子数,控制单棵树复杂度 l2_regularization=1.0, # 正则项,防止 5 人样本过拟合 class_weight="balanced", # 2 人样本最多,必须加权 random_state=42, ), "rf": RandomForestClassifier( n_estimators=500, max_depth=None, class_weight="balanced_subsample", n_jobs=-1, random_state=42, ), # MLP 和 SVM 必须接标准化,否则峰高量纲会压死小特征 "mlp": make_pipeline(StandardScaler(), MLPClassifier(hidden_layer_sizes=(128, 64), alpha=1e-3, max_iter=800, random_state=42)), "svm": make_pipeline(StandardScaler(), SVC(kernel="rbf", C=1.0, gamma="scale", class_weight="balanced", random_state=42)), }

几个参数值得单独说。class_weight="balanced"是必须的,论文里 2 人样本 240 个、3 人 256 个、4 人 176 个、5 人 144 个,看着不算极端,但加上比例类别之后不平衡会放大;不加权的话模型会倾向于把所有样本判成 2 人或 3 人,准确率看着还凑合,但 5 人样本的召回率会掉到 0.6 以下。l2_regularization是 GBDT 在 5 人样本上不崩的关键,5 人样本等位基因重叠严重,特征之间高度共线,正则弱一点就会背训练集。

评估部分,交叉验证一定要按SampleID分组,不能用普通的 KFold:

X = feat.drop(columns=["SampleID", "label"]) y = feat["label"] groups = feat["SampleID"] cv = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=42) for name, model in build_models().items(): scores = cross_val_score(model, X, y, cv=cv, groups=groups, scoring="f1_macro") print(f"{name}: {scores.mean():.4f} ± {scores.std():.4f}")

StratifiedGroupKFold同时保证类别比例和样本不跨折,这是避免「同一张图谱切出来的多个基因座记录同时出现在训练和验证集」的唯一办法。用普通 KFold,你会得到一个漂亮到不真实的 0.99,然后在真实新样本上翻车。

3.3 特征重要性与一个必须警惕的泄漏点

论文给的特征重要性图里,expected_allele_ratio权重超过 0.5,一骑绝尘。这个特征的定义是「等位基因数与理论人数比值」,理论上限是每人每个基因座 2 个等位基因。问题在于:如果分母用的是真实贡献者人数,那这就是拿标签造特征,属于典型的标签泄漏,模型在测试集上会虚高,上线就废。

我在复现时把它改成枚举式写法,彻底避免泄漏:

def add_expected_ratio(feat: pd.DataFrame) -> pd.DataFrame: """不依赖真实人数,而是对每个候选人数 k 各算一列比值""" for k in (2, 3, 4, 5): feat[f"expected_allele_ratio_k{k}"] = ( feat["allele_count_sum"] / (2 * k * 16) # 16 个基因座 ) return feat

这样模型看到的是「假设 k 人时,观测到的等位基因数占理论值的比例」,四个候选人数各一列,让模型自己去选哪个假设最自洽。逻辑上是把标签信息从特征里剥掉,交给模型在训练中学习映射关系,泛化能力会明显好于直接用真实人数归一化。这一改动在原始数据上准确率可能掉一两个点,但在新样本上的稳定性提升远大于这点损失。

4. 比例估计、基因型推断与去噪:问题二到问题四的三段链路

4.1 按人数分桶的比例分类器

问题二在已知人数前提下估计各贡献者的 DNA 比例。论文的思路是把连续比例离散成类别(2 人对应 1:1、1:4、1:9 这类组合),然后针对每种人数单独训一个分类器。这是个很务实的选择——比例回归的误差在小比例上会被放大,分类反而更稳。论文给出的选型是:2 人用 GBDT,3 人用 MLP,4/5 人用随机森林。性能上 2 人 85%,3 人只有 25%,4 人 40%,5 人 45%,这个衰减曲线基本反映了等位基因重叠的真实难度。

代码上是分桶训练,每个桶一个模型文件:

import joblib from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier from sklearn.neural_network import MLPClassifier from sklearn.preprocessing import StandardScaler from sklearn.pipeline import make_pipeline def train_ratio_models(feat: pd.DataFrame, num_col="ContributorNum", label_col="RatioLabel"): """按贡献者人数分桶,每桶独立训练一个比例分类器""" builders = { 2: lambda: HistGradientBoostingClassifier(max_iter=200, learning_rate=0.1, class_weight="balanced", random_state=42), 3: lambda: make_pipeline(StandardScaler(), MLPClassifier(hidden_layer_sizes=(64, 32), alpha=1e-2, max_iter=600, random_state=42)), 4: lambda: RandomForestClassifier(n_estimators=400, class_weight="balanced_subsample", n_jobs=-1, random_state=42), 5: lambda: RandomForestClassifier(n_estimators=400, class_weight="balanced_subsample", n_jobs=-1, random_state=42), } models = {} for k, build in builders.items(): sub = feat[feat[num_col] == k] if sub.empty: continue model = build() model.fit(sub.drop(columns=["SampleID", num_col, label_col]), sub[label_col]) joblib.dump(model, f"model_ratio_k{k}.pkl") models[k] = model return models

参数选择是有讲究的。3 人用 MLP 而不是树模型,是因为 3 人样本的比例结构比 2 人复杂,特征之间的交互更多是连续型的,MLP 在标准化后能拟合这种平滑边界;4/5 人退回随机森林,是因为样本量少、噪声大,树模型的抗噪和鲁棒性此时比拟合能力更重要。alpha=1e-2是 MLP 的 L2 正则,3 人桶样本只有 240 条,正则松了直接过拟合。

4.2 兼容性评分与候选组合搜索

问题三是整份论文里最有工程味道的部分:给定参考库里一批候选个体的基因型,推断混合样本中各个贡献者到底是谁。论文的策略是两步——先算兼容性做粗筛,再对留下的组合做精细匹配。

兼容性评分的定义很直白:候选个体的基因型等位基因与混合样本观测到的等位基因的重合比例。完全重合得满分,部分重合扣分。

from itertools import combinations def compat_score(ref_geno: dict, observed: dict) -> float: """ ref_geno: {marker: set(alleles)},每个基因座 1~2 个等位基因 observed: {marker: set(alleles)},混合样本在该基因座观测到的全部等位基因 """ total = hit = 0 for marker, alleles in ref_geno.items(): obs = observed.get(marker, set()) total += len(alleles) # 分母是该个体的基因座总数 hit += len(alleles & obs) # 完全重合的等位基因数 return hit / total if total else 0.0

这个函数返回 0~1 的分数,1 表示候选个体的所有等位基因都能在混合样本里找到。粗筛时按分数排序取前 N 个候选,再对候选做组合枚举:

def search_combinations(candidates, k, top_n=20): """从候选中选 k 人组合,先按分数粗排再枚举""" scored = sorted(candidates, key=lambda c: c["score"], reverse=True)[:top_n] return [tuple(c["id"] for c in combo) for combo in combinations(scored, k)]

top_n=20是个经验值。取太小会漏掉真正的贡献者,取太大组合数会爆炸——20 选 3 是 1140 组,20 选 5 是 15504 组,再往上算力扛不住。单基因座层面还能再剪枝:某个候选在超过一半的基因座上兼容性为 0,直接踢掉,因为真实贡献者不可能出现这种情况。

4.3 加权最小二乘估计混合比例并打分

对每个候选组合,下一步是用加权最小二乘拟合混合比例,让预测峰高尽量贴近观测峰高。这里构造一个设计矩阵 A,行是观测到的等位基因,列是组合里的每个贡献者,元素是该贡献者在这个等位基因上的拷贝数(0、1、2):

import numpy as np from scipy.optimize import nnls def fit_proportions(alleles, combo_genos, h_obs): """ alleles: 观测到的等位基因列表(顺序固定) combo_genos: [set(alleles), ...],组合里每个贡献者的基因型 h_obs: 与 alleles 对应的观测峰高 """ A = np.zeros((len(alleles), len(combo_genos))) for j, geno in enumerate(combo_genos): for i, a in enumerate(alleles): A[i, j] = 1.0 if a in geno else 0.0 # 非负最小二乘,保证比例不会出现负数 w, residual = nnls(A, h_obs) return w / (w.sum() + 1e-9), residual

nnls相比普通lstsq的关键优势是约束解非负。用普通最小二乘,在等位基因重叠严重的场景下会出现负比例,物理上无法解释,后续归一化还会放大误差。归一化那一步除以w.sum(),把绝对系数转成占比,这就是最终的混合比例估计。

有了比例,综合匹配评分就是两件事相加:解释率(观测峰高被模型解释掉的比例)和峰高分布相关性:

def match_score(w, A, h_obs, alpha=0.5): pred = A @ w pred_n = pred / (pred.sum() + 1e-9) obs_n = h_obs / (h_obs.sum() + 1e-9) corr = np.corrcoef(pred_n, obs_n)[0, 1] if h_obs.size > 2 else 0.0 explain = obs_n[pred_n > 1e-6].sum() # 被预测覆盖的观测峰占比 return alpha * explain + (1 - alpha) * max(corr, 0.0)

alpha=0.5是论文里的默认权重,实际调试时可以调。解释率偏高的组合更可能是真解,但如果两个组合解释率接近,相关性高的那个更可信,因为它把峰高的相对大小也对上了,而不只是"能解释到"。这两项分开看,比只用一项要稳得多。

4.4 小波阈值去噪与伪峰筛除

问题四的定位是给前三问兜底:原始图谱有噪声、有 PCR 扩增产生的拖尾峰,不去掉会污染特征和匹配结果。论文给的是多阶段流程——小波阈值去噪、拖尾峰识别、自适应峰筛选,最后再进匹配打分。效果上,多数样本准确率能到 0.8~1.0。

小波这一步的核心是把同一基因座按片段长度排序的峰高序列当成一维信号,分解后对高频细节系数做软阈值处理:

import pywt import numpy as np def wavelet_denoise(sig: np.ndarray, wavelet="db4", level=2) -> np.ndarray: n = len(sig) if n < 8: # 点数太少,小波没有意义,原样返回 return sig lvl = min(level, pywt.dwt_max_level(n, wavelet)) if lvl < 1: return sig coeffs = pywt.wavedec(sig, wavelet, level=lvl) # 用最高频细节系数的中位绝对偏差估计噪声标准差,对异常值稳健 sigma = np.median(np.abs(coeffs[-1])) / 0.6745 thr = sigma * np.sqrt(2 * np.log(max(n, 2))) # 通用阈值 # 近似层(低频)保留,细节层(高频)软阈值收缩 coeffs[1:] = [pywt.threshold(c, thr, mode="soft") for c in coeffs[1:]] rec = pywt.waverec(coeffs, wavelet)[:n] return np.clip(rec, 0, None) # 峰高不可能为负

参数说明:wavelet="db4"是峰型信号里常用的紧支撑小波,分解层数level=2对应 STR 峰高序列的尺度,层数再高容易把真实峰也削掉。阈值公式sigma * sqrt(2*ln(n))是经典的通用阈值,sigma用中位数绝对偏差估而不是标准差,是因为峰高本身有大量离群值,用标准差会把噪声估高。np.clip那一步必须有,重构后会出现小幅负值,直接当峰高用会把后面的统计特征带偏。

拖尾峰识别靠规则,不靠模型:看峰高比和片段长度差。伪峰通常紧贴主峰、峰高低于主峰的某个比例、且片段长度差很小。筛完之后再过一道自适应评分,把质量分低的峰剔掉,剩下的才进 4.3 的匹配打分。

5. 避坑与排查:跑这份配套代码最容易翻车的五个地方

5.1 特征重要性高得离谱,先怀疑标签泄漏

现象:expected_allele_ratio重要性超过 0.5,测试集准确率 0.96,换一批新样本直接掉到 0.7 以下。原因:分母用了真实贡献者人数,等于把标签编码进了特征。解决:按 3.3 节的写法改成对每个候选人数 k 各算一列比值,或者干脆只用allele_count_sum,让人数信息从标签里学。

5.2 交叉验证分数虚高,实际预测一塌糊涂

现象:普通 KFold 交叉验证 0.98,线上新样本 0.75。原因:同一张图谱切出来的 16 条基因座记录被随机分到了训练集和验证集,模型其实在背样本。解决:用StratifiedGroupKFold,groups传SampleID。这个改动会让交叉验证分数掉几个点,但那才是真实水平。

5.3 5 人样本召回率异常低,混淆矩阵全挤在 2 人

现象:5 人样本召回率 0.6 出头,大量 5 人被预测成 2 人。原因:类别不平衡没处理,加上 5 人样本等位基因重叠严重,特征区分度本来就弱。解决:class_weight="balanced"必加;评估指标换成f1_macro而不是 accuracy;必要时对 5 人桶做上采样或调低 GBDT 的l2_regularization。

5.4 峰高阈值写死 50,换个数据集全崩

现象:换了一批不同试剂盒导出的数据,allele_count全线偏高,人数普遍被高估。原因:height_thr=50.0写死在函数默认值里,新数据的基线噪声水平和老数据不一样。解决:把阈值提到config.py,按数据集的峰高分位数动态取(比如取全样本峰高的 5% 分位),并在特征表旁边记下这次用的阈值。

5.5 小波去噪后峰数变少,真峰被削平

现象:去噪后某些基因座的allele_count从 4 掉到 3,后续推断直接错。原因:阈值取得太大,或者对近似层也做了阈值处理,把低峰高的真实等位基因当成噪声滤掉了。解决:只对coeffs[1:](细节层)做软阈值,近似层保留;去噪前后各统计一次峰数,差异超过 20% 就把阈值调小;点数少于 8 的序列直接跳过,不要硬去噪。

6. 二次开发:把论文模型改成能吃新样本的流水线

这份资源的模型是按赛题数据调的,直接拿去跑新样本,八成会遇到字段对不上、阈值不适用、模型加载路径写死这几个问题。我的做法是先把它压成一条统一流水线,再做增量改造。下面这个骨架可以直接抄:

import joblib from features import load_peaks, sample_features, add_expected_ratio class STRPipeline: def __init__(self, num_model_path="model_num.pkl", ratio_dir="."): self.num_model = joblib.load(num_model_path) # 问题一模型 self.ratio_dir = ratio_dir # 问题二按 k 分桶的模型目录 def run(self, peaks_path, k=None): peaks = load_peaks(peaks_path) feat = add_expected_ratio(sample_features(peaks)) if k is None: # 人数未知,先分类 k = int(self.num_model.predict(feat)[0]) ratio_model = joblib.load(f"{self.ratio_dir}/model_ratio_k{k}.pkl") ratio = ratio_model.predict(feat)[0] # 比例类别 return {"contributor_num": k, "ratio_label": ratio}

改造时我最看重的三个点,按优先级排:

第一,把config.py里的常量全部参数化。峰高阈值、基因座列表、OL 判定规则、参考库路径,这些在不同数据源之间差异最大,写死在代码里等于给自己埋雷。我一般会用一个 YAML 存这批常量,跑之前先--config传进去。

第二,模型持久化要带元信息。joblib.dump只存模型对象,不存训练时的特征列顺序和阈值。我习惯把feature_names、height_thr、训练集人数分布一起打进一个 dict 里存,加载时校验特征列顺序,顺序错了直接报错而不是静默给出错误结果。

第三,验证要分层看。整体 accuracy 会掩盖 5 人样本的崩塌,必须按ContributorNum分层打印classification_report,再单独看混淆矩阵对角线以外的分布。我一般会加一个evaluate.py --by-group的参数,跑完自动输出每个 k 的 P/R/F1。

改造项改之前改之后验证方式
峰高阈值硬编码 50配置文件 / 分位数自适应换数据集后allele_count分布是否漂移
特征列顺序依赖 DataFrame 隐式顺序训练时持久化feature_names加载时断言列名一致
比例模型单文件全量按 k 分桶 + 文件命名规范每个 k 单独跑混淆矩阵
人数特征用真实人数归一化枚举 k 生成多列对比新旧特征的下游准确率

最后说个血泪经验。我最开始复现时图省事,直接把论文里的公式按符号抄进代码,结果符号表里w在问题二表示比例、在问题三表示权重、在问题四又表示小波系数,三处混用导致我调了两天才发现比例估计出来的值一直是错的。从那以后我每次做这种跨章节的复现,都会先建一个symbols.md,把每个字母在每一问里的含义写死,代码里的变量名一律加前缀区分。这个习惯看着笨,但能省掉的调试时间远比写文档的时间多。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询