先差异后GSEA呢还是先ssGSEA后差异呢

# 1.构建表达矩阵 ----------------------------------------------------------------
dir.create("./data")

# 魔幻操作，一键清空
rm(list = ls()) 
options(stringsAsFactors = F)
# BiocManager::install('airway')
# 加载airway数据集并转换为表达矩阵
library(airway,quietly = T)
data(airway)
class(airway)

rawcount <- assay(airway)
colnames(rawcount)

# 查看表达谱
rawcount[1:4,1:4]

# 去除前的基因表达矩阵情况
dim(rawcount)

# 获取分组信息
group_list <- colData(airway)$dex
group_list

# 过滤在至少在75%的样本中都有表达的基因
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)

filter_count <- rawcount[keep,]
filter_count[1:4,1:4]
dim(filter_count)

# 加载edgeR包计算counts per millio(cpm) 表达矩阵,并对结果取log2值
library(edgeR)
express_cpm <- log2(cpm(filter_count)+1)
express_cpm[1:6,1:6]

# 保存表达矩阵和分组结果
save(filter_count,
     express_cpm,
     group_list,
     file = "./data/Step01-airwayData.Rdata") 

 
exprSet <- filter_count
dim(exprSet)
exprSet[1:6,1:6]
table(group_list)

# 加载包
library(DESeq2)

# 第一步，构建DESeq2的DESeq对象
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)

# 第二步，进行差异表达分析
dds2 <- DESeq(dds)

# 提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","trt","untrt"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)

# 去除差异分析结果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2) 
head(DEG_DESeq2)
library(AnnoProbe)
ag=annoGene(rownames(DEG_DESeq2),
         ID_type = 'ENSEMBL',species = 'human'
            )
head(ag)
DEG_DESeq2$ENSEMBL=rownames(DEG_DESeq2)
 
deg_anno=merge(ag,DEG_DESeq2,by='ENSEMBL')
head(deg_anno)
# 保存 
save(deg_anno, file = "./data/Step03-DESeq2_nrDEG.Rdata")

rm(list = ls())
options(stringsAsFactors = F)
load(file = "./data/Step03-DESeq2_nrDEG.Rdata") 

library(clusterProfiler)
library(org.Hs.eg.db)
# 转成ENTREZID
df <- bitr(unique(deg_anno$SYMBOL), fromType = "SYMBOL",
           toType = c( "ENTREZID"),
           OrgDb = org.Hs.eg.db) 
DEG=merge(deg_anno,df,by='SYMBOL')
colnames(DEG)
data_all_sort <- DEG %>%  #排序
  arrange(desc(log2FoldChange))
geneList = data_all_sort$log2FoldChange #把foldchange按照从大到小提取出来
names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(geneList) #排序好的基因序列，而且是entrezeID的形式

R.utils::setOption( "clusterProfiler.download.method",'auto' )
kkgsea <- gseKEGG(geneList     = geneList,
               organism     = 'hsa', 
               minGSSize    = 3,
               maxGSSize    = 5000,
               pvalueCutoff = 1,
               pAdjustMethod = "none" ) #进行gseKEGG富集分析 
tmp = kkgsea@result
save(kkgsea,file = 'gsea_kk.Rdata') 

rm(list = ls())
options(stringsAsFactors = F) 
load(file = "./data/Step01-airwayData.Rdata")
express_cpm[1:6,1:6]
table(group_list)

library(msigdbr)  #install.packages("msigdbr")
library(GSVA) 
library(GSEABase)
library(pheatmap)
library(limma) 

## msigdbr包提取下载  KEGG数据集
KEGG_df_all <-  msigdbr(species = "Homo sapiens", # Homo sapiens or Mus musculus
                        category = "C2",
                        subcategory = "CP:KEGG") 
colnames(KEGG_df_all)
KEGG_df <- dplyr::select(KEGG_df_all,gs_exact_source,gs_exact_source,ensembl_gene)  
kegg_list <- split(KEGG_df$ensembl_gene, KEGG_df$gs_exact_source)  
kegg_list
 
names(kegg_list)
ssgsea_mat <- gsva(expr=express_cpm, 
                 method='ssgsea',
                 gset.idx.list=kegg_list,  
                 verbose=T, 
                 parallel.sz = 4 )#调用所有核，行名变成通路名
ssgsea_mat[1:4,1:4]
library(limma)
g=factor( group_list ,levels = c('untrt','trt'))
design=model.matrix(~g)
fit=lmFit(ssgsea_mat,design)
fit=eBayes(fit)
## 上面是limma包用法的一种方式 
options(digits = 4) #设置全局的数字有效位数为4
#topTable(fit,coef=2,adjust='BH') 
ssgsea_deg = topTable(fit,coef=2,adjust='BH', n=Inf) 
save(ssgsea_deg,file = 'ssgsea_deg.Rdata') 

rm(list = ls())
options(stringsAsFactors = F) 
load(file = 'ssgsea_deg.Rdata')
load(file = 'gsea_kk.Rdata')

deg_gsea = kkgsea@result
gsea_deg = ssgsea_deg

ids = intersect(rownames(deg_gsea),
                rownames(ssgsea_deg)) #186条中有184条，换了个富集源，从40显著增加

# 02构建可视化所需的矩阵（相关性就两行值）
df=data.frame(
  ssgsea_deg=ssgsea_deg[ ids,"logFC"],
  deg_gsea=deg_gsea[ ids,"enrichmentScore"] #就是两列logFC的值（散点图）
)
deg_gsea=deg_gsea[ ids,]
df$Description=deg_gsea$Description  #加名字了的，可视化了的

library(ggstatsplot)
#BiocManager::install("ggpmisc")
library(ggpmisc)
p <- ggplot(df, aes(x=ssgsea_deg,
                    deg_gsea)) +    
  geom_point() +
  labs(x = "ssgsea_deg",
       y = "deg_gsea")+
  scale_color_manual(values=c("blue", "grey","red"))+
  #geom_vline(xintercept = c( -1.5,0,1.5),lty=4,col="#ec183b",lwd=0.8) +
  # geom_hline(yintercept =  c( -1.5,0,1.5),lty=4,col="#18a5ec",lwd=0.8) +
  #xlim(-5,5)+
  ylim(-1,1)+
  theme_ggstatsplot()+
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())
p
p+stat_poly_eq(aes(label=paste(..eq.label..,..adj.rr.label..,..p.value.label..,sep = "~~~~")),formula = y~x,parse=T,size=3.0)