学习小组-D10

完结撒花！开心今天发现了课程里的包更新带来的代码问题～

辛苦花花老师！！！

GSVA

#加载数据和R包
rm(list = ls())
library(Seurat)
library(GSVA)
library(limma)
library(clusterProfiler)
load("scRNA.Rdata")
seu.obj = scRNA
table(Idents(seu.obj))
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
h_df = read.gmt("h.all.v2023.2.Hs.symbols.gmt")[,c(2,1)]
h_list = unstack(h_df)
gsvapar <- gsvaParam(exp, h_list, maxDiff=TRUE)
ES <- gsva(gsvapar)
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 =  AggregateExpression(seu.obj)[[1]]
exp  =  as.matrix(exp)
exp  =  exp[rowSums(exp)>0,] 
exp[1:4,1:4]

h_df = read.gmt("h.all.v2023.2.Hs.symbols.gmt")[,c(2,1)]
h_list = unstack(h_df)
gsvapar <- gsvaParam(exp, h_list, maxDiff=TRUE)
ES <- gsva(gsvapar)
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))

Cellchat

#BiocManager::install("ComplexHeatmap")
if(!require(CellChat))devtools::install_local("CellChat-main/")

rm(list = ls())
if(!require(presto))devtools::install_github('immunogenomics/presto')
library(CellChat)
library(ggplot2)
library(Seurat)
library(ggalluvial)
load("../../day6/sce.Rdata")
table(Idents(sce))

#抽样（实战中不可以
scRNA = subset(scRNA,downsample = 100)
table(Idents(sce))

##CellChatDB.human,CellChatDB.mouse分别是人和小鼠的配受体数据库,根据物种按需修改
str(CellChatDB.human,max.level = 1)
table(CellChatDB.human$interaction$annotation)

##STRING数据库里高等级证据的相互作用关系组成的稀疏矩阵
class(PPI.human)
dim(PPI.human)

#table(as.numeric(PPI.human))

#构建cellchat对象
scRNA$samples = scRNA$orig.ident
cellchat <- createCellChat(scRNA,
                           group.by = "ident",
                           assay = "RNA")

cellchat@DB <- subsetDB(CellChatDB.human, 
                        search = "Secreted Signaling")
#也可以用全部的(cellchat@DB <- CellChatDB)

#细胞通讯分析
# 识别过表达基因
cellchat <- identifyOverExpressedGenes(cellchat)
# 识别配体-受体对
cellchat <- identifyOverExpressedInteractions(cellchat)

## The number of highly variable ligand-receptor pairs used for signaling inference is 739

# 将配体、受体投射到PPI网络
cellchat <- projectData(cellchat, PPI.human)#慢
## 推测细胞通讯网络
cellchat <- computeCommunProb(cellchat) #慢
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

cellchat@netP$pathways
pathways.show <- "GALECTIN" #画某一条具体的通路

#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)  

#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)

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

#heatmap
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
netAnalysis_contribution(cellchat, signaling = pathways.show)

# 热图-展示每一类细胞是什么角色
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)

#气泡图-显示所有的显著的配体-受体对
#可以分开，也可以合到一起
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")
nPatterns = 3
#根据上图选择的，嫌麻烦也可以用默认值5
#传出
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns)
selectK(cellchat, pattern = "incoming")
#nPatterns要根据上面的图来修改
#传入
nPatterns = 3
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns)

#画图展示
# 桑基图
netAnalysis_river(cellchat, pattern = "outgoing")
netAnalysis_river(cellchat, pattern = "incoming")
# 气泡图
netAnalysis_dot(cellchat, pattern = "outgoing",dot.size = 4)
netAnalysis_dot(cellchat, pattern = "incoming",dot.size = 4)