传统KEGG（通路富集分析）和GO（功能富集）分析时，如果富集到的同一通路下，既有上调差异基因，也有下调差异基因，那么这条通路总体的表现形式究竟是怎样？是被抑制还是激活？或者更直观点说，这条通路下的基因表达水平在实验处理后是上升了呢，还是下降了呢?由于没有采用有效的统计学手段去分析某条通路下的差异基因的总体变化趋势，这使得传统的富集分析结果无法回答这些问题。
想要回答这个问题，我们需要GSEA富集方法的结果。GSEA分是根据处理后的差异倍数值对基因进行从大到小排序, 用来表示基因在两组间的表达量变化趋势。排序之后的基因列表其顶部可看做是上调的差异基因，其底部是下调的差异基因。可用于判断某条通路在某组样本中是激活还是抑制！（下面继续给大家展示一个典型的GESA分析案例）
下载GSE174177数据集的表达量矩阵：https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE174177&format=file 下载成功后解压这个GSE174177_RAW.tar 成为文件夹
代码语言：javascript复制setwd('G:/编程/生信菜鸟团学徒练习/作业20')
fs=list.files('GSE174177_RAW', full.names = T)
fs
****读取counts.txt构建counts矩阵，进行样品的重命名和分组
代码语言：javascript复制library(data.table)
tmp= fread(fs[1],data.table = F)
View(tmp)
gid=fread(fs[1],data.table = F)[,1]
head(gid)
rawcount = do.call(cbind,
lapply(fs, function(x){
fread(x,data.table = F)[,7]
}))
rawcount[1:4,1:4]
View(rawcount)
****将新列名作为字符向量传递
代码语言：javascript复制colnames(rawcount) <- c("KD-1", "KD-2", "KD-3","control-1","control-2","control-3")
rownames(rawcount)<- tmp$Geneid
View(rawcount)
****基因ID转换
#由于本次使用的为gencode或ensembl的gtf与cdna文件，因此最后得到的为ensembl_id (gene_id)和 transcript_id，形式为：ENSMUSG00000000001.1 ，而我们下游常用gene symbol进行展示，因此还需要从gtf注释文件中获取ensembl_id 、transcript_id与gene symbol的对应关系文件。#参考此方法获取基因ID转化的对应文件：https://zhuanlan.zhihu.com/p/518137593?utm_id=0
代码语言：javascript复制rawcount=read.csv("rawcount.csv")
rawcount2=rawcount[,-1]
rownames(rawcount2)<- rawcount$X
View(rawcount2)
g2s <- fread('g2s_vm25_gencode.txt',header = F,data.table = F) #载入从gencode的gtf文件中提取的信息文件
colnames(g2s) <- c("geneid","symbol")
View(g2s)
symbol <- g2s[match(rownames(rawcount2),g2s$geneid),"symbol"] #匹配counts行名对应的symbol
table(duplicated(symbol)) #统计重复基因名
****使用aggregate根据symbol列中的相同基因进行合并
代码语言：javascript复制counts <- aggregate(rawcount2, by=list(symbol), FUN=sum)
View(counts)
library(tibble)
counts <- column_to_rownames(counts,'Group.1')
View(counts)
****差异分析
代码语言：javascript复制#加载包
library(DESeq2)
#第一步，构建DESeq2的DESeq对象
group_list = c(rep('KD',3),rep('control',3))
rawcount=counts
colData<-data.frame(row.names = colnames(rawcount),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = rawcount,colData = colData,design = ~ group_list)
#第二步，进行差异表达分析
dds2 <- DESeq(dds)
#提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","KD","control"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)
#去除差异分析结果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2)
DEG_DESeq2['Gapdh',]
****针对以上差异分析结果进行 GSEA分析
代码语言：javascript复制head(DEG_DESeq2)
geneList = DEG_DESeq2$log2FoldChange
names(geneList)= rownames(DEG_DESeq2)
geneList=sort(geneList,decreasing = T)
head(geneList)
library(gplots)
library(ggplot2)
library(clusterProfiler)
library(org.Ms.eg.db)#BiocManager::install('org.Ms.eg.db')
library(GSVA) # BiocManager::install('GSVA')
library(GSEABase)
#install.packages("msigdbr")
library(msigdbr)
all_gene_sets = msigdbr(species = "Mus musculus", # Homo sapiens or Mus musculus
category = "H" )
length(unique(table(all_gene_sets$gs_name)))
egmt <- GSEA(geneList, TERM2GENE= all_gene_sets[,c('gs_name','gene_symbol')] ,
minGSSize = 1,
pvalueCutoff = 0.99,
verbose=FALSE)
head(egmt)
egmt@result
gsea_results_df <- egmt@result
rownames(gsea_results_df)
write.csv(gsea_results_df,file = 'gsea_results_df.csv')
****可视化特定通路
代码语言：javascript复制library(enrichplot)
pdf("gsea_results.pdf", width = 12,height = 8)
gseaplot2(egmt,geneSetID = 'HALLMARK_TNFA_SIGNALING_VIA_NFKB',pvalue_table=T)
gseaplot2(egmt,geneSetID = 'HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION',pvalue_table=T)
gseaplot2(egmt,geneSetID = 'HALLMARK_PANCREAS_BETA_CELLS',pvalue_table=T)
dev.off()
