批量的GSEA及基因表达热图可视化

# 修改镜像并下载BiocManager
options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
getOption("BioC_mirror")
getOption("repos")
BiocManager::install("airway")
BiocManager::install("org.Hs.eg.db")

library(airway)
data(airway) # 将数据集加载出来
class(airway) #查看数据集类型 (发现看不懂输出结果)
dim(airway) # 查看数据集的行列 (输出内容依然不明白)
View(airway)
airway@colData #查看数据集第一列内容

airway$SampleName;airway[[1]]
airway$cell;airway[[2]]

library(airway)
#Biocductor R包为三种：1.功能函数包2.数据包3.注释包（芯片基因之间的转换）
#此为中的一种，为数据包
data(airway)#加载数据
exprSet=assay(airway)#获取表达矩阵，默认airway获取表达矩阵就是assay，没有原因的
colnames(exprSet)#看表达矩阵的列名
dim(exprSet)#查看表达矩阵的维度
View(exprSet)
#设定分组信息
group_list=colData(airway)[,3]#得出分组信息
tmp=data.frame(group_list)#把group_list向量变为数据框tmp
row.names(tmp)=colnames(exprSet)
#把tmp的行名改为exprSet的列名

exprSet=exprSet[apply(exprSet,1,function(x)sum(x>1)>5),]
##分别对数据中每一行的数据进行一个什么运算，1代表行，2代表列

library("DESeq2") 
colData<-data.frame(row.names = colnames(exprSet),group_list=group_list)
#给每个样本进行标签化，设定每个样本的性质特点，分组的前期准备

dds <- DESeqDataSetFromMatrix(countData =exprSet, colData = colData, design = ~group_list)
#countData为表达矩阵，colData样本特点内涵分组信息，design 以某个种（group_list）样本特点设定分组
#keep <- rowSums(counts(dds)) >= 10#筛选表达量之和大于10的基因
#dds <- dds[keep,]
dds <- DESeq(dds)#四步矫正，处理数据
res <- results(dds, contrast=c("group_list","untrt","trt"))#按照"untrt","trt"比较得出差异分析结果
resOrdered<-res[order(res$padj),]#把res排序
head(resOrdered)
DEG=as.data.frame(resOrdered)#把差异分析结果变为数据框
DESeq2_DEG=na.omit(DEG)#删除差异分析中缺少值的结果
View(DESeq2_DEG)

head(DESeq2_DEG)
geneList = DESeq2_DEG$log2FoldChange
names(geneList)= rownames(DESeq2_DEG)

#按logFC值进行排序
geneList=sort(geneList,decreasing = T)
head(geneList)

library(clusterProfiler)
library(org.Hs.eg.db)
library(DOSE)
data(geneList)
head(geneList) #排序好的基因序列，而且是entrezeID的形式

R.utils::setOption( "clusterProfiler.download.method",'auto' )
kkgsea <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa', 
                  minGSSize    = 10,
                  maxGSSize    = 500,
                  pvalueCutoff = 1,
                  pAdjustMethod = "none" ) #进行gseKEGG富集分析 
kkgsea=setReadable(kkgsea,keyType = 'ENTREZID',OrgDb = 'org.Hs.eg.db')
tmp = kkgsea@result[,c(1:5,11)]
save(kkgsea,file = 'G:/编程/生信菜鸟团学徒练习/gsea_kk.Rdata')
View(kkgsea@result)

#可视化前12条通路的分析结果
library(enrichplot)
gseaplot2(kkgsea,geneSetID = head(kkgsea@result$ID,12))#可视化前12条通路
# 按第一个条目做二维码图，并显示p值
#color是enrichment score线的颜色
#base_sizexy轴标题字的大小
#ES_geom是enrichment score线用线或点显示
gseaplot2(kkgsea,5,color="red",pvalue_table = T,title="DNA replication",base_size=10,ES_geom="line")#可视化第5条信号通路
##当然也可以在一张图上展示多个条目
gseaplot2(kkgsea,1:3,color="red",pvalue_table = T,title=kkgsea@result$Description[1],base_size=10,ES_geom="line")
#山脊图
ridgeplot(kkgsea,fill = "pvalue",5)+scale_fill_continuous(type = "viridis")

# 其中 exprSet 是前面的转录组测序后的counts矩阵
# group_list 是矩阵里面的每个样品的分组信息
# up_kegg 是自己挑选好的通路 

library(pheatmap)
#对通路的里面的基因拿去热图可视化
pro = 'up'
print(pro)
dir.create('G:/编程/生信菜鸟团学徒练习/up_kk_gse_heatmap')
savDir <- 'G:/编程/生信菜鸟团学徒练习/up_kk_gse_heatmap/'

lapply(1:nrow(kkgsea), function(i){ 
  i=5 
  cg = strsplit(kkgsea@result$core_enrichment[i],'/')[[1]]
  cg
  length(cg)
  dat=log2(edgeR::cpm(exprSet)+1)
  dat[1:4,1:4] 
  table(group_list)
  pheatmap(dat,show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
  n=t(scale(t(dat))) # 'scale'可以对log-ratio数值进行归一化
  n[n>2]=2 
  n[n< -2]= -2
  #  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(group=group_list)
  rownames(ac)=colnames(n) 
  pheatmap(n,show_colnames =F,show_rownames = F,
           annotation_col=ac)
  
  heatmap.title = paste0('Genes of ', kkgsea@result$Description[i])
  filename =  paste0(savDir, pro,'_',
                     gsub('/','-',kkgsea@result$Description[i]),
                     '.pdf')
  if(length(cg) > 15){plot.height <- length(cg)/5}
  if(length(cg) %in% c(3:15)){plot.height <- length(cg)*0.5}
  if(length(cg) %in% c(1:2)){plot.height <- length(cg)*1.5}
  pheatmap(n,main = heatmap.title,
           show_colnames =F,show_rownames = F,
           annotation_col=ac,
           height = plot.height,
           filename = filename) 
})