q6-批次矫正和差异分析--复现文章Fig3
华哥微信：15623525389 公众号：Bioinformation
此文献标题：融合机器学习与分子对接揭示黄曲霉毒素B1诱导的肝细胞癌分子网络
DOI: 10.1097/JS9.0000000000002455https://journals.lww.com/international-journal-of-surgery/fulltext/2025/07000/integrating_machine_learning_and_molecular_docking.33.aspx
摘要
目的：本研究旨在探究黄曲霉毒素B1（Aflatoxin B1, AFB1）诱导的肝细胞癌（HCC）的分子机制。
方法：研究对多个转录组数据集进行了差异表达分析，以识别与HCC相关的靶基因。结合机器学习算法、网络毒理学与分子对接技术，探索AFB1与靶蛋白之间的结合相互作用。
结果：共识别出48个可能参与AFB1诱导肝癌发生的候选靶基因。随后通过机器学习分析筛选出6个核心调控基因（RND3、PCK1、AURKA、BCAT2、UCK2 和 CCNB1）。其中，RND3 和 PCK1 表现出显著下调，而 AURKA、BCAT2、UCK2 和 CCNB1 明显上调（P<0.05）。分子对接模拟显示AFB1与这些靶蛋白具有高度结合特异性。
结论：本研究表明，AFB1可能通过作用于特定基因与信号通路，促进肝细胞癌的发生。机器学习鉴定出六个核心调控基因，分子对接验证了AFB1与这些关键靶点之间的高度结合亲和力。这些发现为进一步深入揭示AFB1诱导肝癌的机制提供了关键见解。
关键词：黄曲霉毒素B1、肝细胞癌、生物信息学、机器学习、分子对接
引言
肝细胞癌（Hepatocellular carcinoma, HCC）是全球第六大常见恶性肿瘤，也是癌症相关死亡的第三大原因，其复杂的病因机制构成了严峻的公共卫生挑战。HCC的发病机制是多种环境、病毒及代谢因素共同作用的结果，包括乙型/丙型肝炎病毒（HBV/HCV）感染、酒精性肝病、非酒精性脂肪肝病（NAFLD）以及黄曲霉毒素B1（AFB1）暴露。流行病学证据指出，HBV感染与AFB1暴露之间存在显著的协同作用，这种双重病因通过破坏肝脏稳态显著提高HCC风险。值得注意的是，在亚洲和撒哈拉以南非洲的高发地区，这种病因组合造成了全球超过80%的HCC疾病负担，反映出明显的地域性差异。
在分子层面，HCC的发生发展受基因组不稳定性（如DNA修复障碍）、表观遗传调控紊乱和致癌信号通路（尤其是Wnt/β-catenin和PI3K/Akt/mTOR通路）异常激活的共同驱动。AFB1诱导HCC的一个显著分子特征是TP53基因R249S热点突变，该突变源自AFB1-DNA加合物引发的G→T碱基转变突变，明确了肝癌发生中的“基因型–环境型”连接机制。
AFB1是由黄曲霉（Aspergillus flavus）和寄生曲霉（A. parasiticus）产生的真菌毒素，被国际癌症研究机构（IARC）列为I类致癌物，在高温潮湿地区常污染玉米和花生等主食。AFB1的致癌潜能依赖于肝脏CYP450酶（主要是CYP1A2和CYP3A4）介导的生物活化，形成高活性的AFB1-exo-8,9-环氧化物（AFBO）。该代谢产物可与鸟嘌呤残基共价结合，形成AFB1-N7-dG加合物，诱导G→T突变及染色体不稳定性。
慢性AFB1暴露还会通过持续氧化应激、线粒体膜通透性转变及NF-κB介导的促炎信号通路破坏氧化还原平衡，共同促成肿瘤易感微环境。在机制上，AFB1还可通过HBV X蛋白（HBx）抑制核苷酸切除修复（NER）过程，增强DNA损伤积累，同时激活cGAS-STING通路维持慢性炎症与免疫逃逸。此外，AFB1还能通过激活芳香烃受体（AHR）通路，上调PD-L1等免疫抑制分子，促进肿瘤免疫逃逸与HCC进展。
近年来，网络毒理学的兴起为解析复杂毒物与宿主之间的相互作用提供了变革性方法，融合多组学数据、动态网络建模和结构-活性预测。在AFB1研究中，转录组分析已鉴定MTCH1和PPM1D为关键连接毒性与肿瘤发生的枢纽基因，而脂质组学揭示了鞘脂与PI3K的交叉调控；分子对接预测AFB1对Keap1的变构抑制作用可抑制Nrf2通路的抗氧化功能，该机制已在动物实验中得到验证。
尽管已有研究取得进展，但目前多数研究仍聚焦于线性通路，缺乏系统层面的“暴露–反应–效应”网络模型，亦忽略了基因组损伤与免疫微环境重塑之间的动态互作。为填补该空白，我们本研究采用多角度的系统毒理策略：（1）整合多组学数据重构AFB1诱导HCC的分子网络；（2）通过拓扑与功能富集分析结合机器学习算法优选关键枢纽因子；（3）采用集成分子对接策略验证关键靶点与AFB1之间的结合热力学特性。通过揭示AFB1的多尺度致癌网络，本研究旨在为高风险人群的HCC化学预防与精准治疗提供可干预的分子靶点，回应当前该领域的关键研究需求。
华哥微信：15623525389 公众号：Bioinformation
黄曲霉毒素B1（AFB1）靶蛋白的识别
AFB1的分子结构通过PubChem数据库获取（图2A）。利用三个互补数据库（ChEMBL——生物活性小分子数据库、PharmMapper——反向药效团匹配平台、SwissTargetPrediction——基于配体的靶点预测）系统预测了AFB1的潜在生物学靶点。在整合数据并剔除冗余条目后，共鉴定出371个独特的潜在靶点（图2B）。
肝细胞癌相关靶基因的识别
为尽量减少批次效应，我们整合了GSE25097和GSE36376两个数据集，并对基因表达矩阵进行了全面归一化处理。主成分分析（PCA）显示归一化后数据分布得到了优化，归一化数据集表现出更清晰的聚类模式（图3A-B）。
差异表达分析共识别出797个在HCC中显著改变的基因。这些表达变化通过火山图和热图可视化展示（图3C-D）。
随后进行加权基因共表达网络分析（WGCNA）。首先确定最优软阈值功率（β），以确保网络具有无标度特性。系统评估1至20范围内的功率值后发现，β=5 为满足无标度拓扑标准（R² ≥ 0.8）的最小值。在此基础上，我们构建了拓扑重叠矩阵（TOM），并进行层级聚类分析，最终识别出8个不同的共表达模块，每个模块通过颜色编码进行可视化（图3E）。
模块与表型（HCC）之间的关系通过模块特征基因与样本特征之间的相关性分析确定（P < 0.05）（图3F）。将常规差异分析得到的DEGs与WGCNA模块基因进行整合（去重后），最终获得1179个与HCC相关的候选基因（图3G）。
AFB1相关HCC靶点的识别
通过将AFB1预测靶蛋白与HCC相关基因进行交集分析，最终识别出48个可能参与AFB1诱导肝癌发生的关键靶点（图4A-B）。
进一步通过GO与KEGG富集分析对这些靶点进行功能注释（图4C-D），结果揭示了广泛的分子层面信息。GO分析显示这些基因在多种生物过程（如成纤维细胞增殖调控、有机酸代谢过程）、细胞组分（如线粒体基质、纺锤体极）及分子功能（如激酶调控活性、维生素B6结合、抗氧化功能）中显著富集。
KEGG通路分析表明，这些靶点参与多种关键信号通路（如肝细胞癌、前列腺癌通路），包括PI3K-Akt、p53、FoxO等信号级联；同时还涉及代谢调控（如氨基酸生物合成、蛋氨酸代谢）及治疗抗性机制（如药物代谢酶、EGFR抑制剂耐药）。这些结果共同揭示了AFB1相关肝癌中，肿瘤进展通路激活、代谢重编程以及治疗抵抗机制的显著增强。
AFB1诱导肝癌的核心基因识别
通过对48个候选靶点进行全面的机器学习分析，我们构建了113个预测模型，以识别与AFB1相关肝细胞癌（HCC）发病机制相关的核心基因。在所有模型中，Lasso+RF集成模型在训练阶段和验证阶段均表现出最优的准确性（图5A），从而最终筛选出6个关键基因：RND3、AURKA、PCK1、BCAT2、UCK2 和 CCNB1。
ROC曲线分析证实这些核心基因具有良好的诊断潜力（AUC > 0.85，图5B）；其在HCC组织中的差异表达模式通过火山图进一步可视化（图5C）。SHAP可解释性分析揭示了各基因在模型预测中的功能贡献差异：RND3（SHAP值 = 0.158）和PCK1（SHAP值 = 0.120）为最具影响力的预测因子（图5D），AURKA与CCNB1则显示出双向调控作用（图5E）。
值得注意的是，我们识别出几个关键的非线性表达-预测关系：（1）RND3与AURKA的表达呈负相关；（2）PCK1在高表达水平（>6）时表现出预测效应峰值；（3）CCNB1与BCAT2在中等表达水平（3.5–4.5）时具有最佳预测效应（图5F）。进一步的力导向分析（图5G）表明，AURKA（值为3.17，Δ=-0.162）与RND3（值为5.91，Δ=-0.149）是主要的负调控因子，驱动总体预测值 f(x) = -0.0299 低于期望均值 E[f(x)] = 0.55。
AFB1与核心基因相互作用的分子对接验证
为验证AFB1与上述六个核心基因（RND3、AURKA、PCK1、BCAT2、UCK2 和 CCNB1）之间可能存在的结合相互作用，我们进行了全面的分子对接分析。结果显示，AFB1与所有靶蛋白均具有较强的结合亲和力，其结合能普遍低于 -5 kcal/mol，提示其相互作用具有较高的稳定性与自发性。
根据分子对接研究的既有标准：结合能 < 0 kcal/mol 表示存在自发结合能力；结合能 < -5.0 kcal/mol 表明具有良好结合亲和性。AFB1与靶蛋白的对接构象可视化结果（图6）显示所有复合物均具有稳定的对接构型。这些发现在结构层面上进一步支持了机器学习方法识别出的关键HCC靶点与AFB1之间存在直接的分子相互作用。
华哥微信：15623525389 公众号：Bioinformation
讨论
由黄曲霉（Aspergillus flavus）和寄生曲霉（A. parasiticus）产生的遗传毒性致癌物AFB1，已被机制性地证实通过多个相互关联的通路参与肝细胞癌（HCC）的发病过程。在肝脏CYP450酶（尤其是CYP3A4/CYP1A2）介导的代谢激活作用下，AFB1形成的反应性AFB1-8,9-环氧化物代谢物可与DNA结合，诱导特征性的突变损伤，包括标志性的p53 R249S突变和染色体不稳定。这种基因毒性与乙肝病毒（HBV）合并感染具有显著协同效应，后者可通过诱导肝细胞增殖和削弱DNA修复能力进一步加重损伤。
不同物种对AFB1的易感性差异可能与谷胱甘肽S-转移酶（GST）介导的解毒效率有关。人类GST的多态性在一定程度上解释了不同人群中HCC风险的异质性。除了直接的DNA损伤外，我们的研究结果也与当前证据一致，表明AFB1还可通过多种机制促进肝癌发生发展，包括：致癌信号通路（如Ras信号通路）的表观遗传调控失衡、DNA修复基因（如XRCC1、OGG1）表达的抑制，以及通过JAK-STAT信号通路抑制所介导的免疫调节功能紊乱，这些变化最终促成了免疫逃逸。
AFB1的分布呈明显地域性特征，在热带和亚热带地区的粮食中污染严重，其主要通过受污染主食摄入或农业职业暴露导致高风险人群暴露。值得注意的是，AFB1的多重毒性效应还扩展至胎盘毒性与免疫抑制，可能进一步加重某些传染病的临床结局。这些多层次的致病机制强调了在AFB1流行地区开展综合预防干预的必要性。
本研究通过整合多组学数据，结合先进的机器学习与生物信息学技术，系统识别了AFB1暴露与HCC发病之间潜在的分子靶点。我们基于网络毒理学和生信分析，筛选出48个候选基因，并最终确定了6个核心基因（RND3、PCK1、AURKA、BCAT2、UCK2、CCNB1）具有重要作用。差异表达分析显示，在HCC中，RND3与PCK1稳定下调，而AURKA、BCAT2、UCK2、CCNB1则显著上调。机器学习模型进一步证实该基因组合具有较高的诊断价值，SHAP可解释性分析也强调了RND3（SHAP值 = 0.158）和PCK1（0.120）作为关键预测因子的作用。分子对接模拟为这些发现提供了结构学支持，表明AFB1与上述核心基因的蛋白产物之间存在稳定的分子结合，且受特定氨基酸相互作用介导。
上述核心基因在AFB1诱导的HCC发病机制中具有多重生物学作用。RND3（又名RhoE）作为Rho家族小GTP酶的重要成员，参与细胞周期调控、迁移和凋亡等关键生理过程。近期研究发现，RND3在调控肝细胞对AFB1暴露应激反应中具有重要生物学意义。其核心机制可能通过调控AKT信号通路实现，AKT作为细胞生存与增殖的中央通路，其核转位过程可直接控制促增殖基因的表达。特别是在AFB1暴露的病理状态下，RND3可通过动态调节AKT通路的活化状态，显著影响肝细胞的异常增殖–生存平衡，从而为解释AFB1致癌机制提供了新的视角。鉴于RND3在AFB1相关信号网络中的关键地位，未来可通过靶向其表达水平或功能活性，开发出新的干预策略以预防AFB1相关的HCC。进一步阐明RND3与AFB1信号网络之间的分子互作机制，特别是在表观遗传修饰与蛋白互作网络中的细节，将为高效靶向治疗的开发提供关键理论基础。
在AFB1诱导的肝癌机制中，磷酸烯醇式丙酮酸羧激酶1（PCK1）则承担关键的代谢调控角色。作为糖异生通路的核心酶，PCK1在AFB1暴露后表达水平显著下调，这种失调通过“双通路机制”促进HCC发展：一方面直接扰乱肝细胞的脂质合成稳态，另一方面诱导致癌性的代谢重编程。值得注意的是，AFB1同时下调PCK2与肝细胞生长因子（HRG）表达，导致维持细胞正常功能所需的代谢网络全面失衡，而这种多靶点的代谢紊乱与AFB1的肝毒性和致癌风险密切相关。在分子层面，PCK1功能障碍触发的三羧酸（TCA）循环重构，会导致代谢中间产物的异常外泄，产生“代谢关闭”效应，不仅加剧氧化应激损伤，还通过促凋亡信号与p53突变的协同作用，推动肝细胞的恶性转化。其中，具有获得性致癌功能的p53突变已被确认为AFB1暴露的特异性分子标志物。考虑到PCK1在AFB1代谢毒性网络中的关键作用，靶向其活性调控或干预其参与的代谢通路，或将成为开发预防性与靶向治疗HCC的突破口。
AURKA与CCNB1共同参与有丝分裂过程的调控，AURKA促进纺锤体形成与染色体分离，而CCNB1–CDK1复合物则驱动细胞周期G2/M期转换。这两者在多种癌症中均表现为高表达，提示其作为潜在治疗靶点的可能性。代谢相关基因BCAT2与UCK2分别参与支链氨基酸代谢与嘧啶核苷酸合成，可能通过维持肿瘤细胞的能量与生物合成需求而促进其增殖。此外，UCK2还在病毒RNA合成中发挥作用，具有一定的抗病毒治疗潜力。这些基因间的协同效应揭示了细胞周期异常、代谢重编程与肿瘤发生发展的复杂网络，为后续机制研究与靶向干预提供了理论依据与关键靶点。
结论
本研究表明，黄曲霉毒素B1可能通过作用于特定基因与信号通路影响肝细胞癌的发病机制。分子对接模拟显示，AFB1与多个关键靶蛋白具有显著的结合特异性。这些结果为深入研究AFB1诱导肝癌的分子机制提供了重要基础。未来研究应聚焦于AFB1暴露与HCC风险之间的剂量–反应关系，并进一步探索可能缓解AFB1对肝脏健康不良影响的干预措施。
文章的创新性
引入集成机器学习框架，系统识别AFB1致癌的核心基因调控网络。
融合网络毒理学与分子对接技术，从系统层面构建“AFB1-分子靶点-HCC”多尺度机制图谱。
筛选出6个关键致癌基因（RND3、PCK1、AURKA、BCAT2、UCK2、CCNB1），建立高准确性诊断模型。
通过SHAP解释性算法解析黑箱模型，明确各关键基因的预测贡献。
对接模拟实验证实AFB1与核心蛋白之间存在高亲和力结合，形成结构功能双重验证。
文章的技术路线图
整合五个GEO数据库中的肝癌数据集，消除批次效应后进行差异表达分析与WGCNA共表达分析。
预测AFB1潜在靶点并与HCC相关基因交集，获得48个致癌候选靶点。
结合10种机器学习算法构建113个预测模型，通过堆叠集成法筛选出最优模型并识别核心基因。
应用SHAP算法量化每个基因对预测模型的影响力，提升生物学解释性。
采用AutoDock Vina对AFB1与核心蛋白进行分子对接，结合能验证其生物学结合能力。
文章中生信分析的作用
通过limma和WGCNA系统挖掘HCC相关基因，定位潜在致癌模块。
使用ChEMBL、PharmMapper、SwissTargetPrediction预测AFB1靶点，并整合至HCC网络中。
运用clusterProfiler进行GO与KEGG富集分析，明确48个关键靶点的功能通路特征。
搭建高通量机器学习预测框架，量化基因在诊断模型中的关键性。
辅助分子对接分析结果定位高结合能靶点，为靶向干预提供理论依据。
对我们研究的思路启发
利用公共数据库进行多数据集整合分析，可以提高研究的广泛性与可信度。
机器学习在毒理学机制研究中的价值显著，值得在环境致癌因素研究中广泛应用。
生信分析不仅限于差异筛选，还可用于信号网络构建与关键调控节点识别。
SHAP等可解释性算法有助于生信预测结果向临床落地转化，提升模型可信度。
分子对接验证使预测结果更具说服力，未来可作为机制验证的重要补充。
AFB1研究框架可推广至其他化学暴露风险，如污染物、添加剂等，构建“毒物-信号-靶点”新模式。
对涉及化学暴露的流行病学队列研究，建议加强代谢基因（如PCK1）与免疫通路的联动分析。
推动跨学科整合（生信+毒理+结构生物学）将成为复杂疾病致病机制研究的重要趋势。
华哥微信：15623525389 公众号：Bioinformation
为尽量减少批次效应，GEO 数据下载、表达矩阵与注释文件处理，以及探针到基因的转换流程。脚本用于处理 GSE36376（人类肝癌）芯片数据。
# 设置系统语言为英文，有助于在报错时以英文显示，便于查找解决方案。Sys.setenv(LANGUAGE = "en")#设置字符变量默认不转换为因子类型，避免处理数据框时因子带来的问题。options(stringsAsFactors = FALSE)# 清空当前 R 会话的环境变量，确保脚本从干净环境开始运行。rm(list=ls())
# 用于下载和解析GEO数据库数据library(GEOquery)# 数据处理包（管道 %>%、filter、mutate 等）library(dplyr)# 整合多个数据处理包（包含 dplyr、ggplot2、readr 等）library(tidyverse)
##显示当前的工作目录。getwd()setwd("E:/duoz/q4/GSE36376/")
series_matrix.txt.gz 文件（你已事先下载），跳过平台注释信息（getGPL = FALSE）gset = getGEO(filename='GSE36376_series_matrix.txt.gz', destdir=".", getGPL = F)
exp = gset@assayData[["exprs"]]#提取表达矩阵（表达值），并将其转换为标准 data.frame 格式。df.exp = as.data.frame(exp)#检查原始表达矩阵与转换后的数据类型class(exp)class(df.exp)
pdata = gset@phenoData@datatable(pdata$source_name_ch1)table(pdata$`tissue:ch1`)
提取样本的临床信息（phenoData），并查看不同组织类型的样本数量。
source_name_ch1和tissue:ch1是样本注释信息字段。可能包含 "adjacent non-tumor liver"（癌旁组织）与 "liver tumor"（肿瘤组织）。
#将样本信息保存为CSV文件，并重新读取（备用）。write.csv(pdata, file = "GSE36376_pdata.csv")pdata1 = read.csv("GSE36376_pdata.csv")
library(data.table)## 使用 fread 快速读取 GPL 平台注释文件（多线程，高效）。anno = fread("GPL10558-50081.txt", sep = "\t", header = T, data.table = F)View(anno)colnames(anno)gpl_gene = anno[, c(1, 13)]colnames(gpl_gene)##将第二列重命名为 GeneSymbol，便于合并使用colnames(gpl_gene)[2] = "GeneSymbol"
exp.anno = merge(x = gpl_gene, y = df.exp, by.x = 1, by.y = 0)将表达矩阵与注释信息通过“探针 ID”进行合并：
by.x = 1→ 以gpl_gene的第一列（探针 ID）为连接键；by.y = 0→ 表示以df.exp的行名（探针 ID）为连接键。
x = exp.anno$GeneSymbolx[1:50]View(exp.anno)
提取前 50 个基因符号进行检查（是否缺失、是否重复等）。
使用 View 打开整个合并后的数据表。
exp1 = exp.annocolnames(exp1)rownames(exp1) = exp1$GeneSymbol
exp1 设置行名为 GeneSymbol（基因名）——但这一步会报错，因为基因名有重复。exp2.1 = distinct(exp1, GeneSymbol, .keep_all = F)exp2 = distinct(exp1, GeneSymbol, .keep_all = T)
使用 distinct() 去除重复基因名：
exp2.1：保留唯一GeneSymbol，但不保留其他列（结果基本是空）；exp2：保留第一个出现的探针对应的表达量（推荐使用这个）。
rownames(exp2) = exp2$GeneSymbolView(exp2)
## 去除前两列（探针 ID 与基因符号），只保留表达量数据。exp3 = exp2[, -c(1,2)]## 最终将标准化后的表达矩阵（每行为一个基因，每列为一个样本）保存为CSV文件，供下游分析使用write.csv(exp3, file = "GSE36376_exp.csv")
华哥微信：15623525389 公众号：Bioinformation
为尽量减少批次效应，我们整合了GSE25097和GSE36376两个数据集，并对基因表达矩阵进行了全面归一化处理。主成分分析（PCA）显示归一化后数据分布得到了优化，归一化数据集表现出更清晰的聚类模式（图3A-B）。
差异表达分析共识别出797个在HCC中显著改变的基因。这些表达变化通过火山图和热图可视化展示（图3C-D）。
Identification of HCC-related target genes. (A) PCA scatter plot shows distinct separation between GSE22097 and GSE38736 datasets before batch correction, indicating batch effects. (B) PCA scatter plot after batch correction shows the integration of GSE22097 and GSE38736 datasets, indicating reduced batch effects. (C) Volcano plot shows DEGs based on logFC and significance. Red dots are upregulated, green dots are downregulated, and gray dots are nonsignificant. (D) Heatmap shows DEG expression patterns across samples. Red indicates upregulation, blue indicates downregulation.
华哥微信：15623525389 公众号：Bioinformation
Acquisition of disease-related targets
Five HCC transcriptomic datasets (GSE25097, GSE36376, GSE14811, GSE54236, and GSE76427) were curated from the NCBI GEO database. GSE25097 and GSE36376 served as the discovery cohort, while the remaining three datasets (GSE14811, GSE54236, GSE76427) comprised the validation cohort. To mitigate batch effects, a multi-stage normalization pipeline was implemented: Surrogate Variable Analysis (SVA): Latent confounding factors in the discovery cohort were modeled and adjusted using the SVA package. ComBat Harmonization: Residual batch variations were further corrected via parametric empirical Bayes frameworks. Post-correction principal component analysis (PCA) demonstrated improved inter-batch sample clustering in reduced-dimensional space, confirming successful data harmonization. The complete analytical workflow is schematized in Fig 1.
Differential gene expression analysis
Transcriptomic data were analyzed using the limma package. Differentially expressed genes (DEGs) were identified with thresholds of FDR-adjusted P < 0.05and|log2FC| > 0.585(1.5 foldchange). Results were visualized via ggplot2.
华哥微信：15623525389 公众号：Bioinformation
为尽量减少批次效应
https://academic.oup.com/bioinformatics/article/28/6/882/311263?login=false摘要
异质性和潜在变量已被广泛认为是高通量实验中偏差和变异性的主要来源。最广为人知的潜在变异源是批次效应——即样本在不同时间、不同实验批次或由不同人员处理时产生的差异。然而，还有大量其他变量也可能对高通量测量结果产生重大影响。本文介绍了sva包，该包用于识别、估计和去除高通量实验中不需要的变异来源。sva包支持使用sva函数估计替代变量、使用ComBat函数直接调整已知批次效应，以及在预测问题中使用fsva函数调整批次和潜在变量的影响。
1 引言
高通量数据目前已广泛用于分子生物学中，主要用于：(i) 识别与结果相关的基因组特征，以及 (ii) 构建预测模型。但这些目标常常受到高通量数据中的潜在变量或不需要的异质性的影响。批次效应是基因组实验中最常被认识的潜在变量，其影响可能十分严重，甚至完全破坏生物学研究结果（Leek 等，2010）。此外，批次效应并非唯一可能影响研究统计和生物学有效性的潜在变异来源（Leek 和 Storey，2007）。
本文介绍了sva包，该包可用于识别并去除批次效应和其他不希望出现的变异来源。sva包包括两类方法：(i) 通过估计替代变量来去除未知变异来源（Leek 和 Storey，2007, 2008），以及 (ii) 使用ComBat直接移除已知批次效应（Johnson 等，2007）。已有研究表明，去除批次效应并使用替代变量可以减少依赖性、稳定误差率估计并提高结果的可重复性（Leek 和 Storey，2007, 2008；Leek 等，2010）。此外，sva包还包含目前唯一公开的函数fsva，用于在基因组/表观基因组预测问题中识别并去除潜在变量。
2 使用 SVA 包
2.1 数据格式
数据应格式化为一个矩阵，行表示特征（如转录本、基因、蛋白质），列表示样本。需使用model.matrix函数创建两个模型矩阵：一个是“空模型”（null model），包含所有已知的调整变量；另一个是“完整模型”（full model），包含空模型中的变量以及感兴趣的变量（即用于建模的结果变量或表型）。
2.2 使用 sva 函数估计并去除替代变量
sva函数是一个两步流程中的第一步，首先估计替代变量，然后在差异表达分析中将其去除。可以通过将高维数据矩阵（dat）、完整模型（mod）和空模型（mod0）传入sva函数来估计替代变量。输出结果是替代变量矩阵，可将其加入模型矩阵后传入sva包中的f.pvalue函数，计算调整替代变量后的参数F检验P值。
2.3 使用 ComBat 函数去除批次效应
ComBat函数基于经验贝叶斯框架调整已知批次效应（Johnson 等，2007）。此函数直接作用于数据矩阵，需传入未包含批次变量的完整模型矩阵，以及单独传入批次变量（batch）。输出为已去除批次效应的校正数据。后续可以对校正后的数据进行标准分析，也可再使用sva函数去除其他潜在变异来源。
2.4 使用 fsva 进行预测分析
在基因组预测任务中，数据集通常包括训练集和测试集。训练集中已知样本的表型/分类信息，但潜在变异来源未知；测试集则未知表型也未知潜在变量。大多数批次校正方法和替代变量估计方法设计用于群体研究，不能直接应用于预测任务。fsva函数（frozen surrogate variable analysis）允许在训练集、测试集乃至未来采集的单个样本中去除潜在变异，类似于标准化流程（McCall 等，2010）。
使用fsva函数需提供训练集的表达数据（dbdat）、模型矩阵（mod）、训练集生成的sva对象（svaobj），可选地还可传入测试集数据（newdat）。返回结果包括校正后的训练数据（db）和测试数据（new）。若未来有新样本加入，只需将其添加到newdat中即可进行替代变量校正。我们在一项膀胱癌基因表达研究中应用了fsva，结果显示其能提高预测准确性并改善测试样本聚类效果（见补充材料）。
3 讨论
我们介绍了sva包，其中包含广泛使用的ComBat函数，可用于去除批次及其他未建模的变异来源。同时介绍了首个用于预测问题中去除批次效应的函数。该包可通过 Bioconductor 免费获取，并可与常用的差异表达分析软件（如limma）兼容（Smyth, 2004）。
3.1 替代变量与直接调整的对比
sva的目标是尽可能去除所有不希望的变异来源，同时保留与主要变量相关的对比信号。这有助于识别在不同组间稳定变化的特征，屏蔽所有常见的潜在变异来源。
但在某些情形下，潜在变量本身可能是重要的生物学变异来源。如果分析目标是研究某些亚群体中的异质性，那么sva函数可能并不合适。例如，如果认为癌症样本包含两个未知但生物上显著不同的亚群体，而这两个亚群体对表达有显著影响，那么估计出的某些替代变量可能会与这些亚群高度相关（Teschendorff 等，2011）。这种情况无论使用主成分（PCA）、奇异向量（SVD）（Leek 和 Storey，2007, 2008）还是独立成分分析（ICA）（Teschendorff 等，2011）都可能发生。然而，去除与目标表型相关的替代变量可能会导致显著性分析的偏倚和不一致，特别是在这些未知变量与表型相关时（Leek 和 Storey，2007）。因此，是否应去除这些替代变量以改善推断结果，目前仍是未解决的科学问题。
相比之下，直接调整只去除已知批次变量的影响。批次效应是基因组实验中最广为人知的潜在变异来源（Leek 等，2010）。然而，除了批次变量之外，还有许多可能对基因组测量结果造成显著影响的变量，包括环境因素（Gibson，2008）和遗传变异（Brem 等，2002；Schadt 等，2003）。这些变量可能是研究的重点，但在许多研究中，目标是识别基因表达与某种表型的关联，这时环境和遗传变量往往未被测量或建模。如果忽视这些因素，它们也可能像批次效应一样遮蔽信号、降低统计效能，并对生物结论造成偏倚（Leek 和 Storey，2007）。
经验法则：当已知或未知的潜在混杂因素较多时，推荐使用替代变量调整；而当已知批次变量较明确、且生物分组本身存在异质性时，更适合采用直接调整方法。
华哥微信：15623525389 公众号：Bioinformation
1. 代码分析四件套---初始化、加载包、设置路径、加载数据
# 将系统语言设置为英文，方便调试与查看错误信息Sys.setenv(LANGUAGE = "en")# 禁止自动将字符串转为因子（factor）options(stringsAsFactors = FALSE)# 清空当前环境中的所有变量rm(list=ls())
## 2. 加载R包library(GEOquery) # 用于下载和读取GEO数据库中的数据library(data.table) # 提供更快的数据操作功能library(dplyr) # 数据整理语法library(tidyverse) # 数据科学套件，包含ggplot2、dplyr等library(ggplot2) # 可视化绘图包library(limma) # 微阵列和RNA-seq数据分析，含标准化与差异分析等功能library(patchwork) # 用于组合多个ggplot图形library(reshape2) # 数据重塑library(tools) # 提供文件处理功能library(ggpubr) # ggplot2的扩展，增强图表功能library(RColorBrewer) # 提供颜色方案
## 3. 设置工作路径getwd() # 显示当前工作目录setwd("E:/duoz/q4/") # 设置工作目录为数据所在路径
## 4. 读取数据# 读取表达矩阵exp_GSE25097 = read.csv("GSE25097_exp.csv", row.names = 1)# 读取对应的临床/分组信息phe_GSE25097 = read.csv("GSE25097_pdata.csv", row.names = 1)exp_GSE36376 = read.csv("GSE36376/GSE36376_exp.csv", row.names = 1)phe_GSE36376 = read.csv("GSE36376/GSE36376_pdata.csv", row.names = 1)
row.names=1 表示第一列为行名（通常是基因名）exp = exp_GSE25097quantiles = quantile(exp, c(0, 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE)needLog = (quantiles[5] > 100) || ((quantiles[6] - quantiles[1]) > 50 && quantiles[2] > 0)if (needLog) {exp[exp < 0] = 0exp = log2(exp + 1)}## 将小于0的值设为0，再进行 log2(x+1) 转换。
## 使用 limma 包的 normalizeBetweenArrays 函数进行跨样本标准化normalizedData = normalizeBetweenArrays(exp)exp_GSE25097 = normalizedData
## 取两个数据集中共有的基因交集。gene_inter = intersect(rownames(exp_GSE25097), rownames(exp_GSE36376))## 按行名（基因）合并两个数据集，得到合并后的表达矩阵。merge_data = merge(x = exp_GSE25097[gene_inter,], y = exp_GSE36376[gene_inter,], by = 0)rownames(merge_data) = merge_data$Row.namesmerge_data = merge_data[, -1]
library(sva)## 指定每个样本所属的批次；batchType = c(rep("GSE25097", 557), rep("GSE36376", 433))batchType <- as.factor(batchType)##删除NA；merge_data1 = na.omit(merge_data)## 使用 ComBat 方法进行批次效应校正。correctedData = ComBat(merge_data1, batchType, par.prior = TRUE)
setwd("E:/duoz/q6")write.csv(correctedData, file = "correctedData.csv")
华哥微信：15623525389 公众号：Bioinformation
系统解读这个图，并且重点描述什么是pca 以及pca在生信中的作用图 A：Before batch correction
标题：未进行批次矫正
横轴：PC1（主成分1）
纵轴：PC2（主成分2）红色圆点（GSE25097） 和 蓝绿色三角（GSE36376） 完全分离，分别堆积在PC1两侧，表示两个数据集在主成分空间中存在显著差异，说明受到了强烈的批次效应干扰。
图 B：After batch correction
标题：批次矫正后
同样是PC1和PC2主成分轴
此时两个数据集（红、蓝绿）已明显混合在一起，形成一个整体分布，椭圆代表95%置信椭圆区间。
说明批次效应被有效移除，两个数据集的表达谱变得更为一致。
什么是 PCA（主成分分析）？
PCA（Principal Component Analysis） 是一种降维算法，它将高维数据投影到几个主成分轴上：
每个主成分（PC）都是原始变量的线性组合；
第一主成分（PC1）是能解释最多数据变异度的方向；
第二主成分（PC2）是解释次多变异且与PC1正交的方向；
以此类推...
PCA 的本质是为了压缩信息、去噪声和可视化结构。
🔬 PCA 在生物信息学中的作用
在生信分析中，尤其是转录组、单细胞RNA测序、代谢组等高维数据中，PCA 有如下作用：
🧩 总结
这张图通过 PCA 清晰展示了 批次矫正前后两个数据集的融合情况：
批次矫正前：样本严重分离，PCA揭示了强烈批次效应；
批次矫正后：样本在PCA空间中充分混合，说明矫正成功。
华哥微信：15623525389 公众号：Bioinformation
1. 函数理解
huage=function(x,y){h=x+3h2=y+hh3=h2+xprint(h3)}huage(x=2,y=5)https://github.com/Aleksobrad/single-cell-rcc-pipeline/blob/master/singlecell_gex_viper_analysis.R## 941行代码到946行代码for (i in c(1,4,6,12,11)) {w1=huage(x=i,y=5)}i=1 w1=10
huage_pca_plot <- function(mat, batch_info, title, file=NULL){data_pca <- prcomp(t(mat), scale. = TRUE)pc_var <- data_pca$sdev^2 / sum(data_pca$sdev^2) # 方差解释比例pca_df <- data.frame(Sample = rownames(data_pca$x),PC1 = data_pca$x[,1],PC2 = data_pca$x[,2],Batch = factor(batch_info, levels=unique(batch_info)))n_batch <- length(unique(batch_info))mycol <- c("#e1615e","#26a3a3")if(n_batch > 8) mycol <- colorRampPalette(brewer.pal(9,"Set1"))(n_batch)pc1_lab <- paste0("PC1 (", round(pc_var[1]*100, 1),"%)")pc2_lab <- paste0("PC2 (", round(pc_var[2]*100, 1),"%)")p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch)) +geom_point(size=2, alpha=0.90 ) +stat_ellipse(aes(fill=Batch), geom="polygon", alpha=0.14, linetype=2, show.legend=FALSE) +scale_color_manual(values = mycol) +scale_fill_manual(values = mycol) +labs(title=title, x=pc1_lab, y=pc2_lab, color="Batch") +theme_bw(base_size=14) +theme(panel.grid.major = element_line(colour="gray90", linetype=2),legend.title = element_text(face="bold"),legend.background = element_rect(colour="black", fill=NA, size=0.25),plot.title = element_text(size=16, face="bold"),axis.title = element_text(face="bold"))if(!is.null(file)){ggsave(file, p, width=8, height=6)}return(p)}
这段 R 语言代码定义了一个名为 huage_pca_plot() 的函数，用于批次可视化（PCA图）绘制，结合 ggplot2 和 stat_ellipse() 展示不同批次样本在主成分空间（PC1-PC2）中的分布。
✅ 函数定义
huage_pca_plot <- function(mat, batch_info, title, file=NULL){mat：表达矩阵，行为基因、列为样本（需转置）。batch_info：每个样本对应的批次信息（向量），用于分组着色。title：图的标题。file：可选参数，若提供文件路径，将导出图片。
📊 PCA 计算部分
data_pca <- prcomp(t(mat), scale. = TRUE)对输入表达矩阵 转置后 进行 PCA。
scale. = TRUE：表示标准化每个变量（基因），使其具有相同权重。prcomp()返回一个包含主成分载荷、标准差等信息的对象。
pc_var <- data_pca$sdev^2 / sum(data_pca$sdev^2)计算每个主成分的方差贡献比例（即解释的变异百分比）。
pca_df <- data.frame(Sample = rownames(data_pca$x),PC1 = data_pca$x[,1],PC2 = data_pca$x[,2],Batch = factor(batch_info, levels=unique(batch_info)))
构建绘图数据框：
PC1,PC2为前两主成分；Batch为样本所属批次，用于分组。
🎨 颜色设置
n_batch <- length(unique(batch_info))mycol <- c("#e1615e","#26a3a3")if(n_batch > 8) mycol <- colorRampPalette(brewer.pal(9,"Set1"))(n_batch)
默认颜色为红、蓝绿两种，适用于两个批次；
若批次数大于8，则自动调用
Set1调色板生成更多颜色。
🏷️ 标签设置
pc1_lab <- paste0("PC1 (", round(pc_var[1]*100, 1),"%)")pc2_lab <- paste0("PC2 (", round(pc_var[2]*100, 1),"%)")
用来标记 x、y 轴，附带解释变异度百分比。
📈 主图构建
# 用 ggplot2 创建主成分图。p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch)) +
绘制每个样本的点。
stat_ellipse(aes(fill=Batch), geom="polygon", alpha=0.14, linetype=2, show.legend=FALSE)## 添加 95% 置信椭圆轮廓，辅助观察每个批次的分布范围；## 填充透明度为 0.14，线条为虚线scale_color_manual(values = mycol)scale_fill_manual(values = mycol)## 设置点和椭圆的颜色为指定色板labs(title=title, x=pc1_lab, y=pc2_lab, color="Batch")## 设置图标题、轴标题和图例标题。
🧑🎨 图形美化
theme_bw(base_size=14) +theme(panel.grid.major = element_line(colour="gray90", linetype=2),legend.title = element_text(face="bold"),legend.background = element_rect(colour="black", fill=NA, size=0.25),plot.title = element_text(size=16, face="bold"),axis.title = element_text(face="bold"))
使用简洁白底风格（
theme_bw()）；加粗图例标题与轴标题；
设置标题字体大小；
设置背景网格线为灰色虚线等。
💾 图像保存（可选）
if(!is.null(file)){ggsave(file, p, width=8, height=6)}
🔚 最终输出
return(p)## 返回 ggplot 图对象，便于直接在 R 中显示或进一步组合。
华哥微信：15623525389 公众号：Bioinformation
1. 代码分析四件套---初始化、加载包、设置路径、加载数据
Sys.setenv(LANGUAGE = "en") # 设置系统报错信息为英文，有助于搜索解决方案options(stringsAsFactors = FALSE) # 关闭数据框自动将字符串转为因子的行为rm(list=ls()) # 清空当前 R 环境，防止变量冲突library(GEOquery) # 用于从GEO数据库下载数据library(data.table) # 高效数据操作包library(dplyr) # 数据处理（管道式操作）library(tidyverse) # 一整套数据科学工具包library(ggplot2) # 绘图包library(limma) # 差异表达分析核心包library(patchwork) # 多图拼接library(reshape2) # 数据框重塑（melt/cast）library(tools) # 提供文件路径等辅助函数library(ggpubr) # 美化 ggplot 图形library(RColorBrewer) # 提供配色方案getwd() # 显示当前工作目录setwd("E:/duoz/q4/") # 设置工作目录（需手动修改）# 加载两个数据集的表型信息phe_GSE25097 = read.csv("GSE25097_pdata.csv", row.names = 1)phe_GSE36376 = read.csv("GSE36376/GSE36376_pdata.csv", row.names = 1)# 生成肿瘤和癌旁样本分组s1_n = filter(phe_GSE25097, source_name_ch1 == "non_tumor")s1_t = filter(phe_GSE25097, source_name_ch1 == "tumor")s2_n = filter(phe_GSE36376, tissue.ch1 == "adjacent non-tumor liver")s2_t = filter(phe_GSE36376, tissue.ch1 == "liver tumor")
读取表达矩阵 & 排序
setwd("E:/duoz/q6/") # 切换到已矫正表达矩阵目录correctedData = read.csv("correctedData.csv", row.names = 1)# 将样本列按“癌旁在前，肿瘤在后”排序col_index = c(rownames(s1_n), rownames(s2_n), rownames(s1_t), rownames(s2_t))correctedData = correctedData[ , col_index]
构建分组矩阵
group_list = c(rep('Nor', 436), rep('Tumor', 508)) # 创建分组信息group_list = factor(group_list, levels = c("Nor", "Tumor")) # 强制顺序design = model.matrix(~ group_list) # 构建设计矩阵colnames(design) = c("Nor", "Tumor")rownames(design) = colnames(correctedData)
差异分析
fit = lmFit(correctedData, design) # 线性模型拟合fit2 = eBayes(fit) # 经验贝叶斯调整allDiff = topTable(fit2, coef=2, adjust.method = "fdr", number = Inf) # 提取全部差异基因write.csv(allDiff, file = "allDiff.csv") # 保存所有结果
筛选差异表达基因（DEGs）
select.log2FC = abs(allDiff$logFC) > 0.585select.qval = allDiff$adj.P.Val < 0.05select.vec = select.log2FC & select.qvaldeg = allDiff[select.vec, ]write.csv(deg, file = "deg.csv")
华哥微信：15623525389 公众号：Bioinformation
1. 排序差异基因，选取上/下调 top25
ordered_DEGs <- allDiff[order(as.numeric(as.vector(allDiff$logFC))), ]对
allDiff差异分析结果按照 logFC 从小到大排序。as.vector()和as.numeric()排序时不会因因子/字符导致出错（保险做法）。上调的在最后，说明最后是 logFC 最大的基因。
ordered_gene_names <- rownames(ordered_DEGs)total_DEG_count <- length(ordered_gene_names)
提取排序后基因名列表。
记录总共排序了多少个基因（用于索引）。
选出最显著的上调和下调基因（各25个）
display_genes_num = 25selected_gene_set <- ordered_gene_names[c(1:display_genes_num,(total_DEG_count - display_genes_num + 1):total_DEG_count)]
选出 logFC 最小的前 25 个（下调）+ 最大的后 25 个（上调）基因。
用于在热图中展示明显表达变化的代表性基因。
提取这50个基因的表达量子矩阵
heatmap_exp <- correctedData[selected_gene_set, ]从批次矫正后的表达矩阵 correctedData 中取出这些基因的表达数据，构成 heatmap_exp。
重新确认样本列的顺序
col_index = c(rownames(s1_n), rownames(s2_n), rownames(s1_t), rownames(s2_t))重新整理样本顺序：癌旁样本在前，肿瘤样本在后。
顺序分别是 GSE25097 非肿瘤 + GSE36376 非肿瘤 + GSE25097 肿瘤 + GSE36376 肿瘤。
构造列注释，用于热图样本分组展示
annotation_col1 = data.frame(Type = c(rep("Nor", 436), rep("Tumor", 508)),Project = c(rep("GSE25097", 243), # 非瘤243rep("GSE36376", 193), # 非瘤193rep("GSE25097", 268), # 肿瘤268rep("GSE36376", 240) # 肿瘤240))
构建注释信息用于热图顶部标识：
Type表示样本是癌旁还是肿瘤。Project表示来自哪个 GEO 项目。
rownames(annotation_col1) = colnames(heatmap_exp)设置注释行名为表达矩阵的列名（即样本名），保证匹配。
用 pheatmap 绘制最终热图
library(pheatmap)pheatmap::pheatmap(heatmap_exp, # 50个DEG表达矩阵cluster_rows = TRUE, # 行聚类（基因聚类）cluster_cols = FALSE, # 不进行列聚类（样本顺序固定）annotation_col = annotation_col1, # 热图顶部注释show_colnames = FALSE, # 不显示样本名scale = "row", # 每行基因标准化（Z-score）color = colorRampPalette(c("blue", "white", "red"))(100) # 自定义颜色梯度（低-中-高表达）)
华哥微信：15623525389 公众号：Bioinformation
加载必要的绘图包
library(ggpubr) # 提供更优雅的 ggplot 增强函数library(ggplot2) # 核心绘图包library(ggrepel) # 防止文本标签重叠的增强功能包
准备火山图数据框
data <- allDiffdata$significant = "stable"
复制
allDiff（limma 差异分析结果）为新变量data。添加一列
significant，默认值为"stable"，表示非差异表达基因。
标记上调差异表达基因
data$significant[data$logFC >= 0.585 & data$adj.P.Val < 0.05] = "up"将 logFC ≥ 0.585 且 FDR < 0.05 的基因标记为 "up"（上调）
探索更高 fold change 情况
x = data$logFC >= 1.5 & data$adj.P.Val < 0.05table(x)x[1:20]
这是调试代码，用于查看是否有满足更严格阈值（logFC ≥ 1.5）的基因。
并查看这些结果的前20个。
标记下调差异表达基因
data$significant[data$logFC <= -0.585 & data$adj.P.Val < 0.05] = "down"table(data$significant)
将 logFC ≤ -0.585 且 FDR < 0.05 的基因标记为
"down"（下调）。输出各类别的基因数量：
up、down、stable。
调试：log10阈值转换计算（辅助理解）
log10(0.001) # = -3log10(0.00001) # = -5log10(0.0001) # = -4log10(0.05) # ≈ -1.3
这部分是用于判断 -log10(p) 阈值位置，帮助后面横线绘制
绘制火山图（Volcano Plot）
ggplot(data, aes(x = logFC, y = -1 * log10(adj.P.Val))) +xlim(-4, 3) + ylim(0, 310) + # 设置坐标轴范围geom_point(aes(color = significant), size = 0.8) + # 绘制散点图并按上下调着色theme_classic() + # 使用简洁主题（无背景框）scale_color_manual(values = c("green", "grey", "red")) + # 设置颜色：down=green, stable=grey, up=redgeom_hline(yintercept = 1.3, linetype = 4, size = 0.8) + # 横线：p = 0.05 对应 -log10 ≈ 1.3geom_vline(xintercept = c(-0.585, 0.585), linetype = 4, size = 0.8) + # 竖线：logFC 截断值theme(title = element_text(size = 12), text = element_text(size = 12)) +labs(x = "log2 fold change", y = "-log10(p_value)") + # 坐标轴标题theme_bw() # 黑白主题（背景更干净）
华哥微信：15623525389 公众号：Bioinformation
特别告知：版权属于华哥，所有单位或个人未经许可 禁止翻录和传播 违者我们将追究法律责任华哥：15623525389（华哥微信） 公众号：Bioinformation##系统报错改为英文Sys.setenv(LANGUAGE = "en")##禁止转化为因子options(stringsAsFactors = FALSE)##清空环境rm(list=ls())## 加载r包library(GEOquery)library(data.table)library(dplyr)library(tidyverse)#install.packages("ggplot2")library(ggplot2)## BiocManager::install("limma")library(limma)## 加载画图包library(patchwork)library(reshape2)library(tools)library(ggpubr)library(RColorBrewer)##设置路径getwd()###一定要设置自己的路径setwd("E:/duoz/q4/")exp_GSE25097=read.csv("GSE25097_exp.csv",row.names = 1)phe_GSE25097=read.csv("GSE25097_pdata.csv",row.names = 1)exp_GSE36376=read.csv("GSE36376/GSE36376_exp.csv",row.names = 1)phe_GSE36376=read.csv("GSE36376/GSE36376_pdata.csv",row.names = 1)##################################################################################### 数据归一化# 判断是否需要 log2 转换exp=exp_GSE25097quantiles = quantile(exp, c(0, 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE)needLog = (quantiles[5] > 100) || ((quantiles[6] - quantiles[1]) > 50 && quantiles[2] > 0)if (needLog) {exp[exp < 0] = 0 # 小于0的值设为0exp = log2(exp + 1) # 进行 log2 转换}# 对数据进行标准化处理normalizedData = normalizeBetweenArrays(exp)exp_GSE25097=normalizedDataexp=exp_GSE36376quantiles = quantile(exp, c(0, 0.25, 0.5, 0.75, 0.99, 1.0), na.rm = TRUE)needLog = (quantiles[5] > 100) || ((quantiles[6] - quantiles[1]) > 50 && quantiles[2] > 0)if (needLog) {exp[exp < 0] = 0 # 小于0的值设为0exp = log2(exp + 1) # 进行 log2 转换}# 对数据进行标准化处理normalizedData = normalizeBetweenArrays(exp)exp_GSE36376=normalizedData############################################################### 两个数据集合并gene_inter=intersect(rownames(exp_GSE25097),rownames(exp_GSE36376))exp_GSE25097=as.data.frame(exp_GSE25097)exp_GSE36376=as.data.frame(exp_GSE36376)merge_data=merge(x=exp_GSE25097[gene_inter,],y=exp_GSE36376[gene_inter,],by=0)rownames(merge_data)=merge_data$Row.namesmerge_data=merge_data[,-1]#################################################################### 批次去除library(sva)# 对合并后的数据进行批次效应校正?ComBatbatchType=c(rep("GSE25097",557),rep("GSE36376",433))batchType <- as.factor(batchType)merge_data1=na.omit(merge_data)correctedData = ComBat(merge_data1, batchType, par.prior = TRUE)getwd()setwd("E:/duoz/q6")write.csv(correctedData,file = "correctedData.csv")huage_pca_plot <- function(mat, batch_info, title, file=NULL){data_pca <- prcomp(t(mat), scale. = TRUE)pc_var <- data_pca$sdev^2 / sum(data_pca$sdev^2) # 方差解释比例pca_df <- data.frame(Sample = rownames(data_pca$x),PC1 = data_pca$x[,1],PC2 = data_pca$x[,2],Batch = factor(batch_info, levels=unique(batch_info)))n_batch <- length(unique(batch_info))mycol <- c("#e1615e","#26a3a3")if(n_batch > 8) mycol <- colorRampPalette(brewer.pal(9,"Set1"))(n_batch)pc1_lab <- paste0("PC1 (", round(pc_var[1]*100, 1),"%)")pc2_lab <- paste0("PC2 (", round(pc_var[2]*100, 1),"%)")p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch)) +geom_point(size=2, alpha=0.90 ) +stat_ellipse(aes(fill=Batch), geom="polygon", alpha=0.14, linetype=2, show.legend=FALSE) +scale_color_manual(values = mycol) +scale_fill_manual(values = mycol) +labs(title=title, x=pc1_lab, y=pc2_lab, color="Batch") +theme_bw(base_size=14) +theme(panel.grid.major = element_line(colour="gray90", linetype=2),legend.title = element_text(face="bold"),legend.background = element_rect(colour="black", fill=NA, size=0.25),plot.title = element_text(size=16, face="bold"),axis.title = element_text(face="bold"))if(!is.null(file)){ggsave(file, p, width=8, height=6)}return(p)}# 输出PCA PDF# 为PCA分析准备数据mode(correctedData) <- "numeric" # 确保是数值矩阵huage_pca_plot(merge_data1, batchType, "PCA_Before Batch Effect",file.path( "PCA_Before_Batch.pdf"))huage_pca_plot(correctedData, batchType, "PCA_After Batch Effect",file.path( "PCA_After_Batch.pdf"))
特别告知：版权属于华哥，所有单位或个人未经许可 禁止翻录和传播 违者我们将追究法律责任华哥：15623525389（华哥微信） 公众号：Bioinformation##系统报错改为英文Sys.setenv(LANGUAGE = "en")##禁止转化为因子options(stringsAsFactors = FALSE)##清空环境rm(list=ls())## 加载r包library(GEOquery)library(data.table)library(dplyr)library(tidyverse)#install.packages("ggplot2")library(ggplot2)## BiocManager::install("limma")library(limma)## 加载画图包library(patchwork)library(reshape2)library(tools)library(ggpubr)library(RColorBrewer)##设置路径getwd()###一定要设置自己的路径setwd("E:/duoz/q4/")phe_GSE25097=read.csv("GSE25097_pdata.csv",row.names = 1)table(phe_GSE25097$source_name_ch1)s1_n=filter(phe_GSE25097,source_name_ch1=="non_tumor")s1_t=filter(phe_GSE25097,source_name_ch1=="tumor")phe_GSE36376=read.csv("GSE36376/GSE36376_pdata.csv",row.names = 1)table(phe_GSE36376$tissue.ch1)s2_n=filter(phe_GSE36376,tissue.ch1=="adjacent non-tumor liver")s2_t=filter(phe_GSE36376,tissue.ch1=="liver tumor")setwd("E:/duoz/q6/")correctedData=read.csv("correctedData.csv",row.names = 1)### 排序 让所有癌旁在前 所有的肿瘤在后length(c(rownames(s1_n),rownames(s2_n)))length( c(rownames(s1_t),rownames(s2_t)) )col_index=c( rownames(s1_n),rownames(s2_n) , rownames(s1_t),rownames(s2_t) )correctedData=correctedData[ , col_index ]######################################################################################################################################################################################################################################这里我们创建分组信息，意思就是前436个是control，后面508 个是处理length(c(rownames(s1_n),rownames(s2_n)))length( c(rownames(s1_t),rownames(s2_t)) )group_list=c(rep('Nor',436),rep('Tumor',508))##转化为因子 这里可以不用理解 为了后面差异分析的矩阵构建## 强制限定顺序group_list <- factor(group_list, levels = c("Nor", "Tumor"))##构建差异分析的分组矩阵design=model.matrix(~ group_list)colnames(design) <- c("Nor", "Tumor")rownames(design) <- colnames(correctedData)##lmFit()：线性拟合模型构建fit=lmFit(correctedData,design)##eBayes()使用trend=TRUE对标准误差进行经验贝叶斯平滑，计算每个对比中每个基因的moderated t-statistic和log-odds。fit2 <- eBayes(fit)##topTable()给出一个最有可能在给定对比下差异表达的基因列表。allDiff=topTable(fit2,coef=2, adjust.method = "fdr", number = Inf)write.csv(allDiff,file = "allDiff.csv")####################################################################################################################################################################################################################################### Differential gene expression analysis#Transcriptomic data were analyzed using the limma package.#Differentially expressed genes (DEGs) were identified with thresholds of#FDR-adjusted P < 0.05and|log2FC| > 0.585(1.5 foldchange).#Results were visualized via ggplot2.select.log2FC <- abs(allDiff$logFC) > 0.585table(select.log2FC)select.qval <- allDiff$adj.P.Val < 0.05table(select.qval)select.vec= select.log2FC & select.qvaltable(select.vec)deg=allDiff[select.vec,]write.csv(deg,file = "deg.csv")####################################################################################################################################################################################################################################### 画 热图ordered_DEGs <- allDiff[order(as.numeric(as.vector(allDiff$logFC))), ]ordered_gene_names <- rownames(ordered_DEGs)total_DEG_count <- length(ordered_gene_names)### 文章的热图展示了上调25 下调25display_genes_num=25selected_gene_set <- ordered_gene_names[c(1:display_genes_num,(total_DEG_count - display_genes_num + 1):total_DEG_count)]heatmap_exp <- correctedData[selected_gene_set, ]col_index=c( rownames(s1_n),rownames(s2_n) , rownames(s1_t),rownames(s2_t) )##画个热图length( c( rownames(s1_n),rownames(s2_n)))length(c( rownames(s1_t),rownames(s2_t) ))annotation_col1 = data.frame(Type =c(rep("Nor",436),rep("Tumor",508 )) ,Project =c( rep("GSE25097",243), rep("GSE36376",193), rep("GSE25097",268), rep("GSE36376",240) ))rownames(annotation_col1)=colnames(heatmap_exp)#install.packages("pheatmap")library(pheatmap)pheatmap::pheatmap(heatmap_exp, #热图的数据cluster_rows = T,#行聚类cluster_cols =F,#列聚类，可以看出样本之间的区分度annotation_col =annotation_col1,show_colnames=F,scale = "row", #以行来标准化，这个功能很不错color =colorRampPalette(c("blue", "white","red"))(100))####################################################################################################################################################################################################################################### 画 火山图library(ggpubr)library(ggplot2)library(ggrepel)data <- allDiffdata$significant="stable"data$significant[data$logFC >=0.585 & data$adj.P.Val <0.05]="up"#data$significant[ ]="up"x=data$logFC>=1.5 & data$adj.P.Val <0.05table(x)x[1:20]data$significant[data$logFC<= -0.585 & data$adj.P.Val <0.05]="down"table(data$significant)log10(0.001)log10(0.00001)#ggplot( data,aes( x=logFC ,y=-1*log10(adj.P.Val) ) )log10(0.0001)log10(0.05)colnames(data)ggplot(data,aes(x=logFC,y=-1*log10(adj.P.Val)))+xlim(-4,3)+ylim(0,310)+geom_point(aes(color=significant),size=0.8)+theme_classic()+scale_color_manual(values = c("green","grey","red"))+geom_hline(yintercept = 1.3,linetype=4,size=0.8)+geom_vline(xintercept = c(-0.585,0.585),linetype=4,size=0.8)+theme(title=element_text(size = 12),text = element_text(size=12))+labs(x="log2 fold change ",y="-log10(p_value)")+theme_bw()