企业官网建设流程全解析

简介：这是一份整合生物信息学与机器学习方法、预测NAT10相关下游基因的完整项目资源，面向生物信息学研究者、机器学习初学者及关注转录调控机制的科研人员。资源包含从GEO数据获取、清洗、标准化到特征构建的整套预处理流程，并实现SVM、随机森林等分类器训练与交叉验证，同时提供PCA降维、热图绘制、基因互作网络可视化及limma差异分析结果，可直接用于复现和扩展NAT10互作基因的筛选工作。压缩包整体31.5MB，共55个文件，以py脚本、csv表格、txt说明为主，辅以R脚本、xls数据及png结果图，便于对照代码与结论进行学习。目前已有53人浏览学习，适合希望系统掌握多组学数据挖掘与预测建模流程的读者。

1. 用生物信息学与机器学习预测 NAT10 下游基因：这套资源到底能帮你省多少事

NAT10 是核仁里的一个多功能蛋白，参与 rRNA 加工、转录调控和染色质修饰，偏偏它调控哪些下游基因这件事，实验做起来又贵又慢。这个资源包就是用生物信息学与机器学习把这件事的预测环节整个跑通：从 GEO 原始数据下载，到 limma 差异分析、PCA 降维、相关性筛选，再到 SVM 与随机森林建模，最后输出基因互作网络图和候选基因列表。适合正在做基因调控网络、想快速拿到一组可信候选基因做实验验证的从业者，也适合要把这套流程挪到其他基因上的生信初学者。我拆完之后最大的感受是：它不是一个教学 demo，而是一套真跑过、有完整中间产物的实战工程。

2. 数据预处理：把三个 GEO 系列合并成一个表达矩阵，批次效应是第一道坎

2.1 解压之后先看清流水线：从 raw 到 processed 再到 result

这套资源解压后，目录结构其实已经把工作流写脸上了。datasets/raw是下载下来的原始矩阵，datasets/processed是清洗后的中间产物，result里是所有图表和 CSV 输出。脚本分两类：数据合并类（extractGSE206917Data.py、changeGSE82139.py、mergeGSE207002Data.py、mergeDatasets.py）和建模可视化类（PCA.py、SVM.py、randomForest.py、plotNetwork.py等）。

我建议拿到包的第一件事不是急着跑脚本，而是先按这个顺序捋一遍数据流：

# 按依赖顺序执行数据预处理脚本 python extractGSE206917Data.py python changeGSE82139.py python mergeGSE207002Data.py python mergeDatasets.py

第一条命令把 GSE206917 的原始数据抽取成统一格式的矩阵；第二条做 GSE82139 的基因名映射；第三条处理 GSE207002；最后一条把三个矩阵合并。mergeDatasets.py是这条链路的汇合点，它输出的矩阵就是后续所有分析的输入。

合并矩阵时，三个数据集来自不同平台、不同实验室，直接concat一定会出问题。常见做法是先做基因 symbol 统一，再取交集。核心逻辑大概是这样的：

# mergeDatasets.py 的核心逻辑（简化版） import pandas as pd exp1 = pd.read_csv("datasets/processed/GSE206917_symbol.csv", index_col=0) exp2 = pd.read_csv("datasets/processed/GSE82139_symbol.csv", index_col=0) exp3 = pd.read_csv("datasets/processed/GSE207002_symbol.csv", index_col=0) # 取三个矩阵共有的基因集合，避免大量 NaN common_genes = set(exp1.index) & set(exp2.index) & set(exp3.index) common_genes = sorted(common_genes) merged = pd.concat([exp1.loc[common_genes], exp2.loc[common_genes], exp3.loc[common_genes]], axis=1) merged.to_csv("datasets/processed/merged_expression_matrix.csv")

这里的重点是set(exp1.index) & set(exp2.index) & set(exp3.index)。三个数据集里同一个基因的探针名可能完全不同，直接按行拼接会让缺失值爆炸。取交集是最保守的策略，代价是会丢掉部分只在单个数据集里检测到的基因，但对后续机器学习建模来说，一个没有 NaN 的矩阵远比一个基因数全但有大量空洞的矩阵靠谱。

2.2 探针映射与基因 symbol 统一：changeGSE82139.py 到底在干什么

changeGSE82139.py这个脚本名字起得很直白，就是对 GSE82139 做基因名替换。GEO 下载的矩阵通常以探针 ID 为行名，比如AFFX-BioB-5_at这种，而你要做的是把所有平台的探针 ID 统一成ENSG或 gene symbol。

这里有个关键选择：统一成什么格式。如果后续要用limma做差异分析，gene symbol 更直观；如果要对接富集分析，ENSEMBL更好。这套资源里convertGSE82139.R出现得很及时，它多半是用biomaRt或者平台注释包做转换的：

# convertGSE82139.R 的典型写法 library(biomaRt) ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl") probe2gene <- getBM(attributes = c("affy_hg_u133_plus_2", "hgnc_symbol"), filters = "affy_hg_u133_plus_2", values = rownames(expr_matrix), mart = ensembl)

转换完成后，多个探针对应同一个基因的情况很常见。我的建议是保留表达量最高的那个探针，或者直接取平均值，但千万不要直接drop_duplicates留第一个，因为第一个往往不是表达最显著的。这一步做不好，后面差异分析的结果会非常难看。

2.3 批次效应处理：热血上头直接合并的下场就是 PCA 图按数据集分群

三个 GEO 系列合并后，样本间最大的差异往往不是生物学分组，而是来自哪个数据集。这就是批次效应。如果不处理，后续 PCA 的第一主成分会忠实地把三个数据集分成三坨，而你要找的 NAT10 相关信号全被淹没在批次噪声里。

处理批次效应的常见做法有两种：一是用limma包的removeBatchEffect，需要你知道每个样本来自哪个数据集；二是用sva包的ComBat，它能自动估计批次。我个人在这类场景下优先选ComBat，因为它对批次间方差齐性的假设更宽松，对小样本数据更稳。

# 用 sva 包做 ComBat 批次矫正 library(sva) batch <- c(rep(1, ncol(exp1)), rep(2, ncol(exp2)), rep(3, ncol(exp3))) mod <- model.matrix(~group, data = meta) corrected <- ComBat(dat = merged_matrix, batch = batch, mod = mod)

参数上，batch向量的长度必须和矩阵列数完全对齐，这是我见过最多的报错来源。另外mod里放的是你真正关心的分组变量，目的是在矫正批次的同时保留生物学差异。如果mod里放了 NAT10 的表达量分组，那矫正后的矩阵就已经带上了你要研究的信号。

数据预处理做完，输出物应该是一个行是基因、列是样本、批次效应已矫正的表达矩阵。判断这一步是否做好的标准很简单：画出 PCA 图，看样本是否按生物学分组聚集，而不是按数据集聚集。

3. 差异表达与 PCA 降维：先筛到几十个候选基因再进机器学习

3.1 用 limma 做差异表达分析：阈值的设定决定候选基因数量

DEA.R是这个资源里最标准的差异分析脚本，用的是limma的经典贝叶斯流程。它的输出all.limmaOut.csv是所有基因的完整结果表，significant_genes_from_limma.csv是经过阈值过滤后的显著基因列表。

# DEA.R 的核心流程 library(limma) design <- model.matrix(~0 + factor(group), data = meta) colnames(design) <- levels(factor(meta$group)) contrast.matrix <- makeContrasts(high_vs_low = group_high - group_low, levels = design) fit <- lmFit(expr_matrix, design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) all_results <- topTable(fit2, coef = "high_vs_low", number = Inf, sort.by = "none")

topTable里的number = Inf表示输出全部基因，sort.by = "none"保持原始顺序。拿到all.limmaOut.csv之后，筛选阈值一般这么定：log2FC绝对值大于 1，adj.P.Val小于 0.05。注意这里用的是调整后 p 值，不是原始 p 值，因为差异分析做了几千次检验，不用 FDR 校正的话假阳性会失控。

这里有个血泪经验：adj.P.Val的阈值直接决定了你后面机器学习输入特征的数量。定 0.05 可能筛出几百个基因，定 0.01 可能只剩几十个。我的习惯是先看all.limmaOut.csv里adj.P.Val的分布，找一个自然拐点，而不是死磕某个固定值。

3.2 PCA 降维与载荷解读：pca_with_NAT10.png 在看什么

PCA.py做的事情是对表达矩阵做主成分分析，输出两个文件：pca_with_NAT10.png和pca_loadings.csv。前者是样本在主成分空间的分布图，后者是每个基因在每个主成分上的载荷，也就是贡献度。

# PCA.py 的核心逻辑 from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(expr_matrix.T) # 样本为行，基因为列 pca = PCA(n_components=5) scores = pca.fit_transform(X_scaled) loadings = pd.DataFrame(pca.components_.T, index=expr_matrix.index, columns=[f"PC{i+1}" for i in range(5)]) loadings.to_csv("result/pca_loadings.csv")

这里为什么要对矩阵转置再标准化？因为 sklearn 的PCA默认样本是行，而表达矩阵的常规格式是基因在行、样本在列，不转置的话算出来的不是样本主成分而是基因主成分。StandardScaler这一步也很关键——表达量范围从几十到几万，不标准化的话，高表达基因会主导主成分方向，低表达但可能有生物学意义的基因直接被淹没。

在pca_with_NAT10.png这张图上，我一般会做两件事：一是看样本是否按高低表达 NAT10 分组分开，二是看 NAT10 本身在载荷图中的位置。如果 NAT10 在 PC1 或 PC2 的载荷绝对值很大，说明这个基因对样本分群的贡献显著，后续把它当作核心特征是有依据的。

3.3 相关性筛选：NAT10_high_correlation_genes.csv 的由来

除了差异表达，这套资源还走了一条平行路径——直接计算每个基因与 NAT10 表达量的相关性。这个思路很朴素：如果一个基因的表达趋势跟着 NAT10 走，不管正相关还是负相关，它都可能是 NAT10 的下游靶点或共调控基因。

# 计算每个基因与 NAT10 的 Pearson 相关系数 nat10_expr = expr_matrix.loc["NAT10"] correlations = expr_matrix.apply(lambda row: row.corr(nat10_expr), axis=1) correlations = correlations.sort_values(ascending=False) # 取相关系数绝对值大于阈值的基因 selected = correlations[abs(correlations) > 0.6] selected.to_csv("result/NAT10_high_correlation_genes.csv")

阈值定多少要看数据分布。0.6 是相对严格的标准，如果筛完一个都没有，就放松到 0.5，再不行就 0.4。我在这类项目里一般会画一个相关系数直方图看一眼分布再定阈值，而不是直接拍一个数。

到这里，资源包里已经有了两个维度的候选基因列表：significant_genes_from_limma.csv是差异表达显著基因，NAT10_high_correlation_genes.csv是与 NAT10 表达高度相关的基因。这两组基因取交集，就得到了selected_genes_common_to_both_criteria.csv——这是后续机器学习建模的输入特征。取交集的好处是同时满足"在高低表达组间有差异"和"与 NAT10 表达趋势一致"两个条件，比单一准则稳健得多。

4. 机器学习建模：SVM 与随机森林的选型、参数与结果合并

4.1 为什么是 SVM 和随机森林，而不是深度学习

NAT10 下游基因预测这个任务有几个鲜明特征：样本量小（几十到一两百）、特征维度高（候选基因通常几十到几百）、标签是二分类（NAT10 高表达 vs 低表达，或者患病 vs 对照）。这种小而高维的数据，深度学习极易过拟合，你辛辛苦苦搭的网络大概率在训练集上准确率 99%，验证集上打回原形。

SVM 和随机森林在这类任务里是经过大量验证的经典选择。SVM 擅长处理高维小样本，RBF 核可以把特征映射到更高维空间找分界面；随机森林天然带特征重要性输出，对异常值和噪声的鲁棒性好，还能告诉你哪些基因在分类决策中最关键。这两个模型互补性很强：SVM 看整体分类边界，随机森林看单个特征的贡献，最后把两者的特征重要性合并，比单模型可信得多。

4.2 randomForest.py 参数解读：n_estimators 和 max_features 是主要旋钮

randomForest.py的核心训练代码是标准的 sklearn 写法：

# randomForest.py 核心训练与特征重要性提取 from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score, StratifiedKFold X = feature_matrix # 行是样本，列是 selected_genes_common_to_both_criteria.csv 里的基因 y = labels # NAT10 高表达 vs 低表达 rf = RandomForestClassifier( n_estimators=500, max_features="sqrt", max_depth=10, min_samples_leaf=2, random_state=42, n_jobs=-1 ) cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) scores = cross_val_score(rf, X, y, cv=cv, scoring="accuracy") print(f"CV accuracy: {scores.mean():.3f} ± {scores.std():.3f}") # 特征重要性直接取模型的属性 importance = pd.DataFrame({ "gene": X.columns, "importance": rf.feature_importances_ }).sort_values("importance", ascending=False) importance.to_csv("result/feature_importances.csv", index=False)

几个参数在实际跑的时候最影响结果。n_estimators是树的数量，500 棵树在这个规模的数据上有富余，再往上加收益很小但训练时间线性增加。max_features="sqrt"是分类问题的默认推荐值，它让每棵树只随机看一部分特征，增加树之间的多样性，降低过拟合。max_depth=10是我在类似项目里的习惯值，数据集小、特征多，限制树深能有效防止单棵树学过头。min_samples_leaf=2保证叶子节点至少有两个样本，进一步抑制过拟合。

交叉验证用StratifiedKFold而不是普通KFold，是因为二分类标签可能不平衡，分层抽样能保证每一折里两类样本的比例和整体一致。random_state=42固定住随机种子，让结果可复现——这是机器学习项目里最容易忽略的一步，不设随机种子的话，每次跑出来的特征重要性排名都不同，你根本没法判断哪些基因是稳定的。

4.3 SVM.py 参数解读：C 和 gamma 要网格搜索，别用默认值

SVM.py用的是 RBF 核的支持向量机。RBF 核有两个核心参数：C是误分类惩罚系数，越大越不容忍错误分类，越小越追求间隔最大化；gamma决定单个样本的影响力半径，越大决策边界越复杂。

# SVM.py 核心训练与网格搜索 from sklearn.svm import SVC from sklearn.model_selection import GridSearchCV param_grid = { "C": [0.1, 1, 10, 100], "gamma": [0.001, 0.01, 0.1, 1] } svm = SVC(kernel="rbf", probability=True, random_state=42) grid = GridSearchCV(svm, param_grid, cv=5, scoring="accuracy", n_jobs=-1) grid.fit(X, y) best_svm = grid.best_estimator_ print(f"Best params: {grid.best_params_}") # sklearn 的 SVC 没有直接的特征重要性，用系数或 permutation importance from sklearn.inspection import permutation_importance perm = permutation_importance(best_svm, X, y, n_repeats=10, random_state=42)

SVM 的特征重要性是个玄学问题：线性核可以直接取coef_的绝对值，RBF 核没有天然的特征重要性，只能用permutation_importance做置换检验。它的原理是把一个特征的值随机打乱，看模型准确率掉多少，掉得越多说明这个特征越重要。svm_feature_importances.csv就是这一步的输出。

网格搜索的param_grid里C和gamma各取 4 个值，共 16 组参数组合。这个范围不一定够，我一般会先跑一遍看最优参数落在网格边缘还是内部，如果落在边缘（比如C=100最优），就把网格往那个方向扩展再搜一遍。直接拿默认参数跑 SVM 在这个场景下基本等于摆烂，RBF 核的默认gamma是1/n_features，特征多的时候决策边界会过于复杂，几乎必然过拟合。

4.4 mergeMLResults.py：把两个模型的结果合起来才有说服力

mergeMLResults.py做的事情是把feature_importances.csv和svm_feature_importances.csv合到一起，再结合前面的差异表达和相关性的分数，生成combined_gene_scores.csv。这个综合分数大概的逻辑是把每个基因在不同维度的排名归一化后加权求和。

# mergeMLResults.py 的合并思路 rf_imp = pd.read_csv("result/feature_importances.csv") svm_imp = pd.read_csv("result/svm_feature_importances.csv") limma_res = pd.read_csv("result/significant_genes_from_limma.csv") corr_res = pd.read_csv("result/NAT10_high_correlation_genes.csv") # 每个维度都转成排名，然后归一化到 0-1 def rank_to_score(series): return series.rank(pct=True) combined = pd.DataFrame({"gene": rf_imp["gene"]}) combined["rf_score"] = rank_to_score(rf_imp["importance"]) combined["svm_score"] = rank_to_score(svm_imp["importance"]) combined["limma_score"] = rank_to_score(limma_res.set_index("gene")["adj.P.Val"], ascending=False) combined["corr_score"] = rank_to_score(corr_res.set_index("gene")["correlation"].abs()) # 最终分数取加权平均 combined["final_score"] = combined[["rf_score", "svm_score", "limma_score", "corr_score"]].mean(axis=1) combined = combined.sort_values("final_score", ascending=False) combined.to_csv("result/combined_gene_scores.csv", index=False)

合并排名的思路我觉得是这个资源里最值得学的部分。不同模型的分数量纲不同，随机森林的重要性是 0 到 1 的概率值，limma 的 p 值是从 0 到 1 但越小越好，相关系数是 -1 到 1，直接相加等于把不同尺度的东西硬凑在一起。转成百分位排名后，每个维度都用相对位置说话，这才有可比性。

最终的selected_genes_common_to_both_criteria.csv是交集准则的产物，top_genes_in_selected_criteria.csv是按照综合分数排在前面的基因。我个人更信赖后者，因为它在多维度之间做了权衡，而不是只看单一准则。拿到这个列表之后，才算真正完成了"预测 NAT10 相关下游基因"这个目标的第一阶段。

5. 避坑与排查：跑这套 NAT10 预测脚本时最常翻车的五个点

5.1 探针注释后基因匹配率不到一半

现象：changeGSE82139.py跑完后，处理好的矩阵里基因数量从原始探针数几万个骤减到几千个，和另外两个数据集取交集后只剩一两千个基因。

原因：最常见的是平台注释版本不一致。GSE82139 的平台可能是 GPL570，用的注释文件基于旧版 RefSeq，很多探针在当前版本的biomaRt里已经查不到对应基因名。另外，部分探针本身设计时就落在非编码区或基因间区，永远不可能注释到 gene symbol。

解决：不要死磕 100% 注释率。先查 GPL 平台的软注释文件（GPL 表格），找到对应探针与基因的映射关系。对多探针映射同一基因的情况，按表达量取最大值或中位数；对注释不到的探针直接丢弃。匹配率在 60% 以上就属于正常水平，非要追求 90% 以上反而可能引入错误的映射。

5.2 PCA 图上样本不按分组分群，而是按数据集分群

现象：pca_with_NAT10.png画出来，样本点聚成三四坨，仔细一看每一坨恰好对应一个 GEO 系列，而不是按照 NAT10 高表达和低表达分开。

原因：批次效应没处理。三个数据集来自不同实验室、不同测序批次，技术差异远大于你要找的生物学差异。mergeDatasets.py如果只做了简单的矩阵拼接，PCA 第一主成分承载的主要是批次信息。

解决：回到预处理阶段，用ComBat或removeBatchEffect做批次矫正。注意ComBat的mod参数要包含分组信息，否则它会连你要找的生物学差异一起抹掉。做完之后重新跑一遍PCA.py，看样本是否按分组分群。这一步是整套流程里最容易返工的地方，建议在数据合并完成后第一时间画图确认，不要等跑完机器学习再回头查。

5.3 SVM 准确率接近 100%，反而慌了

现象：SVM.py网格搜索后，交叉验证准确率高达 0.99 甚至 1.0，看着非常漂亮。

原因：极大概率是数据泄露。最常见的有两种：一是特征选择在交叉验证之前完成，也就是说你在全部样本上先筛了显著基因，再用这些基因去做交叉验证，导致验证集的信息提前泄露给了训练过程；二是样本标签和特征矩阵对不齐，比如按基因排序后索引错位，模型实际学到了样本批次信息而不是生物学信号。

解决：检查特征选择是否在交叉验证的循环内部。正确做法是每一折交叉验证里，先在该折的训练集上做差异分析或相关性筛选，再对验证集应用同样的筛选标准。如果资源里的脚本没这么做，你需要手动改一下流程。另一个快速检查：把 N 个随机基因放进模型，看准确率是不是依然很高。如果随机基因也能达到 0.9 以上，几乎可以断定是数据泄露或标签错位。

5.4 randomForest.py 报错特征维度不一致

现象：randomForest.py运行时抛出ValueError: Number of features of the model must match the input，或者训练完成后特征重要性列表和基因名单对不上。

原因：selected_genes_common_to_both_criteria.csv里的基因集合和表达矩阵的行名没有对齐。常见情况是合并矩阵时取了三数据集交集，但后面筛选时混入了不在交集里的基因；或者 CSV 文件里有重复基因名，导致merge时行数不对。

解决：在进入建模前加一个校验步骤。用set比较一下特征基因列表和矩阵行名的差异，打印出缺失的基因名前 20 个，一般一眼就能看出问题。另外，处理重复基因名时用df.index.is_unique检查一下，不唯一就按表达量去重。

5.5 热图和网络图的输出不符合预期：要么全一色，要么全是孤立点

现象：plotHeapMap.py画出的热图颜色几乎没有变化，所有格子挤在一个色阶区间里；plotNetwork.py画出的网络图只有孤立的几个点，没有连边。

原因：热图颜色单调通常是数据没做标准化，表达量分布偏态严重，大部分数值集中在低区间；网络图全是孤立点通常是因为相关性阈值设得太高，在我见过的案例里，基因表达相关系数绝对值能超过 0.8 的本来就屈指可数，阈值定 0.8 以上等于直接把连边全砍掉。

解决：热图画之前先对表达矩阵做 z-score 标准化或者取 log2 变换，把分布拉正再画；网络图先看相关系数的整体分布，至少让排名前 5% 的基因连边数不为零。顺便提一句，这个脚本文件名plotHeapMap.py拼错了，正确拼写是 Heatmap，但这不影响运行，只是读代码的时候别被带偏。

6. 从候选基因列表到一张能讲故事的图：网络可视化和下游验证的落地技巧

机器学习模型输出的是一串基因名和分数排名的 CSV，但要让别人信服，你需要一张直观的图。plotNetwork.py做的事情是读入combined_gene_scores.csv和相关性矩阵，以基因为节点、相关性为连边、综合分数为节点大小，画出gene_interaction_network_based_on_scores.png。跑这个脚本之前，我的习惯是先把连边阈值确定下来：排序后取相关系数绝对值前 5% 的基因对作为连边。阈值太低图会乱成一团毫无信息量，阈值太高全是孤立点，这个 5% 的经验值在大多数表达谱数据上都适用。

拿到候选基因列表后，下一步的落脚点是生物学验证。把selected_genes_common_to_both_criteria.csv里的基因符号整理成一份 txt，上传到 DAVID 或 Enrichr 这类在线富集工具做 GO 和 KEGG 富集分析，重点关注核仁、rRNA 加工、RNA 结合这些 NAT10 已知功能相关的通路。如果富集到的通路恰好落在这些类别里，说明你的模型预测结果和已知生物学一致，可信度大幅提升。实验中可以用免疫共沉淀或 ChIP-seq 验证模型预测的基因确实与 NAT10 有物理互作。

最后说一个我自己的强制习惯：在拿到这套流程准备跑新数据集的时候，我会在进入建模前强制走一遍 sanity check——确认特征矩阵的行名和标签顺序完全对齐，确认交叉验证在特征选择之后，再确认一遍 PCA 图里样本没按批次分群。这套资源本身已经帮你把坑踩了大半，但每个新数据集的批次效应、注释质量和样本量都不一样，你不重新验证一遍，就很难分清结果是生物学信号还是技术噪声。希望这一套流程梳理下来，能帮你在 NAT10 下游基因预测的路上少绕几个弯。

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