简介这份文档是2000年全国大学生数学建模竞赛DNA序列分类赛题的完整解答资料面向参加数学建模竞赛的学生、生物信息学入门者以及需要模式识别案例的读者。资源包内仅含1个doc文件约228KB集中呈现了从问题重述、模型假设到建模求解的全过程。文档以20个已知类别的人工DNA序列为学习样本先统计1字符串、2字符串、3字符串的出现频率构成含41个变量的基本特征集再用主成分分析法提取4个特征以降维去噪最后采用Fisher线性判别法完成分类并给出20个人工序列与182个自然序列的具体分类结果。读者可从中获得完整的特征提取思路、主成分分析与判别分析的应用范例以及检验模型效率的方法适合作为数学建模竞赛训练与模式识别学习的参考。目前已有216人学习。1. 从 41 维频率表到 4 维主成分这份 2000 年赛题文档为什么今天还值得拆如果你手头只有一份 Word 文档里面塞满了频率表、Fisher 判别公式和一段 Fortran 代码你大概会犹豫要不要打开它。这份《DNA序列分类2000年数学建模竞赛题》就是这种资源——它不是教程是一份完整的赛题解答存档包含 20 个人工序列和 182 个自然序列的分类全过程。核心思路很清晰把 A、T、C、G 四种碱基的 1 字符串、2 字符串、3 字符串出现频率全部统计出来拼成 41 维特征向量再用主成分分析压到 4 维最后用 Fisher 线性判别做二分类。适合正在准备数学建模竞赛、想复现经典分类流程、或者需要一套“频率特征 PCA Fisher”完整代码骨架的人。文档里连每个样本的 AT 含量、64 种密码子压缩成 20 类的映射关系都列了表省去了自己查遗传密码表的功夫。2. 特征工程把 41 维频率向量拆成可复现的统计步骤2.1 为什么选 1/2/3 字符串频率而不是 k-mer 谱DNA 序列分类在生物信息学里常见做法是用 k-mer 频率作为特征但 k 取多少、怎么压缩直接决定后续分类器能不能跑通。这份文档选了 1、2、3 三种长度1 字符串给出 A、T、C、G 的全局比例2 字符串捕捉相邻碱基的偏好性3 字符串对应密码子直接关联氨基酸编码。关键决策在于 3 字符串没有直接用 64 维而是按图 1 的密码子简并表压缩成 20 类对应 20 种氨基酸。这样特征维度从 4166484 降到 4162040再加一个 AT 含量正好 41 维。文档里明确写了假设三联组的起始位置不影响分类64 种密码子压缩为 20 组也不影响结果。这两个假设是整条流水线能跑通的前提复现时不要跳过。2.2 滚动窗口统计的代码实现统计频率时用的是“滚动”算法比如序列 ATTCG2 字符串取 AT、TT、TC、CG 共 4 个3 字符串取 ATT、TTC、TCG 共 3 个。下面用 Python 重写这段逻辑输入是纯文本序列输出是 41 维特征向量。import numpy as np from collections import Counter # 密码子简并表20 类每类对应若干 3 字符串 codon_groups { b1: [aaa,ata], b2: [aca,aga], b3: [cac,ctc], b4: [ccc,cgc], b5: [gag,gtg], b6: [gcg,ggg], b7: [tat,ttt], b8: [tct,tgt], b9: [aac,caa,atc,cta], b10: [aag,gaa,atg,gta], b11: [aat,taa,att,tta], b12: [acc,cca,agc,cga], b13: [acg,gac,ctg,gtc], b14: [act,tca,agt,tga], b15: [cag,gac,ctt,ttc], b16: [cat,tac,ctt,ttc], b17: [ccg,gcc,cgg,ggc], b18: [cct,tcc,cgt,tgc], b19: [gat,tag,gtt,ttg], b20: [gct,tcg,ggt,tgg] } def kmer_freq(seq, k): 滚动窗口统计 k-mer 频率返回归一化后的字典 seq seq.lower().replace(\n, ).replace( , ) kmers [seq[i:ik] for i in range(len(seq)-k1)] cnt Counter(kmers) total sum(cnt.values()) return {kmer: cnt[kmer]/total*100 for kmer in cnt} def extract_features(seq): 输入一条 DNA 序列输出 41 维特征向量 seq seq.lower().replace(\n, ).replace( , ) feats [] # 1 字符串A, C, T, G 频率 f1 kmer_freq(seq, 1) for base in [a, c, t, g]: feats.append(f1.get(base, 0.0)) # AT 含量 feats.append(f1.get(a, 0.0) f1.get(t, 0.0)) # 2 字符串16 种 f2 kmer_freq(seq, 2) for b1 in actg: for b2 in actg: feats.append(f2.get(b1b2, 0.0)) # 3 字符串压缩成 20 类 f3 kmer_freq(seq, 3) for group, codons in codon_groups.items(): freq sum(f3.get(c, 0.0) for c in codons) feats.append(freq) return np.array(feats)这段代码里kmer_freq用Counter统计所有 k-mer 出现次数再除以总数得到百分比。extract_features按 1 字符串、AT、2 字符串、3 字符串分组的顺序拼接最终长度是 41162041。注意 3 字符串部分遍历codon_groups的 20 个键把每组内所有密码子的频率加总这样即使某个密码子没出现对应组频率也是 0不会缺维度。参数方面k只取 1、2、3不要改成 4 或 5因为文档的 41 维定义就是基于这三个长度。如果序列长度小于 3kmer_freq会返回空字典后续.get全部取 0特征向量全零这种短序列在分类时会被判为无法分类文档里也允许“无法分类的不写入”。2.3 主成分分析降维从 41 维到 4 维的阈值选择41 维特征直接送进 Fisher 判别也能算但文档明确说样本数只有 20模式识别一般要求样本数至少是变量数的 3 倍否则结果不可靠。所以先用 PCA 降维。文档计算后前 4 个主成分累计贡献率达到 96%于是取 4 维。复现时用sklearn.decomposition.PCA即可但要注意PCA 的n_components不要写死 4而是根据累计贡献率阈值动态选。from sklearn.decomposition import PCA def pca_reduce(X, var_threshold0.96): X: shape (n_samples, 41) 返回降维后的特征和 PCA 对象 pca PCA() pca.fit(X) cum_var np.cumsum(pca.explained_variance_ratio_) n_components np.searchsorted(cum_var, var_threshold) 1 print(f累计贡献率达到 {var_threshold:.0%} 需要 {n_components} 个主成分) X_reduced pca.transform(X)[:, :n_components] return X_reduced, pcavar_threshold默认 0.96和文档一致。searchsorted找到第一个累计贡献率超过阈值的位置加 1 是因为索引从 0 开始。打印出来的n_components应该等于 4如果数据预处理不同导致结果有偏差可以微调阈值到 0.95 或 0.97。注意 PCA 之前要不要标准化文档没有提标准化因为频率特征本身量纲一致都是百分比所以直接做 PCA 没问题。但如果你的序列长度差异极大建议先做 z-score 标准化否则方差大的特征会主导主成分方向。3. Fisher 线性判别从两类均值到分类门槛的完整推导3.1 判别函数的矩阵形式与参数估计Fisher 线性判别的目标是找一个投影方向让两类样本投影后的均值差最大、类内方差和最小。文档给出的判别函数是U(x) (X̄₁ - X̄₂)ᵀ (Σ₁ Σ₂)⁻¹ X其中 X̄₁、X̄₂ 是两类样本的均值向量Σ₁、Σ₂ 是两类样本的协方差矩阵。分类门槛值 U₀ 取两类均值中点处的判别值U₀ U(α·X̄₁ (1-α)·X̄₂)α1/2因为两类样本数相等各 10 个所以 α 取 1/2 合理。如果两类样本数不等α 应该按样本比例调整否则门槛会偏向样本多的那一类。下面用 Python 实现 Fisher 判别训练和预测。def fisher_train(X, y): X: shape (n_samples, n_features) y: 标签数组0 或 1 返回判别向量 w 和门槛值 u0 X0 X[y 0] X1 X[y 1] mean0 np.mean(X0, axis0) mean1 np.mean(X1, axis0) cov0 np.cov(X0, rowvarFalse) cov1 np.cov(X1, rowvarFalse) # 防止协方差矩阵奇异加一个小正则项 reg 1e-6 * np.eye(X.shape[1]) w np.linalg.solve(cov0 cov1 reg, mean0 - mean1) # 门槛值两类均值中点处的投影值 u0 w (0.5 * mean0 0.5 * mean1) return w, u0 def fisher_predict(X, w, u0): 返回预测标签投影值大于 u0 为类 0否则为类 1 scores X w return np.where(scores u0, 0, 1)fisher_train里np.linalg.solve解线性方程组 (Σ₁Σ₂)w (X̄₁-X̄₂)比直接求逆更稳定。加reg是因为 4 维特征下协方差矩阵可能接近奇异尤其样本少的时候。u0用两类均值中点计算对应文档的 α1/2。预测时scores u0判为类 0否则类 1这个方向取决于w的符号如果发现预测结果和标签反了把改成即可。参数方面reg的大小可以调整1e-6 是保守值如果数据量级很大可以适当放大到 1e-4。3.2 留一法交叉验证20 个样本怎么验证模型效率文档用了留一法每次取出一个学习样本用其余 19 个训练预测取出的那个结果 20 个学习样本预报成功率 100%。同时对未知的 20 个人工序列做预报发现除取出 4、15、17、20 号样本外其余 16 次预报结果完全一致都是 22、23、25、27、29、34、35、36、37 为 A 类。这说明模型对训练集扰动不敏感稳定性好。复现时可以用sklearn.model_selection.LeaveOneOut快速实现。from sklearn.model_selection import LeaveOneOut def loo_validate(X, y): loo LeaveOneOut() correct 0 for train_idx, test_idx in loo.split(X): w, u0 fisher_train(X[train_idx], y[train_idx]) pred fisher_predict(X[test_idx], w, u0) if pred[0] y[test_idx][0]: correct 1 print(f留一法正确率{correct}/{len(y)} {correct/len(y):.1%}) return correct / len(y)LeaveOneOut对每个样本做一次训练和预测20 个样本就是 20 次循环。fisher_train内部对协方差加正则避免某次训练时协方差奇异导致报错。如果正确率不是 100%检查特征提取是否和文档一致尤其是 3 字符串压缩表有没有抄错。文档里 b15 和 b16 的密码子列表有重复ctt、ttc 同时出现在两组这是原文的笔误还是有意为之复现时建议按标准遗传密码表修正否则特征会有冗余。3.3 对 182 个自然序列的分类结果解读文档最后给出了 182 个自然序列的分类A 类 142 个B 类 40 个。B 类序号是 1、4、8、10、27、29、32、41、43、48、54、63、70、72、75、76、81、86、90、92、102、110、116、119、126、131、144、150、157、159、160、161、162、163、164、165、166、169、170、182。这个结果可以直接作为基准如果你用自己的代码跑出来 B 类数量差很多大概率是特征提取或 PCA 维度不对。注意自然序列比人工序列长很多滚动窗口统计时计算量会大一些但 182 条序列用 Python 跑也就几秒钟。如果序列里有非 ATCG 字符比如 N 表示未知碱基kmer_freq会把它当成独立字符统计导致频率分母变大、特征偏移。常见做法是先把 N 替换成 A 或直接剔除含 N 的窗口文档没有提这一点但实际数据里经常遇到。4. 避坑与排查复现这份 2000 年赛题时最容易翻车的五个地方4.1 现象PCA 累计贡献率到不了 96%只能到 90% 左右原因特征提取时 3 字符串压缩表用错了或者 2 字符串统计时没有用滚动窗口而是按不重叠切分。文档明确说“滚动”算法ATTCG 的 2 字符串是 AT、TT、TC、CG 共 4 个不是 AT、TC 两个。如果按不重叠切分2 字符串频率会严重失真PCA 的方差结构改变累计贡献率自然下降。解决检查kmer_freq里range(len(seq)-k1)这个循环确保窗口每次移动 1 位。另外确认 3 字符串压缩表是否和文档表 3 下方的 b1~b20 定义一致尤其注意 b15 和 b16 的重复项建议按标准密码子表重新整理。4.2 现象Fisher 判别预测时所有样本都被判为同一类原因协方差矩阵奇异np.linalg.solve解出来的w全是零或者极大值导致投影值全部落在门槛同一侧。20 个样本、4 维特征协方差矩阵接近奇异是常态。解决在fisher_train里加正则项reg 1e-6 * np.eye(n_features)如果还不行就增大到 1e-4。另一个办法是先对特征做标准化让每个维度的方差为 1这样协方差矩阵的条件数会改善。文档里没有提正则因为当年可能用 Fortran 直接算的数值精度和现在不同。4.3 现象留一法正确率不是 100%总有 1~2 个样本分错原因标签编码反了。文档说 1~10 为 A 类11~20 为 B 类但 Fisher 判别里类 0 和类 1 的对应关系取决于w的符号。如果fisher_predict里scores u0判为类 0而你的标签里 A 类编码为 1就会全部反掉。解决统一标签编码比如 A 类为 0、B 类为 1然后在fisher_train里用mean0 - mean1计算w预测时scores u0判为类 0。如果发现正确率是 0% 而不是 100%直接把预测结果取反即可。4.4 现象对 182 个自然序列分类时B 类数量远多于 40 个原因自然序列长度差异大频率特征没有做长度归一化。虽然kmer_freq返回的是百分比但短序列的频率波动更大PCA 降维后短序列容易落在 B 类区域。文档里的自然序列“都较长”但具体长度分布没有给出如果某条序列特别短特征会异常。解决在extract_features开头加一个长度过滤比如if len(seq) 30: return None分类时跳过这些序列。或者对频率特征做对数变换压缩动态范围。文档没有提这一点但实际数据里短序列是常见的干扰源。4.5 现象Fortran 代码附录里的READ*,LINE读不进去数据原因文档附录一的 Fortran 代码是片段CHARACTER*121 LINE(40)声明了 40 行、每行 121 字符但READ*,LINE需要外部文件按固定格式提供输入。如果你直接复制这段代码编译运行没有对应的输入文件会卡在读取阶段。解决不用纠结 Fortran 代码用 Python 重写更现实。如果一定要跑 Fortran把序列按每行 121 字符、共 40 行的格式写入t1.dat然后OPEN(5, FILEt1.dat, STATUSOLD)读取。注意文档里READ*,LINE是从标准输入读不是从文件读需要改成READ(5,*) LINE才能读文件。5. 进阶技巧用 sklearn 的 LDA 替代手写 Fisher 并做分类边界可视化手写 Fisher 判别能帮你理解原理但实际项目里直接用sklearn.discriminant_analysis.LinearDiscriminantAnalysis更省事而且它内部做了数值优化不容易遇到奇异矩阵问题。下面用 LDA 重跑一遍完整流程并画出 4 维特征降到 2 维后的分类边界方便你直观检查两类是否线性可分。import matplotlib.pyplot as plt from sklearn.discriminant_analysis import LinearDiscriminantAnalysis from sklearn.decomposition import PCA # 假设 X 是 20x41 的特征矩阵y 是 20 个标签A0, B1 # 先用 PCA 降到 4 维 pca PCA(n_components4) X_pca pca.fit_transform(X) # 用 LDA 训练 lda LinearDiscriminantAnalysis() lda.fit(X_pca, y) # 为了可视化再用 PCA 把 4 维降到 2 维 pca2 PCA(n_components2) X_2d pca2.fit_transform(X_pca) # 生成网格点预测分类边界 x_min, x_max X_2d[:, 0].min() - 1, X_2d[:, 0].max() 1 y_min, y_max X_2d[:, 1].min() - 1, X_2d[:, 1].max() 1 xx, yy np.meshgrid(np.linspace(x_min, x_max, 200), np.linspace(y_min, y_max, 200)) # 把网格点映射回 4 维空间需要逆变换这里简化处理 # 直接用 LDA 在 2 维上的投影方向画线 # 更严谨的做法是用 pca2.inverse_transform 再送进 lda grid_2d np.c_[xx.ravel(), yy.ravel()] grid_4d pca2.inverse_transform(grid_2d) Z lda.predict(grid_4d) Z Z.reshape(xx.shape) plt.contourf(xx, yy, Z, alpha0.3, cmapcoolwarm) plt.scatter(X_2d[y0, 0], X_2d[y0, 1], cblue, labelA 类) plt.scatter(X_2d[y1, 0], X_2d[y1, 1], cred, labelB 类) plt.xlabel(PC1) plt.ylabel(PC2) plt.legend() plt.title(Fisher/LDA 分类边界PCA 降维后) plt.show()这段代码先用 PCA 把 41 维降到 4 维再用 LDA 训练最后为了画图又降到 2 维。pca2.inverse_transform把 2 维网格点映射回 4 维再送进lda.predict得到分类边界。注意inverse_transform只能近似还原因为 PCA 降维是有损的所以边界图仅供参考不能替代原始 4 维空间的准确率评估。如果你发现边界图里两类点混在一起但留一法正确率又是 100%说明 2 维可视化丢失了区分信息以 4 维空间的预测为准。另一个进阶技巧是调整 PCA 的whiten参数。PCA(n_components4, whitenTrue)会让每个主成分的方差为 1这样 Fisher 判别对特征尺度不敏感协方差矩阵更稳定。文档里没有用白化但如果你遇到fisher_train解不稳定可以试试打开白化。参数whitenTrue之后pca.transform返回的特征已经除以了奇异值后续 Fisher 的reg可以调小到 1e-8。最后说一个我自己的习惯每次复现这类老赛题我都会先把文档里的结果表抄下来当基准然后跑自己的代码逐项对比。比如 20 个人工序列的 A 类序号是 22、23、25、27、29、34、35、36、37如果我的代码跑出来多一个 30 或者少一个 34我就知道特征提取或者分类门槛有问题。这种“对答案”的方式比盲目调参快得多。从那以后我每次拿到带结果的赛题文档都强制走一遍“抄基准 → 跑代码 → 逐项对比”的流程省得在玄学调参上浪费时间。希望帮到你。本文还有配套的精品资源点击获取
