生命大数据与R语言-5-差异基因火山图
XiaYun
2025-03-26 ~ 2025-04-10
😊一、安装并载入包
setwd("/Volumes/MacData/1_SynologyDrive/1_Job/03-课程-公选课-生命大数据与R语言——数据驱动的生命科学探索/文件/")
# 安装 BiocManager（如果尚未安装）
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager", quiet = TRUE)
}
# 静默检查并安装 CRAN 包
cran_pkgs <- c("ggplot2", "ggrepel", "dplyr")
for (pkg in cran_pkgs) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg, quiet = TRUE)
}
suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}
# 静默检查是否存在DESeq2，不存在就安装，存在则加载，这种组合在自动化脚本或者 R Markdown 文档中非常有用，尤其是当你希望加载多个包时，而不想让启动消息干扰输出内容。
if (!requireNamespace("DESeq2", quietly = TRUE)) {
BiocManager::install("DESeq2", quiet = TRUE, ask = FALSE)
}
suppressPackageStartupMessages(library(DESeq2))😊二、读取数据，过滤数据
#这里的75仅仅是我项目中的一个count的序号的标记，代码中还是使用count，便于函数的书写
Count75 <- readRDS("Count75.RDS")
print(Count75[1:5,1:5])## ENSEMBL SRR23104045 SRR23104046 SRR23104047 SRR23104048
## 1 ENSMUSG00000102628 0 0 0 0
## 2 ENSMUSG00000100595 0 0 1 0
## 3 ENSMUSG00000097426 0 0 0 0
## 4 ENSMUSG00000104478 0 0 0 0
## 5 ENSMUSG00000104385 0 0 1 0Meta75 <- readRDS("Meta75.RDS")
print(Meta75[1:5,1:2])## Run Group
## 1 SRR23104057 Control
## 2 SRR23104058 Control
## 3 SRR23104059 Control
## 4 SRR23104060 Control
## 5 SRR23104053 ICB#geneid_df1.mus <- readRDS("~/1_home_gaoy/ICBS_r/xy/geneid_df1.mus.RDS")
geneid_df1.mus <- readRDS("geneid_df1.mus.RDS")
count <- Count75
count.use <- dplyr::left_join(geneid_df1.mus,count,by="ENSEMBL")
rownames(count.use) <- count.use$SYMBOL
meta <- Meta75
meta5<-meta[meta$Group %in% c("ICB", "Combination"),]
#根据列名取需要的count列，例如count.use[,"SRR23104053"]
count5<-as.data.frame(count.use[,meta5$Run])
#判断是否一致，避免出错，一些同学会漏了meta5中的“5”，导致False
identical(colnames(count5),meta5$Run)## [1] TRUEmeta <- meta5
#count5 <- count5[rowSums(count5) >= 20, ]
CountFilter=function(count,meta){
Group<-meta$Group[!duplicated(meta$Group)]
#unique(meta$Group)
n1<-nrow(meta[meta$Group==Group[1],])
# length(which(meta$Group==Group[1]))
#n1 <- meta %>% dplyr::filter(., Group == Group[1]) %>% dim() %>% .[1]
n2<-nrow(meta[meta$Group==Group[2],])
n<-ifelse(n1<n2,n1,n2)
count = count[rowSums(count >=10) >= n,]
return(count)
}
count5<-CountFilter(count5,meta5)😊三、求取差异基因
library(DESeq2)
count <- count5
treatment <- "Combination"
Control <- "ICB"
dds <- DESeqDataSetFromMatrix(count, colData=meta, design= ~ Group)## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factorsdds <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingres <- results(dds, contrast=c("Group",treatment,Control))
res$SYMBOL<-row.names(res)
res.df <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)
res.df <- res.df[!is.na(res.df$padj), ]
# library (ggplot2)
# library(ggrepel)
# library(dplyr)
#(res.df$padj < 0.05) & abs(res.df$log2FoldChange) >=1
res.df$change <- ifelse((res.df$padj < 0.05) & (abs(res.df$log2FoldChange) >= 1),
ifelse(res.df$log2FoldChange> 1 ,'Up','Down'),
'Stable')分步解释：
ifelse(条件, 为真时的值, 为假时的值) 是 R 的向量化条件判断函数。
外层判断：
(res.df$padj < 0.05) & (abs(res.df$log2FoldChange) >= 1)👉 这是判断基因是否差异显著且表达变化足够大的标准
• 统计学显著性（FDR校正后）:
padj < 0.05：• 表达倍数差异至少是 2 倍（log2FoldChange ≥ 1 或 ≤ -1）:
abs(log2FoldChange) >= 1✅ 如果这两个条件都满足，就进入下一步判断。
内层判断：
ifelse(res.df$log2FoldChange > 1, 'Up', 'Down')在差异显著的前提下：
• 如果 log2FoldChange > 1 → 表达上调，标注为 “Up”
• 否则（即 ≤ -1）→ 表达下调，标注为 “Down”
否则（不显著或变化太小）：表示无显著表达变化，标注为 “Stable”
总结:
“Up” → 表达上调（显著＆倍数大于2）
“Down” → 表达下调（显著&倍数小于1/2）
“Stable” → 不显著或倍数变化小
😊四、绘制火山图
😞1-基础的差异基因火山图
ggplot( #设置数据
res.df,
aes(x = log2FoldChange,
y = -log10(padj),
colour=change)) +
geom_point() 
😔2-修正颜色-颜色为基本的红灰蓝
colorlist=c("red","grey","blue")
ggplot( #设置数据
res.df,
aes(x = log2FoldChange,
y = -log10(padj),
colour=change)) +
geom_point() +
scale_color_manual(values=colorlist,na.translate = FALSE )
😟3-增加辅助线
colorlist=c("red","grey","blue")
ggplot( #设置数据
res.df,
aes(x = log2FoldChange,
y = -log10(padj),
colour=change)) +
geom_point() +
scale_color_manual(values=colorlist,na.translate = FALSE )+
geom_vline(xintercept=c(-1,1),lty=3,col="grey",lwd=1) +
geom_hline(yintercept = -log10(0.05),lty=3,col="grey",lwd=1)
😕4-修改点的大小、透明度
colorlist=c("red","grey","blue")
ggplot( #设置数据
res.df,
aes(x = log2FoldChange,
y = -log10(padj),
colour=change)) +
geom_point(alpha=0.6, size=3) +
scale_color_manual(values=colorlist,na.translate = FALSE )+
geom_vline(xintercept=c(-1,1),lty=3,col="grey",lwd=1) +
geom_hline(yintercept = -log10(0.05),lty=3,col="grey",lwd=1)
🙁5-美化点的颜色，使用添加灰色系的红色和蓝色
colorlist=c("#0072B5", "grey", "#BC3C28")
ggplot( #设置数据
res.df,
aes(x = log2FoldChange,
y = -log10(padj),
colour=change)) +
geom_point(alpha=0.6, size=3) +
scale_color_manual(values=colorlist,na.translate = FALSE )+
geom_vline(xintercept=c(-1,1),lty=3,col="grey",lwd=1) +
geom_hline(yintercept = -log10(0.05),lty=3,col="grey",lwd=1)
😊6-进一步美化
# 自定义差异基因颜色标注
# https://colorhunt.co/
# 提取 res.df$change（表示基因是 Up / Down / Stable）中的 唯一值，并去除 NA，以确定当前结果中到底有哪些类型的变化。
uniq_change <- res.df$change[!duplicated(res.df$change)]
uniq_change <- uniq_change[!is.na(uniq_change)] # 去除 NA
# 如果三种状态都有（Up、Down、Stable），就给三种颜色：
# Down: #0072B5 蓝
# Stable: grey 灰
# Up: #BC3C28 红
# 如果只有 Up 和 Stable：灰 + 红
# 如果只有 Down 和 Stable：蓝 + 灰
if (length(uniq_change) == 3) {
colorlist <- c("#0072B5", "grey", "#BC3C28")
} else if ("Up" %in% uniq_change) {
colorlist <- c("grey", "#BC3C28")
} else {
colorlist <- c("#0072B5", "grey")
}
p <- ggplot( #设置数据
res.df,
aes(x = log2FoldChange,
y = -log10(padj),
colour=change)) +
geom_point(alpha=0.6, size=3) +
scale_color_manual(values=colorlist,na.translate = FALSE )+
# "#0072B5","grey","#BC3C28"
# "#546de5", "#d2dae2","#ff4757"
geom_vline(xintercept=c(-1,1),lty=3,col="grey",lwd=1) +
geom_hline(yintercept = -log10(0.05),lty=3,col="grey",lwd=1) +
# 坐标轴
labs(x="log2FoldChange",
y="-log10 (Adjusted P-Value)")+
theme_bw()+
theme(legend.position='bottom',
legend.title = element_blank(),
panel.background = element_blank(),
text=element_text(family="sans"),
axis.title = element_text(size = 14),
legend.text = element_text(size = 14),
axis.text = element_text(size = 10))
print(p)
😃😄😁😆🤩……
补充知识点：
1. %in% 是什么？
%in% 是一个 逻辑运算符，用于判断一个向量中的元素是否存在于另一个向量中。 x %in% y 表示：判断 x 中的每个元素是否“在” y 中，如果在，返回 TRUE，否则返回 FALSE。
genes <- c("TP53", "BRCA1", "EGFR", "MYC")
target_genes <- c("TP53", "EGFR")
genes %in% target_genes## [1] TRUE FALSE TRUE FALSE意思是：
-“TP53” 在 target_genes 里 → TRUE
-“BRCA1” 不在 → FALSE
-“EGFR” 在 → TRUE
-“MYC” 不在 → FALSE
总结
x %in% y 判断x中每个元素是否出现在y中
!x %in% y 判断x中哪些元素不在y中
df[df$col %in% vec, ] 筛选数据框中某留属于某集合的行
df[df$Gene %in% target_genes, ] 取出 df 中所有 Gene 列属于目标基因集 的行。
2. which的用法
- 向量中查找满足条件的位置
x <- c(2, 4, 6, 8, 10)
which(x > 5)## [1] 3 4 5返回：3 4 5 （表示 x 的第3、4、5个元素 > 5）
- 布尔向量
bool_vec <- c(TRUE, FALSE, TRUE, FALSE)
which(bool_vec)## [1] 1 3返回：1 3 （即 TRUE 的位置）
3.自定义函数CountFilter意义讲解
3.1 整体逻辑回顾
函数目的是对表达矩阵 count 进行过滤，保留在两个分组中有足够表达量的基因。
关键变量是：
-count：一个基因 × 样本 的表达矩阵。
-meta$Group：每个样本的分组信息。3.2 什么是 n？
函数中这几句是为了找出两个组中较小的样本数：
Group <- meta$Group[!duplicated(meta$Group)]
n1 <- nrow(meta[meta$Group == Group[1], ])
n2 <- nrow(meta[meta$Group == Group[2], ])
n <- ifelse(n1 < n2, n1, n2) # 取两组中较小的那个样本数量比如：
如果 Group1 有 5 个样本，Group2 有 4 个样本，那么 n = 4。进一步举例说明：
# 创建一个包含两组的meta数据框
meta0 <- data.frame(
SampleID = c("S1", "S2", "S3", "S4", "S5"),
Group = c("Control", "Control", "Control", "Case", "Case")
)
# 查看生成的meta数据框
print(meta0)## SampleID Group
## 1 S1 Control
## 2 S2 Control
## 3 S3 Control
## 4 S4 Case
## 5 S5 CaseGroup <- meta0$Group[!duplicated(meta0$Group)]
#这行的作用是获取不重复的 Group 值。结果是：
# [1] "Control" "Case"
#更清晰：
unique(meta0$Group)## [1] "Control" "Case"# n1 是 Control 组的样本数（3）
# n2 是 Case 组的样本数（2）
n1 <- nrow(meta0[meta0$Group == Group[1], ])
n2 <- nrow(meta0[meta0$Group == Group[2], ])
#这行是取较小的样本数。因为 2 < 3，所以结果是：# [1] 2
n <- ifelse(n1 < n2, n1, n2)
#
# 总结一下这个代码做了什么：
# 确定两组名称（比如 Control 和 Case）
# 计算每组的样本数
# 取两个组中较小的那个样本数（例如 2）
# 通常用于 **等量抽样分析3.3. 关键句解读：
count = count[rowSums(count >= 10) >= n, ]这句话的意思是：
-对于 count 矩阵的每一行（每个基因），看它在所有样本中有多少个表达量 >=10。
-rowSums(count >= 10) 计算的是每一行中表达量 ≥10 的样本个数。
-只有当这个数量 >= n（即在最小那组样本数中全部或更多样本中表达量都≥10）时，才保留这个基因。3.4. 为什么和 n 有关系？
因为如果一个基因在最小组的所有样本中都没有较高表达（比如都<10），就不能说明它在两个组中都有代表性表达，因此过滤掉。
举个例子：
-假设 n = 4（两组中较小的一组有 4 个样本）。
-那么一个基因在所有样本中，至少有 4 个样本表达量 ≥10，才能保留下来。
-如果它只在 3 个样本中表达量 ≥10，就说明它不够“普遍表达”。