STR混合图谱贡献者人数判别:特征工程与GBDT实践
简介:面向2025深圳杯数学建模竞赛D题的完整分析论文,以单份PDF形式呈现,适合参赛选手、指导教师及数学建模学习者用于复盘赛题、理解完整建模流程。包内共1个PDF文件,大小约2.71MB,围绕问题重述、问题分析、模型假设与符号定义展开,并依次给出问题一至问题三的模型建立与求解。正文重点覆盖特征工程设计、分类模型结构与数学表达、模型总体性能对比及不同贡献者数量下的表现,问题三还涉及贡献者兼容性筛选、混合比例优化估计与组合匹配评分函数等内容。论文同时含可运行代码与相关数据说明,便于读者对照模型结构、评估指标和求解结果进行复现与查漏补缺。已有330人学习,适合需要系统掌握赛题方案、模型评估与优化思路的读者参考。
1. 混合图谱里有几个人:为什么数峰个数必然翻车
拿一份混合 STR 图谱问“这里面有几个人的 DNA”,最直觉的做法是数峰:一个位点最多两个等位基因,出现六个峰就该是三个人。这套推理在模拟数据上几乎必然翻车,因为三个人的理论峰数上限是 6 个,实际观测里常常只落到 3 到 4 个——共有等位基因把位置占了,重叠直接把信息吃掉。真正让判别器跑起来的,不是“有多少个峰”这个绝对量,而是观测峰数相对于理论上限的比值结构。STR 混合图谱的贡献者人数判别,本质上是一次对“重叠程度”的建模,而不是一次计数。这条线索值得从头拆到底,因为后面所有关于比例估计、基因型还原、去噪的步骤,全都建立在这个判别结果是否可靠之上。
2. 特征工程:一个比值特征凭什么压过几十个统计量
2.1 四类特征各自在描述什么
从图谱到特征向量,这份分析文档拆了四个维度:等位基因结构特征、峰高统计特征、片段长度分布特征、基因座结构逻辑特征。这不是为了凑数量,四类特征分别对应四种不同的物理机制。
等位基因结构特征描述的是“有没有可能这么多人”。总等位基因数、非 OL 等位基因数(OL 指无法与标准梯度比对的片段)、等位基因数占理论上限的比例,这三者构成一个下界约束:如果某个位点观测到 5 个峰,那贡献者人数至少是 3,不可能是 2。这类特征提供的是硬约束,判别力最强也最不容易被噪声污染。
峰高统计特征描述的是“各贡献者的浓度结构”。峰高的均值、标准差、变异系数刻画离散程度,极值统计和降序排列的前几高峰比值刻画的是主次峰关系。一个人贡献 80% DNA 的两人混合样本,最高峰和次高峰能差一个数量级;两个人各贡献一半,峰高就趋于平坦。top1_top2_ratio这类比值特征直接把这个直觉编码成了数字。
片段长度分布特征描述的是“个体间差异有多大”。Size 的均值、标准差、极差反映的是扩增片段长度的离散程度,当混合个体在这些位点上差异显著时,分布会更宽。这类特征的判别力弱于前两类,但在人数较多、峰高被打散的情况下能提供补充证据。
基因座结构逻辑特征是另一回事:16 个核心标记的 one-hot 编码、OL 等位基因比例、是否存在 k 个等位基因的布尔阈值。它们提供的是上下文,让模型知道“当前这个位点本身的扩增特性是什么样的”,因为不同基因座的杂峰率和扩增效率本来就不一样。
2.2 expected_allele_ratio 权重超过 0.5 的背后
特征重要性排序里,expected_allele_ratio以超过 0.5 的权重遥遥领先,把allele_count、height_count、height_min、height_top4、size_range全都甩在后面。这个结果乍看很反直觉:一个标量比值,凭什么比一整组峰高统计还有用?
原因在于它把“绝对计数”转换成了“归一化后的重叠率”。allele_count的问题是有天花板效应——2 人样本峰数在 4 到 5 之间,3 人样本在 4 到 6 之间,两个区间大面积重叠,单看计数分不开。但除以理论上限之后,2 人样本的比值集中在 1.0 附近(峰数已经饱和),3 人样本落在 0.6 到 0.8,4 人、5 人继续往下掉。重叠越严重,比值越低,这个单调关系比原始计数稳定得多。
提示:这个特征的构造方式需要多问一句。如果分母用的是“真实人数 × 2”,那训练阶段的比值直接由标签算出来,属于标签泄漏,0.96 的准确率会虚高。更稳妥的实现是把人数候选值 k=2..5 分别代入,生成一组
allele_count / (2k)的辅助特征,让模型自己在候选空间里选,推理阶段才可复现。
这也解释了为什么 GBDT 在这个特征空间上远好于其他模型:树模型天然擅长处理这种“分段单调 + 阈值切分”的关系,一次分裂就能把比值小于 0.75 的样本划到 4 人以上那一侧,不需要拟合任何连续函数。
2.3 特征向量的构造代码
单个基因座的特征提取逻辑大致是这样,注意expected_ratio按候选人数展开而不是直接用标签:
import numpy as np def locus_features(peaks, marker_name, core_markers, k_candidates=(2, 3, 4, 5)): """ peaks: [(allele, height, size), ...] 单个基因座的有效峰 k_candidates: 贡献者人数的候选取值,用于展开比值特征 """ h = np.array([p[1] for p in peaks], dtype=float) s = np.array([p[2] for p in peaks], dtype=float) h_desc = np.sort(h)[::-1] feat = {} # 结构特征:绝对量与相对量 feat["allele_count"] = len(h) feat["height_count"] = int((h > 50).sum()) # 有效峰阈值可按仪器灵敏度调 for k in k_candidates: # 分母用候选人数而不是真实人数,推理阶段可复现 feat[f"allele_ratio_k{k}"] = len(h) / (2.0 * k) # 峰高统计:离散度 + 主次峰关系 feat["height_mean"] = h.mean() feat["height_std"] = h.std(ddof=0) feat["height_cv"] = h.std(ddof=0) / (h.mean() + 1e-9) feat["height_min"] = h.min() feat["height_top4"] = h_desc[:4].sum() / (h.sum() + 1e-9) # 前四高峰能量占比 feat["top1_top2_ratio"] = h_desc[0] / (h_desc[1] + 1e-9) if len(h) > 1 else np.nan # 片段长度分布 feat["size_mean"] = s.mean() feat["size_std"] = s.std(ddof=0) feat["size_range"] = s.max() - s.min() # 基因座上下文 feat["is_core_marker"] = int(marker_name in core_markers) return feat几个参数值得说明。height > 50这个阈值不是随便定的,它对应仪器的荧光检测下限,低于这个值的峰大概率是基线噪声而不是真实等位基因,阈值调高会漏掉低比例贡献者的峰,调低会把噪声当信号。height_top4用能量占比而不是绝对高度,是为了消除不同样本上样量差异带来的量纲影响。size_range在单个基因座上通常意义有限,真正起作用的是跨位点聚合后的方差,实际实现里往往按样本再算一层统计量。
把 16 个核心位点的特征拼起来,一个样本就是一个几百维的向量。维度不算高,但特征之间的相关性很强,这也是后面模型选型要考虑的因素。
3. 分类器选型:GBDT 的 0.96 与 MLP 的 0.36 差在哪
3.1 四种模型的结构差异决定了什么
四个候选模型在多分类任务上的表现差距大到不像在同一个数据集上跑出来的:GBDT 准确率 0.9608,MLP 0.7549,随机森林 0.7255,SVM 0.5931。这个排序不是偶然,它精确对应了特征空间的性质。
GBDT 是加法模型,每轮拟合上一轮的残差,最终输出是若干棵回归树的叠加。它的核心优势在于每次分裂都在找“让损失下降最多的那个阈值切分点”,而这份数据里最有判别力的特征恰好是阈值型结构——allele_ratio小于某个值就大概率是 4 人以上,top1_top2_ratio超过某个值就说明存在主导贡献者。树模型对这种分段常数关系几乎不需要学习成本,一次分裂就到位。
随机森林是 Bagging 思路,多棵树投票,每棵树在随机子空间上生长。它的方差控制得好,但偏差压不下去——单棵树的切分粒度有限,投票又会把边界模糊化。0.7255 的准确率说明它对中等复杂度样本还行,对 4 人、5 人这种需要精细切分的场景就力不从心了。
SVM 用 RBF 核做一对多分类,理论上能拟合非线性边界,但它对特征尺度极其敏感,而且一对多策略在类别不平衡时会严重偏向多数类。0.5931 的准确率加上 4 人样本 F1 只有 0.3095,就是典型的“多数类吃掉一切”。
MLP 的情况最有意思。训练集准确率 0.7549 看着比随机森林还高一点,但交叉验证准确率掉到 0.3603——这个 0.36 甚至低于随机猜测(5 分类的基线是 0.2,但类别不平衡下实际基线更低)。这不是模型不行,是它在当前特征维度上直接过拟合了。
3.2 交叉验证暴露的过拟合
单次train_test_split的评估会骗人。MLP 在测试集上 0.75、在 5 折交叉验证上 0.36,差了一倍多,说明它对训练集的划分方式极度敏感。原因有两个层面。
一是特征维度相对于样本量偏高。几百维特征、几百个样本,MLP 的参数量轻松超过样本量,它有能力把训练集里的每一个样本都记住,包括那些由噪声决定的边界。二是特征之间的强相关性让损失曲面出现大量平坦方向,不同的初始化会收敛到泛化能力完全不同的解。
GBDT 在这两点上天然免疫。树模型的容量由树的深度和数量控制,剪枝和最小叶节点样本数就是正则化手段;而且它每次只在一维特征上做切分,不需要同时估计多个特征的联合分布,对特征相关性不敏感。交叉验证 0.9584 ± 0.0136 的波动幅度,说明它对数据划分基本无感。
注意:看到“测试集准确率很高、交叉验证很低”这种组合,第一反应应该是查数据泄漏和样本量,而不是调超参。这份分析里 MLP 的 0.36 就是最直接的信号。
3.3 分人数拆解:3 人以上为什么会塌
把准确率按贡献者人数拆开看,问题结构立刻清楚:
| 贡献者人数 | GBDT F1 | 随机森林 F1 | MLP F1 | SVM F1 |
|---|---|---|---|---|
| 2 人 | 0.9916 | 0.8189 | 0.9027 | 0.8403 |
| 3 人 | 0.9466 | 0.6957 | 0.7786 | 0.6515 |
| 4 人 | 0.9663 | 0.6316 | 0.5977 | 0.3095 |
| 5 人 | 0.9275 | 0.7164 | 0.6753 | 0.4110 |
GBDT 在四个类别上都保持在 0.92 以上,而其他三个模型从 3 人开始就集体下滑。随机森林在 4 人样本上掉到 0.6316,SVM 直接崩到 0.3095。
根源在等位基因重叠的数学性质。2 人混合时,两个个体在每个位点最多贡献 4 个不同等位基因,观测到的峰数接近上限,allele_ratio接近 1,特征区分度极高。到 4 人、5 人时,理论峰数上限是 8 到 10,但实际观测往往只有 5 到 6 个——大量等位基因在人群中共有,重叠把可观测的信号压到了一个很窄的区间里。这时候不同人数的特征分布大面积交叠,模型必须依赖峰高的细微结构来分辨,而峰高本身受 PCR 扩增偏倚影响,噪声很大。
GBDT 能扛住这个场景,靠的是它能把多个弱特征组合起来做联合判断:单看allele_ratio分不开 4 人和 5 人,但加上height_cv、height_top4、non_ol_allele_count的组合切分,就能在交叠区域里划出还不错的边界。这也是为什么特征重要性里除了比值特征,height_min、size_range这些看起来边缘的统计量还能排进前十五——它们在边界样本上提供了决定性的补充信息。
4. 比例估计与基因型推断:离散化、最小二乘与组合搜索怎么接上
4.1 把连续比例拆成离散类别,图的是什么
问题二要在已知人数前提下估计各贡献者的 DNA 比例。最直接的想法是回归,输出一组和为 1 的连续值。这份分析选了另一条路:把比例离散化成类别,2 人样本定义 1:1、1:4、1:9 三类,3 人样本六类,4 人八类,5 人七类,然后对每个人数单独训练一个分类器。
这个取舍值得琢磨。回归看起来更精细,但混合比例本身是个有歧义的目标——“1:4:1”和“1:1:4”归一化之后重排一下就是同一个比例,标签的定义需要先做规范化(通常按降序排列),否则模型学到的是一堆等价类。离散化之后,每个类别对应一个明确的生物学含义(谁是主要贡献者、次要贡献者占多少),标签一致性容易保证。
代价是精度损失。真实案件里的比例是连续分布的,1:3.7 和 1:4 被强行归到同一类,模型学到的边界是人为设定的。这个损失在 2 人样本上还能接受(准确率 85%),到 3 人样本直接掉到 25%——六类里的每一类都要求三个比例同时判对,误差累积得非常快。
4.2 峰高线性假设与最小二乘反解
所有比例估计方法的地基是一个假设:峰高与 DNA 含量近似线性相关。用公式说,某个位点上第 i 个贡献者贡献的等位基因峰高,正比于他的混合比例乘以该等位基因的扩增效率。
在这个假设下,比例估计可以写成最小二乘问题。设 A 是“贡献者 × 等位基因”的指示矩阵(某个贡献者携带某等位基因则对应位置为 1),y 是观测到的峰高向量,w 是待求的比例向量,那么求解min ||Aw - y||²且w ≥ 0就得到比例。用非负最小二乘(NNLS)实现:
import numpy as np from scipy.optimize import nnls def estimate_ratio(genotype_matrix, observed_heights, max_contributors=5): """ genotype_matrix: shape (n_alleles, n_candidates),候选贡献者的等位基因指示 observed_heights: shape (n_alleles,),观测峰高 返回归一化后的比例向量 """ # 加一列全 1 作为基线项,吸收 PCR 偏倚带来的整体增益 A = np.hstack([genotype_matrix, np.ones((genotype_matrix.shape[0], 1))]) w, residual = nnls(A, observed_heights) ratio = w[:-1] if ratio.sum() < 1e-9: return np.full(ratio.shape, 1.0 / len(ratio)), residual return ratio / ratio.sum(), residual # 残差越小,组合越可信这段代码的关键在两点。加基线列是为了处理不同样本的扩增效率差异,否则所有峰高会整体偏高或偏低,比例虽然对但残差没意义。归一化是必须的,因为 NNLS 输出的是绝对幅度,只有相对大小才有生物学含义。
残差是这个方法最有价值的副产品。它直接衡量“这组候选者能不能解释观测到的峰高”,残差大说明候选组合选错了,或者样本里还有没被纳入的贡献者。问题三的组合搜索就是靠这个残差排序的。
4.3 兼容性筛选与组合评分的两段式设计
问题三要在已知候选个体基因型数据库的条件下,推断混合样本里有哪几个人。这是一个组合搜索问题:候选池有 N 个人,要选出 k 个人,组合数是 C(N, k),N 上百的时候直接枚举不现实。
这份分析用了两段式:先用兼容性评分粗筛,再用匹配评分精排。
兼容性评分很朴素——计算候选个体的基因型与混合样本观测等位基因的重合比例。完全重合得满分,部分重合扣分。这一步的计算量是 O(N × 位点数),非常快,可以把候选池从上百人压到几十人。它的问题是只用了“有没有”的信息,完全忽略了峰高,所以只能做粗筛。
精排阶段用的是组合匹配评分函数,把两个量加权组合:等位基因解释率(这组人的等位基因能覆盖观测峰的多少)和峰高分布相关性(用这组人反解出的预测峰高与观测峰高的相关程度)。得分为两者加权,最高分的组合胜出。
这个设计的好处是把计算量分布到了两个阶段:粗筛便宜但粗糙,精排贵但只在小候选集上跑。代价是粗筛可能把正确答案筛掉——如果某个真实贡献者的基因型因为突变或检出缺失导致兼容性评分偏低,他在第一阶段就被淘汰了,后面再精确也救不回来。
提示:兼容性筛选的阈值不能设太死。实际实现里通常保留 top-M 而不是卡死分数阈值,M 取候选池的 10% 到 20%,宁可多留一些让精排去判断。
5. 去噪与排查:小波阈值调错一步,准确率反而往下掉
5.1 小波软阈值去噪在峰高信号上做了什么
混合图谱里的噪声有两类:仪器基线噪声(高频、低幅度、随机分布)和 PCR 扩增伪峰(拖尾峰,紧邻主峰、高度比主峰低一个量级)。小波阈值去噪处理的是第一类。
小波分解把峰高信号拆成多个频带,近似分量保留整体轮廓,细节分量承载高频波动。噪声几乎全部落在最高频的细节分量里,而真实的等位基因峰是有限宽度的脉冲,能量集中在低频。软阈值处理就是对每个细节系数做收缩:绝对值小于阈值的直接置零,大于阈值的向零方向压缩一个阈值量。
软阈值而不是硬阈值,是因为硬阈值会在系数上产生跳变,重构回来的信号出现振铃,等位基因峰的边缘被削出假的小峰——这在后续的峰高统计里会被当成真实信号,直接污染比例估计。
import numpy as np import pywt def wavelet_denoise(signal, wavelet="db4", level=3): """ signal: 单个基因座沿片段长度方向的峰高序列 level: 分解层数,按信号长度取 log2(len)-2 较稳 """ coeffs = pywt.wavedec(signal, wavelet, level=level) # 用最高频细节系数的中位数估计噪声标准差,比用全局方差稳健 sigma = np.median(np.abs(coeffs[-1])) / 0.6745 # 通用阈值:随信号长度对数增长 thr = sigma * np.sqrt(2 * np.log(len(signal))) # 第 0 层是近似分量,不动;其余层做软阈值 coeffs[1:] = [pywt.threshold(c, thr, mode="soft") for c in coeffs[1:]] return pywt.waverec(coeffs, wavelet)[:len(signal)]db4小波的支撑长度适中,对脉冲型信号的重构误差小。level=3是个经验值,层数太少噪声压不干净,层数太多会把宽峰的低频成分也当成噪声削掉。0.6745是中位数绝对偏差到标准差的换算系数,用中位数而不是均值估噪声,是因为细节系数里混着少量真实峰边缘的大值,用均值会被拉高。
这个阈值的作用是“宁可放过一点噪声,也不要削掉低比例贡献者的峰”。低比例贡献者的峰高只有主峰的百分之几,阈值稍微大一点,这些峰就被压到噪声水平以下,后面的人数判别和比例估计全部失真。
5.2 拖尾峰识别规则与自适应峰值筛选
拖尾峰(stutter)是 PCR 扩增过程中聚合酶滑动产生的伪峰,特征很明确:紧邻真实峰、片段长度短一个重复单元、高度通常是主峰的 5% 到 15%。
识别规则就是两条组合:峰高比值落在 [0.03, 0.20] 区间,且与相邻主峰的片段长度差等于一个重复单元的长度。两条同时满足才判为拖尾峰,只满足一条不动——因为真实的低比例贡献者峰也可能落在同样的高度比区间,只靠高度比会误杀。
筛掉拖尾峰之后还有一步自适应峰值筛选:对候选峰做质量评分,综合峰高绝对幅度、峰形对称性、与相邻峰的间距,低于阈值的剔除。这一步的作用是给不同样本自适应调整灵敏度,因为不同样本的上样量差异很大,固定阈值不可行。
这个流程跑完之后,贡献者识别的准确率在多数样本上能到 0.8 到 1.0,极端样本上也优于只做预处理的方案。
5.3 四条血泪踩坑记录
现象一:交叉验证准确率 0.36,训练集 0.75,怀疑模型写错了。原因:单次train_test_split的随机性在样本量只有几百的时候影响巨大,某一折恰好抽到了特征分布不同的样本,评估结果就完全失真。MLP 这类高方差模型对划分方式尤其敏感。 解决:强制走 StratifiedKFold 至少 5 折,看均值和标准差。标准差超过 0.05 就先别调模型,回去查特征和样本量。
现象二:3 人样本的比例分类准确率只有 25%,六个类别等于瞎猜。原因:两个叠加问题。一是类别极度不平衡,1:4:1占了 27%,剩下五类各占百分之十几,模型倾向于全predict多数类;二是比例语义本身有歧义,1:4:1和1:1:4在重排后是同一类,标签生成逻辑如果不统一,同一样本会被打上不同标签。 解决:标签生成前先做降序规范化,并且用class_weight="balanced"或者分层采样保证每类样本量可比。25% 这个数字已经低于 1/6 的均匀基线,说明模型学到的东西是负的。
现象三:加了小波去噪之后,整体准确率反而从 0.72 掉到 0.65。原因:阈值取得太大,把低比例贡献者的次峰一起压掉了。这些人本来就是最难分辨的样本,去噪把他们仅有的信号抹平,模型只能靠猜。 解决:把阈值系数从 1.0 降到 0.6 左右试,或者改成按位点自适应——每个位点单独估噪声水平,而不是全样本共用一个阈值。验证方法是看去噪前后“峰数量”的分布变化,如果平均峰数下降超过 15%,基本可以确定削过头了。
现象四:训练时特征维度 320,推理时报维度不匹配。原因:基因座的 one-hot 编码是按训练集里出现过的基因座名生成的,推理数据里如果出现了训练集没见过的标记名,编码器会直接报错或者静默丢弃。 解决:编码器固定 16 个核心标记的白名单,遇到白名单外的标记统一归到“其他”类,绝不动态扩展维度。这个坑在跨批次数据上几乎必踩一次。
现象五:同一份数据跑两次 GBDT,准确率差了 3 个百分点。原因:没有固定随机种子,而且特征里混进了行号或者样本 ID 这类无意义但因为数据排序而带信息的列。 解决:全局random_state固定,特征构造阶段显式排除 ID 类字段。3 个百分点的波动在 0.96 这个量级上不算小,足以影响模型选型的判断。
6. 这份模型真正值钱的取舍在哪
整条链路里最脆弱的不是模型,是六个假设中的第二个:峰高与 DNA 比例近似线性相关。这份分析里所有比例估计、组合评分、残差排序,全部建立在这个线性关系上。真实样本里 PCR 扩增对不同长度片段的效率差异能到 20% 以上,长片段扩增效率低,峰高被系统性压低。这个偏差在 2 人样本上还能靠基线列吸收一部分,到 4 人、5 人时不同贡献者的片段长度分布交错,偏差就没法用单一系数校正了。
第二个值得盯的取舍是离散化的粒度。比例分成三类还是十类,直接决定了问题二的准确率天花板。三类的 85% 里有多少是因为类别少而虚高,很难拆开看。更合理的评估是同时报一个连续误差指标,比如估计比例与真实比例的 MAE,这样即使分类准确率低,也能看出模型是不是“差得不远”。
第三个是组合搜索的搜索空间。兼容性粗筛把候选池压到几十人,精度上去了,但如果真实贡献者在粗筛阶段被误杀,后面所有步骤都是白费。这类方法的失败模式是“静默错误”——它不会报错,只会给出一个看起来合理但完全错误的答案。实际使用里应该保留 top-M 的完整排序,人工复核排名靠前但得分接近的几组。
验证这套流程是否可信,我一般会做两件事:一是构造已知答案的合成样本,把噪声、拖尾峰、比例偏移逐项加上去,看准确率是在哪一步开始掉的;二是拿掉去噪模块跑一遍,如果准确率没变化,说明去噪根本没起作用,或者阈值根本没设对。
从那以后我每次碰这类多阶段流水线,都会强制走一遍“逐模块消融”——把每个模块单独关掉,看指标掉多少。掉得多的模块才值得调,掉不掉的模块要么是冗余,要么是根本没生效。希望帮到你。
本文还有配套的精品资源,点击获取