代码语言：javascript复制前些天的教程：直接为CellPhoneDB创建一个独立的conda环境，以及：把Seurat对象里面表达量矩阵和细胞表型信息输出给CellPhoneDB做细胞通讯，给大家演示了如何对pbmc3k单细胞数据集做细胞通讯，并且在：CellPhoneDB的单细胞通讯结果的理解 给大家演示了细胞通讯结果的多个txt文件的含义，前面的代码是：
conda activate cellphonedb
cellphonedb method statistical_analysis test_meta.txt test_counts.txt
--counts-data hgnc_symbol --threads 4
# 必须要保证当前路径下面有前面的步骤输出的out文件夹哦
cellphonedb plot dot_plot
cellphonedb plot heatmap_plot cellphonedb_meta.txt
# 做完后为了跟前面的区分，我把 out文件夹，修改名字为 out-by-symbol 文件夹啦
把 out文件夹，修改名字为 out-by-symbol文件夹，其结构如下所示：
代码语言：javascript复制ls -lh out-by-symbol|cut -d" " -f 7-
1.5K 2 12 09:43 count_network.txt
64K 2 12 09:30 deconvoluted.txt
11K 2 12 09:43 heatmap_count.pdf
11K 2 12 09:43 heatmap_log_count.pdf
113B 2 12 09:43 interaction_count.txt
149K 2 12 09:30 means.txt
803K 2 12 09:43 plot.pdf
128K 2 12 09:30 pvalues.txt
61K 2 12 09:30 significant_means.txt
如果你对这些文件的理解还不够，继续看 ：CellPhoneDB的单细胞通讯结果的理解
接下来我们对这些文件进行高度定制化的可视化！
我们肉眼看的：
代码语言：javascript复制 all_sum
Memory_CD4_T 101
B 105
CD14_Mono 157
NK 171
CD8_T 100
Naive_CD4_T 68
FCGR3A_Mono 241
DC 180
Platelet 60
同样是单核细胞， CD14_Mono 和 FCGR3A_Mono 的总通讯数量不一样，那就看看它们跟两种CD4的T细胞（Memory_CD4_T和Naive_CD4_T）的通讯情况吧。
主要是因为细胞亚群太多了，如果全部可视化出来，也很难给大家讲清楚。
代码语言：javascript复制
rm(list = ls())
getwd()
setwd('pbmc3k-CellPhoneDB/')
mypvals <- read.table("./out-by-symbol/pvalues.txt",header = T,sep = "t",stringsAsFactors = F)
mymeans <- read.table("./out-by-symbol/means.txt",header = T,sep = "t",stringsAsFactors = F)
# 我们重点看 CD14_Mono 和 FCGR3A_Mono ，以及 Memory_CD4_T和Naive_CD4_T 的通讯情况。
kp = grepl(pattern = "Mono", colnames(mypvals)) & grepl(pattern = "CD4_T", colnames(mypvals))
table(kp)
pos = (1:ncol(mypvals))[kp]
choose_pvalues <- mypvals[,c(c(1,5,6,8,9),pos )]
choose_means <- mymeans[,c(c(1,5,6,8,9),pos)]
logi <- apply(choose_pvalues[,5:ncol(choose_pvalues)]<0.05, 1, sum)
# 只保留具有细胞特异性的一些相互作用对
choose_pvalues <- choose_pvalues[logi>=1,]
# 去掉空值
logi1 <- choose_pvalues$gene_a != ""
logi2 <- choose_pvalues$gene_b != ""
logi <- logi1 & logi2
choose_pvalues <- choose_pvalues[logi,]
# 同样的条件保留choose_means
choose_means <- choose_means[choose_means$id_cp_interaction %in% choose_pvalues$id_cp_interaction,]
得到了如下所示的两个变量，choose_means和 choose_pvalues ， 他们前面的5列是一样的，但是后面只有8列，并不是 4X4的16列，说明很多单细胞亚群的两两组合之间其实并没有统计学显著的细胞通讯受体和配体对哦。
代码语言：javascript复制> head(choose_means,1)
id_cp_interaction gene_a gene_b receptor_a receptor_b CD14_Mono.Memory_CD4_T CD14_Mono.Naive_CD4_T
40 CPI-SS0703338F5 LGALS9 CD44 False True 0.862 0.803
FCGR3A_Mono.Memory_CD4_T FCGR3A_Mono.Naive_CD4_T Memory_CD4_T.CD14_Mono Memory_CD4_T.FCGR3A_Mono
40 0.949 0.89 0.552 0.71
Naive_CD4_T.CD14_Mono Naive_CD4_T.FCGR3A_Mono
40 0.513 0.67
> head(choose_pvalues,1)
id_cp_interaction gene_a gene_b receptor_a receptor_b CD14_Mono.Memory_CD4_T CD14_Mono.Naive_CD4_T
40 CPI-SS0703338F5 LGALS9 CD44 False True 0 0
FCGR3A_Mono.Memory_CD4_T FCGR3A_Mono.Naive_CD4_T Memory_CD4_T.CD14_Mono Memory_CD4_T.FCGR3A_Mono
40 0 0 0.888 0
Naive_CD4_T.CD14_Mono Naive_CD4_T.FCGR3A_Mono
40 1 1
可以看到之前是 302行， 是90多列，也就是说9个细胞亚群，每个都是自己跟包括自己在内的9个单细胞亚群进行通讯分析结果都存在，81种两两组合都存在 ，但是现在只剩下8个单细胞亚群组合啦，而且仅仅是19个受体配体对啦 ：
代码语言：javascript复制> dim(mypvals)
[1] 302 92
> dim(mymeans)
[1] 302 92
> dim(choose_means)
[1] 19 13
> dim(choose_pvalues)
[1] 19 13
接下来对这19个受体配体对，在8个单细胞亚群组合的结果进行简单的可视化，虽然说简单，但是代码也有点长啊 ：
代码语言：javascript复制
# 将choose_pvalues和choose_means数据宽转长
library(tidyverse)
meansdf <- choose_means %>% reshape2::melt()
meansdf <- data.frame(interacting_pair = paste0(meansdf$gene_a,"_",meansdf$gene_b),
CC = meansdf$variable,
means = meansdf$value)
pvalsdf <- choose_pvalues %>% reshape2::melt()
pvalsdf <- data.frame(interacting_pair = paste0(pvalsdf$gene_a,"_",pvalsdf$gene_b),
CC = pvalsdf$variable,
pvals = pvalsdf$value)
# 合并p值和mean文件
pvalsdf$joinlab<- paste0(pvalsdf$interacting_pair,"_",pvalsdf$CC)
meansdf$joinlab<- paste0(meansdf$interacting_pair,"_",meansdf$CC)
pldf <- merge(pvalsdf,meansdf,by = "joinlab")
# dotplot可视化
summary((filter(pldf,means >0))$means)
head(pldf)
pcc = pldf%>% filter(means >0) %>%
ggplot(aes(CC.x,interacting_pair.x) )
geom_point(aes(color=means,size=-log10(pvals 0.0001)) )
scale_size_continuous(range = c(1,3))
scale_color_gradient2(high="red",mid = "yellow",low ="darkblue",midpoint = 1 )
theme_bw()
# scale_color_manual(values = rainbow(100))
theme(axis.text.x = element_text(angle = -45,hjust = -0.1,vjust = 0.8))
pcc
如下所示：
可以看到，从mono到cd4方向的受体配体对，跟从 cd4到mono的，基本上是完全不一样的哦。这19个受体配体对具体信息如下所示
代码语言：javascript复制> unique(pldf$interacting_pair.x)
[1] "ALOX5_ALOX5AP" "ANXA1_FPR1" "C5AR1_RPS19" "CD2_CD58" "CD28_CD86" "CD47_SIRPG"
[7] "CD52_SIGLEC10" "CD74_MIF" "CD99_PILRA" "HLA-F_LILRB1" "HLA-F_LILRB2" "LAMP1_VSTM1"
[13] "LGALS9_CD44" "LGALS9_CD47" "LGALS9_SORL1" "MIF_TNFRSF14" "PLXNB2_SEMA4D" "SCGB3A1_MARCO"
[19] "SELL_SELPLG"
而且，肉眼上看，是 CD74_MIF 和 C5AR1_RPS19 这两个配对比较明显， 让我们看看他们的表达量情况，代码如下所示：
代码语言：javascript复制
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("sce")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
sce = sce[,Idents(sce) %in% c('Naive CD4 T', 'Memory CD4 T' , 'CD14 Mono','FCGR3A Mono')]
DimPlot(sce,label = T,repel = T)
cg = unique(unlist(str_split(unique(pldf[,2]),'_')) )
library(ggplot2)
th= theme(axis.text.x = element_text(angle = -45,hjust = -0.1,vjust = 0.8))
pdopt <- DotPlot(sce, features = cg,
assay='RNA' ) coord_flip() th
library(patchwork)
pdopt DimPlot(sce,label = T,repel = T)
pdopt pcc
出结果如下所示：
是不是很简单，前面的CD74_MIF 和 C5AR1_RPS19 这两个配对情况，在mono到cd4的方向比较明显，所以它们就是CD74和C5AR1在mono高表达，而MIF和RPS19在cd4的t细胞里面高表达。
所谓的细胞通讯，就是差异分析从一个基因单位到两个基因
我们差异分析很容易找到不同单细胞亚群各自的高表达量基因，但是这些基因在每个单细胞亚群是独立的列表，所谓的单细胞通讯，就是找到那些不同单细胞亚群的特异性基因的两两组合，而且是具有受体配体的生物学背景知识的那些。
