免费代码+过程｜复现IF 50+期刊（Figure 1）！

pythonic生物人2025年10月16日约 7 分钟

使用R语言实现一下Nature medcine中的Figure 1，数据代码来自文章s41591-023-02790-x。
数据自己下载，替换下文代码中的路径，例如r_data/PDAC2_data_results/data/proteomics/20230412_PDAC_PRO_exp.rds为PDAC2_data_results/data/proteomics/20230412_PDAC_PRO_exp.rds：http://www.genetictargets.com/PDAC2BOOKLET/PDAC2_data_results.zip。
复现图和原图存在部分差异，没给全代码或者使用了AI等工具做了二次处理，不过多解释，只是学习下。

✅Figure 1 b图

该图为3d PCA图，该主成分分析图基于4,787种可分析蛋白质的表达数据。

# 加载必要的R包
library(scatterplot3d)  # 用于绘制3D散点图[ 
library(tidyverse)     # 数据清洗和整理工具集

# 定义颜色方案，用于区分不同的样本类型
color.bin <- c("#00599F","#D01910")

#----------------------------------------------------------------------------------
#  Step 1: 导入数据
#----------------------------------------------------------------------------------
# 读取蛋白质组学表达数据（4787个蛋白 × 281个样本）
pro  <- readRDS("r_data/PDAC2_data_results/data/proteomics/20230412_PDAC_PRO_exp.rds") 
# 读取样本元数据
meta <- readRDS("r_data/PDAC2_data_results/data/proteomics/20230412_PDAC_PRO_meta.rds")

#  Step 2: 数据预处理
#----------------------------------------------------------------------------------
# 执行主成分分析（PCA），不对数据进行标准化（scale = F）
res.pca.comp <- prcomp(pro, scale = F)
#----------------------------------------------------------------------------------

# 提取前10个主成分，并准备绘图数据
plot.data <- as.data.frame(res.pca.comp$rotation[, 1:10])
plot.data <- plot.data %>% 
              mutate(ID=rownames(plot.data),                    # 保留样本ID
                     Type=meta$Type,                            # 添加样本类型信息
                     TypeColor=color.bin[as.numeric(as.factor(Type))])  # 根据类型分配颜色

#  Step 3: 画图
#----------------------------------------------------------------------------------
# 创建3D散点图，使用PC2、PC1、PC3三个主成分 
scatterplot3d(x = plot.data$PC2,     # x轴：第二主成分
              y = plot.data$PC1,     # y轴：第一主成分（注意：通常PC1在前）
              z = plot.data$PC3,     # z轴：第三主成分
              color = plot.data$TypeColor,  # 按样本类型着色
              pch = 16,              # 点形状：实心圆
              cex.symbols = 1,       # 点大小
              scale.y = 0.7,         # y轴缩放比例 
              angle = 45,            # 视角角度 
              xlab = "PC2",          # x轴标签
              ylab = "PC1",          # y轴标签  
              zlab = "PC3",          # z轴标签
              main="3D Scatter Plot of proteomics",  # 图形标题
              col.axis = "#444444",  # 坐标轴颜色
              col.grid = "#CCCCCC")  # 网格线颜色

# 添加图例[2](@ref)
legend("bottom", legend = levels(as.factor(meta$Type)),
      col =  color.bin,  pch = 16,
      inset = -0.15,     # 图例位置调整（负值表示向绘图外移动）
      xpd = TRUE,        # 允许在绘图区域外绘制
      horiz = TRUE)      # 水平排列图例项

图中每个点代表一个样本，红色圆点表示肿瘤组织样本，蓝色圆点表示其配对的癌旁正常组织样本。主成分1和主成分2下方的百分比数值，分别表示各主成分所能解释的数据集总体方差的比例。
分析结果显示，肿瘤样本与癌旁正常样本在主成分空间中呈现出明显的分离趋势，这表明两组样本在整体的蛋白质表达谱上存在显著差异。

✅Figure 1 c图

火山图，火山图展示了肿瘤组织与癌旁正常组织（TAC）之间的差异表达蛋白质情况。

# 加载必要的R包
library(limma)       # 差异表达分析
library(tidyverse)   # 数据清洗和整理工具集
library(ggpubr)      # 基于ggplot2的增强绘图
library(ggthemes)    # 提供额外的图形主题

#----------------------------------------------------------------------------------
#  Step 1: 导入数据
#----------------------------------------------------------------------------------
# 读取蛋白质组学表达数据（4787个蛋白 × 281个样本）
exp  <- readRDS("r_data/PDAC2_data_results/data/proteomics/20230412_PDAC_PRO_exp.rds") 
# 读取样本元数据
meta <- readRDS("r_data/PDAC2_data_results/data/proteomics/20230412_PDAC_PRO_meta.rds")

#----------------------------------------------------------------------------------
#  Step 2: 数据预处理和差异分析
#----------------------------------------------------------------------------------
# 准备元数据：创建对比因子变量
meta <- meta %>% mutate(contrast = as.factor(Type)) 

# 构建设计矩阵：使用无截距模型（~ 0 + contrast）
design <- model.matrix(~ 0 + contrast, data = meta)

# 使用limma进行线性模型拟合
fit <- lmFit(exp, design)

# 定义对比矩阵：比较Tumor vs Normal
contrast <- makeContrasts(Tumor_Normal = contrastT - contrastN, levels = design)

# 应用对比到拟合模型
fits <- contrasts.fit(fit, contrast)

# 应用经验贝叶斯平滑来稳定方差估计
ebFit <- eBayes(fits)

# 提取差异分析结果（使用FDR校正，返回所有基因）
limma.res <- topTable(ebFit, coef = "Tumor_Normal", adjust.method = 'fdr', number = Inf)

## 结果后处理
limma.res <- limma.res %>% 
  filter(!is.na(adj.P.Val)) %>%     # 移除P值为NA的行
  mutate(logP = -log10(adj.P.Val)) %>%   # 计算-log10(adj.P.Val)用于火山图
  mutate(tag = "Tumor -vs- Normal") %>%   # 添加对比标签
  mutate(Gene = ID)                      # 确保有基因名列

# 设置差异表达筛选阈值：FC:1.5 (logFC: ±0.58)，adj.p:0.05
# log2(1.5) ≈ 0.58
limma.res <- limma.res %>% mutate(group = case_when(
  (adj.P.Val < 0.05 & logFC > 0.58) ~ "up",      # 显著上调
  (adj.P.Val < 0.05 & logFC < -0.58) ~ "down",   # 显著下调
  .default = "not sig"# 不显著
))

#----------------------------------------------------------------------------------
#  Step 3: 绘制火山图
#----------------------------------------------------------------------------------
# 确保分组因子的水平顺序正确
limma.res <- limma.res %>% mutate(group = factor(group, levels = c("up", "down", "not sig")))

# 创建图例标签：显示筛选标准和上下调基因数量
my_label <- paste0("FC:1.5 ; AdjP:0.05 ; ",
"Up:", table(limma.res$group)[1], " ; ",
"Down:", table(limma.res$group)[2])

# 使用ggscatter绘制火山图
p <- ggscatter(limma.res,
               x = "logFC", y = "logP",          # x轴：log2倍变化，y轴：-log10(adj.P.Val)
               color = "group", size = 2,         # 按分组着色，点大小=2
               main = paste0("Tumor -vs- TAC"),   # 图形标题
               xlab = "log2FoldChange",          # x轴标签
               ylab = "-log10(adjusted P.value)", # y轴标签
               palette = c("#D01910", "#00599F", "#CCCCCC"),  # 颜色方案：上调红，下调蓝，不显著灰
               ylim = c(-1, 70), xlim = c(-8, 8)) +           # 坐标轴范围

# 应用基本主题
  theme_base() +

# 添加显著性阈值线
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "#222222") +
  geom_vline(xintercept = 0.58, linetype = "dashed", color = "#222222") +    # 上调阈值
  geom_vline(xintercept = -0.58, linetype = "dashed", color = "#222222") +   # 下调阈值

# 添加包含统计信息的副标题
  labs(subtitle = my_label) +

# 设置透明背景
  theme(plot.background = element_blank())

# 显示图形
p

与TACs相比，共有1,213种蛋白质在肿瘤中显著上调，864种蛋白质显著下调。
火山图的解读看之前的文章👉复现Nature正刊上3张图，换数据即可用！

✅Figure 1 d图

富集分析气泡图，该图展示了基于显著上调/下调蛋白质所识别的富集通路。

# 加载数据清洗和整理工具集
library(tidyverse)

# 导入富集分析数据
# 读取Excel文件，从第2个sheet的第2行开始读取（跳过标题行）
plot.data <- read.xlsx("r_data/PDAC2_data_results/data/Extended Data Table 3.xlsx", 
                       sheet = 2, startRow = 2)

# 定义要展示的特定通路/术语列表
# 包含GO术语（GO:开头）和KEGG通路（hsa:开头）的重要生物学通路
plot.pathway <- c("GO:0006730~one-carbon metabolic process","GO:0006888~ER to Golgi vesicle-mediated transport","hsa00020:Citrate cycle (TCA cycle)","hsa00071:Fatty acid degradation","hsa04062:Chemokine signaling pathway","hsa04066:HIF-1 signaling pathway","hsa04151:PI3K-Akt signaling pathway","hsa04512:ECM-receptor interaction","hsa04610:Complement and coagulation cascades","hsa04621:NOD-like receptor signaling pathway","hsa04666:Fc gamma R-mediated phagocytosis")

# 数据预处理和筛选
plot.data <- plot.data %>% 
  filter(Term %in% plot.pathway) %>%  # 筛选出指定的通路术语
  mutate(LogFDR = -log10(FDR))        # 计算FDR的负对数，用于统计显著性可视化

# 定义颜色方案
color.bin <- c("#00599F","#D01910")

# 使用ggscatter绘制气泡图
p <- ggscatter(plot.data, 
               x = "LogFDR",              # x轴：-log10(FDR)，表示富集显著性
               y = "Fold.Enrichment",     # y轴：富集倍数，表示富集程度
               color = "Type",            # 按类型着色（如上下调相关通路）
               main = "Enrichment of tumor/TAC protein",  # 图形标题
               size = "Ratio",            # 点大小由Ratio（富集基因比例）决定
               shape = 16,                # 点形状：实心圆
               label = plot.data$Term,    # 添加通路术语标签
               palette = color.bin) +     # 使用自定义颜色方案
  theme_base()                            # 应用基本主题

# 图形美化优化
p <- p + 
  scale_size(range = c(4, 20)) +         # 设置点大小范围（最小4，最大20）
  scale_x_continuous(limit = c(-10, 40)) + # 设置x轴范围从-10到40
  theme(plot.background = element_blank()) # 设置透明背景

# 显示图形
p