企业官网建设流程全解析

简介：这份PDF为《基于机器学习方法的胃癌分型标志基因提取》论文原文，发表于《中国生物医学工程学报》2009年，面向生物信息学、医学数据挖掘方向的研究者与学生，可作为专业指导与参考文献使用。压缩包内仅含1个PDF文件，大小约722KB，内容覆盖论文全文，已有117人学习浏览。研究以33例中国人的胃癌Oligo基因芯片数据为基础，采用SAM、PLS与BD-SFS结合的多步骤降维方法，从21378个基因中筛选出20个弥漫型与肠型胃癌的区分特征基因；随后用SVM分类器达到89.43%准确率，用层次聚类进一步验证到93.94%。文章还分析了所选标志基因的生物学意义，指出大部分基因与人类恶性肿瘤的诊断和分型密切相关，可为胃癌早期检测与靶向治疗提供潜在标记物。读者可借此掌握高维基因表达谱中特征选择、分类建模与聚类验证的完整流程，也可将方法迁移至其他肿瘤分型或标志基因筛选任务。

1. 高维小样本的经典示范：怎么拿两万个基因筛出20个胃癌分型标志基因

机器学习在生信上最典型的场景之一，就是拿两万个基因的表达值去回答一个临床问题：这个胃癌样本到底是弥漫型还是肠型。这篇论文给的答案很干脆——33 个样本、21378 维基因，先做三层降维，最后锁定了 20 个基因。用这 20 个基因做 SVM 分类，准确率 89.45%；做层次聚类，准确率 93.94%，只有 2 例肠型样本被分到弥漫型那边。它适合两类人：一类是刚拿到芯片或 RNA-seq 数据、样本量却只有几十例的从业者，想知道在 p 远大于 n 时怎么选特征；另一类是想把一套成熟的生信分析流程搬到自己数据集上的人，尤其是 SAM、PLS-VIP、SFS 这套组合拳，每一步的坑和参数思路都值得拆开看。

2. 数据长什么样：33例样本、21378个基因、0.6%缺失，先解决三个预处理问题

2.1 Lauren分型的底子和Oligo芯片的数据形态

这篇论文用的是 Lauren 分型，把胃癌分成肠型和弥漫型两类。样本构成是 13 例弥漫型加 20 例肠型，一共 33 例，全部来自中国人，由北京市肿瘤防治研究所提供。芯片平台是 Oligo 基因芯片，每一个样本扫出来的原始信号经过 GenePix Pro 处理、Lowess 归一化之后，形成一张 21378 × 33 的表达矩阵——行是基因，列是样本。这个矩阵在机器学习里属于标准的“高维小样本”：特征数远超样本数。

直接拿 21378 个基因去训分类器，不是不能跑，而是结果基本不可信。支持向量机、逻辑回归这类模型在高维空间里很容易找到一个完全分开训练集的超平面，但那是在记住噪声而不是学习规律。这也是为什么论文要把特征选择放在模型之前，而不是把全部基因丢进 SVM 里。特征选择解决的不只是计算量问题，更是泛化能力问题。

2.2 缺失值：为什么用最近邻法而不是均值填补

原文说原始数据有 0.6% 的缺失，采用最近邻法填补，k=10。0.6% 听起来不多，但放在 21378 维的基因向量里，意味着矩阵里有相当数量的基因在某几个样本上是空值。而 SAM、PLS 这些方法都要求输入矩阵是完整的，缺失值必须提前处理。

常见做法是用均值或中位数填补，简单但对后续分析会有隐性影响：均值填补会把该样本的表达值向整体中心拉，人为缩小基因在两类样本间的方差，SAM 的差异检验就更容易漏掉真阳性。最近邻法不一样，它用表达模式最接近的 10 个样本来估计缺失值，保留的是局部表达结构。k=10 这个数字不算大，说明作者希望填充值尽量受局部邻居影响；如果 k 太大，填充结果会趋向全局均值，失去意义。

2.3 Lowess归一化的顺序问题

先归一化还是先补缺失，这个顺序不同工具链有不同习惯。论文的流程是：芯片扫描 → 图像信号转数字信号 → Lowess 归一化 → 构建 21378 × 33 表达矩阵，然后才是 KNN 补缺失。这样做的合理性在于：Lowess 处理的是荧光强度依赖的系统误差，如果不先做，不同芯片间的基线可能不在一个量级，KNN 计算样本距离时会被整体偏移干扰。

我自己处理这类数据时会坚持“先归一化、后填补”。先归一化保证所有样本在同一尺度下，再算缺失值才有意义；反过来先补缺失再归一化也能跑，但补缺失时用的距离可能被芯片间批次效应污染。论文里这点没有细节展开，但背后逻辑是清楚的。

2.4 复现前先核对这张表

拿到任何数据集，第一步永远是核对样本标签和缺失情况。尤其要注意论文里提到“两组采用相同的 20 例正常样本作为共同参照”——这 20 例是芯片杂交时的参照池，不是加入模型的正常样本。真正进入机器学习的样本只有 33 个，别把参照样本混进标签里，这个错我在实际项目里见过不止一次。

3. 第一道闸口SAM：用置换检验控制FDR，从21378个基因里筛出1642个差异基因

3.1 SAM为什么比普通t检验更适合同类数据

SAM，全称 Significance Analysis of Microarrays，核心思路是在 t 检验的基础上加了一个稳定的“小常数”。单基因的 t 检验会面临一个问题：表达量本身很低的基因，在两组间只要有微小差异，t 值就会虚高，看起来显著，其实是噪声。SAM 用“差值 ÷ (标准差 + s0)”这种形式替换纯 t 值，s0 的作用是把低表达基因的高方差压住。这样筛选出来的差异基因更倾向于表达量稳定、差异幅度真实的基因。

另一个关键机制是置换检验。只做 33 个样本、21378 个基因的差异检验，每个基因都算一个 p 值，直接按 p<0.05 切会得到大量假阳性。SAM 把样本标签随机打乱，重复 100 次，用随机状态下能产生多少“显著基因”来估计误判率 FDR。论文里 permutation 次数设为 100，在 2009 年的计算条件下是常规选择；现在跑的话可以加倒 1000，FDR 估计会更稳定，但这不影响主线流程。

3.2 delta=0.859和FDR=4.81%怎么配合

SAM 的 delta 是一个滑块，delta 越大，被判定为显著的基因越少，FDR 越小。论文最终取值 delta=0.859，对应 FDR=4.81%。配合的另一个条件是基因表达改变最小倍数设为 2 倍，也就是 fold change ≥ 2。这两个条件同时生效后，从 21378 个基因里筛出 1642 个差异表达基因，其中上调 612 个，下调 1030 个。

实际操作中 delta 不是拍脑袋定的。SAM 软件会输出一张 delta 与 FDR 的对照表，通常的做法是：先看 FDR 在 5% 附近对应的 delta 是多少，再回头看这个 delta 下选出的基因数量是否平衡。如果筛出来一两百个，后面的 PLS 就没有多少压缩空间；如果筛出来上万个，说明 delta 太松。论文卡在 1642 个，数量适中，既过滤了大部分噪声，又保留了足够的候选基因，这个体量对下一步 PLS-VIP 是友好的。

3.3 现在复现SAM的代码路径

当年跑 SAM 用的是独立软件或 Excel 插件，现在复现这个流程，我一般用 R 的 samr 包。核心调用方式如下：

library(samr) # 表达矩阵：行=基因，列=样本；先做quantile归一化 expr <- readRDS("expr_21378x33.rds") # 21378行，33列 X <- normalize.quantiles(as.matrix(expr)) # 标签：前13例弥漫型，后20例肠型 y <- c(rep(1, 13), rep(2, 20)) sam_data <- list( x = X, y = y, geneid = rownames(X), genenames = rownames(X), logged2 = TRUE ) # 非配对两组比较，置换次数按论文设为100 samfit <- SAM( x = X, y = y, resp.type = "Two class unpaired", nperms = 100, logged2 = TRUE ) # delta设为0.859，对应FDR约4.81% sig_table <- samr.compute.sig.table( samfit, delta = 0.859, data = sam_data, R = 100 ) print(sig_table)

代码里的 SAM() 会返回拟合对象，samr.compute.sig.table() 用来在指定 delta 下提取显著基因表。logged2=TRUE 表示表达值已经做过 log2 变换，如果手头是原始荧光强度，需要先变换。normalize.quantiles() 是 preprocessCore 包里的函数，作用是对所有样本做分位数归一化，保证芯片间的表达分布一致。

这里要留意：论文的筛选条件是“2倍表达改变 + FDR<5%”组合，samr 输出的是基于 delta 的显著基因列表，Fold Change 筛选需要额外处理。实操时可以先看 sig_table，再根据基因的平均表达值或样本均值差做一次 FC≥2 的过滤，两边取交集。

3.4 筛完以后怎么检查结果

1642 个基因拿到手，先别急着往下送。我会检查三件事：第一，已知的胃癌相关基因在不在这个集合里，比如 P53、c-erbB2，如果在，说明筛出的基因方向是合理的；第二，上调与下调基因的比例，论文里是 612:1030，下调多于上调，这个比例在肿瘤组织 vs 正常组织里很常见，但如果你自己的数据里全是上调，要警惕可能是归一化或标签方向出了问题；第三，把 1642 个基因的表达谱做一次快速热图，看两组样本是否能大致分开。如果热图上完全分不开，说明 SAM 的 delta 选松了，回去调大 delta。

4. PLS-VIP与BD-SFS接力：候选基因从1642个压到139个，再定成20个

4.1 PLS一次主成分就解释了77.8%的类别信息

1642 个基因对 SVM 来说还是太多。论文的第二道降维是偏最小二乘（PLS），具体用的是 VIP 系数。PLS 和 PCA 的区别在于：PCA 只看自变量 x 的结构，不加区分类别；PLS 在抽取成分时同时考虑 x 与类别 y 的相关性，所以它在监督降维任务里比 PCA 更合适。

论文里 PLS 只取了一个主成分，就能解释 77.8% 的因变量 y 和 64.2% 的自变量信息。这是一个很关键的数字——一个成分已经够了。很多人在复现时习惯把 n_components 调大，觉得成分越多越好，实际上在小样本场景下，后续成分往往是在拟合噪声，反而干扰 VIP 排序。遇到这种情况，我的习惯是先用一个主成分跑一遍，看解释率，如果 y 的解释率已经超过 70%，就没必要加第二个成分。

4.2 VIP>1.5而不是VIP>1，是为了给后面的搜索留余地

VIP>1 一般认为基因对分类有正面贡献。论文在 PLS 取一个主成分时，VIP>1 的基因有 675 个，VIP>2 的只有 9 个，最大 VIP 系数是 2.159。如果卡在 VIP>1，剩下 675 个基因对 BD-SFS 来说还是太多；卡在 VIP>2 又太少，可能丢掉有协同作用的基因。最终选 VIP>1.5，剩下 139 个，这是一个很务实的折中。

这里有一个可以借鉴的判断方式：把 VIP 从高到低排序，看数量分布曲线。如果大部基因集中在 VIP 1~1.5 之间，说明分类信号分散在很多基因上，取 1.5 能保留信号最强的头部，又不会把候选集缩得太死。如果曲线很平滑、没有明显拐点，我会把 1.5 当作起点，向下调或者向上调各试一轮，看最终 SVM 准确率的变化，而不是死守固定阈值。

4.3 BD-SFS：巴氏距离排序加上SVM准确率做栅栏

第三道降维是基于巴氏距离的顺序前向搜索。逻辑分两段：

第一段，对 139 个基因计算巴氏距离（Bhattacharyya distance）。这个指标同时考虑基因在两组样本中的均值和方差，距离越大，说明这个基因单独区分两类样本的能力越强。按距离从大到小排序，得到基因搜索顺序。

第二段，按顺序做前向搜索。第一个基因先进入候选集，用 SVM 分类准确率作为准则 C1；再加入第二个基因，计算 C2。只有 C2 > C1，第二个基因才留下，否则丢弃。以此类推，直到遍历完所有按巴氏距离排序的候选基因。这样选出来的 20 个基因，每一个都对分类准确率有增量贡献，而不是单纯堆特征。

这里实际上是一个 filter 加 wrapper 的混合策略：巴氏距离负责初排序，SVM 准确率负责终审。单靠 filter 会漏掉组合效应强但单基因区分度弱的基因；单靠 wrapper 从 139 个基因里穷举子集计算量又太大。论文把两者接起来，是典型的次优搜索思路，工程上非常务实。

4.4 SVM的RBF核、独立测试集和200次重复

SVM 的核选择用 RBF 径向基核，这是处理非线性可分数据的常用选择。训练集和测试集的划分方式是独立测试集：训练集从弥漫型里随机抽 9 例、肠型里抽 13 例，共 22 例；测试集是剩下的弥漫型 4 例加肠型 7 例，共 11 例。

33 例样本的随机划分会有相当大的偶然性，论文的处理方式是重复 200 次，每次重新随机抽样，最后取 200 次选出的特征基因的并集作为最终特征基因集，取平均分类准确率作为最终指标，得到 89.45%。这个做法的关键点在于“并集”——不是取 200 次里出现次数最多的 20 个，而是把所有被选过的基因合并。这样得到的是分类器依赖的完整基因集合，而不是某一次随机划分下的偶然结果。

4.5 复现PLS-VIP与SFS的代码骨架

PLS 部分可以用 scikit-learn 的 PLSRegression，VIP 系数需要手工算：

import numpy as np from sklearn.cross_decomposition import PLSRegression from sklearn.svm import SVC # 输入：SAM筛出的 1642 x 33 表达矩阵 X = expr_1642.T # 33 行样本，1642 列基因 y = np.array([0]*13 + [1]*20) # PLS 只取 1 个主成分 pls = PLSRegression(n_components=1, scale=True) pls.fit(X, y) # 手工计算 VIP 系数 def vip_score(pls, X, y): t = pls.x_scores_ # 样本得分 w = pls.x_weights_ # 基因权重 q = pls.y_loadings_ # y 载荷 p = X.shape[1] # 基因数 ssy = np.sum(y**2) ssy_h = np.sum(q**2) * np.sum(t**2) vips = np.sqrt(p * ssy_h * w.ravel()**2 / ssy) return vips vip = vip_score(pls, X, y) keep_idx = np.where(vip > 1.5)[0] # 139 个基因

BD-SFS 的搜索循环骨架如下：

# 对 139 个基因按巴氏距离降序排列 genes_sorted = sorted(genes_139, key=bhattacharyya_distance, reverse=True) selected = [] best_acc = 0.0 for g in genes_sorted: trial = selected + [g] acc = evaluate_svm(X[:, trial], y) # RBF核，独立22/11划分 if acc > best_acc: selected = trial best_acc = acc

evaluate_svm() 里可以封装 SVM 分类流程：训练集 22 例、测试集 11 例，RBF 核。注意这一步的 SVM 参数不能直接用默认值，gamma 和 C 要通过小范围网格搜索确定。论文没有给出具体数值，复现时常见做法是在训练集内部用交叉验证选参，再在独立测试集上评估。

5. 实战避坑：小样本特征选择最容易翻车的五个细节

5.1 坑一：把SAM和PLS跑完才切训练/测试集

现象：复现流程时先把全部 33 个样本放进 SAM 和 PLS 做特征选择，选完 20 个基因后再切训练集和测试集。结果分类准确率奇高，接近 98%，但换到验证集或新数据上直接崩。

原因：提前用全量数据做特征选择，测试集的信息已经通过特征选择过程泄露进了模型。SVM 看到的特征是在包含测试样本的情况下选出来的，相当于考试前先看了答案。

解决：严格按论文的流程来——特征选择的依据是 SAM、PLS 的统计量和巴氏距离，这些在数据生成后是固定的，不会因为划分方式改变。但 SVM 的参数（RBF 的 gamma、C）必须在训练集内部调，独立测试集只能最后碰一次。如果要用交叉验证做整体评估，特征选择要嵌入每一折内部，重新在训练折上跑一遍 SAM+PLS+SFS，再评估测试折。这一步会慢很多，但结果可信。

5.2 坑二：delta和VIP阈值来回试，试出一个漂亮数就拍板

现象：调 delta 从 0.8 到 0.9，VIP 阈值从 1.2 到 1.8，最后选了一组让 SVM 准确率达到 91% 的参数，比论文的 89.45% 还高，觉得复现成功了。

原因：这是典型的过拟合到超参数。小样本场景下，阈值微调就能显著影响入选基因，而 SVM 在小特征集上准确率波动很大。高准确率很可能只是阈值组合碰巧适合当前的 33 个样本。

解决：每个阈值组合下重复 200 次随机抽样，观察平均准确率的方差。论文的关键不是那个 89.45% 数字，而是 200 次重复取平均的做法。我一般会记录不同阈值组合下准确率的均值和标准差，选均值高且标准差小的组合，而不是只看单次结果。如果某个阈值下准确率均值 91% 但标准差 8%，另一个阈值下 89% 但标准差 3%，后者显然更稳。

5.3 坑三：表达矩阵行和列放反，SAM和sklearn的输入全乱套

现象：R 的 samr 包报错说样本数与标签长度不一致，或者 PLS 跑完发现 VIP 值长得像样本数而不是基因数。

原因：SAM 要求矩阵行是基因、列是样本；scikit-learn 要求 X 矩阵行是样本、列是特征。两套库的矩阵方向正好相反，中间转置漏了一次，后续全废。

解决：在代码开头写死接口函数，专门做格式转换。R 侧按“行基因、列样本”组织，转给 Python 前用 t() 转置，并在变量名里标注清楚，比如 X_samples_genes 表示行样本、列基因；X_genes_samples 表示行基因、列样本。每次调用模型前检查 X.shape 的第一个维度是否等于样本数 33，第二个维度是否等于当前基因数。这一步能省掉 80% 的对接事故。

5.4 坑四：SVM的RBF核参数没调，直接默认参数跑SFS

现象：按论文流程跑完，SVM 准确率只有 60% 多，或者 SFS 在第一轮就停住，只选出两三个基因。

原因：RBF 核的 gamma 参数决定单个样本的影响半径，sklearn 默认 gamma 是 1/n_features。当特征维度从 139 降到 20 时，默认 gamma 变化很大，不加调整直接跑 SFS，SVM 的分类行为完全不符合预期。

解决：在 SFS 内部包一个参数搜索。候选基因数量变化时，gamma 每次都要重调。常见做法是用训练集做 5 折交叉验证，在 2^(-5) 到 2^5 的指数网格里搜 gamma 和 C。SFS 循环会很慢，可以把参数网格变粗，只搜几个典型值，比如 C ∈ {0.1, 1, 10}，gamma ∈ {0.001, 0.01, 0.1}。对于小样本场景，粗网格往往比精网格更稳定，细调的参数绑在 33 个样本上没有意义。

5.5 坑五：200次随机抽样选出的基因集不稳定

现象：重复跑 200 次，每次选出的基因集合不一样，两次试验之间重合的基因只有 10~12 个。担心是算法不稳定，怀疑复现错了。

原因：样本量太小，随机划分训练/测试集对 SFS 的搜索路径影响很大。某些基因在一种划分下有用，在另一种划分下没增量贡献，这是高维小样本的固有属性，不是代码 bug。

解决：论文用并集是有道理的——取 200 次选出的所有基因，保留的是任何一次划分下都有用的基因，交集反而会丢掉互补基因。实际操作中我可以进一步：记录每个基因在 200 次里被选中的频次，按频次排序列出 TOP 30，观察哪几个基因上榜率超过 80%。那些高频率基因才是真正稳定的分类信号。如果最终结论需要提供给生物学验证，建议从高频率基因里挑核心候选集，用 KEGG 或功能注释去补证据。

6. 把20个特征基因送到生物学验证：三种值得照搬的验证路径

6.1 路径一：层次聚类看无监督表现

论文用 Cluster3.0 做层次聚类，基于 20 个特征基因的表达数据，33 例样本中只有 2 例肠型胃癌错分到弥漫型一侧，聚类准确率 93.94%。

用 Python 复现这个验证不复杂：

from scipy.cluster.hierarchy import linkage, dendrogram from scipy.spatial.distance import pdist # expr20：20个特征基因 x 33个样本 Z = linkage(pdist(expr20.T, metric="correlation"), method="average") dendrogram(Z, labels=sample_names)

correlation 距离对应 Cluster3.0 里常用的相似性度量，average 对应平均连接法。看树状图时重点不是“分成几簇”，而是两簇的标签构成是否与 Lauren 分型一致。如果样本没有预先标记，单靠这 20 个基因也能把绝大多数样本分成两组，说明这些基因本身携带分类结构——这是对 SVM 结果的重要交叉验证，不依赖监督信息。

6.2 路径二：通路富集时注意基因名映射

论文用 KEGG、GenMAPP 等数据库分析，发现 20 个基因中有 8 个参与了 32 条信号转导通路的节点，典型代表是 WNT16、MLL3、LAMA2、CPB2。现在复现这条路径，我一般用 R 的 clusterProfiler，但要注意基因名映射问题。

Oligo 芯片时代的探针注释是以当时的 RefSeq 或 Unigene ID 为准，映射到今天的标准基因名时经常发生断链。比如论文里的 MLL3，在 HGNC 的标准命名里已经改成了 KMT2C，如果你拿着旧名字去查注释，很可能查不到。做通路富集前，先用最新注释文件把 20 个基因统一映射到 GeneSymbol，再跑 KEGG 富集，否则分析结果会打对折。

映射时还要检查反义链和内参基因。有些探针同时比对到多个转录本，富集分析会把同一基因的不同转录本算成多个条目，导致通路信号虚高。

6.3 路径三：三张表交叉验证才算闭环

我复现这类论文，习惯把结果整理成三张表交叉核对。第一张是 SVM 分类表现：平均准确率、每次划分的方差、20 个基因在训练集上的独立表现；第二张是层次聚类结果与 SVM 结果的标签一致性；第三张是基因集合在 KEGG、GenMAPP 等通路数据库命中情况。

如果 SVM 准确率低于论文的 89.45%，差距不大且标准差也在合理范围，说明特征基因是可复现的；如果聚类准确率从 93.94% 掉到 70% 以下，说明选出的基因可能没有独立的分类结构，要从 SAM 的 delta 和 PLS 的 VIP 阈值往回查；如果富集通路全是代谢大类、没有肿瘤相关通路，重点检查基因名映射是否出错。

从那以后，我每次做完一轮特征筛选，都会强制走一遍这三条验证路径：先看无监督聚类是否稳，再做基因名映射和通路富集，最后把三张表放在一起判断这组基因值不值得往实验方向推进。希望帮到你。

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