简介这套资源面向生物信息学与机器学习交叉领域的研究者聚焦与NAT10相关的下游基因预测任务完整覆盖从数据获取、清洗预处理、特征构建到模型训练与可视化的分析流程。压缩包共55个文件包含Python与R脚本、CSV结果表、TXT说明文档以及PNG分析图整体大小约31.5MB目录按数据处理、差异分析、模型评估和结果展示分层组织。项目以公共基因表达数据为基础整合支持向量机、随机森林、主成分降维与limma差异分析等方法输出特征重要性排序、基因筛选共识和互作网络图等结果。无论是想复现完整基因预测流程还是学习多组学数据整合建模思路都能从脚本注释和中间结果中受益。资源目前已有53人学习下载适合具备一定编程基础、希望深入理解NAT10调控网络的入门及进阶学习者。1. 为什么单独挑出“NAT10 下游基因”来做机器学习预测拿到“基于生物信息学及机器学习预测与 NAT10 相关的下游基.zip”这个压缩包多数人的第一反应是解压后跑一遍差异表达把显著基因直接当答案。但真正常见的场景是按 NAT10 表达量高低分组后top 差异基因列出来全是核糖体蛋白和增殖标志基因换一个数据集立刻翻车。原因不在于数据质量而在于“下游基因”本身没有标准定义——你把 NAT10 当因变量还是自变量选什么混杂因素做控制用什么模型从两万个基因里挑特征每一步都影响最终列表。这篇笔记把能稳定复现的完整流程讲清楚从 zip 包整理、差异表达筛选、偏相关分析到机器学习选特征与验证适合正在做 ac4C 修饰、肿瘤转录组挖掘或想用机器学习替代人工挑基因的从业者。2. 拿到 zip 包之后解压、文件辨认与环境准备2.1 zip 包里通常装了什么先列清单再动手这份压缩包没有固定的文件清单不同作者打包习惯差异很大但做 NAT10 下游预测最常见的内容是三块表达矩阵、样本临床信息、以及作者跑过的中间结果差异表达表或相关性表。先不要急着解压第一步用unzip -l看目录结构确认没有混入奇怪的可执行文件。unzip -l NAT10_predict.zip # -l 只列出内容不实际解压适合先看结构 unzip -n NAT10_predict.zip -d NAT10_project # -n 表示不覆盖已存在的文件防止覆盖你之前改过的脚本或注释文件 unzip -t NAT10_predict.zip # -t 测试压缩包完整性如果输出 No errors detected 再继续-l参数在 Windows 下可能显示中文文件名乱码这是 zip 编码问题而不是文件损坏。常见做法是先用python -c import zipfile; zipfile.ZipFile(NAT10_predict.zip).namelist()在 Python 环境里看一遍真实文件名再决定用unzip -O gbk还是直接手动改名。如果解压时报extra bytes at beginning or within zipfile大概率是 zip 伪加密或文件头被改动过用zip -FF input.zip --out output.zip修复后通常能读出来。zip 伪加密在热词里经常被提到它的特征是解压时要求密码但文件其实没有加密这种情况下别直接丢给密码工具先尝试修复文件头。解压后我最关心的是表达矩阵的形态因为后面所有分析都建立在它上面。常见做法是用 pandas 读进来确认行是基因、列是样本还是反过来。很多翻车现场都是在这里开始的作者给的表达矩阵是基因在行但样本 ID 里混合了正常组织和肿瘤组织或者矩阵里有 NaN 没有处理或者某个基因名重复出现多次。这些都是后续分析的隐藏雷区。import pandas as pd expr pd.read_csv(NAT10_project/expression_matrix.csv, index_col0) print(expr.shape) print(expr.iloc[:4, :4]) print(expr.index.duplicated().sum()) # 检查重复基因名如果基因名是 Symbol如 TP53要留意同名的非编码 RNA如果索引是 Entrez ID建模前最好统一注释版本否则外部验证时基因对不上。第一列如果是探针 IDGEO 数据常见的格式还需要先用平台注释文件完成探针到基因名的映射这一步偷懒会导致后面机器学习特征名全是探针编号交叉验证都不知道哪些基因在互相替代。2.2 环境依赖最小可运行的那组包NAT10 下游预测的全流程不需要 GPU也不需要复杂环境一个 Python 3.8 环境足够。用到的包和各自承担的角色如下import pandas as pd import numpy as np from scipy import stats from statsmodels.stats.multitest import multipletests from sklearn.model_selection import train_test_split, GroupKFold, StratifiedKFold from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LassoCV, LogisticRegression from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import roc_auc_score import shap如果你想把差异表达分析也统一在 Python 里做statsmodels的 t 检验加上multipletests做 BH 矫正可以应付但正式项目我仍建议用 R 的limma跑差异表达后面第 3 章会给完整代码。机器学习部分 scikit-learn 够用shap用来做可解释性输出LassoCV用于特征压缩。提示不要在这时候升级numpy和scikit-learn到最新版老代码里shap的TreeExplainer与新版scikit-learn偶尔不兼容。如果你用的是 conda建议单独为这个项目建一个新环境避免把日常开发环境搞坏。3. 差异表达与偏相关筛选候选基因的两条腿3.1 普通相关分析为什么不可靠用偏相关控制增殖混杂NAT10 是 mRNA ac4C 修饰的关键酶高表达样本通常伴随更强的细胞增殖和翻译活性。如果直接用 Pearson 或 Spearman 相关找与 NAT10 表达共变的基因捞出来的往往是 MKI67、PCNA、核糖体蛋白这类增殖标志物。这些基因不能说和 NAT10 无关但它们作为“下游候选基因”没有操作价值因为敲低 NAT10 后这些基因的表达变化是间接效应。常见做法是计算偏相关——在控制增殖标志基因表达的前提下看某个基因与 NAT10 的独立关联。from scipy import stats # 假设 expr 是基因 x 样本的 DataFrameNAT10 在 index 中 nat10_expr expr.loc[NAT10].values proliferation_markers [MKI67, PCNA, TOP2A] marker_expr expr.loc[proliferation_markers].T.values partial_pvals {} for gene in expr.index: if gene in proliferation_markers or gene NAT10: continue target expr.loc[gene].values # 先把 NAT10 和目标基因分别对增殖标志物回归取残差 X_confound np.column_stack([np.ones(len(target)), marker_expr]) beta_nat10, _, _, _ np.linalg.lstsq(X_confound, nat10_expr, rcondNone) beta_gene, _, _, _ np.linalg.lstsq(X_confound, target, rcondNone) resid_nat10 nat10_expr - X_confound beta_nat10 resid_gene target - X_confound beta_gene r, p stats.spearmanr(resid_nat10, resid_gene) partial_pvals[gene] p逻辑解释把 NAT10 和目标基因分别对增殖标志物做线性回归所得到的残差代表“去除增殖信号后仍留下的变化”残差之间的 Spearman 相关就是偏相关。np.linalg.lstsq用的是最小二乘解残差归一化后乘上向量长度才是方差这里只用于后续相关计算不需要手动除自由度。参数上选择 log2(TPM1) 表达量做计算比使用原始 count 更稳因为 count 数据的离散度受测序深度影响大偏相关对异方差很敏感。p值需要用 BH 方法做多重检验矫正直接把multipletests(partial_pvals.values(), methodfdr_bh)的结果合并回 DataFrame。筛选阈值我一般取 FDR 0.05 且 |偏相关系数| 0.3。这个阈值不要定死样本量 200 和样本量 1000 时同一阈值代表的显著性完全不同。3.2 差异表达筛选的正规姿势R 的 limma 与 Python 的近似实现差异表达分析在转录组领域默认用 Rlimma的线性模型加经验贝叶斯方法对中小样本量的方差估计比 Python 里简单的 t 检验稳定得多。按 NAT10 表达量中位数分高低组就得到一个最简单的二分组设计library(limma) expr - readRDS(NAT10_project/expr_matrix.rds) nat10 - expr[NAT10, ] group - ifelse(nat10 median(nat10), High, Low) design - model.matrix(~ group) fit - lmFit(expr, design) fit - eBayes(fit) res - topTable(fit, coef 2, number Inf, p.value 0.05, lfc 1) write.csv(res, NAT10_DEGs.csv)coef 2表示取groupLow这一项的系数lfc 1对应 log2 差异倍数阈值为 1也就是表达量相差 2 倍。topTable内部的p.value已经过 BH 矫正不需要额外做 FDR。这里很容易踩的坑是model.matrix(~ group)默认把第一个水平设为基线R 按字母序会以 High 为基线得到的系数方向可能和你的直觉相反。分析前先levels(factor(group))确认分组顺序。如果你坚持纯 Python 环境跑可以用scipy.stats.ttest_ind加 BH 矫正近似替代但要注意样本量不均衡时 Welch 校正和 Student t 的选择。from scipy import stats from statsmodels.stats.multitest import multipletests high expr.loc[:, expr.loc[NAT10] np.median(expr.loc[NAT10])] low expr.loc[:, expr.loc[NAT10] np.median(expr.loc[NAT10])] pvals [] logFCs [] for gene in expr.index: hv high.loc[gene].values lv low.loc[gene].values t, p stats.ttest_ind(hv, lv, equal_varFalse) pvals.append(p) logFCs.append(np.log2(np.mean(hv) / np.mean(lv))) fdr multipletests(pvals, methodfdr_bh)[1] deg_df pd.DataFrame({logFC: logFCs, PValue: pvals, FDR: fdr}, indexexpr.index) deg_df deg_df[(deg_df[FDR] 0.05) (deg_df[logFC].abs() 1)]equal_varFalse对应 Welch t 检验在高低组样本量相差大时比普通 t 检验稳健。logFC的计算用均值比在表达值为 0 的基因上会报inf这类基因直接删掉不要填 0 再继续算否则机器学习特征里会引入无意义的常数项。3.3 候选基因不是越多越好控制规模再进模型把偏相关筛选出的基因和差异表达基因取交集你大概会得到几百到几千个基因。这些基因直接作为机器学习特征样本量只有几百模型极容易过拟合。我一般会在取交集后再做一次硬性过滤优先保留偏相关显著的基因并统计差异表达方向与偏相关方向一致性。所谓方向一致性是看一个基因在 NAT10 高表达组里上调其偏相关是否为正反过来成立。这个过滤会去掉大部分“差异表达但方向不稳定”的基因把候选集压到 50~300 之间。如果交集后仍超过 300就按偏相关 p 值排序取前 300。这个数字不是金标准但它保证了后续 Lasso 和随机森林的特征维度在可控范围且不损失主要信号。4. 机器学习建模从候选基因到预测模型4.1 标签与特征定义把预测问题改成二分类机器学习的标签我们没得选就是要预测“NAT10 高表达 vs 低表达”这个二分类状态。特征是 3.3 节筛选出的候选基因表达值。这里有个容易犯的逻辑错误尽量不要把 NAT10 自身放进特征模型会直接学成 NAT10 预测 NAT10得到的特征重要性全是噪音。标签按中位数切分虽然损失了部分信息但在样本量只有几百的转录组数据里回归任务的误差通常比分类更大业界最常用的就是分组后分类。candidate_genes list(deg_intersect_df.index[:300]) # 前 300 个候选基因 X expr.loc[candidate_genes].T # 样本 x 特征 y (expr.loc[NAT10] np.median(expr.loc[NAT10])).astype(int) X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, stratifyy, random_state42 )stratifyy保证训练集和测试集里高低组的比例与全集一致类别不平衡的时候尤其重要。random_state42固定切分结果方便排查为什么某次改动后指标出现微小波动。4.2 Lasso 特征压缩与随机森林对照Lasso 在线性模型基础上加 L1 正则会把不重要特征的系数压缩为 0非常适合基因表达这种高维小样本场景。LassoCV通过交叉验证自动选择正则强度不需要手工调alpha。from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LassoCV scaler StandardScaler() X_train_s scaler.fit_transform(X_train) X_test_s scaler.transform(X_test) lasso LassoCV(cv5, random_state42, max_iter100000).fit(X_train_s, y_train) selected_lasso X.columns[lasso.coef_ ! 0] print(Lasso 选中的基因数:, len(selected_lasso)) print(测试集 AUC:, roc_auc_score(y_test, lasso.predict(X_test_s)))standard_scaler必须在训练集上fit测试集只transform如果对测试集重新fit_transform就是数据泄漏。max_iter100000是为了防止某些基因数值尺度差异过大导致坐标下降不收敛Lasso 默认的 1000 次迭代在表达矩阵上偶尔会报收敛警告。随机森林用作对照看非线性交互能否进一步提升判别力from sklearn.ensemble import RandomForestClassifier rf RandomForestClassifier( n_estimators500, max_depthNone, min_samples_leaf5, random_state42, n_jobs-1, ) rf.fit(X_train_s, y_train) print(RF 测试集 AUC:, roc_auc_score(y_test, rf.predict_proba(X_test_s)[:, 1]))min_samples_leaf5是控制过拟合的关键参数默认的 1 在基因表达数据上会学到大量只覆盖两三个样本的规则。max_depthNone让树自由生长但靠min_samples_leaf兜底。这个组合在多次实践里比手工限制max_depth10表现更稳定。如果 Lasso 和随机森林给出的特征列表交集很小不要慌张。Lasso 是线性系数随机森林能捕捉非线性关系和交互效应两个模型从不同角度找信号。最后用于下游验证的基因建议以两者交集优先单模型独有的基因作为备选。4.3 SHAP 解释特征重要性的另一种视角随机森林的黑匣子问题靠 SHAP 解决它计算每个特征对每个样本预测结果的边际贡献。TreeExplainer只支持树模型所以要用随机森林而不是 Lasso。explainer shap.TreeExplainer(rf) shap_values explainer.shap_values(X_test_s) shap.summary_plot(shap_values[1], X_test_s, feature_namesX_test.columns, max_display20)shap_values[1]取的是“正类”的贡献值也就是 NAT10 高表达组。summary_plot的横轴是 SHAP 值颜色从蓝到红表示特征值从低到高一眼能看出哪些基因高表达时把样本推向 NAT10 高表达组。SHAP 和差异表达 top 基因对不上是很常见的现象。差异表达看的是单变量边际效应SHAP 看的是多变量条件效应。一个基因单独看差异不大但与其他几个基因组合在一起时对判别 NAT10 状态有独特贡献这在生物学上是真实的信号而不是冲突。我建议把 SHAP 排名前 20 的基因作为湿实验验证的优先对象因为它们经过了多变量控制比单纯差异表达表的重复性更好。5. 避免踩坑与常见问题排查从数据泄漏到假信号5.1 现象相关分析 top 基因全是核糖体蛋白换了数据集就消失原因NAT10 参与 rRNA 加工和核糖体生物合成天然与核糖体蛋白表达共变同时细胞增殖状态会同时驱动 NAT10 和核糖体相关基因上调这种混杂在没有控制变量时无法被相关分析区分。解决用第 3 章的偏相关方法至少控制 MKI67、PCNA 两个增殖标志物。另一个补充做法是把 NAT10 表达量自身放入回归模型当自变量看某个基因的残差是否仍与 NAT10 显著相关——这本质上就是偏相关。如果做完控制后核糖体蛋白仍然显著再考虑它们确实是 NAT10 的直接调控下游而不是直接淘汰。关键判断标准是换一个独立数据集后这些基因是否重复出现。5.2 现象随机切分数据时 AUC 0.98按患者切分后掉到 0.65原因同一个患者的多个组织样本被同时分到了训练集和测试集模型记住了患者特有的表达模式而不是 NAT10 相关的通用规律。这在 TCGA 多组织样本数据里尤其严重。解决用GroupKFold按患者 ID 分组交叉验证而不是普通的train_test_split。from sklearn.model_selection import GroupKFold groups patient_id # 与 X 对齐的患者 ID 列表 gkf GroupKFold(n_splits5) aucs [] for train_idx, val_idx in gkf.split(X, y, groups): X_tr, X_va X.iloc[train_idx], X.iloc[val_idx] y_tr, y_va y.iloc[train_idx], y.iloc[val_idx] scaler StandardScaler().fit(X_tr) lasso LassoCV(cv5, random_state42).fit(scaler.transform(X_tr), y_tr) auc roc_auc_score(y_va, lasso.predict(scaler.transform(X_va))) aucs.append(auc) print(GroupKFold AUC:, np.mean(aucs))GroupKFold不像StratifiedGroupKFold那样同时保证每个折内的类别比例和患者分组当每个患者只有一两个样本且类别分布极端失衡时部分折里可能全是高表达样本。遇到这种情况就用StratifiedGroupKFold或者自己写一层按患者分组再分层的逻辑。比较 AUC 时要看 5 折的均值与方差单一折的 AUC 没有参考价值。5.3 现象Lasso 选出的基因和差异表达表对不上原因差异表达是单变量分析Lasso 是多变量联合分析。差异表达显著的基因可能两两高度相关Lasso 只需保留其中一个代表就能达到同等判别力其他相关基因系数被压缩为 0。这不是 bug而是正则化模型的正常行为。解决不要把差异表达表当作“真值”把 Lasso 的筛选结果当作“浓缩特征”。报告结果时说明筛选策略是差异表达预筛 Lasso 压缩而不是直接用差异表达列表当预测特征。另外Lasso 选基因有随机性建议用 100 次自助采样跑 Lasso保留重复出现频次超过 80% 的基因作为稳定特征。这比单次 Lasso 的结果更有说服力也更容易过审。5.4 现象zip 解压后文件路径带中文或空格脚本一跑就报错原因Windows 下生成的 zip 包在 Linux 解压后文件名编码错乱或者路径名含空格导致read_csv默认分隔符识别错误。解决用 Python 的zipfile手动解压并重命名文件。import zipfile with zipfile.ZipFile(NAT10_predict.zip) as zf: for info in zf.infolist(): # 规范化文件名替换空格与特殊字符 new_name info.filename.replace( , _).replace(#, ) info.filename new_name zf.extract(info, pathNAT10_project/) import pandas as pd expr pd.read_csv(NAT10_project/expression_matrix.csv, index_col0, sepNone, enginepython)sepNone配合enginepython会让 pandas 自动检测分隔符在逗号、制表符混用的文件里很实用但会牺牲一部分读取速度。文件只有几 MB 时可以忽略性能问题。读取时如果报Expected 5 fields in line 3, saw 7说明分隔符被手动指定错了。6. 让预测结果真正落地独立验证与可解释性检查6.1 在独立数据集上评估模型的泛化能力本地交叉验证只能证明流程没写错不能证明基因列表可靠。最硬核的验证方式是拿一个完全独立的数据集比如 GEO 上的其他癌症数据集或另外一批细胞系转录组跑同样流程。需要注意三点基因名要统一到同一种 ID 类型NAT10 表达量必须重新计算中位数分组而不能沿用训练集的分组阈值训练好的scaler只能transform不能重新fit。external pd.read_csv(GSE_external_expr.csv, index_col0) common_genes [g for g in selected_lasso if g in external.index] X_ext external.loc[common_genes].T y_ext (external.loc[NAT10] np.median(external.loc[NAT10])).astype(int) X_ext_s scaler.transform(X_ext) # 复用训练集的 scaler ext_auc roc_auc_score(y_ext, lasso.predict(X_ext_s)) print(外部验证 AUC:, ext_auc)如果外部验证 AUC 接近 0.5不要急着怪数据。先检查训练集和外部数据的样本类型是否一致组织、细胞系、不同平台再检查是否做过批次矫正最后检查基因名匹配数量是否过少。外部数据集样本量不足 50 时AUC 的置信区间会非常宽单看一个点估计容易得出错误结论。6.2 一组常用的参数参考与取舍依据环节参数常用取值调整依据表达定量log2(TPM1)统一转换直接对 count 做回归会放大高表达基因的离散度差异表达lfc / FDR1 / 0.05候选基因太多或太少时先调 lfc偏相关阈值相关系数 / FDR0.3 / 0.05样本量 300 可放宽到 0.25LassoCVcv / max_iter5 / 100000样本量 100 时 cv 调 3 更稳随机森林min_samples_leaf5特征维度 1000 时调大到 10独立验证基因 ID 匹配数80%匹配太少说明注释版本不一致这组参数是我在这些项目里的默认起点不是金标准。真正稳定的做法是跑参数敏感性分析每个参数上下浮动 50%看最终基因列表前 20 名的重合率。重合率超过 40% 说明结论稳固低于 20% 则说明数据信号太弱强行选出的特征没有生物学意义。6.3 我的验证习惯与最后提醒每次建模完成后我不会只盯着 AUC而是把模型选出的基因做一次 GSEA 富集看是否落到 RNA 修饰、细胞周期、p53 通路这类与 NAT10 生物学功能可能相关的通路上。如果富集结果全是细胞因子或免疫球蛋白我大概率会重新检查样本组成因为这类信号通常来自基质细胞污染而不是 NAT10 本身的调控。做湿实验前的最后一个习惯是把特征基因按 SHAP 值排序取前 10 个验证表达变化方向而不是看差异表达倍数最大的 10 个。差异表达倍数大的基因受动态范围影响大SHAP 排名更能反映在多变量条件下的独特贡献。这个区别我花过不少时间才理解最初直接用全部差异表达基因做预测换数据集后 AUC 掉到 0.5 附近后来才意识到标签定义和多变量控制才是真正的核心。希望这篇笔记能帮你在 NAT10 下游预测上少走这些弯路。本文还有配套的精品资源点击获取
基于机器学习与生物信息学的NAT10下游基因预测全流程解析
简介这套资源面向生物信息学与机器学习交叉领域的研究者聚焦与NAT10相关的下游基因预测任务完整覆盖从数据获取、清洗预处理、特征构建到模型训练与可视化的分析流程。压缩包共55个文件包含Python与R脚本、CSV结果表、TXT说明文档以及PNG分析图整体大小约31.5MB目录按数据处理、差异分析、模型评估和结果展示分层组织。项目以公共基因表达数据为基础整合支持向量机、随机森林、主成分降维与limma差异分析等方法输出特征重要性排序、基因筛选共识和互作网络图等结果。无论是想复现完整基因预测流程还是学习多组学数据整合建模思路都能从脚本注释和中间结果中受益。资源目前已有53人学习下载适合具备一定编程基础、希望深入理解NAT10调控网络的入门及进阶学习者。1. 为什么单独挑出“NAT10 下游基因”来做机器学习预测拿到“基于生物信息学及机器学习预测与 NAT10 相关的下游基.zip”这个压缩包多数人的第一反应是解压后跑一遍差异表达把显著基因直接当答案。但真正常见的场景是按 NAT10 表达量高低分组后top 差异基因列出来全是核糖体蛋白和增殖标志基因换一个数据集立刻翻车。原因不在于数据质量而在于“下游基因”本身没有标准定义——你把 NAT10 当因变量还是自变量选什么混杂因素做控制用什么模型从两万个基因里挑特征每一步都影响最终列表。这篇笔记把能稳定复现的完整流程讲清楚从 zip 包整理、差异表达筛选、偏相关分析到机器学习选特征与验证适合正在做 ac4C 修饰、肿瘤转录组挖掘或想用机器学习替代人工挑基因的从业者。2. 拿到 zip 包之后解压、文件辨认与环境准备2.1 zip 包里通常装了什么先列清单再动手这份压缩包没有固定的文件清单不同作者打包习惯差异很大但做 NAT10 下游预测最常见的内容是三块表达矩阵、样本临床信息、以及作者跑过的中间结果差异表达表或相关性表。先不要急着解压第一步用unzip -l看目录结构确认没有混入奇怪的可执行文件。unzip -l NAT10_predict.zip # -l 只列出内容不实际解压适合先看结构 unzip -n NAT10_predict.zip -d NAT10_project # -n 表示不覆盖已存在的文件防止覆盖你之前改过的脚本或注释文件 unzip -t NAT10_predict.zip # -t 测试压缩包完整性如果输出 No errors detected 再继续-l参数在 Windows 下可能显示中文文件名乱码这是 zip 编码问题而不是文件损坏。常见做法是先用python -c import zipfile; zipfile.ZipFile(NAT10_predict.zip).namelist()在 Python 环境里看一遍真实文件名再决定用unzip -O gbk还是直接手动改名。如果解压时报extra bytes at beginning or within zipfile大概率是 zip 伪加密或文件头被改动过用zip -FF input.zip --out output.zip修复后通常能读出来。zip 伪加密在热词里经常被提到它的特征是解压时要求密码但文件其实没有加密这种情况下别直接丢给密码工具先尝试修复文件头。解压后我最关心的是表达矩阵的形态因为后面所有分析都建立在它上面。常见做法是用 pandas 读进来确认行是基因、列是样本还是反过来。很多翻车现场都是在这里开始的作者给的表达矩阵是基因在行但样本 ID 里混合了正常组织和肿瘤组织或者矩阵里有 NaN 没有处理或者某个基因名重复出现多次。这些都是后续分析的隐藏雷区。import pandas as pd expr pd.read_csv(NAT10_project/expression_matrix.csv, index_col0) print(expr.shape) print(expr.iloc[:4, :4]) print(expr.index.duplicated().sum()) # 检查重复基因名如果基因名是 Symbol如 TP53要留意同名的非编码 RNA如果索引是 Entrez ID建模前最好统一注释版本否则外部验证时基因对不上。第一列如果是探针 IDGEO 数据常见的格式还需要先用平台注释文件完成探针到基因名的映射这一步偷懒会导致后面机器学习特征名全是探针编号交叉验证都不知道哪些基因在互相替代。2.2 环境依赖最小可运行的那组包NAT10 下游预测的全流程不需要 GPU也不需要复杂环境一个 Python 3.8 环境足够。用到的包和各自承担的角色如下import pandas as pd import numpy as np from scipy import stats from statsmodels.stats.multitest import multipletests from sklearn.model_selection import train_test_split, GroupKFold, StratifiedKFold from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LassoCV, LogisticRegression from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import roc_auc_score import shap如果你想把差异表达分析也统一在 Python 里做statsmodels的 t 检验加上multipletests做 BH 矫正可以应付但正式项目我仍建议用 R 的limma跑差异表达后面第 3 章会给完整代码。机器学习部分 scikit-learn 够用shap用来做可解释性输出LassoCV用于特征压缩。提示不要在这时候升级numpy和scikit-learn到最新版老代码里shap的TreeExplainer与新版scikit-learn偶尔不兼容。如果你用的是 conda建议单独为这个项目建一个新环境避免把日常开发环境搞坏。3. 差异表达与偏相关筛选候选基因的两条腿3.1 普通相关分析为什么不可靠用偏相关控制增殖混杂NAT10 是 mRNA ac4C 修饰的关键酶高表达样本通常伴随更强的细胞增殖和翻译活性。如果直接用 Pearson 或 Spearman 相关找与 NAT10 表达共变的基因捞出来的往往是 MKI67、PCNA、核糖体蛋白这类增殖标志物。这些基因不能说和 NAT10 无关但它们作为“下游候选基因”没有操作价值因为敲低 NAT10 后这些基因的表达变化是间接效应。常见做法是计算偏相关——在控制增殖标志基因表达的前提下看某个基因与 NAT10 的独立关联。from scipy import stats # 假设 expr 是基因 x 样本的 DataFrameNAT10 在 index 中 nat10_expr expr.loc[NAT10].values proliferation_markers [MKI67, PCNA, TOP2A] marker_expr expr.loc[proliferation_markers].T.values partial_pvals {} for gene in expr.index: if gene in proliferation_markers or gene NAT10: continue target expr.loc[gene].values # 先把 NAT10 和目标基因分别对增殖标志物回归取残差 X_confound np.column_stack([np.ones(len(target)), marker_expr]) beta_nat10, _, _, _ np.linalg.lstsq(X_confound, nat10_expr, rcondNone) beta_gene, _, _, _ np.linalg.lstsq(X_confound, target, rcondNone) resid_nat10 nat10_expr - X_confound beta_nat10 resid_gene target - X_confound beta_gene r, p stats.spearmanr(resid_nat10, resid_gene) partial_pvals[gene] p逻辑解释把 NAT10 和目标基因分别对增殖标志物做线性回归所得到的残差代表“去除增殖信号后仍留下的变化”残差之间的 Spearman 相关就是偏相关。np.linalg.lstsq用的是最小二乘解残差归一化后乘上向量长度才是方差这里只用于后续相关计算不需要手动除自由度。参数上选择 log2(TPM1) 表达量做计算比使用原始 count 更稳因为 count 数据的离散度受测序深度影响大偏相关对异方差很敏感。p值需要用 BH 方法做多重检验矫正直接把multipletests(partial_pvals.values(), methodfdr_bh)的结果合并回 DataFrame。筛选阈值我一般取 FDR 0.05 且 |偏相关系数| 0.3。这个阈值不要定死样本量 200 和样本量 1000 时同一阈值代表的显著性完全不同。3.2 差异表达筛选的正规姿势R 的 limma 与 Python 的近似实现差异表达分析在转录组领域默认用 Rlimma的线性模型加经验贝叶斯方法对中小样本量的方差估计比 Python 里简单的 t 检验稳定得多。按 NAT10 表达量中位数分高低组就得到一个最简单的二分组设计library(limma) expr - readRDS(NAT10_project/expr_matrix.rds) nat10 - expr[NAT10, ] group - ifelse(nat10 median(nat10), High, Low) design - model.matrix(~ group) fit - lmFit(expr, design) fit - eBayes(fit) res - topTable(fit, coef 2, number Inf, p.value 0.05, lfc 1) write.csv(res, NAT10_DEGs.csv)coef 2表示取groupLow这一项的系数lfc 1对应 log2 差异倍数阈值为 1也就是表达量相差 2 倍。topTable内部的p.value已经过 BH 矫正不需要额外做 FDR。这里很容易踩的坑是model.matrix(~ group)默认把第一个水平设为基线R 按字母序会以 High 为基线得到的系数方向可能和你的直觉相反。分析前先levels(factor(group))确认分组顺序。如果你坚持纯 Python 环境跑可以用scipy.stats.ttest_ind加 BH 矫正近似替代但要注意样本量不均衡时 Welch 校正和 Student t 的选择。from scipy import stats from statsmodels.stats.multitest import multipletests high expr.loc[:, expr.loc[NAT10] np.median(expr.loc[NAT10])] low expr.loc[:, expr.loc[NAT10] np.median(expr.loc[NAT10])] pvals [] logFCs [] for gene in expr.index: hv high.loc[gene].values lv low.loc[gene].values t, p stats.ttest_ind(hv, lv, equal_varFalse) pvals.append(p) logFCs.append(np.log2(np.mean(hv) / np.mean(lv))) fdr multipletests(pvals, methodfdr_bh)[1] deg_df pd.DataFrame({logFC: logFCs, PValue: pvals, FDR: fdr}, indexexpr.index) deg_df deg_df[(deg_df[FDR] 0.05) (deg_df[logFC].abs() 1)]equal_varFalse对应 Welch t 检验在高低组样本量相差大时比普通 t 检验稳健。logFC的计算用均值比在表达值为 0 的基因上会报inf这类基因直接删掉不要填 0 再继续算否则机器学习特征里会引入无意义的常数项。3.3 候选基因不是越多越好控制规模再进模型把偏相关筛选出的基因和差异表达基因取交集你大概会得到几百到几千个基因。这些基因直接作为机器学习特征样本量只有几百模型极容易过拟合。我一般会在取交集后再做一次硬性过滤优先保留偏相关显著的基因并统计差异表达方向与偏相关方向一致性。所谓方向一致性是看一个基因在 NAT10 高表达组里上调其偏相关是否为正反过来成立。这个过滤会去掉大部分“差异表达但方向不稳定”的基因把候选集压到 50~300 之间。如果交集后仍超过 300就按偏相关 p 值排序取前 300。这个数字不是金标准但它保证了后续 Lasso 和随机森林的特征维度在可控范围且不损失主要信号。4. 机器学习建模从候选基因到预测模型4.1 标签与特征定义把预测问题改成二分类机器学习的标签我们没得选就是要预测“NAT10 高表达 vs 低表达”这个二分类状态。特征是 3.3 节筛选出的候选基因表达值。这里有个容易犯的逻辑错误尽量不要把 NAT10 自身放进特征模型会直接学成 NAT10 预测 NAT10得到的特征重要性全是噪音。标签按中位数切分虽然损失了部分信息但在样本量只有几百的转录组数据里回归任务的误差通常比分类更大业界最常用的就是分组后分类。candidate_genes list(deg_intersect_df.index[:300]) # 前 300 个候选基因 X expr.loc[candidate_genes].T # 样本 x 特征 y (expr.loc[NAT10] np.median(expr.loc[NAT10])).astype(int) X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, stratifyy, random_state42 )stratifyy保证训练集和测试集里高低组的比例与全集一致类别不平衡的时候尤其重要。random_state42固定切分结果方便排查为什么某次改动后指标出现微小波动。4.2 Lasso 特征压缩与随机森林对照Lasso 在线性模型基础上加 L1 正则会把不重要特征的系数压缩为 0非常适合基因表达这种高维小样本场景。LassoCV通过交叉验证自动选择正则强度不需要手工调alpha。from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LassoCV scaler StandardScaler() X_train_s scaler.fit_transform(X_train) X_test_s scaler.transform(X_test) lasso LassoCV(cv5, random_state42, max_iter100000).fit(X_train_s, y_train) selected_lasso X.columns[lasso.coef_ ! 0] print(Lasso 选中的基因数:, len(selected_lasso)) print(测试集 AUC:, roc_auc_score(y_test, lasso.predict(X_test_s)))standard_scaler必须在训练集上fit测试集只transform如果对测试集重新fit_transform就是数据泄漏。max_iter100000是为了防止某些基因数值尺度差异过大导致坐标下降不收敛Lasso 默认的 1000 次迭代在表达矩阵上偶尔会报收敛警告。随机森林用作对照看非线性交互能否进一步提升判别力from sklearn.ensemble import RandomForestClassifier rf RandomForestClassifier( n_estimators500, max_depthNone, min_samples_leaf5, random_state42, n_jobs-1, ) rf.fit(X_train_s, y_train) print(RF 测试集 AUC:, roc_auc_score(y_test, rf.predict_proba(X_test_s)[:, 1]))min_samples_leaf5是控制过拟合的关键参数默认的 1 在基因表达数据上会学到大量只覆盖两三个样本的规则。max_depthNone让树自由生长但靠min_samples_leaf兜底。这个组合在多次实践里比手工限制max_depth10表现更稳定。如果 Lasso 和随机森林给出的特征列表交集很小不要慌张。Lasso 是线性系数随机森林能捕捉非线性关系和交互效应两个模型从不同角度找信号。最后用于下游验证的基因建议以两者交集优先单模型独有的基因作为备选。4.3 SHAP 解释特征重要性的另一种视角随机森林的黑匣子问题靠 SHAP 解决它计算每个特征对每个样本预测结果的边际贡献。TreeExplainer只支持树模型所以要用随机森林而不是 Lasso。explainer shap.TreeExplainer(rf) shap_values explainer.shap_values(X_test_s) shap.summary_plot(shap_values[1], X_test_s, feature_namesX_test.columns, max_display20)shap_values[1]取的是“正类”的贡献值也就是 NAT10 高表达组。summary_plot的横轴是 SHAP 值颜色从蓝到红表示特征值从低到高一眼能看出哪些基因高表达时把样本推向 NAT10 高表达组。SHAP 和差异表达 top 基因对不上是很常见的现象。差异表达看的是单变量边际效应SHAP 看的是多变量条件效应。一个基因单独看差异不大但与其他几个基因组合在一起时对判别 NAT10 状态有独特贡献这在生物学上是真实的信号而不是冲突。我建议把 SHAP 排名前 20 的基因作为湿实验验证的优先对象因为它们经过了多变量控制比单纯差异表达表的重复性更好。5. 避免踩坑与常见问题排查从数据泄漏到假信号5.1 现象相关分析 top 基因全是核糖体蛋白换了数据集就消失原因NAT10 参与 rRNA 加工和核糖体生物合成天然与核糖体蛋白表达共变同时细胞增殖状态会同时驱动 NAT10 和核糖体相关基因上调这种混杂在没有控制变量时无法被相关分析区分。解决用第 3 章的偏相关方法至少控制 MKI67、PCNA 两个增殖标志物。另一个补充做法是把 NAT10 表达量自身放入回归模型当自变量看某个基因的残差是否仍与 NAT10 显著相关——这本质上就是偏相关。如果做完控制后核糖体蛋白仍然显著再考虑它们确实是 NAT10 的直接调控下游而不是直接淘汰。关键判断标准是换一个独立数据集后这些基因是否重复出现。5.2 现象随机切分数据时 AUC 0.98按患者切分后掉到 0.65原因同一个患者的多个组织样本被同时分到了训练集和测试集模型记住了患者特有的表达模式而不是 NAT10 相关的通用规律。这在 TCGA 多组织样本数据里尤其严重。解决用GroupKFold按患者 ID 分组交叉验证而不是普通的train_test_split。from sklearn.model_selection import GroupKFold groups patient_id # 与 X 对齐的患者 ID 列表 gkf GroupKFold(n_splits5) aucs [] for train_idx, val_idx in gkf.split(X, y, groups): X_tr, X_va X.iloc[train_idx], X.iloc[val_idx] y_tr, y_va y.iloc[train_idx], y.iloc[val_idx] scaler StandardScaler().fit(X_tr) lasso LassoCV(cv5, random_state42).fit(scaler.transform(X_tr), y_tr) auc roc_auc_score(y_va, lasso.predict(scaler.transform(X_va))) aucs.append(auc) print(GroupKFold AUC:, np.mean(aucs))GroupKFold不像StratifiedGroupKFold那样同时保证每个折内的类别比例和患者分组当每个患者只有一两个样本且类别分布极端失衡时部分折里可能全是高表达样本。遇到这种情况就用StratifiedGroupKFold或者自己写一层按患者分组再分层的逻辑。比较 AUC 时要看 5 折的均值与方差单一折的 AUC 没有参考价值。5.3 现象Lasso 选出的基因和差异表达表对不上原因差异表达是单变量分析Lasso 是多变量联合分析。差异表达显著的基因可能两两高度相关Lasso 只需保留其中一个代表就能达到同等判别力其他相关基因系数被压缩为 0。这不是 bug而是正则化模型的正常行为。解决不要把差异表达表当作“真值”把 Lasso 的筛选结果当作“浓缩特征”。报告结果时说明筛选策略是差异表达预筛 Lasso 压缩而不是直接用差异表达列表当预测特征。另外Lasso 选基因有随机性建议用 100 次自助采样跑 Lasso保留重复出现频次超过 80% 的基因作为稳定特征。这比单次 Lasso 的结果更有说服力也更容易过审。5.4 现象zip 解压后文件路径带中文或空格脚本一跑就报错原因Windows 下生成的 zip 包在 Linux 解压后文件名编码错乱或者路径名含空格导致read_csv默认分隔符识别错误。解决用 Python 的zipfile手动解压并重命名文件。import zipfile with zipfile.ZipFile(NAT10_predict.zip) as zf: for info in zf.infolist(): # 规范化文件名替换空格与特殊字符 new_name info.filename.replace( , _).replace(#, ) info.filename new_name zf.extract(info, pathNAT10_project/) import pandas as pd expr pd.read_csv(NAT10_project/expression_matrix.csv, index_col0, sepNone, enginepython)sepNone配合enginepython会让 pandas 自动检测分隔符在逗号、制表符混用的文件里很实用但会牺牲一部分读取速度。文件只有几 MB 时可以忽略性能问题。读取时如果报Expected 5 fields in line 3, saw 7说明分隔符被手动指定错了。6. 让预测结果真正落地独立验证与可解释性检查6.1 在独立数据集上评估模型的泛化能力本地交叉验证只能证明流程没写错不能证明基因列表可靠。最硬核的验证方式是拿一个完全独立的数据集比如 GEO 上的其他癌症数据集或另外一批细胞系转录组跑同样流程。需要注意三点基因名要统一到同一种 ID 类型NAT10 表达量必须重新计算中位数分组而不能沿用训练集的分组阈值训练好的scaler只能transform不能重新fit。external pd.read_csv(GSE_external_expr.csv, index_col0) common_genes [g for g in selected_lasso if g in external.index] X_ext external.loc[common_genes].T y_ext (external.loc[NAT10] np.median(external.loc[NAT10])).astype(int) X_ext_s scaler.transform(X_ext) # 复用训练集的 scaler ext_auc roc_auc_score(y_ext, lasso.predict(X_ext_s)) print(外部验证 AUC:, ext_auc)如果外部验证 AUC 接近 0.5不要急着怪数据。先检查训练集和外部数据的样本类型是否一致组织、细胞系、不同平台再检查是否做过批次矫正最后检查基因名匹配数量是否过少。外部数据集样本量不足 50 时AUC 的置信区间会非常宽单看一个点估计容易得出错误结论。6.2 一组常用的参数参考与取舍依据环节参数常用取值调整依据表达定量log2(TPM1)统一转换直接对 count 做回归会放大高表达基因的离散度差异表达lfc / FDR1 / 0.05候选基因太多或太少时先调 lfc偏相关阈值相关系数 / FDR0.3 / 0.05样本量 300 可放宽到 0.25LassoCVcv / max_iter5 / 100000样本量 100 时 cv 调 3 更稳随机森林min_samples_leaf5特征维度 1000 时调大到 10独立验证基因 ID 匹配数80%匹配太少说明注释版本不一致这组参数是我在这些项目里的默认起点不是金标准。真正稳定的做法是跑参数敏感性分析每个参数上下浮动 50%看最终基因列表前 20 名的重合率。重合率超过 40% 说明结论稳固低于 20% 则说明数据信号太弱强行选出的特征没有生物学意义。6.3 我的验证习惯与最后提醒每次建模完成后我不会只盯着 AUC而是把模型选出的基因做一次 GSEA 富集看是否落到 RNA 修饰、细胞周期、p53 通路这类与 NAT10 生物学功能可能相关的通路上。如果富集结果全是细胞因子或免疫球蛋白我大概率会重新检查样本组成因为这类信号通常来自基质细胞污染而不是 NAT10 本身的调控。做湿实验前的最后一个习惯是把特征基因按 SHAP 值排序取前 10 个验证表达变化方向而不是看差异表达倍数最大的 10 个。差异表达倍数大的基因受动态范围影响大SHAP 排名更能反映在多变量条件下的独特贡献。这个区别我花过不少时间才理解最初直接用全部差异表达基因做预测换数据集后 AUC 掉到 0.5 附近后来才意识到标签定义和多变量控制才是真正的核心。希望这篇笔记能帮你在 NAT10 下游预测上少走这些弯路。本文还有配套的精品资源点击获取
/float()与正则提取的对比测试)
