1. 从一道赛题说起STR混合样本识别到底难在哪法医DNA鉴定里有个绕不开的经典难题犯罪现场提取到的生物样本往往不是单一来源。两个人搏斗后留下的血迹、多人接触过的物品表面、混合了受害者和嫌疑人细胞的微量检材——这些样本在STR分型图谱上表现为多个个体等位基因的叠加。实验室拿到一份混合样本需要回答的核心问题是这里面到底有几个人的DNA各自贡献了多少能不能把每个人的基因型拆出来2025年东三省D题深圳杯把这个法医物证学的真实痛点搬上了数学建模赛场。题目要求参赛队伍构建一套智能识别系统从混合STR样本的峰高、峰面积、等位基因位置等电泳图谱数据出发推断混合样本的贡献者数量、混合比例并尽可能还原每个贡献者的基因型。这不是一道纯统计题也不是一道纯优化题它需要你把信号处理、概率建模、组合优化、机器学习几套工具揉在一起用。我拿到这道题的第一反应是单模型绝对做不出来。原因很简单——混合样本的STR数据有三个层面的不确定性叠加在一起。第一层是贡献者数量未知你事先不知道是2人、3人还是4人混合第二层是混合比例未知各人的DNA贡献量可能从1:20到1:1不等第三层是等位基因 dropout 和 drop-in低模板量样本会出现等位基因丢失高模板量又可能出现stutter峰干扰。这三层不确定性互相耦合用单一模型去解要么假设太强导致泛化差要么搜索空间太大导致算不出来。所以我的整体思路是多模型融合用概率模型处理峰高与基因型的映射关系用组合优化搜索最优贡献者组合用机器学习做贡献者数量的预判和噪声过滤最后用贝叶斯推理做基因型的后验修正。下面我把这套系统的完整搭建过程拆开讲包括每一步为什么这么选、参数怎么定、实测中遇到什么问题。2. 数据预处理从电泳图谱到可建模的数值矩阵2.1 原始数据的结构长什么样法医STR分型数据通常以.fsa或.hid格式从测序仪导出经过分析软件如GeneMapper处理后得到一张峰表。这张表的核心字段包括Marker基因座名称、Allele等位基因编号、Height峰高RFU、Area峰面积、Size片段长度bp。一个典型的STR复合扩增体系包含15-24个基因座每个基因座在混合样本中可能出现2-8个峰。我拿到的训练数据是CSV格式的峰表每个样本一行记录列是各基因座的各等位基因峰高。这里第一个坑就出现了不同基因座的等位基因编号范围不同比如TH01的等位基因是3-14而D21S11可以到30以上。如果直接把这些编号当数值特征喂给模型模型会误以为D21S11的等位基因30和TH01的14有数值上的远近关系实际上它们只是标签。处理方式对每个基因座单独做one-hot编码把等位基因编号映射为离散的类别变量。具体来说对基因座$l$假设其等位基因集合为$A_l {a_1, a_2, ..., a_{K_l}}$则每个等位基因对应一个$K_l$维的one-hot向量。这样每个样本在每个基因座上表示为一个$K_l$维的峰高向量$\mathbf{h}_l \in \mathbb{R}^{K_l}$。2.2 峰高归一化与stutter峰剔除STR图谱里有个必须处理的干扰stutter峰。PCR扩增时模板链在重复序列处滑链产生比真实等位基因短一个重复单位的峰通常出现在真实峰的前一个位置高度约为真实峰的5%-15%。如果不剔除模型会把stutter峰当成一个独立的等位基因导致贡献者数量被高估。我的做法是对每个基因座先找出所有峰然后对每个峰检查其右侧片段长度大一个重复单位是否存在更高的峰。如果存在且当前峰高与右侧峰高的比值小于阈值$\tau_{stutter}$通常取0.15则判定为stutter峰并剔除。这个阈值不是拍脑袋定的我参考了法医遗传学文献中不同基因座的stutter比例范围对每个基因座单独设定阈值。实测下来用固定阈值0.15会误删一些真实的低峰改用基因座特异性阈值后召回率提升了约8%。注意stutter剔除的阈值不能设得太激进。我试过用0.10结果在一些低模板量样本中把真实的次要贡献者等位基因也删掉了导致后续贡献者数量估计偏低。2.3 峰高到等位基因计数的转换原始峰高是RFU相对荧光单位数值范围从几十到几万不等。直接拿RFU做建模不同样本之间的尺度差异会淹没真实信号。我采用样本内归一化对每个样本将所有峰的RFU除以该样本所有峰RFU的中位数得到相对峰高。这样处理的好处是保留了样本内各峰之间的相对比例关系同时消除了样本间的绝对亮度差异。归一化后的数据矩阵记为$\mathbf{H} \in \mathbb{R}^{N \times M}$其中$N$是样本数$M$是所有基因座所有等位基因的总数。这个矩阵就是后续所有模型的输入。3. 贡献者数量估计先用无监督聚类定个调3.1 为什么不用简单的峰计数最朴素的想法是数一数每个基因座最多有几个峰取最大值除以2就是贡献者数量。这个思路在理想情况下成立——每个贡献者在每个基因座最多贡献2个等位基因杂合子2个峰纯合子1个峰$k$个贡献者最多产生$2k$个峰。但实际数据里峰的数量受三个因素影响等位基因重叠不同贡献者携带相同等位基因时峰合并、dropout低模板量导致某些等位基因没扩增出来、stutter残留。所以直接数峰会严重低估或高估。我试过在一个包含200个样本的训练集上直接用最大峰数除以2准确率只有52%左右基本上跟瞎猜差不多。3.2 基于峰数分布的混合模型更靠谱的做法是把每个基因座的峰数看作一个随机变量它的分布取决于贡献者数量$k$和混合比例。我构建了一个混合比例感知的峰数期望模型对于$k$个贡献者假设各人的等位基因在人群中的频率服从 Hardy-Weinberg 平衡混合比例为$\boldsymbol{\pi} (\pi_1, ..., \pi_k)$则某个等位基因在混合样本中出现的概率为$1 - \prod_{i1}^{k}(1 - p_i)^{g_i}$其中$p_i$是该等位基因在人群中的频率$g_i$是贡献者$i$携带该等位基因的拷贝数0, 1, 2。基于这个概率模型可以计算每个基因座峰数的期望和方差然后用负二项分布拟合实际观测到的峰数分布。对每个候选的$k$值2, 3, 4计算观测峰数在该$k$下的似然取似然最大的$k$作为估计值。实测下来这个方法在混合比例不太极端最小贡献者占比10%的样本上准确率能到78%左右。但在极端比例如1:20下次要贡献者的峰经常被噪声淹没峰数分布和主要贡献者的单人样本几乎一样模型会误判为2人混合。3.3 用XGBoost做贡献者数量的分类修正为了处理极端比例的情况我引入了一个有监督的分类器做修正。特征工程是关键我构造了以下几类特征峰数特征每个基因座的峰数、峰数的均值/方差/最大值峰高比特征每个基因座中最高峰与最低峰的比例、次高峰与最高峰的比例峰高分布特征峰高的偏度、峰度、熵基因座间一致性特征不同基因座之间峰数的一致性如果某个基因座峰数明显多于其他可能是stutter残留用这些特征训练一个XGBoost多分类器类别为2、3、4人混合在验证集上准确率提升到了86%。更重要的是对于极端比例的样本分类器的召回率比纯概率模型高了近20个百分点。实操心得XGBoost的scale_pos_weight参数在这里很关键。因为训练集中2人混合的样本占大多数约70%3人和4人混合样本较少如果不做类别平衡模型会倾向于预测2人。我用了compute_sample_weight(balanced)来自动调整权重效果比手动设权重稳定。4. 混合比例与基因型推断贝叶斯网络与组合优化的接力4.1 问题形式化从峰高到基因型的概率映射确定了贡献者数量$k$之后下一步是推断每个贡献者的基因型和混合比例。形式化地设第$i$个贡献者在基因座$l$的基因型为$g_{il} \in {(a,b): a,b \in A_l, a \leq b}$混合比例为$\pi_i$$\sum_i \pi_i 1$。观测到的峰高向量$\mathbf{h}_l$由各贡献者的等位基因峰高叠加而成并受到噪声干扰。我采用的概率模型是对于等位基因$a$其期望峰高为 $$\mu_{al} \sum_{i1}^{k} \pi_i \cdot c \cdot n_{ial}$$ 其中$n_{ial}$是贡献者$i$在基因座$l$携带等位基因$a$的拷贝数0, 1, 2$c$是比例常数与总模板量有关。观测峰高$h_{al}$服从以$\mu_{al}$为均值的某种分布我试过正态分布和伽马分布伽马分布对低峰高的拟合更好。4.2 用MCMC采样做后验推断直接计算后验概率$P({g_{il}}, \boldsymbol{\pi} | \mathbf{H})$是不可行的因为基因型的组合空间太大。我采用马尔可夫链蒙特卡洛MCMC方法具体是Metropolis-Hastings算法初始化随机生成一组满足$k$个贡献者的基因型组合和混合比例提议步骤随机选择一个贡献者在其某个基因座上随机改变一个等位基因或者微调混合比例接受/拒绝计算新状态的对数后验概率按Metropolis准则接受或拒绝重复直到链收敛这里有个关键技巧混合比例的提议分布不能太宽。我一开始用均匀分布提议链的接受率极低5%因为混合比例的微小变化会导致似然剧烈波动。后来改用以当前值为中心的Beta分布提议接受率提升到了30%左右收敛速度明显加快。MCMC跑完之后取后验样本的众数作为基因型估计取后验均值作为混合比例估计。实测下来对于2人混合样本主要贡献者的基因型准确率能到95%以上次要贡献者在比例15%时准确率约85%。4.3 组合优化做基因型的精修MCMC的缺点是计算量大而且在高维情况下可能陷入局部最优。为了提升精度我在MCMC之后加了一步组合优化精修固定混合比例把基因型推断转化为一个二次分配问题用模拟退火求解。具体来说对每个基因座枚举所有可能的$k$人基因型组合在等位基因集合较小的时候可行计算每个组合的对数似然取最大的。如果等位基因集合太大就用分支定界剪枝。这一步的计算量比MCMC小得多但能把MCMC中一些明显的错误修正过来。踩坑记录模拟退火的初始温度设得太高会导致搜索时间过长设得太低又容易陷入局部最优。我最后用的是自适应降温策略初始温度设为当前似然标准差的10倍降温系数0.95每个温度下迭代100次。这个参数组合在测试集上表现最稳定。5. 多模型融合策略投票、堆叠与不确定性传播5.1 三个基模型的输出怎么合并到这一步我手里有三个模型的输出模型A基于峰数分布的负二项模型输出贡献者数量$k_A$和粗略的混合比例模型BXGBoost分类器输出贡献者数量$k_B$和各类别概率模型CMCMC模拟退火输出基因型组合和混合比例这三个模型的输出维度不同不能简单平均。我的融合策略是分层融合第一层贡献者数量的融合。用模型B的概率输出作为先验模型A的似然作为证据做贝叶斯模型选择 $$P(k | \text{data}) \propto P(k) \cdot P(\text{data} | k)$$ 其中$P(k)$来自XGBoost的softmax输出$P(\text{data} | k)$来自负二项模型的似然。取后验概率最大的$k$作为最终估计。第二层基因型的融合。以模型C的MCMC后验样本为基础用模型A和B的输出做重要性重加权。具体来说如果模型B对某个$k$的置信度很高而模型C在该$k$下的基因型似然也高则提升该基因型组合的权重。5.2 不确定性传播与置信区间多模型融合的一个额外好处是能做不确定性量化。单一模型给出的基因型估计是一个点估计而融合模型可以给出每个等位基因的后验概率。我定义了一个等位基因置信度指标 $$\text{Conf}(a) \frac{\text{后验样本中包含等位基因}a\text{的样本数}}{\text{总后验样本数}}$$这个指标在实际应用中很有价值法医鉴定报告里低置信度的等位基因需要标注不确定不能作为定罪依据。我在测试集上统计了一下置信度0.9的等位基因实际正确率约97%置信度在0.7-0.9之间的正确率约82%低于0.7的正确率不到60%。5.3 融合模型的实测表现在包含500个混合样本的测试集上2人混合300个3人混合150个4人混合50个融合模型的整体表现如下指标模型A模型B模型C融合模型贡献者数量准确率78%86%-91%主要贡献者基因型准确率--95%97%次要贡献者基因型准确率--85%89%混合比例MAE0.08-0.050.04单样本推理时间0.1s0.05s12s13s融合模型在准确率上全面优于单模型代价是推理时间增加了约1秒主要是MCMC的开销。在实际法医鉴定场景中这个时间代价完全可以接受。6. 代码实现中的关键细节与性能优化6.1 MCMC的向量化加速MCMC是整套系统中最耗时的部分。我最初用纯Python循环实现单个样本要跑30秒以上。后来做了三处优化第一似然计算的向量化。把峰高矩阵和基因型矩阵都用NumPy数组表示用矩阵乘法一次性计算所有等位基因的期望峰高避免Python循环。第二多链并行。用multiprocessing开4条链同时跑每条链跑5000步最后合并后验样本。这样不仅速度快了4倍还能用Gelman-Rubin统计量诊断收敛性。第三自适应步长。在burn-in阶段动态调整提议分布的宽度使接受率维持在25%-40%之间。这个技巧来自贝叶斯统计的经典实践效果立竿见影。优化后的MCMC单样本推理时间降到了3秒左右。6.2 XGBoost的特征工程代码片段import numpy as np import xgboost as xgb from scipy.stats import skew, kurtosis, entropy def extract_features(peak_matrix): peak_matrix: shape (num_loci, max_alleles), 每个基因座的峰高向量 features [] for locus in peak_matrix: peaks locus[locus 0] # 只取非零峰 if len(peaks) 0: features.extend([0, 0, 0, 0, 0, 0]) continue # 峰数 n_peaks len(peaks) # 峰高比 ratio_max_min peaks.max() / (peaks.min() 1e-6) ratio_second_max np.sort(peaks)[-2] / peaks.max() if n_peaks 1 else 0 # 分布特征 sk skew(peaks) if n_peaks 2 else 0 ku kurtosis(peaks) if n_peaks 3 else 0 ent entropy(peaks / peaks.sum()) features.extend([n_peaks, ratio_max_min, ratio_second_max, sk, ku, ent]) # 基因座间一致性 n_peaks_per_locus [np.sum(locus 0) for locus in peak_matrix] consistency np.std(n_peaks_per_locus) / (np.mean(n_peaks_per_locus) 1e-6) features.append(consistency) return np.array(features) # 训练XGBoost分类器 model xgb.XGBClassifier( n_estimators300, max_depth6, learning_rate0.05, subsample0.8, colsample_bytree0.8, objectivemulti:softprob, num_class3 # 2人、3人、4人 )6.3 模拟退火的参数调优记录模拟退火的关键参数是初始温度$T_0$、降温系数$\alpha$和每个温度下的迭代次数$L$。我做了网格搜索$T_0$$\alpha$$L$基因型准确率平均耗时1000.905082%1.2s1000.9510087%2.8s5000.9510089%4.5s5000.9820090%8.2s10000.9510089%5.1s最终选了$T_0500$、$\alpha0.95$、$L100$这组参数在准确率和耗时之间取得了较好的平衡。7. 这道题给数学建模竞赛的几点启示7.1 多模型融合不是万能药但在这类问题上几乎是必选项混合STR样本识别涉及离散变量基因型和连续变量混合比例的联合推断还叠加了模型选择贡献者数量的问题。单一模型很难同时处理好这三个层面。我的经验是用无监督方法做粗估计用有监督方法做修正用贝叶斯方法做精细推断三者各司其职最后融合。但融合也有代价计算复杂度上升、调参难度加大、模型解释性下降。如果题目对推理时间有硬性要求可能需要简化融合策略比如只用两个模型。7.2 领域知识是特征工程的源泉XGBoost分类器的特征里峰高比和基因座间一致性这两个特征对准确率的贡献最大。这两个特征不是拍脑袋想出来的而是来自法医STR分型的领域知识混合比例越极端峰高比越大stutter残留越多基因座间峰数的一致性越差。如果不懂这些背景很难构造出有效的特征。所以做数学建模题尤其是应用题先花时间读背景文献比急着写代码重要得多。我在这道题上花了整整一天读法医遗传学的综述才把问题的物理机制搞清楚。7.3 不确定性量化是加分项很多参赛队伍只给出点估计但法医鉴定场景对置信度的要求极高。我在论文里专门用了一节讲不确定性量化包括等位基因置信度、混合比例的置信区间、贡献者数量的后验概率分布。这部分内容在评审中很受青睐因为它体现了对问题本质的理解——法医鉴定不是猜一个答案而是给出一个带有置信度的判断。7.4 代码可复现性决定论文的可信度数学建模竞赛的论文里方法描述再漂亮如果代码跑不出来评审也会打折扣。我的做法是所有随机过程固定随机种子所有参数在论文中明确列出提供完整的运行环境说明。这样评审如果想复现可以直接跑通。最后分享一个小技巧在论文附录里放一个最小可运行示例用10行代码展示核心模型的调用方式。这比放几百行完整代码更有效评审一眼就能看懂你的方法怎么用。8. 后续可以继续深挖的方向这套系统在2-3人混合样本上表现不错但4人以上混合的准确率明显下降。一个可能的方向是引入深度学习用卷积神经网络处理STR图谱的原始信号自动学习峰形特征可能比手工特征更有效。另一个方向是迁移学习用大规模单人样本预训练一个基因型编码器再在小规模混合样本上微调缓解标注数据不足的问题。混合比例极端如1:50以上的情况仍然是难点。我试过用加权似然给低峰高更高的权重但效果不稳定。可能需要专门针对低模板量样本设计一套独立的推断流程。法医STR混合样本识别这个方向数学建模能做的还有很多。这道题的价值在于它把一个真实的法医物证学问题抽象成了可计算的模型而多模型融合的思路在类似的多源信号解混问题中都有借鉴意义。