GEO—芯片GSE3292 _pd 中无法找到分组信息—火山图中添加目标基因

原创

sheldor没耳朵

发布于 2024-10-07 09:56:19

690

发布于 2024-10-07 09:56:19

文章被收录于专栏：数据挖掘

GEO—芯片GSE3292 _pd 中无法找到分组信息—火山图中添加目标基因

最近做了一个需求，这里简单介绍下实现步骤

要求是在这张火山图中显示出目标基因（ELF3和CCNE2）的具体位置

1 数据导入

采用的数据集是GSE3292，2005年发表的基于芯片的转录组测序数据，按照常规方法导入即可

rm(list = ls())
options(timeout = 100000) 
options(scipen = 20)#不要以科学计数法表示

#传统下载方式
library(GEOquery)
#eSet = getGEO("GSE3292", destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#试试geoChina,只能下载2019年前的表达芯片数据
library(AnnoProbe)
eSet = geoChina("GSE3292") #选择性代替第8行
eSet = eSet[[1]] 
class(eSet)

#(1)提取表达矩阵exp
exp <- exprs(eSet)
#表达矩阵行列数，正常是几万行，列数=样本数，
dim(exp)
#二个要检查的地方
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
#可能要修改的地方
#exp = log2(exp+1) #需要log才log，不需要log要注释掉这一句
#第三个要检查的地方
boxplot(exp,las = 2) #看是否有异常样本

#(2)提取临床信息
pd <- pData(eSet)
#⭐多分组中提取两分组的代码示例，二分组不需要

#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
  s = intersect(rownames(pd),colnames(exp))
  exp = exp[,s]
  pd = pd[s,]
}

#(4)提取芯片平台编号
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")

2 分组信息

这是这个数据集比较特殊的地方，也是我写这篇帖子的原因。可以看到这个数据集pd中是不包含分组信息（HPV阳性和阴性）的。

在网页中找到分组信息，如下。可以看到分组信息对应的ID号是pd表格中title列中内容的后面的数字。所以我们要额外处理pd表格把对应的分组信息加到对应的GSM数据集后。

我采用的方法是直接复制上表内容，形成sup.tsv，然后读取到R中，按照ID值从小到大排列。处理pd的title列，将"UNC HNSCC01-0394"、 "UNC HNSCC02-0387"等的“-”去掉，再按照title列内容后面的数字，如010394、020387等进行从小到大排列。

可以看到现在sup.tsv和pd的行顺序是一致了，然后直接将sup.tsv的HPV列加到pd中即可。

值得注意的是原始pd的行顺序是不能改变的（上一步#(3)让exp列名与pd的行名顺序完全一致），因此我们要提前设置变量记录原始的pd行顺序，pd加完HPV列后，还需要按照这个变量重新排列，恢复原始顺序。

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
library(data.table)
tmp <- fread("sup.tsv",data.table = F)
tmp_sorted <- tmp[order(tmp$ID), ]
#记录原始pd的行顺序
pd$row_order <- 1:nrow(pd)
pd$title <- gsub("-", "", pd$title)
# 提取 title 列后面的数字并转为数值型
pd$numeric_part <- as.numeric(gsub(".*?(\\d+)$", "\\1", pd$title))
# 按照提取的数字排序数据框
pd_sorted <- pd[order(pd$numeric_part), ]
# 可选：删除临时列
pd_sorted$numeric_part <- NULL
pd_sorted$HPV <- tmp_sorted$HPV

pd_sorted <- pd_sorted[order(pd_sorted$row_order), ]
k = str_detect(pd_sorted$HPV,"Negative");table(k) #不在title就在pd的其他列
Group = ifelse(k,"Negative","Positive")

# 需要把Group转换成因子，并设置参考水平，指定levels
#对照组在前，处理组在后
Group = factor(Group,levels = c("Negative","Positive"))
Group

#⭐检查自己得到的分组是否正确
data.frame(pd_sorted$HPV,Group)
#2.探针注释的获取-----------------
#四种方法，方法1里找不到就从方法2找，以此类推。
#方法1 BioconductorR包(最常用)
#⭐要操作的地方
library(tinyarray)
gpl_number 
pkg_all[pkg_all$gpl==gpl_number,2]
#用上面找到的前缀替换下面所有的hgu133plus2，共5处
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db",ask = F,update = F)
library(hgu133plus2.db)
ls("package:hgu133plus2.db") 
ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
# 方法2 下载并读取GPL网页的表格文件，按列取子集
#library(tinyarray)
#get_gpl_txt(gpl_number) #获取表格文件的下载链接
# 接下来是复制网址去浏览器下载、放在工作目录下、读取、提取探针id和基因symbol(没有现成的需要拆分和转换)
# 注意:最终的数据ids只能有两列，第一列列名是probe_id,第二列列名是symbol,且都是字符型，否则后面代码要报错
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
save(exp,Group,ids,file = "step2output.Rdata")

3 火山图绘制

原图中并没有标出差异基因logFC的阈值，采用limma包分析表达矩阵，可以看到ELF3和CCNE2的logFC的值分别为0.48和0.74。若采用常规的logFC设置为1的话，则ELF3和CCNE2算不上差异基因。因此原火山图并非标注的是差异基因而是满足P<0.05的上调和下调基因，即logFC设置的是0。

rm(list = ls()) 
load(file = "step2output.Rdata")
#差异分析
library(limma)
design = model.matrix(~Group)
fit = lmFit(exp,design)
fit = eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)

#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg = mutate(deg,probe_id = rownames(deg))
#2.加上探针注释
ids = distinct(ids,symbol,.keep_all = T)
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
#⭐检查
nrow(deg) #如果行数为0就是你找的探针注释是错的。

#3.加change列,标记上下调基因
#阈值，可按需修改
logFC_t = 0
p_t = 0.05

k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)



#火山图
library(ggplot2)
library(ggrepel)
highlight_genes <- c("CCNE2", "ELF3")
ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) +
  geom_point(size=1.5, aes(color=change)) +
  scale_color_manual(values=c("#4CA8F8", "#4d4d4d","red"))+
  scale_y_continuous(limits = c(0, 14), breaks = seq(0, 14, by = 2))+
  labs(x = "log2(fold change)", y = "-log10(P.value)", 
       title = "Volcano plot\nGSE3292: Gene expression signature of HPV\nin head and neck squamous.\npositive vs negative, Padj<0.05")+
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme(
    panel.background = element_rect(fill = "white"),  # 设置背景为白色
    panel.grid.major = element_blank(),  # 去掉主网格线
    panel.grid.minor = element_blank(),  # 去掉次网格线
    axis.line = element_line(colour = "black"),  # 坐标轴线设置为黑色
    axis.title = element_text(size=12),
    axis.text = element_text(size=10),
    plot.title = element_text(hjust = 0.5, size=16, face="bold"),
    plot.title.position = "plot",
    legend.position = "none" 
  ) +
  geom_text_repel(data = subset(deg, symbol %in% highlight_genes), 
                  aes(label = symbol), 
                  vjust = -1, hjust = 1, size = 4, 
                  color = "black", 
                  nudge_y = 0.3, # 标签偏移
                  arrow = arrow(length = unit(0.01, "npc"),  # 调整短线长度
                                type = "closed",             # 设置封闭的箭头
                                angle = 15))  
ggsave("Figure-1.tif",height = 7.5,width = 8)
write.csv(deg, file = "deg_data.csv", row.names = FALSE)

注意

设置目标基因箭头指向的话，应该用ggrepel这个包，在绘图代码中添加geom_text_repel参数；

色号的确定，可是直接使用ishot截图工具，指针指向哪里，就会显示哪里的色号；

如何在图中添加p值阈值的水平线，和logFC阈值的竖直线；

geom_hline(yintercept = -log10(p_t), lty=4, col="black", linewidth=0.8):
#geom_hline() 绘制了一条水平线，用来标识P值的阈值。
#yintercept = -log10(p_t) 指定了这条线的y值。p_t 应该是P值的阈值，例如0.05。
#lty=4 指定了线型，4 表示虚线；col="black" 设置线的颜色为黑色，linewidth=0.8 控制线的粗细。

#geom_vline(xintercept = c(-logFC_t, logFC_t), lty=4, col="black", linewidth=0.8)

最后的效果如图