代码语言：javascript复制我们在前些天的教程：直接为CellPhoneDB创建一个独立的conda环境，以及：把Seurat对象里面表达量矩阵和细胞表型信息输出给CellPhoneDB做细胞通讯，给大家演示了如何对pbmc3k单细胞数据集做细胞通讯，得到了如下所示的结果：
669 deconvoluted.txt
300 means.txt
300 pvalues.txt
300 significant_means.txt
就是4个文本文件，其中 means.txt 和 pvalues.txt以及， significant_means.txt都是300行，而且都是90多列，说明这3个文件是相互关联的。下面我们来一一解读：
其中 means.txt 如下所示：
代码语言：javascript复制 head -1 means.txt|tr 't' 'n'
id_cp_interaction
interacting_pair
partner_a
partner_b
gene_a
gene_b
secreted
receptor_a
receptor_b
annotation_strategy
is_integrin
B|B
B|CD14_Mono
而pvalues.txt 如下所示：
代码语言：javascript复制 head -1 pvalues.txt |tr 't' 'n'
id_cp_interaction
interacting_pair
partner_a
partner_b
gene_a
gene_b
secreted
receptor_a
receptor_b
annotation_strategy
is_integrin
B|B
B|CD14_Mono
最后的 significant_means.txt 文件 ：
代码语言：javascript复制head -1 significant_means.txt |tr 't' 'n'
id_cp_interaction
interacting_pair
partner_a
partner_b
gene_a
gene_b
secreted
receptor_a
receptor_b
annotation_strategy
is_integrin
rank
B|B
B|CD14_Mono
我们前面的接近3千个细胞，归类如下；
代码语言：javascript复制cut -f 2 test_meta.txt |sort |uniq -c
344 B
480 CD14_Mono
271 CD8_T
32 DC
162 FCGR3A_Mono
483 Memory_CD4_T
155 NK
697 Naive_CD4_T
14 Platelet
也就是说9个细胞亚群，每个都是自己跟包括自己在内的9个单细胞亚群进行通讯分析。也就是说81种结果，占81 列。而means.txt 和 pvalues.txt以及， significant_means.txt 都是300行，而且都是90多列。
但是我稍微看了看前面的十多行内容，发现基因仍然是ensembl数据库的ID，如下所示：
代码语言：javascript复制id_cp_interaction interacting_pair partner_a partner_b gene_a gene_b secreted receptor_a receptor_b annotation_strategy
CPI-SS0B84DAE3D PVR_CD96 simple:P15151 simple:P40200 ENSG00000073008 ENSG00000153283 True True True curated
CPI-SS0A8627ED6 PVR_CD226 simple:P15151 simple:Q15762 ENSG00000073008 ENSG00000150637 True True True curated
CPI-SS00561BBD7 PVR_TIGIT simple:P15151 simple:Q495A1 ENSG00000073008 ENSG00000181847 True True False curated
CPI-SS0889E281D CADM1_CADM1 simple:Q9BY67 simple:Q9BY67 ENSG00000182985 ENSG00000182985 False False False curated
CPI-SS0C9696AAA CRTAM_CADM1 simple:O95727 simple:Q9BY67 ENSG00000109943 ENSG00000182985 False False False curated
CPI-SC02B111AB0 COL6A1_a10b1 complex simple:P12109 complex:a10b1 complex ENSG00000142156 True False False curated
CPI-SC0CE851065 COL6A2_a10b1 complex simple:P12110 complex:a10b1 complex ENSG00000142173 True False False curated
CPI-SC0017F527E COL18A1_a10b1 complex simple:P39060 complex:a10b1 complex ENSG00000182871 True False False curated
CPI-SC00C0800EB COL4A4_a10b1 complex simple:P53420 complex:a10b1 complex ENSG00000081052 True False False curated
所以我调整了前面的代码，这次直接输出表达量矩阵就是基因名字的，不需要ensembl数据库的ID ：
代码语言：javascript复制library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("sce")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
DimPlot(sce,label = T,repel = T)
test_counts <- as.data.frame( sce@assays$RNA@data)
test_counts[1:4,1:4]
table(sce$seurat_annotations )
test_meta <- data.frame(Cell = rownames(sce@meta.data),
cell_type = sce$seurat_annotations )
head(test_meta)
test_meta$cell_type=gsub(' ','_',test_meta$cell_type)
test_meta$cell_type=gsub('\ ','',test_meta$cell_type)
table(test_meta$cell_type)
length(unique(test_meta$Cell))
identical(colnames(test_counts),test_meta$Cell)
test_counts=cbind(rownames(test_counts),test_counts)
colnames(test_counts)[1]='Gene'
test_counts[1:4,1:4]
write.table(test_counts, "test_counts.txt", row.names=F, sep='t',quote = F)
write.table(test_meta, "test_meta.txt", row.names=F, sep='t',quote = F)
运行CellPhoneDB的代码如下所示：
代码语言：javascript复制conda activate cellphonedb
cellphonedb method statistical_analysis test_meta.txt test_counts.txt --counts-data hgnc_symbol --threads 4
运行仍然是很快，仍然是不到十分钟解决战斗 ：
代码语言：javascript复制[ ][APP][12/02/22-09:21:59][WARNING] Latest local available version is `v2.0.0`, using it
[ ][APP][12/02/22-09:21:59][WARNING] User selected downloaded database `v2.0.0` is available, using it
[ ][CORE][12/02/22-09:21:59][INFO] Initializing SqlAlchemy CellPhoneDB Core
[ ][CORE][12/02/22-09:21:59][INFO] Using custom database at /Users/jmzeng/.cpdb/releases/v2.0.0/cellphone.db
[ ][APP][12/02/22-09:21:59][INFO] Launching Method cpdb_statistical_analysis_local_method_launcher
[ ][APP][12/02/22-09:21:59][INFO] Launching Method _set_paths
[ ][APP][12/02/22-09:21:59][INFO] Launching Method _load_meta_counts
[ ][APP][12/02/22-09:22:04][INFO] Launching Method _check_counts_data
[ ][CORE][12/02/22-09:22:04][INFO] Launching Method cpdb_statistical_analysis_launcher
[ ][CORE][12/02/22-09:22:04][INFO] Launching Method _counts_validations
[ ][CORE][12/02/22-09:22:04][INFO] Launching Method get_interactions_genes_complex
[ ][CORE][12/02/22-09:22:04][INFO] [Cluster Statistical Analysis] Threshold:0.1 Iterations:1000 Debug-seed:-1 Threads:4 Precision:3
[ ][CORE][12/02/22-09:22:05][INFO] Running Real Analysis
[ ][CORE][12/02/22-09:22:05][INFO] Running Statistical Analysis
[ ][CORE][12/02/22-09:30:20][INFO] Building Pvalues result
[ ][CORE][12/02/22-09:30:20][INFO] Building results
有意思的是这个时候的结果稍微有一点点不同：
代码语言：javascript复制 wc -l out/*
673 out/deconvoluted.txt
303 out/means.txt
303 out/pvalues.txt
303 out/significant_means.txt
应该是前面的ID转换过程的一些基因损失造成的。
这个时候，当然是可以使用CellPhoneDB自带的一些统计绘图函数：
代码语言：javascript复制conda activate cellphonedb
# 必须要保证当前路径下面有前面的步骤输出的out文件夹哦
cellphonedb plot dot_plot
cellphonedb plot heatmap_plot cellphonedb_meta.txt
# 做完后为了跟前面的区分，我把 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
可以看到是多了两个热图，一个气泡图，其中热图是展现 其计算好的 count_network.txt 里面的内容，如下所示：
计算好的 count_network.txt 里面的内容
而气泡图就是展现前面的 means.txt 和 pvalues.txt 以及， significant_means.txt 。
这个最新计算好的 count_network.txt 里面的内容，就是9个细胞亚群，每个都是自己跟包括自己在内的9个单细胞亚群进行通讯分析后的合格的通讯（受体和配体对）数量，也就是81行，再加上一个表头，就是82行啦 ：
代码语言：javascript复制SOURCE TARGET count
Memory_CD4_T Memory_CD4_T 3
Memory_CD4_T B 11
Memory_CD4_T CD14_Mono 13
Memory_CD4_T NK 17
Memory_CD4_T CD8_T 7
Memory_CD4_T Naive_CD4_T 4
Memory_CD4_T FCGR3A_Mono 23
Memory_CD4_T DC 20
Memory_CD4_T Platelet 3
这个表格其实很容易自己使用代码，读取 pvalues.txt 内容后做统计。
代码语言：javascript复制mypvals <- read.table("./out-by-symbol/pvalues.txt",header = T,sep = "t",stringsAsFactors = F)
n = which(colnames(mypvals) == "Memory_CD4_T.Memory_CD4_T")
choose_p = mypvals[,c(1,5,6,8,9,n)]
table(choose_p$Memory_CD4_T.Memory_CD4_T < 0.05)
# 可以看到，就是 Memory_CD4_T Memory_CD4_T 3 的内容
# FALSE TRUE
# 299 3
如果需要做全部的，就是9个细胞亚群，每个都是自己跟包括自己在内的9个单细胞亚群进行通讯分析后的合格的通讯（受体和配体对）数量，也就是81次，这样的统计即可。
现在我们理解了这些文本文件，接下来就是如何读取它们，并且进行合适的可视化啦。
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
代码语言：javascript复制We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。
