完结撒花!开心今天发现了课程里的包更新带来的代码问题~
辛苦花花老师!!!
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)