这篇文章中规中矩，应该是药企推出的免疫治疗药物最后做了个科研探索治疗机制。分析得比较浅，不过人家主要看疗效~
1.下载数据并QC聚类分群
首先下载数据 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE225689
前面的QC和降维聚类流程就不放了。 很简单的注释
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(SingleR)
library(celldex)
library(singleseqgset)
library(devtools)
getwd()
dir.create("3-cell")
setwd('3-cell/')
sce.all=readRDS( "../2-harmony/sce.all_int.rds")
sce.all
#文章中给的信息
#tumour cells (expressing EpCAM, PGR and TFF3—all markers of ECs21),
#immune cells (marked by CD45 (PTPRC) expression),
#cancer-associated fibroblasts (CAFs, marked by αSMA (ACTA2)) and
#endothelial cells (PECAM1 positive)
library(ggplot2)
genes_to_check = c('PTPRC', # immune
'EPCAM', 'KRT19', 'PROM1', 'ALDH1A1', # tumor
'PECAM1','VWF', #endothelial
'ACTA2','FGF7' ) #CAFs
library(stringr)
p_paper_markers <- DotPlot(sce.all, features = genes_to_check,
assay='RNA' ) coord_flip()
p_paper_markers
ggsave(plot=p_paper_markers,
filename="check_paper_marker_by_seurat_cluster.pdf",width = 12)
p_umap=DimPlot(sce.all, reduction = "umap",
group.by = "RNA_snn_res.0.1",label = T)
library(patchwork)
p_paper_markers p_umap
#ggsave('markers_umap.pdf',width = 15)
DimPlot(sce.all, reduction = "umap",split.by = 'orig.ident',
group.by = "RNA_snn_res.0.8",label = T)
#ggsave('orig.ident_umap.pdf',width = 45)
2.注释亚群
代码语言：javascript复制# 需要自行看图，定细胞亚群：
celltype=data.frame(ClusterID=0:9 ,
celltype= 0:9)
#定义细胞亚群
celltype[celltype$ClusterID %in% c( 1,3,4,9),2]='Immune'
celltype[celltype$ClusterID %in% c( 5),2]='Endothelial'
celltype[celltype$ClusterID %in% c( 7),2]='CAFs'
celltype[celltype$ClusterID %in% c( 0,6,8 ),2]='Tumor'
celltype[celltype$ClusterID %in% c( 2 ),2]='Unknown'
head(celltype)
celltype
table(celltype$celltype)
sce.all@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
sce.all@meta.data[which(sce.all@meta.data$RNA_snn_res.0.1 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
table(sce.all@meta.data$celltype)
th=theme(axis.text.x = element_text(angle = 45,
vjust = 0.5, hjust=0.5))
library(patchwork)
p_all_markers=DotPlot(sce.all, features = genes_to_check,
assay='RNA' ,group.by = 'celltype' ) coord_flip() th
p_umap=DimPlot(sce.all, reduction = "umap", group.by = "celltype",label = T)
p_all_markers p_umap
#ggsave('markers_umap_by_celltype.pdf',width = 12,height = 8)
sce.all
table(Idents(sce.all))
Idents(sce.all)=sce.all$celltype
table(Idents(sce.all))
挺有趣的，看来每次还是得看一下所有的marker表达，不然出现unknown就尴尬。
col =c("#3176B7","#F78000","#3FA116","#CE2820","#9265C1",
"#885649","#DD76C5","#BBBE00","#41BED1")
DimPlot(sce.all, reduction = "umap",split.by = 'orig.ident',
group.by = "celltype",label = T,cols = col,label.box=T)
很明显可以看出来，在C3D1时Tumor细胞减少了很多
Unkown这一群竟然增加了挺多，来看看这一群到底是什么吧
代码语言：javascript复制###看看2这一细胞亚群是什么
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A',
'CD19', 'CD79A', 'MS4A1' ,
'IGHG1', 'MZB1', 'SDC1',
'CD68', 'CD163', 'CD14',
'TPSAB1' , 'TPSB2', # mast cells,
'RCVRN','FPR1' , 'ITGAM' ,
'C1QA', 'C1QB', # mac
'S100A9', 'S100A8', 'MMP19',# monocyte
'FCGR3A',
'LAMP3', 'IDO1','IDO2',## DC3
'CD1E','CD1C', # DC2
'KLRB1','NCR1', # NK
'FGF7','MME', 'ACTA2', ## fibo
'DCN', 'LUM', 'GSN' , ## mouse PDAC fibo
'MKI67' , 'TOP2A',
'PECAM1', 'VWF', ## endo
'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' )
p <- DotPlot(sce.all, features = genes_to_check,
assay='RNA' ,group.by = 'RNA_snn_res.0.1' ) coord_flip()
p
看起来这一群是plasma
3.featureplot&heatmap
代码语言：javascript复制library(grid)
library(ggplot2)
library(gridExtra)
##基因在细胞中的表达
##EPCAM，PGR,TFF3
#PTPRC
#PECAM1
#ACTA2
tumor_cellmarker <- c('EPCAM','PGR','TFF3')
immune_cellmarker <- c("PTPRC")
endo_cellmarker <- c("PECAM1")
CAFs_cellmarker <- c("ACTA2")
FeaturePlot(sce.all,features = tumor_cellmarker,cols = c("lightgrey" ,"#DE1F1F"),ncol=3,raster=FALSE)
FeaturePlot(sce.all,features = immune_cellmarker,cols = c("lightgrey" ,"#DE1F1F"),raster=FALSE)
FeaturePlot(sce.all,features = endo_cellmarker,cols = c("lightgrey" ,"#DE1F1F"),raster=FALSE)
FeaturePlot(sce.all,features = CAFs_cellmarker,cols = c("lightgrey" ,"#DE1F1F"),raster=FALSE)
plot1 <- FeaturePlot(sce.all,features = CAFs_cellmarker,cols = c("lightgrey" ,"#DE1F1F"),raster=FALSE)
# 添加标题
title <- textGrob("CAFs", gp = gpar(fontsize = 20, fontface = "bold"))
# 将标题放在图片上方
plot_with_title <- grid.arrange(top = title, plot1)
#2.热图
Idents(sce.all)
sce.markers <- FindAllMarkers(object = sce.all, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
DT::datatable(sce.markers)
#pro='markers'
#write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
library(dplyr)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_log2FC)
#为了防止数据量太大不好出图，这里在每个亚群提取出来100个
DoHeatmap(subset(sce.all,downsample=100),top10$gene,size=3) scale_fill_gradientn(colors=c("#94C4E1","white","red"))
看看每个群里差异表达的基因
4.治疗前后细胞亚群比例
代码语言：javascript复制rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(celldex)
library(singleseqgset)
library(devtools)
getwd()
dir.create("4-prop")
setwd('4-prop/')
sce.all=readRDS( "3-cell/sce_celltype.rds")
sce.all
##堆积柱状图
library(tidyr)
library(reshape2)
tb=table(sce.all$group, sce.all$celltype)
head(tb)
library (gplots)
balloonplot(tb)
看看C1D1和C3D1组中细胞亚群的比例
bar_data <- as.data.frame(tb)
bar_per <- bar_data %>%
group_by(Var1) %>%
mutate(sum(Freq)) %>%
mutate(percent = Freq / `sum(Freq)`)
head(bar_per)
#write.csv(bar_per,file = "celltype_by_group_percent.csv")
col =c("#3176B7","#F78000","#3FA116","#CE2820","#9265C1",
"#885649","#DD76C5","#7F7F7F","#BBBE00","#41BED1")
colnames(bar_per)
p1 = ggplot(bar_per, aes(x = Freq, y = Var1))
geom_bar(aes(fill = Var2) , stat = "identity") coord_flip()
theme(axis.ticks = element_line(linetype = "blank"),
legend.position = "top",
panel.grid.minor = element_line(colour = NA,linetype = "blank"),
panel.background = element_rect(fill = NA),
plot.background = element_rect(colour = NA))
labs(y = " ", fill = NULL) labs(x = 'Total cell number')
scale_fill_manual(values=col)
theme_classic()
theme(plot.title = element_text(size=12,hjust=0.5))
p1
p2 = ggplot(bar_per, aes(x = percent, y = Var1))
geom_bar(aes(fill = Var2) , stat = "identity") coord_flip()
theme(axis.ticks = element_line(linetype = "blank"),
legend.position = "top",
panel.grid.minor = element_line(colour = NA,linetype = "blank"),
panel.background = element_rect(fill = NA),
plot.background = element_rect(colour = NA))
labs(y = " ", fill = NULL) labs(x = 'Relative proportion(%)')
scale_fill_manual(values=col)
theme_classic()
theme(plot.title = element_text(size=12,hjust=0.5))
p1 p2
治疗后肿瘤细胞明显减少，plasma增加。
5.EMT_score评分
先找到基因集
代码语言：javascript复制rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(celldex)
library(singleseqgset)
library(devtools)
getwd()
sce.all=readRDS( "3-cell/sce_celltype.rds")
sce.all
# BiocManager::install("msigdbr")
library(msigdbr)
#https://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION
?msigdbr
#EMT评分基因集在hallmark里
m = msigdbr(species = "Homo sapiens", category = "H")
m_EMT = m[m$gs_name == 'HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION',]
#评分计算
#BiocManager::install("UCell")
library(UCell)
EMT <- list(m_EMT$gene_symbol)#将基因整成list
names(EMT)[1] <- 'EMT'
DefaultAssay(sce.all) <- 'RNA'
EMT_score <- AddModuleScore_UCell(sce.all,
features=EMT,
name="Hallmark_EMT_score")
EMT_score$EMTHallmark_EMT_score
library(viridis)
FeaturePlot(EMT_score,features = "EMTHallmark_EMT_score",
order = T,cols = viridis(256), split.by = 'orig.ident')
#提取数据
library(ggpubr)
df<- FetchData(EMT_score,vars = c("orig.ident","EMTHallmark_EMT_score"))
#设置比较组
my_comparisons <- list(c("C1D1", "C3D1"))
colnames(df)
colnames(df)[2] = "Hallmark_EMT_score"
#小提琴图
ggplot(df,aes(x=orig.ident,y=Hallmark_EMT_score,fill=orig.ident))
geom_violin(color='black',size=1) #小提琴
theme_classic()
theme(text = element_text(size=10, colour = "black"))
theme(plot.title = element_text(hjust = 0.5, size = 15),
axis.text.x = element_text(colour = "black", size = 12),
axis.text.y = element_text(colour = "black", size = 10),
axis.title.y = element_text(color = 'black', size = 12),
axis.line = element_line(size = 1))
theme(legend.position="none")
xlab(NULL)
geom_boxplot(width=0.1, fill="white", outlier.shape = NA) #箱线图
stat_compare_means(method="t.test",hide.ns = F,
comparisons =c(my_comparisons),
label="p.signif",
bracket.size=0.8,
size=6)#添加显著性比较
6.肿瘤细胞重聚类分群
细胞分亚群重新降维聚类
代码语言：javascript复制## create subsets ----
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
getwd()
sce.all=readRDS( "3-cell/sce_celltype.rds")
table(Idents(sce.all))
immune <- subset(sce.all, idents = "Immune")
epithelial <- subset(sce.all, idents = "Tumor")
endothelial <- subset(sce.all, idents = "Endothelial")
CAFs <- subset(sce.all, idents = "CAFs")
saveRDS(immune, "immune.rds")
saveRDS(epithelial, "epithelial.rds")
saveRDS(endothelial, "endothelial.rds")
saveRDS(CAFs, "CAFs.rds")
####肿瘤细胞重聚类####
sce=readRDS( "epithelial.rds")
DefaultAssay(sce) <- "RNA"
table(sce$orig.ident)
table(sce$celltype)
table(sce$group,sce$celltype)
sce=CreateSeuratObject(counts = sce@assays$RNA@counts,
meta.data = sce@meta.data)
sce <- NormalizeData(sce, normalization.method = "LogNormalize",
scale.factor = 1e4)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce,
selection.method = "vst", nfeatures = 2000)
sce <- ScaleData(sce)
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
ElbowPlot(sce)
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.1)
table(sce@meta.data$RNA_snn_res.0.1)
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
sce <- RunUMAP(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "umap",label=T)
head(sce@meta.data)
DimPlot(sce,reduction = "umap",label=T,group.by = 'orig.ident')
mycolors <-c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3',
'#E95C59', '#E59CC4', '#AB3282', '#BD956A',
'#9FA3A8', '#E0D4CA', '#C5DEBA', '#F7F398',
'#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
'#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
'#968175')
library(harmony)
library(patchwork)
library(RColorBrewer)
library(scales)
col4<-colorRampPalette(brewer.pal(8,'Set2'))(12)
DimPlot(sce, reduction = 'umap',pt.size = 2,cols = mycolors,label = T,group.by = "RNA_snn_res.0.1",label.size = 6)
#####治疗前后之间的差异
DimPlot(sce, reduction = 'umap',pt.size = 1,cols = col4,label = T,split.by = 'orig.ident',label.size = 6,raster=FALSE)
#ggsave(filename = 'umap_sce_recluster_treatment.pdf',width = 10,height = 6.5)
跟文中趋势一致，基本上治疗后每个细胞群都降低了一些。
###看看比例
library(tidyr)
library(reshape2)
tb=table(sce$group, sce$RNA_snn_res.0.1)
head(tb)
library (gplots)
balloonplot(tb)
bar_data <- as.data.frame(tb)
bar_per <- bar_data %>%
group_by(Var1) %>%
mutate(sum(Freq)) %>%
mutate(percent = Freq / `sum(Freq)`)
head(bar_per)
#write.csv(bar_per,file = "celltype_by_group_percent.csv")
col =c("#3176B7","#F78000","#3FA116","#CE2820","#9265C1",
"#885649","#DD76C5","#7F7F7F","#BBBE00","#41BED1")
colnames(bar_per)
p1 = ggplot(bar_per, aes(x = Freq, y = Var1))
geom_bar(aes(fill = Var2) , stat = "identity") coord_flip()
theme(axis.ticks = element_line(linetype = "blank"),
legend.position = "top",
panel.grid.minor = element_line(colour = NA,linetype = "blank"),
panel.background = element_rect(fill = NA),
plot.background = element_rect(colour = NA))
labs(y = " ", fill = NULL) labs(x = 'Total cell number')
scale_fill_manual(values=col)
theme_classic()
theme(plot.title = element_text(size=12,hjust=0.5))
p1
p2 = ggplot(bar_per, aes(x = percent, y = Var1))
geom_bar(aes(fill = Var2) , stat = "identity") coord_flip()
theme(axis.ticks = element_line(linetype = "blank"),
legend.position = "top",
panel.grid.minor = element_line(colour = NA,linetype = "blank"),
panel.background = element_rect(fill = NA),
plot.background = element_rect(colour = NA))
labs(y = " ", fill = NULL) labs(x = 'Relative proportion(%)')
scale_fill_manual(values=col)
theme_classic()
theme(plot.title = element_text(size=12,hjust=0.5))
p1 p2
可以看到0这一个亚群，降低的特别明显。
##EMT评分
library(msigdbr)
#EMT评分基因集在hallmark里
m = msigdbr(species = "Homo sapiens", category = "H")
m_EMT = m[m$gs_name == 'HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION',]
#评分计算
library(UCell)
EMT <- list(m_EMT$gene_symbol)#将基因整成list
names(EMT)[1] <- 'EMT'
DefaultAssay(sce.all) <- 'RNA'
EMT_score <- AddModuleScore_UCell(sce,
features=EMT,
name="Hallmark_EMT_score")
EMT_score$EMTHallmark_EMT_score
library(viridis)
FeaturePlot(EMT_score,features = "EMTHallmark_EMT_score",
order = T,cols = viridis(256))
##看看另一种评分方式
#AddModuleScore不是Ucell哦
Inscore <- AddModuleScore(sce,
features = EMT,
ctrl = 100,
name = "Hallmark_EMT_score")
colnames(Inscore@meta.data)
VlnPlot(Inscore,features = 'Hallmark_EMT_score1',
pt.size = 0, adjust = 2,group.by = "orig.ident")
###密度图
library(UCell)
library(irGSEA)
library(Nebulosa)
library(decoupleR)
library(AUCell)
library(irGSEA)
install.packages('Nebulosa')
BiocManager::install('Nebulosa')
BiocManager::install('decoupleR')
BiocManager::install('AUCell')
if (!requireNamespace("irGSEA", quietly = TRUE)) {
devtools::install_github("chuiqin/irGSEA")
}
?
sce.final <- irGSEA.score(object = sce, assay = "RNA",
slot = "data", msigdb = T, species = "Homo sapiens",
category = "H", geneid = "symbol",
method = c("AUCell", "UCell", "singscore", "ssgsea"), kcdf = 'Gaussian')
irGSEA.density.scatterplot(object = sce.final,
method = "singscore",
show.geneset = "HALLMARK-EPITHELIAL-MESENCHYMAL-TRANSITION",
reduction = "umap")
这里是使用 singscore的评分方式。 所以密度图真的比FeaturePlot更为直观一点，颜色更为分明。 看出主要是0这一群EMT特征更强，治疗后0这一群显著减少，寿命治疗后抑制了EMT。
7.免疫细胞重聚类分群
下面就是免疫细胞的细分亚群
代码语言：javascript复制rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
getwd()
sce=readRDS( "immune.rds")
DefaultAssay(sce) <- "RNA"
table(sce$orig.ident)
table(sce$celltype)
table(sce$group,sce$celltype)
sce=CreateSeuratObject(counts = sce@assays$RNA@counts,
meta.data = sce@meta.data)
sce <- NormalizeData(sce, normalization.method = "LogNormalize",
scale.factor = 1e4)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce,
selection.method = "vst", nfeatures = 2000)
sce <- ScaleData(sce)
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)
ElbowPlot(sce)
sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.5)
table(sce@meta.data$RNA_snn_res.0.5)
set.seed(123)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)
sce <- RunUMAP(object = sce, dims = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "umap",label=T)
head(sce@meta.data)
DimPlot(sce,reduction = "umap",label=T,group.by = 'orig.ident')
mycolors <-c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3',
'#E95C59', '#E59CC4', '#AB3282', '#BD956A',
'#9FA3A8', '#E0D4CA', '#C5DEBA', '#F7F398',
'#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
'#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
'#968175')
library(harmony)
library(patchwork)
library(RColorBrewer)
library(scales)
col4<-colorRampPalette(brewer.pal(8,'Set2'))(13)
DimPlot(sce, reduction = 'umap',pt.size = 2,cols = mycolors,label = T,group.by = "RNA_snn_res.0.5",label.size = 6)
#####治疗前后之间的差异
DimPlot(sce, reduction = 'umap',pt.size = 1,cols = col4,label = T,split.by = 'orig.ident',label.size = 6,raster=FALSE)
#ggsave(filename = 'umap_sce_recluster_treatment.pdf',width = 10,height = 6.5)
使用singleR注释
代码语言：javascript复制library(SingleR)
library(celldex)
library(Seurat)
library(dplyr)
library(stringr)
library(pheatmap)
library(ReactomeGSA)
library(monocle)
library(limma)
library(ggplot2)
library(msigdbr)
library(singleseqgset)
library(devtools)
#install_github("arc85/singleseqgset")
# singleR注释
#hpca.se <- HumanPrimaryCellAtlasData()
#save(hpca.se,file = 'hpca.RData')
load('hpca.RData')
unique(hpca.se$label.main)
unique(hpca.se$label.fine)
unique(bpe.se$label.main)
unique(bpe.se$label.fine)
#bpe.se <- BlueprintEncodeData()
#save(bpe.se,file = 'bpe.RData')
load('bpe.RData')
str(sce)
anno <- SingleR(sce@assays$RNA$data,
ref = list(BP=bpe.se,HPCA=hpca.se),
labels = list(bpe.se$label.fine,hpca.se$label.main),
clusters = sce@meta.data$seurat_clusters
)
plotScoreHeatmap(anno,clusters = anno@rownames,show_colnames = T)
table(anno$labels)
celltype = data.frame(ClusterID=rownames(anno),
celltype=anno$labels,
stringsAsFactors = F)
sce@meta.data$singleR = celltype[match(sce@meta.data$seurat_clusters,celltype$ClusterID),'celltype']
DimPlot(sce, reduction = "umap",
group.by = "singleR",label = T,cols = col,label.box=T)
8.细胞间相互作用
有一点不理解的是，本来前面一直运行的好好的，突然上一步变成了Assay5，于是卡住。难道是在不知不觉中由于更新或安装某个包，Seurat就变成了 v5？ 首先装一下Y叔的Seurat4_4.4.0
代码语言：javascript复制rm(list=ls())
install.packages("Seurat4_4.4.0.tar.gz", repos = NULL, type = "source")
library(Seurat4)
options(stringsAsFactors = F)
library(SeuratObject)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
getwd()
sce_immue=readRDS( "immune_sce_celltype.rds")
sce_tumor=readRDS( "epithelial.rds")
library(CellChat)
library(tidyverse)
library(ggalluvial)
table(sce_tumor@meta.data$celltype)
table(sce_immue$orig.ident)
sce_tumor$singleR = 'tumor'
table(sce_tumor$singleR)
merge_sce <- merge(sce_immue, y = sce_tumor, add.cell.ids = c("immue", "tumor"), project = "CD",merge.data = TRUE)
table(merge_sce$singleR)
table(Idents(merge_sce))
Idents(merge_sce) = merge_sce$singleR
#将这些细胞类型提取出来
table(merge_sce$singleR)
con1 <- merge_sce$singleR %in% c("CD8 Tem","NK_cell","Tregs" ,"tumor")
con2 <- merge_sce$singleR %in% c("DC","Macrophages","Monocytes","tumor")
con3 <- merge_sce$orig.ident %in% c("C1D1")
con4 <- merge_sce$orig.ident %in% c("C3D1")
#C1D1-T-tumor
#C1D1-APC-tumor
#C3D1-T-tumor
#C3D1-APC-tumor
sce.T_C1D1 <- merge_sce[,con1&con3]
sce.T_C3D1 <- merge_sce[,con1&con4]
sce.APC_C1D1 <- merge_sce[,con2&con3]
sce.APC_C3D1 <- merge_sce[,con2&con4]
##超过三个了，循环一下
conditions = c("sce.T_C1D1" , "sce.T_C3D1", "sce.APC_C1D1", "sce.APC_C3D1")
cellchat_object = list()
for (i in 1:length(conditions)) {
# 提取input数据并添加生成meta信息
input_data <- GetAssayData(get(conditions[i]), assay = "RNA", slot = "data")
labels <- Idents(get(conditions[i]))
meta_data <- data.frame(group = labels, row.names = names(labels))
# 建立CellChat对象
cellchat_object[[i]] <- createCellChat(object = input_data, meta = meta_data, group.by = "group")
}
#设置配体-受体相互作用数据库
CellChatDB <- CellChatDB.human
CellChatDB.use <- CellChatDB
#这部分的处理包括：
#预处理表达数据以进行细胞间通讯分析；
#计算通信概率并推断通信网络；
#推断信号通路水平的细胞间通讯；
#计算聚合的cell-cell通信网络。
cc.fun <- function(cellchat){
levels(cellchat@idents)
#cellchat对象中增加DB插槽
cellchat@DB <- CellChatDB.use
#预处理表达数据以进行细胞间通讯分析
# subset the expression data of signaling genes for saving computation cost
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
# project gene expression data onto PPI network (optional)
cellchat <- projectData(cellchat, PPI.mouse) #PPI.human PPI.mouse
#计算通信概率并推断通信网络
#Compute the communication probability and iC3D1er cellular communication network
cellchat <- computeCommunProb(cellchat)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
#这里还是不要过滤了，不然有些细胞亚群细胞数量过少被过滤掉之后会造成数据不一致而无法进行后续分析
#cellchat <- filterCommunication(cellchat, min.cells = 10)
#推断信号通路水平的细胞间通讯
cellchat <- computeCommunProbPathway(cellchat)
#计算聚合的cell-cell通信网络
cellchat <- aggregateNet(cellchat)
# Compute the network centrality scores
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
# the slot 'netP' means the iC3D1erred intercellular communication network of signaling pathways
}
cellc = list()
#样本处理
for (i in 1:length(conditions)) {
cellchat <- cellchat_object[[i]]
cellc[[i]] <- cc.fun(cellchat)
}
conditions
object.list <- list(C1D1_T = cellc[[1]], C3D1_T = cellc[[2]])
cellchat1 <- mergeCellChat(object.list, add.names = names(object.list))
#比较相互作用的连接数和连接强度
gg1 <- compareInteractions(cellchat1, show.legend = F, group = c(1,2), measure = "count")
gg2 <- compareInteractions(cellchat1, show.legend = F, group = c(1,2), measure = "weight")
gg1 gg2
很明显治疗后tumor和T细胞间的连接和强度都比治疗前高。
# 比较不同细胞亚群间的连接数和连接强度的差异-网络图
par(mfrow = c(1,2), xpd=TRUE)
cellchat1@meta$datasets
g3 = netVisual_diffInteraction(cellchat1, weight.scale = T, measure = "count", color.edge = c("#b2182b", "#2166ac"))
g4 = netVisual_diffInteraction(cellchat1,weight.scale = T, measure = "weight", color.edge = c("#b2182b", "#2166ac"))
这里的对比顺序，呈现红色线条和蓝色线条的顺序，需要看看Levels
weight.max <- getMaxWeight(object.list, attribute = c("idents","count"))
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
netVisual_circle(object.list[[i]]@net$count, weight.scale = T, label.edge= F, edge.weight.max = weight.max[2], edge.width.max = 12, title.name = paste0("Number of interactions - ", names(object.list)[i]))
}
可以看出C3D1相比于C1D1，tumor和免疫细胞之间的连接强度大了很多
文中还显示在MHC-I和MHC-II在C3D1中增加
代码语言：javascript复制gg7 <- rankNet(cellchat1, mode = "comparison", stacked = T, do.stat = TRUE)
gg8 <- rankNet(cellchat1, mode = "comparison", stacked = F, do.stat = TRUE)
gg7 gg8
接下是tumor和APC之间的相互作用
代码语言：javascript复制####再看看另外两组
conditions
object.list <- list(C1D1_APC = cellc[[3]], C3D1_APC = cellc[[4]])
cellchat <- mergeCellChat(object.list, add.names = names(object.list))
#比较相互作用的连接数和连接强度
gg11 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "count")
gg22 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
gg11 gg22
# 比较不同细胞亚群间的连接数和连接强度的差异-网络图
par(mfrow = c(1,2), xpd=TRUE)
cellchat@meta$datasets
g33 = netVisual_diffInteraction(cellchat, weight.scale = T, measure = "count", color.edge = c("#b2182b", "#2166ac"))
g44 = netVisual_diffInteraction(cellchat,weight.scale = T, measure = "weight", color.edge = c("#b2182b", "#2166ac"))
治疗后tumor和抗原递呈细胞间的连接和强度都比治疗前高。
可以看出一点治疗后，tumor和DC之间的连接变多，但是DC和单核细胞之间更明显，所以并没有看到文中说治疗后tumor-monocyte转变为tumor-DC的趋势。
