简介这份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% 的缺失采用最近邻法填补k10。0.6% 听起来不多但放在 21378 维的基因向量里意味着矩阵里有相当数量的基因在某几个样本上是空值。而 SAM、PLS 这些方法都要求输入矩阵是完整的缺失值必须提前处理。常见做法是用均值或中位数填补简单但对后续分析会有隐性影响均值填补会把该样本的表达值向整体中心拉人为缩小基因在两类样本间的方差SAM 的差异检验就更容易漏掉真阳性。最近邻法不一样它用表达模式最接近的 10 个样本来估计缺失值保留的是局部表达结构。k10 这个数字不算大说明作者希望填充值尽量受局部邻居影响如果 k 太大填充结果会趋向全局均值失去意义。2.3 Lowess归一化的顺序问题先归一化还是先补缺失这个顺序不同工具链有不同习惯。论文的流程是芯片扫描 → 图像信号转数字信号 → Lowess 归一化 → 构建 21378 × 33 表达矩阵然后才是 KNN 补缺失。这样做的合理性在于Lowess 处理的是荧光强度依赖的系统误差如果不先做不同芯片间的基线可能不在一个量级KNN 计算样本距离时会被整体偏移干扰。我自己处理这类数据时会坚持“先归一化、后填补”。先归一化保证所有样本在同一尺度下再算缺失值才有意义反过来先补缺失再归一化也能跑但补缺失时用的距离可能被芯片间批次效应污染。论文里这点没有细节展开但背后逻辑是清楚的。2.4 复现前先核对这张表数据项数值/方式样本总数33 例弥漫型 13 肠型 20基因向量21378 个缺失比例约 0.6%缺失填补最近邻法k10归一化Lowess拿到任何数据集第一步永远是核对样本标签和缺失情况。尤其要注意论文里提到“两组采用相同的 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 值直接按 p0.05 切会得到大量假阳性。SAM 把样本标签随机打乱重复 100 次用随机状态下能产生多少“显著基因”来估计误判率 FDR。论文里 permutation 次数设为 100在 2009 年的计算条件下是常规选择现在跑的话可以加倒 1000FDR 估计会更稳定但这不影响主线流程。3.2 delta0.859和FDR4.81%怎么配合SAM 的 delta 是一个滑块delta 越大被判定为显著的基因越少FDR 越小。论文最终取值 delta0.859对应 FDR4.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 下提取显著基因表。logged2TRUE 表示表达值已经做过 log2 变换如果手头是原始荧光强度需要先变换。normalize.quantiles() 是 preprocessCore 包里的函数作用是对所有样本做分位数归一化保证芯片间的表达分布一致。这里要留意论文的筛选条件是“2倍表达改变 FDR5%”组合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 VIP1.5而不是VIP1是为了给后面的搜索留余地VIP1 一般认为基因对分类有正面贡献。论文在 PLS 取一个主成分时VIP1 的基因有 675 个VIP2 的只有 9 个最大 VIP 系数是 2.159。如果卡在 VIP1剩下 675 个基因对 BD-SFS 来说还是太多卡在 VIP2 又太少可能丢掉有协同作用的基因。最终选 VIP1.5剩下 139 个这是一个很务实的折中。VIP 阈值基因数量VIP 1675VIP 1.5139论文采用VIP 29最大 VIP 值2.159这里有一个可以借鉴的判断方式把 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 例。样本分组弥漫型肠型合计训练集91322测试集471133 例样本的随机划分会有相当大的偶然性论文的处理方式是重复 200 次每次重新随机抽样最后取 200 次选出的特征基因的并集作为最终特征基因集取平均分类准确率作为最终指标得到 89.45%。这个做法的关键点在于“并集”——不是取 200 次里出现次数最多的 20 个而是把所有被选过的基因合并。这样得到的是分类器依赖的完整基因集合而不是某一次随机划分下的偶然结果。4.5 复现PLS-VIP与SFS的代码骨架PLS 部分可以用 scikit-learn 的 PLSRegressionVIP 系数需要手工算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_components1, scaleTrue) 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, keybhattacharyya_distance, reverseTrue) 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 accevaluate_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必须在训练集内部调独立测试集只能最后碰一次。如果要用交叉验证做整体评估特征选择要嵌入每一折内部重新在训练折上跑一遍 SAMPLSSFS再评估测试折。这一步会慢很多但结果可信。5.2 坑二delta和VIP阈值来回试试出一个漂亮数就拍板现象调 delta 从 0.8 到 0.9VIP 阈值从 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 变化很大不加调整直接跑 SFSSVM 的分类行为完全不符合预期。解决在 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 # expr2020个特征基因 x 33个样本 Z linkage(pdist(expr20.T, metriccorrelation), methodaverage) dendrogram(Z, labelssample_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 阈值往回查如果富集通路全是代谢大类、没有肿瘤相关通路重点检查基因名映射是否出错。从那以后我每次做完一轮特征筛选都会强制走一遍这三条验证路径先看无监督聚类是否稳再做基因名映射和通路富集最后把三张表放在一起判断这组基因值不值得往实验方向推进。希望帮到你。本文还有配套的精品资源点击获取