数据下载及处理

页面介绍

GEO（Gene Expression Omnibus）是 NCBI 中公共高通量基因数据共享平台

Link：Home - GEO - NCBI

如下是进入搜索的一个关于甲状腺癌的文件

GEO Accession viewer

Platforms：（GPLxxxx）平台编号对应的检测技术体系，包含了探针序列与基因的对应关系，是后续进行探针 ID 转换为基因 Symbol/Ensembl ID 的必需文件，直接影响下游功能富集分析的准确性。
Samples：（GSMxxxx）样本列表，可以查看每一个样本信息
Relations：关联信息，可找到该项目的其他数据
Download family：
- SOFT formatted family file(s) 、 MINiML formatted family file(s)：包含原始信号和所有元数据，信息非常全，但格式比较复杂，适合进阶分析。
- Series Matrix File(s)：GEO 整理好的表达矩阵，TXT 格式，直接就能用。里面的行是基因或探针，列是样本，每个单元格是对应的表达量。
Supplementary file：作者上传的补充材料，未经过GEO统一处理的文件，是作者对原数据的个性化解读加工。

Platform里面的GPL注释表和文件是注释信息，探针ID 到基因。而Download family里面的文件是表达量数据，只能确定探针中的表达值。

需要Download full table…

GEO数据处理

下载Download family处的Series Matrix File(s)文件和Platform里的full table文件。

安装BiocManager包，再用BiocManager安装其他包

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82


BiocManager::install("GEOquery") # 下载GEOquery包
library(GEOquery)
library(tidyverse)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis') # 设置工作路径；注意只能是“/”
gset<-getGEO("GSE299988",destdir=".",AnnotGPL=F,getGPL=F)
# getGEO(GEO,GSEMatrix = TRUE,AnnotGPL = FALSE,destdir = ".",getGPL = TRUE)，这是getGEO()函数的使用格式。
# GSEMatrix指是否下载整理好的 Series Matrix 表达矩阵
# AnnotGPL指是否同时下载平台注释表（GPL）
# destdir指数据保存到本地的路径（默认为"."）
# getGPL指是否自动下载对应的 GPL 平台文件。

class(gset) # 读取的为list类型
gset[[1]] # 查看信息，GSE只有一个平台就是1，可能有列表中其他平台，要查看就是更改成gset[[2]]等等
exp<-exprs(gset[[1]])# 提取基因表达矩阵数据，由于没有载入注释文件，矩阵的行名还是探针名

# 1、标准化
boxplot(exp,outline=F,las=2) # 绘制表达量分布箱线图，快速检查数据的整体质量和标准化效果。在同一条直线上则是已经标准化了
# boxplot(x, outline, las) x为表达矩阵，outline为是否显示异常值点，las为控制x轴标签方向：1 = 水平，2 = 垂直
library(limma)
exp<-normalizeBetweenArrays(exp) # 手动标准化，再去查看箱图

# 2、归一化
range(exp) # 查看range范围，数值过大需要归一化，log2处理
exp<-log2(exp+1)

# 3、基因注释，GPL导入
gpl<-read.table("GPL21185-21174.txt",header=T,fill=T,sep="\t",comment.char="#",stringsAsFactors=F,quote='')
# header=T列名保留
# fill=T文件中某些行的列数不一致时，自动用空值填充缺失的列
# sep="\t"指定列之间的分隔符是制表符（Tab）
# comment.char="#"忽略文件中以 # 开头的注释行
# stringsAsFactors=F把字符型数据（如探针 ID、基因名）保留为字符串，而不是自动转换成因子
# quote=''告诉 R 文件中没有引号包裹的字符串，避免把基因名或 ID 中的特殊字符误判为引号
colnames(gpl) # 查看列名
ids<-gpl[,c("ID","GENE_SYMBOL")] # 选取查看的两列
colnames(ids)<-c("probe_id","symbol")
library(stringr)
ids$symbol<-trimws(str_split(ids$symbol,"//",simplify=T)[,3])# 如果是第3个是基因，就选择第3列
# stringr包是tidyverse生态里专门处理字符串的工具包
# str_split()是按指定分隔符拆分字符串
# simplify=T默认返回列表，设为TRUE后会直接返回矩阵
# trimws()去除提取出的第 3 个基因名两端的空白字符，确保基因名格式干净
ids<-ids[ids$symbol!='',] # 删除了空白探针
ids<-ids[ids$symbol!='---',]
ids<-ids[ids$symbol!='--- ',] # 具体情况具体格式删除
ids<-ids[ids$probe_id %in% rownames(exp),] # 保留ids中与exp匹配的行
exp2<-exp[ids$probe_id,] # 选取exp中与ids相互匹配的行
table(rownames(exp2)==ids$probe_id)
exp3<-cbind(ids,exp2) # 合并两者为数据框

# 4、基因去重（一个基因会有多个探针）
exp3<-aggregate( . ~symbol,exp3,max) # 合并同一基因的多探针表达量
# . ~ symbol去重的核心逻辑，“.”代表除了symbol列之外的所有列，~symbol代表按照symbol列进行分组
# max指取每个基因对应探针的最大表达量作为代表值，也可以选择mean平均
rownames(exp3)<-exp3$symbol
exp3<-exp3[,-c(1,2)] # 再删去合并的前两列

# 5、获取临床数据
library(Biobase) # pData是Biobase里面的函数
rt1<-pData(gset[[1]]) # 先前exprs函数是获取的gset中基因表达矩阵
# pData函数用来从 ExpressionSet 对象中提取 phenotype data（表型数据），也就是每个样本的临床信息，比如分组、性别、年龄、疾病分期等。
write.csv(rt1,file='GSE299988_cli.csv',row.names=T)
#（1）疾病分组
colnames(rt1)
group_list<-ifelse(rt1$`tissue types:ch1`=='Normal','Normal','Tumor')
# 将Normal与Tumor细胞分开
# 更多组别的分组：group_list<-ifelse(str_detect(rt1$`title`,'normal_adjace'),'Normal',ifelse(str_detect(rt1$`title`,'positiveLN_RAIrefractive'),'P_tumor','N_tumor'))
# str_detect可以模糊字符选择
rt1$group<-group_list
rt2<-rt1 %>% dplyr::select("group")
# %>%管道操作符，把左边的rt1作为右边函数的输入
# dplyr::select()用于从数据框中选择特定的列
# （2）分组排列
rt3<-rt2 %>% arrange(rt2) # 升序排列，Normal在前，Tumor在后
# rt3<-rt2 %>% arrange(desc(group)) Tumor在前，Normal在后

# 6、数据保存
# exp3与rt3的列行一一对应
exp3<-exp3[,rownames(rt3)]
identical(rownames(rt3),colnames(exp3)) # 返回为T确认一一对应
save(rt3,exp3,file="GSE299988.rda") # 保存数据
# 之后可以直接load("GSE299988.rda")进行使用

GEO非symbol基因名的处理

常见标识类型： Gene Symbol直接对应基因名的标识；Entrez ID是NCBI Entrez 数据库的基因编号；RefSeq Accession是NCBI RefSeq 数据库的转录本编号，如 NM_000546；ENSEMBL ID对应Ensembl 数据库的基因 / 转录本编号，如 ENSG00000141510；GenBank Accession（GB ACC）对应GenBank 数据库的核酸序列编号，如 U94788。

如GSE200539的GPL文件中不包含symbol文件

GEO Accession viewer

Platforms进入的网页下可能有的没有Download full table…，则需要下载table下方的Download family中的SOFT formatted family file(s)，下载为soft.gz文件

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34


library(GEOquery)
library(tidyverse)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE200539')
# 中间部分与上述GEO数据处理部分一致
# 基因注释，soft导入（没有Download full table…的情况下）
GPL<-getGEO(filename="GPL23159_family.soft.gz",destdir=".")
gpl<-GPL@dataTable@table

# 由于这个注释文件没有Gene symbol，有NM_序列，因此提取这一列
ids<-gpl[,c("ID","SPOT_ID.1")] # 选取查看的两列
colnames(ids)<-c("probe_id","symbol") # 重命名
library(stringr)
ids$symbol<-trimws(str_split(ids$symbol,"//",simplify=T)[,1]) # 选取//分割的第一列NM_序列
# 下载两个包：BiocManager::install(c("org.Hs.eg.db", "AnnotationDbi"))
# org.Hs.eg.db包是人类基因注释数据库包，Hs 代表 Homo sapiens（人），eg 代表 Entrez Gene。包含了人类基因的各种标识符（如 Entrez Gene ID、Gene Symbol、Ensembl ID 等）之间的映射关系
# AnnotationDbi包是所有注释数据库包的核心接口包
library(org.Hs.eg.db)
library(AnnotationDbi)
symbol_ids<-mapIds(x=org.Hs.eg.db,keys=ids$symbol,keytype="REFSEQ",column="SYMBOL",multivals="first")
# org.Hs.eg.db指定要使用的注释数据库，这里是人类基因注释库
# keytype告诉输入的类型
# column指定想要转换后的目标 ID 类型
# multivals="first" 处理 “一个输入 ID 对应多个输出 ID” 的情况。当一个 REFSEQ ID 对应多个基因符号时，只取第一个结果，避免返回列表，保证输出是一个向量
symbol_ids
ids$symbol<-symbol_ids # NM_序列转换为基因名
ids<-na.omit(ids) # 删除NA

# 注意：exp行名默认是字符，如果ids$symbol是纯数字，应当转变为字符格式ids$probe_id<-as.character(ids$probe_id)

exp2<-exp[ids$probe_id,] # 选取exp中与ids相互匹配的行
table(rownames(exp2)==ids$probe_id)
exp3<-cbind(ids,exp2) # 合并两者为数据框

# 之后与之前GEO数据处理部分一致

GEO的Affymetrix原始数据整理

上述的GEO处理基于Download family中的Series Matrix File(s)下载后处理，这是已经预处理过的表达矩阵数据。通常是经过背景校正、归一化和探针注释后的表格。无法验证或调整预处理过程（如归一化方法）。

因此有时选取Supplementary file(s)里面的RAW原始文件，是原始的 CEL 文件压缩包，包含芯片扫描得到的原始荧光信号数据。

如上述的GSE200539的Supplementary file里面的GSE200539_RAW.tar

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45


setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE200539/Affymetrix原始数据处理')
library(tidyverse)
# 下载：BiocManager::install("oligo")
# Bioconductor 生态里专门用来处理 Affymetrix 芯片原始数据的 R 包
library(oligo)
# 1、数据导入
FileName<-list.files(path='./GSE200539_RAW') # 注意提前解压
FileName
sample_id<-str_sub(FileName,1,10) # 重命名，此处选取前10个字符，即命名为如GSM6037488的样式
rt<-data.frame(FileName=FileName,row.names=sample_id,group=rep(c('CTL','IDMF','DRF','MMF'),each=3))
# 分组数据参考网站上的Sample的分组，共4组，每组3个样本
cel_files<-list.celfiles('./GSE200539_RAW',full.names=T,listGzipped=T)
# list.celfiles()是oligo包的核心函数，用于搜索指定目录下的所有 .CEL文件
# './GSE200539_RAW'当前目录下的GSE200539_RAW文件夹
# full.names=T返回完整的文件路径（包含文件夹名）
# listGzipped=T支持识别并列出被gzip压缩的 .CEL.gz文件
cel_files
cel_data<-read.celfiles(cel_files) # 通过cel_files获取的所有文件路径导入并读取

# 2、质控处理
# 下载：BiocManager::install('arrayQualityMetrics')
# arrayQualityMetrics包是Bioconductor中专门用于芯片数据质量控制的工具包，能自动生成详细的质控报告，包含多种图表来评估数据质量。
library(arrayQualityMetrics)
# 版本问题手动下载吧……然后根据指示下一堆依赖包
library(GEOquery)
expr<-exprs(cel_data) #读取数据为原始表达量的矩阵
colnames(expr)<-sample_id # 列名重命名
eset<-ExpressionSet(assayData=expr,phenoData=AnnotatedDataFrame(rt))
# ExpressionSet()是Biobase包的核心函数，用来整合表达数据、样本信息和注释信息
# AnnotatedDataFrame()把rt数据框转换成 ExpressionSet 要求的注释数据框格式
# 此处是转换为ExpressionSet用于质控分析
arrayQualityMetrics(eset,outdir='GSE200539_QulityMetrics',force=T,do.logtransform=T,intgroup='group')
# force=T如果输出文件夹已经存在，强制覆盖它，避免报错
# do.logtransform=T对表达数据进行对数转换，这是芯片数据分析的常规预处理步骤
# intgroup='group'指定用于分组的表型信息列（比如 rt 数据框里的 group 列）

#rma标准化,对原始 CEL 文件数据进行背景校正、归一化和汇总
eset_rma<-oligo::rma(cel_data)
class(eset_rma) # "ExpressionSet"
exp<-exprs(eset_rma)
colnames(exp)<-str_sub(colnames(exp),1,10) # 保留列名前面10个字符
boxplot(exp,outline=F,notch=T,las=2) # 箱图看结果

# 这里的exp与之前章节得到的直接读取的exp相同
# 之后部分是基因注释，与先前相同

GEO高通量数据

先前选取的都是GEO中的Expression profiling by array（芯片表达谱），通过荧光强度判断表达量。

而GEO中Expression profiling by high throughput sequencing（高通量测序表达谱，通常指 RNA-seq）将样本 RNA 反转录建库后，直接进行大规模测序，通过比对基因组计数 reads 数。

芯片表达谱是对RNA与探针结合的荧光进行测定，反映相对值。而高通量测序是对RNA基因片段数量进行测定，反映绝对值。

高通量数据除了可以富集分析外，还可做可变剪接、融合基因、新转录本发现。

选取GSE309139高通量数据为例，下载Series Matrix File(s)与GSE309139_raw_counts.txt.gz，原始的counts数据需解压

GEO Accession viewer

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62


setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE309139')
library(tidyverse)
exp<-data.table::fread("GSE309139_raw_counts.txt",header=T)
# 用data.table::fread()高效读取RNA-seq原始计数文件
# header=T 确保文件第一行（样本ID）作为列名
table(is.na(exp)) # 查看是否有缺失值
exp[is.na(exp)]<-0
table(is.na(exp)) # 查看是否有缺失值
exp<-column_to_rownames(exp,var='Gene')
# 把基因名那一列挪到最前面成为行名

# 1、设置分组信息
sample_id<-colnames(exp)
sample_id
group<-rep(c('L','H','LL','HL'),each=6) # 具体分组看具体数据
group
rt<-data.frame(row.names=sample_id,group=group)
identical(rownames(rt),colnames(exp))

# 2、DESeq2分析
# BiocManager::install("DESeq2")
library(DESeq2)
rt$group<-factor(rt$group) # 变成因子的格式
rt$group
# 构建DESeq2对象
dds<-DESeqDataSetFromMatrix(countData=exp,colData=rt,design=~group)
# design=~group根据group分组来比较差异
# dds把原始计数数据、样本分组信息、实验设计逻辑打包在一起，后续所有分析（标准化、差异检验）都基于这个对象

# 低表达基因过滤（自定义，可选）
filter<-rowSums(counts(dds))>1 # 筛选每行加起来大于1的基因
dds<-dds[filter,]
# 差异表达分析
dds<-DESeq(dds)
# 格式为c("因子名"，“实验组”，“对照组”)
res<-results(dds,contrast=c("group","H","L"))
# 提取差异表达结果
res1<-as.data.frame(res)
# log2FoldChange差异倍数的对数值：H 组表达量 / L 组表达量 的 log₂ 结果
# lfcSE是log2FC 的标准误，代表差异倍数的统计波动程度（值越小，log2FC 越可靠）
# pvalue越小越拒绝原假设H0，越显著性差异
# padj是校准后的P值，控制多检验的假阳性
# 结果筛选padj<0.05 & abs(log2FoldChange)>=1，选出差异表达
degs<-filter(res1,padj<0.05 & abs(log2FoldChange)>=1)
view(degs)

# 绘制离散图
plotDispEsts(dds)
# 离散图横轴是基因经过标准化后的平均表达量，纵轴离散度，衡量基因表达量的变异程度
# P值直方图，看P值显著基因多少
hist(res$padj)

# 3、count数据vst转化，这里得到exp1的数据就不是整数了，可以理解为之前章节的exp，之后就是一样的标准化操作等等
vst<-vst(dds,blind=F)
# vst()将RNA-seq的原始计数数据转换为适合可视化和下游分析的数值
# blind=F变换时考虑实验设计信息，不会忽略分组差异，这样能更真实地保留生物学差异
exp1<-assay(vst)
range(exp1)
exp1<-as.data.frame(exp1)
save(exp,exp1,rt,file='GSE309139.rda')

# 注意，这里得到exp1之后记得看看需要标准化，归一化那一系列操作，之后再进行热图绘制。

差异分析

差异分析：芯片数据——array limma ；高通量——count DESeq2

limma芯片数据差异分析

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/limma差异分析')
library(limma)
library(pheatmap)
library(tidyverse)
library(ggpubr)

# 数据读取，直接使用之前第一部分处理好的数据
load('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE299988/GSE299988.rda')
rt<-rt3
exp<-exp3
str(exp) # 数据全是字符型：$ GSM9051410: chr  "8.561352698" "5.812654031" "2.853925376" "13.89898258" ...
exp<-exp %>% mutate(across(everything(),as.numeric))
# 把数据框 exp 里所有列的数据类型都转换成数值型
str(exp) # 把数据转换成了数字型：$ GSM9051410: num  8.56 5.81 2.85 13.9 8.65 ...
identical(rownames(rt),colnames(exp))
# 确认行名列名一致

火山图

接上述limma差异分析

如下是下述代码的DEG结果数据框图：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69


group_list<-factor(rt$group,levels=c("Normal","Tumor"))
# 把临床数据的分组弄成因子形式
design<-model.matrix(~group_list)
# model.matrix()：把因子型的分组信息转化为统计模型可识别的数值矩阵
# 基于上面的 group_list，生成的设计矩阵会以 Normal 作为截距（基线），用一列数值来表示 Tumor 组与基线的差异
design
# 线性模型拟合
fit<-lmFit(exp,design)
# 贝叶斯检验
fit2<-eBayes(fit)
deg<-topTable(fit2,coef=2,number=Inf)
# 是差异表达基因结果数据框，包含了所有基因的差异分析结果
DEG=na.omit(deg) # 去除可能的NA值
write.table(DEG,file='DEG1.txt',sep="\t",row.names=T,col.names=NA,quote=F)
# quote=F：不在字符串两端添加引号，避免后续读取时出现多余的引号干扰

# 通常logFC设置为1，adj.P.Val<0.05，下述代码进行显著水平与正负相关性的区分
logFC_cutoff<-1
type1=(DEG$adj.P.Val<0.05)&(DEG$logFC < -logFC_cutoff)
type2=(DEG$adj.P.Val<0.05)&(DEG$logFC > logFC_cutoff)
DEG$change=ifelse(type1,"DOWN",ifelse(type2,"UP","NOT"))
write.table(DEG,file='DEG2.txt',sep="\t",row.names=T,col.names=NA,quote=F)
table(DEG$change) # 统计上述3中情况："DOWN","UP","NOT"的个数

library(ggplot2)
library(cowplot)
library(ggrepel)
DEG<-mutate(DEG,Gene_symbol=row.names(DEG)) # 把Gene_symbol弄到最后一列
UPDEG<-subset(DEG,change=="UP")
UPDEG_5<-top_n(x=UPDEG,n=-5,wt=P.Value)
# n = -5：- 表示按升序取最小的 5 个，+ 表示按降序取最大的 5 个
# wt = P.Value：按 P.Value 列来排序，P 值越小代表差异越显著
# 由此选取最小的5组就是最小的5个P，最显著的5个
DOWNDEG<-subset(DEG,change=="DOWN")
DOWNDEG_5<-top_n(x=DOWNDEG,n=-5,wt=P.Value)
p <- ggplot(data = DEG,
            aes(x = logFC,
                y = -log10(adj.P.Val))) +
  geom_point(alpha = 0.5, size = 4.5,
             aes(color = change)) +
  ylab("-log10(adj.P.Val)") +
  xlab("log2(Fold Change)") +
  scale_color_manual(values = c("#1F77B4", "grey", "#FF7F0E")) +
  geom_vline(xintercept = c(-1, 1), lty = 5, col = "black", linewidth = 0.8) +
  geom_hline(yintercept = -log10(0.05), lty = 4, col = "black", linewidth = 0.8) +
  theme_half_open() +
  # 标记最显著的5个上调基因
  geom_text_repel(data = UPDEG_5, 
                  aes(label = Gene_symbol), 
                  vjust = 1.5,
                  box.padding = unit(0.35, "lines"),
                  point.padding = unit(0.35, "lines"),
                  segment.color = 'grey50',
                  nudge_y = 0.1,
                  direction = "y") +
  # 标记最显著的5个下调基因
  geom_text_repel(data = DOWNDEG_5, 
                  aes(label = Gene_symbol), 
                  vjust = 1.5,
                  box.padding = unit(0.35, "lines"),
                  point.padding = unit(0.35, "lines"),
                  segment.color = 'grey50',
                  nudge_y = 0.1,
                  direction = "y")
p # 查看火山图
# 保存为pdf
pdf(file="volcano.pdf",width=10,height=6)
print(p)
dev.off()

热图

接上述火山图

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55


library(pheatmap)
identical(rownames(rt),colnames(exp))
DEG<-read.table("DEG2.txt",sep="\t",check.names=F,stringsAsFactor=F,header=T,row.names=1)
table(DEG$change) 
diff<-rownames(DEG)[DEG$change != "NOT"] # 选出仅有显著上调或下调的基因
exp_diff<-exp[diff,]
range(exp_diff)

# 补：从DEG中选出目标基因来做热图
DEG_target <- DEG[grepl(pattern = "CASP1|ZDHHC|NLRP3|IL1B|IL18|GSD", x = DEG$name, ignore.case = TRUE), ]
target<-exp[rownames(DEG_target),]
# 作图
# 绘制差异基因表达热图
map<-pheatmap(
  exp_diff,
  annotation_col = rt, # 列注释（样本分组信息）
  color = colorRampPalette(c("blue", "white", "red"))(100), # 颜色渐变（100个梯度）
  cluster_cols = F, # 不对样本（列）进行聚类
  cluster_rows = T,  # 对基因（行）进行聚类
  show_colnames = F, # 不显示样本名（避免拥挤）
  show_rownames = T, # 显示基因名
  border_color = NA, # 去掉单元格边框
  scale = 'row', # 按行标准化（让每个基因的表达在不同样本间可比较）
  treeheight_row = 20, # 行聚类树的高度（可选）
  fontsize_row = 8, # 行名字体大小（可选）
  fontsize_col = 8 # 列名字体大小（可选）
)
pdf(file='heatmap.pdf',width=6,height=20)
print(map)
dev.off()

# 特定展示（可选）
library(tidyverse)
diff=DEG[DEG$change != "NOT",]
up<-diff %>% top_n(25,logFC)
dw<-diff %>% top_n(-25,logFC)
all<-c(rownames(up),rownames(dw))
exp_diff2<-exp[all,]
map2<-pheatmap(
  exp_diff2,
  annotation_col = rt, # 列注释（样本分组信息）
  color = colorRampPalette(c("blue", "white", "red"))(100), # 颜色渐变（100个梯度）
  cluster_cols = F, # 不对样本（列）进行聚类
  cluster_rows = T,  # 对基因（行）进行聚类
  show_colnames = F, # 不显示样本名（避免拥挤）
  show_rownames = T, # 显示基因名
  border_color = NA, # 去掉单元格边框
  scale = 'row', # 按行标准化（让每个基因的表达在不同样本间可比较）
  treeheight_row = 20, # 行聚类树的高度（可选）
  fontsize_row = 8, # 行名字体大小（可选）
  fontsize_col = 8 # 列名字体大小（可选）
)
pdf(file='heatmap2.pdf',width=6,height=8)
print(map2)
dev.off()

高通量数据分析

先前GEO高通量数据处理时运行得到degs那一步，仅在得到res1差异结果后增加一步，将其进行txt文件导出，作为DEG结果，如下是得到degs的代码：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46


setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE309139')
library(tidyverse)
exp<-data.table::fread("GSE309139_raw_counts.txt",header=T)
# 用data.table::fread()高效读取RNA-seq原始计数文件
# header=T 确保文件第一行（样本ID）作为列名
table(is.na(exp)) # 查看是否有缺失值
exp[is.na(exp)]<-0
table(is.na(exp)) # 查看是否有缺失值
exp<-column_to_rownames(exp,var='Gene')
# 把基因名那一列挪到最前面成为行名

# 1、设置分组信息
sample_id<-colnames(exp)
sample_id
group<-rep(c('L','H','LL','HL'),each=6) # 具体分组看具体数据
group
rt<-data.frame(row.names=sample_id,group=group)
identical(rownames(rt),colnames(exp))

# 2、DESeq2分析
# BiocManager::install("DESeq2")
library(DESeq2)
rt$group<-factor(rt$group) # 变成因子的格式
rt$group
# 构建DESeq2对象
dds<-DESeqDataSetFromMatrix(countData=exp,colData=rt,design=~group)
# design=~group根据group分组来比较差异
# dds把原始计数数据、样本分组信息、实验设计逻辑打包在一起，后续所有分析（标准化、差异检验）都基于这个对象

# 低表达基因过滤（自定义，可选）
filter<-rowSums(counts(dds))>1 # 筛选每行加起来大于1的基因
dds<-dds[filter,]
# 差异表达分析
dds<-DESeq(dds)
# 格式为c("因子名"，“实验组”，“对照组”)
res<-results(dds,contrast=c("group","H","L"))
# 提取差异表达结果
res1<-as.data.frame(res)
# log2FoldChange差异倍数的对数值：H 组表达量 / L 组表达量 的 log₂ 结果
# lfcSE是log2FC 的标准误，代表差异倍数的统计波动程度（值越小，log2FC 越可靠）
# pvalue越小越拒绝原假设H0，越显著性差异
# padj是校准后的P值，控制多检验的假阳性
# 结果筛选padj<0.05 & abs(log2FoldChange)>=1，选出差异表达
write.table(res1,file='DEG1.txt',sep='\t',row.names=T,col.names=NA,quote=F)
degs<-filter(res1,padj<0.05 & abs(log2FoldChange)>=1)
view(degs)

这是最后选出的degs表格样式：

根据筛选出来的表格样式，使用导出的DEG1.txt进行火山图和热图绘制，与之前limma部分一致，注意纵轴名的相应更改即可。

如下是limma后进行的火山图与热图分析结果：

富集分析

在先前差异分析得到的差异表格基础上进行，先前差异分析已经注释好的DOWN、UP或NOT的数据，即DEG2.txt

GO 富集：关注基因的功能属性（这个基因干什么、在哪里、怎么参与过程）

KEGG 富集：关注基因所在的通路 / 代谢网络（这些基因一起参与哪条通路、哪条代谢链）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33


library(tidyverse)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/富集分析')
DEG<-read.table("DEG2.txt",sep="\t",check.names=F,stringsAsFactors=F,header=T,row.names=1)
table(DEG$change)
diff<-rownames(DEG)[DEG$change!='NOT']

# 将基因ID从Symbol转换为Entrez ID（富集分析主要用Entrez ID 或 Ensembl ID来索引基因，因此要进行转换）
gene_entrez <- bitr(diff, fromType ="SYMBOL", toType ="ENTREZID", OrgDb = org.Hs.eg.db)
# OrgDb = org.Hs.eg.db：使用人类（Homo sapiens）的注释数据库
cat("所有基因ID转换", nrow(gene_entrez),"/", length(diff)) # 看转换成功多少个
untrans <- setdiff(diff, gene_entrez$SYMBOL) # 两部分的差，得出哪些转换成功
untrans 
isValid <- "C9orf92" %in% keys(org.Hs.eg.db, keytype = "SYMBOL")
isValid # 查看T还是F，判断没有被转换的基因是否被收录在org.Hs.eg.db里面

# 使用mapIds函数转换
library(AnnotationDbi)
entrez_ids <- mapIds(
  x = org.Hs.eg.db,
  keys = untrans,       # 输入的基因别名列表
  keytype = "ALIAS",    # 明确指定输入类型为别名（ALIAS）
  column = "ENTREZID",  # 转换为ENTREZID作为中间步骤
  multiVals = "first"   # 多个匹配时取第一个
)
entrez_df <- data.frame(ENTREZID = entrez_ids, check.names = FALSE)
entrez_df <- na.omit(entrez_df)
entrez_df <- entrez_df %>% rownames_to_column(var = 'SYMBOL')
#合并
gene_entrez2 <- rbind(gene_entrez,entrez_df)

GO富集

上述完成ID转换后得到gene_entrez2，再进行GO富集

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21


GO <- enrichGO(gene = gene_entrez2$ENTREZID,
               OrgDb = org.Hs.eg.db,
               keyType ="ENTREZID",
               ont ="ALL", # "BP", "MF", "CC"或"ALL"
               pAdjustMethod ="BH",
               pvalueCutoff =0.05,
               qvalueCutoff =0.05,readable = T)
# GO 分为三个分支：BP（Biological Process，生物学过程）；MF（Molecular Function，分子功能）；CC（Cellular Component，细胞组分）；ALL 表示同时分析三个分支
# pAdjustMethod = "BH"：多重检验校正方法，BH（Benjamini-Hochberg）是最常用的方法，用于控制假阳性率
# readable = T：将结果中的Entrez ID自动转换为更易读的基因Symbol，方便解读
GO_res <- GO@result
write.table(GO_res,file="GO_res.txt",sep="\t",quote=F,row.names = F)
#简单作图
barplot(GO, showCategory = 20)
dotplot(GO, showCategory = 20)
#分类展示，按 GO 的三个分支（BP/MF/CC）分别展示最显著的富集条目
barplot(GO, showCategory = 5,,split="ONTOLOGY") +
  facet_grid(ONTOLOGY~., scale='free')
# showCategory = 5：每个分类下只展示最显著的5个GO条目
# split="ONTOLOGY"：按ONTOLOGY（即BP/MF/CC三个分支）对结果进行分组
# facet_grid(ONTOLOGY~.)：将分组后的结果分面展示，生成三个独立的子图

KEGG分析

也是紧接上述完成ID转换后得到gene_entrez2进行

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27


KEGG <- enrichKEGG(gene = gene_entrez2$ENTREZID,
                   organism     = 'hsa', # 物种：hsa = Homo sapiens（人类）
                   pvalueCutoff = 0.05, 
                   qvalueCutoff =0.05) 
KEGG_res <- KEGG@result
write.table(KEGG_res,file="KEGG_res.txt",sep="\t",quote=F,row.names = F)
#简单作图
barplot(KEGG, showCategory = 20)
dotplot(KEGG, showCategory = 20)
#通路图可视化
#BiocManager::install('pathview')
library(pathview)
DEG <- rownames_to_column(DEG,var = 'SYMBOL')
pgenes <- inner_join(DEG,gene_entrez2,by = 'SYMBOL')
pgenes <- dplyr::select(pgenes,ENTREZID,P.Value)

pvalue <- -log10(pgenes$P.Value)
gene_names <- as.character(pgenes$ENTREZID)

gene_data <- pvalue
names(gene_data) <- gene_names
#选择兴趣的通路
hsa04978 <- pathview(gene.data  = gene_data,
                     pathway.id = "hsa04978", #通路名字
                     species    = "hsa", #代表人类
                     out.suffix = "pathview",  
                     limit      = list(gene= 10, cpd= 1)) #对基因数量限制，cpd化合物数量

GSEA分析

GSEA分析（Gene Set Enrichment Analysis，基因集富集分析），与差异富集分析不同。

GSEA 不是看 “差异显著的基因” 富集到哪些通路 / 功能，而是看整个基因列表（不分差异与否） 中，某一 “基因集”（比如某条 KEGG 通路的所有基因）是否整体偏向于列表的顶端或底端，以此判断该基因集是否被显著激活 / 抑制。

GSEA官网：GSEA

如下是进入页面后的Database页面，右侧的 Human Collections 就是给人类研究用的 10 类基因集

H：特征基因集，把多个相似基因集合并，代表最核心、最清晰的生物学状态 / 过程，快速抓最关键的通路，结果干净、易解释。

选取H的右侧Gene Symbols进行下载，下载得到gmt文件

常规GSEA

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58


library(tidyverse)
library(limma)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/富集分析/GSEA分析')#设置工作路径
load('GSE299988.rda')
rt <- rt3
exp <- exp3
str(exp) 
exp <- exp %>% mutate(across(everything(), as.numeric))
str(exp)
#确认行名列名一致#确认行名列名一致
identical(rownames(rt),colnames(exp))
#组间差异分析
group_list <- factor(rt$group,levels = c("Normal","Tumor"))
group_list
design <- model.matrix(~group_list)
#比较矩阵命名
design
#线性模型拟合
fit <- lmFit(exp,design)
#贝叶斯检验
fit2 <- eBayes(fit)
deg <- topTable(fit2, coef = 2, number = Inf)
DEG = na.omit(deg) 
DEG <- rownames_to_column(DEG,var = 'Gene')

# 至此其实和差异分析的处理相类似

#整理输入数据
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
#转换ENTREZID（如果在GSEA中下载的是ENTREZID则需要转换，下载的Symbol则不需要）
# genelist <- bitr(DEG$Gene, fromType="SYMBOL",
#                  toType="ENTREZID", OrgDb='org.Hs.eg.db')
# DEG <- inner_join(DEG,genelist,by=c("Gene"="SYMBOL"))
geneList <- DEG[,2]
names(geneList) <- as.character(DEG[,'Gene'])
#排序
geneList <- sort(geneList, decreasing = TRUE) # sort默认是从小到大排列
geneList

#读取gmt文件
gmt <- read.gmt('h.all.v2026.1.Hs.symbols.gmt')
GSEA <- GSEA(geneList,
             TERM2GENE = gmt, 
             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             minGSSize = 20,
             maxGSSize = 500) 

res <- GSEA@result
write.csv(res,file = 'H_GSEA.csv')
# 选择一条通路进行展示
p <- gseaplot2(GSEA,geneSetID = 'HALLMARK_ALLOGRAFT_REJECTION',title = 'HALLMARK_ALLOGRAFT_REJECTION',color = 'red')
p
# 选择多条通路进行展示
p1 <- gseaplot2(GSEA,geneSetID = 1:3,pvalue_table = T, subplots = 1:2)
p1  

如下是富集分析得到的结果res：

enrichmentScore (ES)：原始富集分数

NES (Normalized ES)：标准化富集分数，大于0表示上调

rank：该基因集在所有结果中的排名

leading_edge：核心富集基因的标签。tags：核心基因占该基因集的比例；list：占全基因列表的比例；signal：富集信号强度

core_enrichment：核心富集基因列表

以下是分析结果图：

单基因GSEA

单基因GSEA分析可以明确一个基因表达的变化时通路的变化情况

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75


library(tidyverse)
library(limma)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/富集分析/GSEA分析')#设置工作路径
load('GSE299988.rda')
rt <- rt3
exp <- exp3
str(exp) 
exp <- exp %>% mutate(across(everything(), as.numeric))
str(exp)
#确认行名列名一致#确认行名列名一致
identical(rownames(rt),colnames(exp))
#提取 自定义基因,整理高低表达组
rt1 <- as.data.frame(t(exp['PTPRC',]))
# 选取的基因表达如果是大于这一组中的中位数，则认为是上调，完成单个基因的高表达或低表达区分
rt1$group <- ifelse(rt1$PTPRC > median(rt1$PTPRC),'up','down')
# 依照group进行排列，中间的'PTPRC'没用
rt1 <- arrange(rt1,'PTPRC',group)
identical(rownames(rt1),colnames(exp))
# 让exp的列名与rt1的行名一致
exp <- exp[,rownames(rt1)]
identical(rownames(rt1),colnames(exp))

# 组间差异分析
group_list <- factor(rt1$group,levels = c("down","up"))
group_list
design <- model.matrix(~group_list)
# 比较矩阵命名
design
# 线性模型拟合
fit <- lmFit(exp,design)
# 贝叶斯检验
fit2 <- eBayes(fit)
deg <- topTable(fit2, coef = 2, number = Inf)
DEG = na.omit(deg) 
DEG <- rownames_to_column(DEG,var = 'Gene')
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
# 转换ENTREZID
# genelist <- bitr(DEG$Gene, fromType="SYMBOL",
#                  toType="ENTREZID", OrgDb='org.Hs.eg.db')
# DEG <- inner_join(DEG,genelist,by=c("Gene"="SYMBOL"))
geneList <- DEG[,2]
names(geneList) <- as.character(DEG[,'Gene'])
# 排序
geneList <- sort(geneList, decreasing = TRUE)
geneList

# 读取gmt文件，去官网下载的C2基因集
gmt <- read.gmt('c2.cp.kegg_medicus.v2026.1.Hs.symbols.gmt')
GSEA <- GSEA(geneList,
             TERM2GENE = gmt,
             pvalueCutoff = 0.05, #改成1获得全部结果
             pAdjustMethod = "BH",
             minGSSize = 20,
             maxGSSize = 500) 
res <- GSEA@result
write.csv(res,file = 'KEGG_GSEA.csv')
p <- gseaplot2(GSEA,geneSetID = 'KEGG_MEDICUS_REFERENCE_COPI_VESICLE_FORMATION',
               title = 'KEGG_MEDICUS_REFERENCE_COPI_VESICLE_FORMATION',color = 'red')
p
p1 <- gseaplot2(GSEA,geneSetID = 1:3)
p1  

p2 <- gseaplot2(GSEA,
                geneSetID = c('KEGG_MEDICUS_REFERENCE_WNT_SIGNALING_PATHWAY',
                              'KEGG_MEDICUS_REFERENCE_ITGA_B_FAK_CDC42_SIGNALING_PATHWAY'))
p2

#3.是不是可以自创基因集？更改文件内容即可
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
gmt <- read.gmt('GOBP_MITOPHAGY.v2025.1.Hs.gmt')
# GO 生物学过程（GOBP）下的「线粒体自噬」功能模块

如下是单基因GSEA分析的作图结果：结果曲线为正即为该通路上调