day 10 GSVA和CellChat

rm(list = ls())
library(Seurat)
library(GSVA)
library(clusterProfiler)
load("../day7/scRNA.Rdata")
seu.obj = scRNA

#文献和教程一般使用平均值来做GSVA
exp  =  AverageExpression(seu.obj)[[1]] #exp =  AggregateExpression(seu.obj)[[1]]
exp  =  as.matrix(exp)
exp  =  exp[rowSums(exp)>0,] ;exp[1:4,1:4]

## GSVA 官网下载gmt文件 并读取
h_df = read.gmt("h.all.v2023.2.Hs.symbols.gmt")[,c(2,1)]
h_list = unstack(h_df)
ES  =  gsva(exp, h_list)
#ES = gsva(gsvaParam(exp,h_list,maxDiff = T)) #针对R4.4.0
ES[1:4,1:4]

#可视化
library(pheatmap)
pheatmap(ES, scale = "row",angle_col = "45",
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50))

#添加分组信息
seu.obj$group = ifelse(seu.obj$orig.ident %in% paste0("sample",1:3),"treat","control")

exp  =  AverageExpression(seu.obj,group.by = c("ident","group"))[[1]]
exp  =  as.matrix(exp)
exp  =  exp[rowSums(exp)>0,] 

h_df = read.gmt("h.all.v2023.2.Hs.symbols.gmt")[,c(2,1)]
h_list = unstack(h_df)
ES  =  gsva(exp, h_list);ES[1:4,1:4]

#可视化
library(pheatmap)
pheatmap(ES, scale = "row",angle_col = "45",
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50),
         cluster_cols = F,
         gaps_col = seq(2,ncol(ES)-1,2)) 
         #是第二行开始空格，按两个两个空n到整个ES

rm(list = ls())
if(!require(CellChat))devtools::install_local("CellChat-master/")
if(!require(presto))devtools::install_github('immunogenomics/presto')
library(CellChat)
library(ggplot2)
library(Seurat)
library(ggalluvial)

load("../day5-6/sce.Rdata")
scRNA = sce
DimPlot(scRNA)
scRNA = subset(scRNA,downsample = 100) #抽样节省计算

str(CellChatDB.human,max.level = 1) #mouse就换
table(CellChatDB.human$interaction$annotation)
scRNA$samples = scRNA$orig.ident
cellchat <- createCellChat(scRNA,
                           group.by = "ident",
                           assay = "RNA") 

cellchat@DB <- subsetDB(CellChatDB.human, 
                        search = "Secreted Signaling") 

cellchat <- subsetData(cellchat) 
dim(cellchat@data.signaling) 

cellchat <- identifyOverExpressedGenes(cellchat) # 识别过表达基因
cellchat <- identifyOverExpressedInteractions(cellchat) # 识别配体-受体对
cellchat <- projectData(cellchat, PPI.human)#慢 将配体、受体投射到PPI网络
cellchat <- smoothData(cellchat, adj = PPI.human) #v2版本的
cellchat <- computeCommunProb(cellchat) #慢 推测细胞通讯网络
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

cellchat@netP$pathways
pathways.show <- "GALECTIN"
#1 hierarchy plot
groupSize <- as.numeric(table(cellchat@idents)) 
vertex.receiver = seq(1,nlevels(scRNA)/2);vertex.receiver
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "hierarchy", vertex.receiver = vertex.receiver, vertex.weight  = groupSize)  

#2 circle plot
par(mfrow = c(1,1), xpd=TRUE,mar = c(2, 2, 2, 2))
netVisual_aggregate(cellchat, signaling = pathways.show, 
                    layout = "circle", 
                    vertex.receiver = vertex.receiver,
                    vertex.weight  = groupSize)

#3 chord plot
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "chord", vertex.receiver = vertex.receiver, vertex.weight  = groupSize)

#4 heatmap
#纵坐标是发射端，横坐标是接收端，有颜色代表横纵坐标所指的两类细胞之间有通讯，颜色深浅代表通讯概率。右侧和上方的条形图是该行/列通讯概率之和
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds") 
#5 计算配体-受体对信号网络的贡献度
netAnalysis_contribution(cellchat, signaling = pathways.show) 



#6 热图展示 分析细胞在信号网络中角色：发送者、接收者、调解者和影响者。
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP") # the slot 'netP' means the inferred intercellular communication network of signaling pathways
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 12, height = 5, font.size = 10)


#7 气泡图
#可以分开，也可以合到一起 p<0.05的才会被画出来，颜色仍然是通讯概率，圈的大小是按照p值，p值越小圈越大
netVisual_bubble(cellchat, sources.use = 1, 
                 targets.use = 1:nlevels(scRNA), 
                 remove.isolate = FALSE) 

#从第一类细胞到全部细胞
netVisual_bubble(cellchat, sources.use = 1:nlevels(scRNA), 
                 targets.use = 1:nlevels(scRNA), 
                 remove.isolate = FALSE) #从全部细胞到全部细胞

library(NMF)
selectK(cellchat, pattern = "outgoing")

#选择合适的细胞通讯模式数量
#这两个指标都是评估聚类稳定性的，二者都突然下降的值对应的横坐标就是合适的聚类数，2到3突然下降就选2,3到4出现了下降就选3
nPatterns = 9 #根据上图选择的#identifyCommunicationPatterns函数识别通讯模式，并画出热图
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns) #传出
nPatterns = 6
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns) #传入

# 桑基图 用桑基图和气泡图展示每种细胞传入/传出的信号都是属于哪些通路的
p1 <- netAnalysis_river(cellchat, pattern = "outgoing")
p2 <- netAnalysis_river(cellchat, pattern = "incoming")
p1+p2
# 气泡图 气泡图的颜色是按照细胞类型来分配，大小按照每个通路对每个细胞类型的贡献程度分配
p3 <- netAnalysis_dot(cellchat, pattern = "outgoing",dot.size = 4)
p4 <- netAnalysis_dot(cellchat, pattern = "incoming",dot.size = 4)
p3+p4