我们生信马拉松学员在学完转录组入门后，开始阅读一些文献做复现，他看到一个文献做GSEA分析，发现这个GSEA结果图很明显不是用的全部基因做的，问：只用统计学显著的差异基因做gsea也可以吗？这个文献的标题是《Obstructive Sleep Apnea Syndrome Exacerbates NASH Progression via Selective Autophagy-Mediated Eepd1 Degradation》，于 2024年9月发表在一个高分杂志Adv Sci (Weinh)（IF=14.3）上。
文献中的GSEA结果图如下：
图片底部的基因排序数只到了1200左右，如果是全部的基因应该在2w左右，中间的那个灰色区域的排序也可以看到明显的断层。
首先，我们来问一下人工智能大模型kimi：
https://kimi.moonshot.cn/ 问：gsea分析之前筛选统计学显著的差异基因有影响吗？
在进行GSEA（基因集富集分析）之前，是否需要筛选统计学显著的差异基因，以及这一筛选过程对分析结果的影响，是许多研究者关心的问题。根据最新的研究和分析工具的使用经验，以下是一些关键点：
1. 差异基因筛选对GSEA结果的影响
差异基因筛选方法的选择：在进行GSEA之前，是否需要筛选差异基因以及如何筛选，可能会对结果产生一定的影响。然而，研究表明，即使使用不同的差异基因筛选方法（如基于t统计量、对数变化率LogFC或结合LogFC和p值），GSEA的最终结果在生物学解释上仍然具有一致性。 生物学解释的一致性：尽管在差异基因的排序上可能会因筛选方法的不同而有所变化，但这些变化对GSEA的整体结果影响较小。换句话说，无论使用哪种预排序方法，生物学信号的解释通常保持一致。 2. GSEA分析的灵活性
无需严格筛选差异基因：GSEA的一个优势在于它可以直接使用基因的排序列表进行分析，而无需预先筛选出统计学显著的差异基因。这意味着GSEA可以利用完整的基因表达数据来评估基因集的富集情况，从而提供更全面的生物学信息。 单样本与成对样本方法的比较：研究还表明，单样本GSEA方法（如ssGSEA）能够提供更细致的生物学信号，尤其是在评估个体样本的基因集富集情况时。这与传统的成对样本GSEA方法相比，可能更能揭示样本间的生物学差异。 3. 建议
综合考虑统计学显著性和生物学意义：在进行GSEA分析时，虽然差异基因的筛选不是必须的，但研究者仍需注意平衡统计学显著性和生物学意义。过度依赖统计学阈值可能会忽略一些具有生物学意义的基因集。 使用可视化工具辅助分析：为了更好地解释GSEA结果，建议使用可视化工具（如ROC曲线、分布图等）来评估基因集在不同样本组中的表现，而不仅仅是依赖于统计学指标。 综上所述，虽然在GSEA分析之前筛选统计学显著的差异基因不是必须的，但这一过程可能会对结果的解释产生一定影响。研究者应根据具体的研究目标和数据特点，灵活选择分析方法，并综合考虑统计学显著性和生物学意义。
可以看到 kimi是不建议我们使用筛选后的基因做GSEA的，下面我们比较一下看看 全部基因 与 只用部分基因做对结果的影响吧！
看看 全部基因 与 只用部分基因做对结果的影响
由于 上述文献并没有提供数据可以供我们分析，这里我们使用来自另外一个文献的数据进行简单的探索：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE224354。
GSE224354_read_count_AT2.txt.gz
GSE224354_read_count_Macrophage.txt.gz
数据情况如下：
从具有AT2细胞特异性Cpt1a基因敲除（Cpt1a-KO）的小鼠以及对照组（Sftpc-CreERt2 +/-）小鼠的肺部分离出AT2细胞，分离时间点包括基线水平以及气管内给予LPS后24小时。为了探究AT2细胞中受损的长链脂肪酸氧化（LCFAO）是否调节肺泡巨噬细胞在急性肺损伤（ALI）中的反应，我们还在气管内给予LPS后4小时从小鼠中分离出肺泡巨噬细胞，分别来自Cpt1a-KO小鼠和对照小鼠，以及在基线水平从对照小鼠中分离出的肺泡巨噬细胞。将2只小鼠的肺泡巨噬细胞合并为一个样本进行转录组分析。
1、简单处理数据
读取 GSE224354_read_count_Macrophage.txt.gz
rm(list = ls())#清空当前的工作环境
options(stringsAsFactors = F)#不以因子变量读取
options(scipen = 20)#不以科学计数法显示
library(data.table)
library(tinyarray)
data <- data.table::fread("GSE224354_read_count_Macrophage.txt.gz", data.table = F)
data <- data[!duplicated(data$GeneID),]
mat <- data[,c( 2:ncol(data))]
rownames(mat) <- data[,1]
mat[1:4,1:4]
# 过滤低表达
keep_feature <- rowSums (mat > 1) > 1 ;table(keep_feature)
ensembl_matrix <- mat[keep_feature, ]
rownames(ensembl_matrix) <- rownames(mat)[keep_feature]
rownames(ensembl_matrix) <- gsub('[.][0-9]*','',rownames(ensembl_matrix))
ensembl_matrix[1:4,1:4]
转换 gene symbol
library(AnnoProbe)
head(rownames(ensembl_matrix))
ids=annoGene(rownames(ensembl_matrix),'ENSEMBL','mouse')
head(ids)
tail(sort(table(ids$biotypes)))
ids=ids[ids$biotypes=='protein_coding',]
ids=ids[!duplicated(ids$SYMBOL),]
ids=ids[!duplicated(ids$ENSEMBL),]
symbol_matrix= ensembl_matrix[match(ids$ENSEMBL,rownames(ensembl_matrix)),]
rownames(symbol_matrix) = ids$SYMBOL
colnames(symbol_matrix)
symbol_matrix[1:5,1:5]
获取样本分组信息：
gset = getGEO("GSE224354",filename = 'GSE224354_series_matrix.txt.gz', destdir = '.', getGPL = F)
pd = pData(gset )
colnames(pd)
pd$description
colnames(symbol_matrix)
pd=pd[match(colnames(symbol_matrix),pd$description), ]
pd$title
group_list=ifelse(grepl('baseline',pd$title ), 'baseline',pd$title )
group_list=ifelse(grepl('control',group_list), 'control',group_list )
group_list=ifelse(grepl('Cpt1a-KO',group_list), 'Cpt1a-KO',group_list )
group_list
table(group_list)
save(symbol_matrix,group_list,file = 'symbol_matrix.Rdata')
dat=log(edgeR::cpm(symbol_matrix)+1)
save(dat,group_list,file = 'step1-output.Rdata')
source('./codes/step2-check.R')
2、差异分析
这里三分组，我们就用 Cpt1a-KO vs baseline 做差异吧：
rm(list = ls()) ## 魔幻操作，一键清空~
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
getOption('timeout')
options(timeout=10000)
load('./symbol_matrix.Rdata')
symbol_matrix[1:4,1:4]
group_list
table(group_list)
gse_number <- 'Cpt1a-KO_vs_baseline'
gse_number
dir.create(gse_number)
setwd( gse_number )
getwd()
kp=group_list %in% c('Cpt1a-KO','baseline');table(kp)
symbol_matrix=symbol_matrix[,kp]
group_list=group_list[kp]
table(group_list)
group_list=ifelse(group_list=='baseline','control','case')
group_list = factor(group_list,levels = c('control','case' ))
save(symbol_matrix,group_list,file = 'symbol_matrix.Rdata')
source('../scripts/step2-qc-counts.R')
source('../scripts/step3-deg-deseq2.R')
source('../scripts/step3-deg-edgeR.R')
source('../scripts/step3-deg-limma-voom.R')
source('../scripts/step4-qc-for-deg.R') # 需要调整logFC和p值的阈值
3、GSEA 分析
使用全部基因做：
rm(list = ls())
library(clusterProfiler)
library(org.Mm.eg.db)
library(clusterProfiler)
load( file = 'DEG_deseq2.Rdata' )
colnames(DEG_deseq2)
nrDEG=DEG_deseq2[,c("log2FoldChange", "padj")]
head(nrDEG)
colnames(nrDEG)=c('logFC','padj')
gene <- bitr(rownames(nrDEG), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
head(gene)
head(nrDEG)
nrDEG = nrDEG[rownames(nrDEG) %in% gene$SYMBOL,]
nrDEG$ENTREZID = gene$ENTREZID[match(rownames(nrDEG) , gene$SYMBOL)]
# https://www.ncbi.nlm.nih.gov/gene/?term=SPARCL1
geneList=nrDEG$logFC
names(geneList)=nrDEG$ENTREZID
geneList=sort(geneList,decreasing = T)
head(geneList)
# 一本书，很多数据库，很多可视化
# ～～～这里需要替换物种～～～
# mmu
str(geneList)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'mmu',#按需替换
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.99,
verbose = FALSE)
tmp=kk_gse@result
kk=DOSE::setReadable(kk_gse, OrgDb='org.Mm.eg.db', keyType='ENTREZID')#按需替换
kk@result$Description=gsub(' - Mus musculus \\(house mouse\\)', '',kk@result$Description )
head(kk@result$Description)
tmp=kk@result
pro='test'
write.csv(kk@result,file = file.path(paste0(pro,'_kegg.gsea.csv')))
save(kk,file = file.path( 'gsea_kk.Rdata'))
同样的代码 ，使用差异基因拿到差异进对应的GSEA结果：
两次结果对比
rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
load('差异挑选后gsea/gsea_kk.Rdata')
kk_after = kk
after = kk@result[,c(2,3,4,6)]
load('原始的矩阵gsea/gsea_kk.Rdata')
kk_before = kk
before = kk@result[,c(3,4,6)]
ids=intersect(rownames(after),rownames(before));ids
tmp=cbind(before[ids,],after[ids,])
plot(tmp[,2],tmp[,6])
tmp$df = tmp[,2] - tmp[,6]
colnames(before)
df= data.frame(
before = before[ids,'enrichmentScore'],
after = after[ids,'enrichmentScore']
)
rownames(df) = ids
df$diff = df$before - df$after
cid='mmu04728'
p1=enrichplot::gseaplot2(kk_after, geneSetID = cid,title = paste0(kk_after[cid,"Description"],"：DEG")) ;p1
p2=enrichplot::gseaplot2(kk_before, geneSetID = cid,title = paste0(kk_before[cid,"Description"],"：All Gene") );p2
p1
p2
library(ggpubr)
p=ggscatter(df, x = "before", y = "after",
color = "black", shape = 21, size = 1, # Points color, shape and size
add = "reg.line", # Add regressin line
add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.sep = "\n")
)
p=p+geom_hline(yintercept = -1, linetype = "dashed", color = "blue") +
geom_hline(yintercept = 0, linetype = "dotted", color = "red") +
geom_hline(yintercept = 1, linetype = "longdash", color = "green") +
geom_vline(xintercept = -1, linetype = "dashed", color = "blue") +
geom_vline(xintercept = 0, linetype = "dotted", color = "red") +
geom_vline(xintercept = 1, linetype = "longdash", color = "green")
p
ggsave('cor-of-two-enrichmentScore.pdf')
结果如下：
大部分通路得到的打分结果趋势一致（我觉得这里跟使用FC值作为排序指标有关系，如果是其他的排序指标呢？），但也有一些通路在两次分析中得到的结果截然相反，如通路 mmu04728：
更多操作，可以观看生信技能树的视频号：
友情宣传：
