单细胞转录组高级分析三：细胞通讯分析

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨!

细胞通讯分析原理

细胞通讯研究领域涵盖的内容很广，如上图所示包括通讯方式、功能、信号分子以及各种途径的机制。细胞之间通讯的介质有很多，例如钙离子、脂质、多肽、蛋白、外泌体以及电信号等。利用单细胞转录组数据分析的细胞通讯，仅限于蛋白质配体-受体复合物介导的细胞间通讯。其分析的基础是基因表达数据和配体-受体数据库信息，例如转录组数据表明A、B细胞分别表达了基因α和β，通过数据库查询α和β是配体-受体关系，则认为A-B通过α-β途径进行了通讯。

细胞通讯分析工具

目前利用单细胞转录组数据分析细胞通讯的工具有CellPhoneDB，celltalker和iTALK。CellPhoneDB相关的文章2018年发表在Nature上，升级后又发了一篇Nature Protocols，在三者之中最权威。iTALK的可视化效果最好，相关文献2019年发表在bioRxiv。celltalker由本专题演示数据来源文章的作者开发，相关文献2020年发表于Immunity。iTALK和celltalker都只有R语言版本，CellPhoneDB可以使用网页版和python版。

可视化结果对比 

CellPhoneDB

celltalker

iTALK

CellPhoneDB网页版

CellPhoneDB有网页版和python单机版，建议有linux基础的朋友使用python版，官方教程：https://github.com/Teichlab/cellphonedb. 由于本教程只针对非专业生信人员，所以只介绍网页版的使用。

CellPhoneDB的优势 

CellPhoneDB的独特之处不是它的算法，而是他们开发的细胞通讯配体-受体数据库。他们的数据库不仅包含公共资源注释的受体和配体，还有人工挑选的参与细胞通讯的特定蛋白质家族。他们的数据库还考虑了配体和受体的亚基结构，这是大多数据库和研究没有涉及的，这点对于准确研究多亚基复合物介导的细胞通讯很关键。CellPhoneDB数据库概况如下图所示：

CellPhoneDB官网 

登录官网：https://www.cellphonedb.org/explore-sc-rna-seq

网页打开后如下图所示，注意红色标记区域，这些是需要按情况更改的。统计迭代次数最好改成1000

分析数据准备 

细胞通讯肯定只能是样本内的，因此我们需要按样本提取数据子集。

library(Seurat)
library(tidyverse)
dir.create("CellTalk")
scRNA <- readRDS("scRNA.rds")
#提取HNSCC肿瘤样本HNC01TIL
sp1 <- scRNA[,str_detect(colnames(scRNA),'HNC01TIL')]
sp1_counts <- as.matrix(sp1@assays$RNA@data)
sp1_counts <- data.frame(Gene=rownames(sp1_counts), sp1_counts)
sp1_meta <- data.frame(Cell=rownames([e-mail]), cell_type=[e-mail]$celltype_Monaco)
write.table(sp1_counts, "CellTalk/sp1_counts.txt", row.names=F, sep='\t')
write.table(sp1_meta, "CellTalk/sp1_meta.txt", row.names=F, sep='\t')

正常情况下将sp1_counts.txt和sp1_meta.txt文件提交到服务器就可以分析了，不巧的是我遇到CellPhoneDB服务器忙碌而暂停服务了，因此我用python版本跑了一下，可视化结果如下：

部分细胞通讯关系的气泡图

细胞通讯关系热图

细胞通讯关系网络图

celltalker的安装使用

celltalker分析原理 

celltalker的分析同样依赖配体-受体配对信息，它采纳的数据来自Ramilowski et al (Nature Communications, 2015)的研究成果，并且允许用户自己添加配对信息。celltalker分析细胞通讯有它自己的特点：

• celltalker可能更倾向于分析样本之间具有差异的细胞通讯关系，我没有找到分析单个样本细胞通讯的教程，其构建celltalker对象的函数也要求输入分组信息。
  

• 为了保证分析的可靠性，它要求一个分组要有几个重复样本。

• celltalker认定细胞进行通讯的前提是配体和受体的表达值在通讯的细胞之间具有一致性。
  

 安装celltalker 

library(devtools)
install_github("arc85/celltalker")

实例演示 

本专题来源于GSE139324数据集的10个样本，我们提取其中的肿瘤浸润免疫细胞和正常扁桃体免疫细胞对比分析。

library(Seurat)
library(tidyverse)
library(celltalker)
scRNA <- readRDS("scRNA.rds")
##调整scRNA的meta.data，方便满足celltalker的数据要求
cell_meta = [e-mail][,1:9]
names(cell_meta)[which(names(cell_meta)=='orig.ident')] <- "sample.id"
names(cell_meta)[which(names(cell_meta)=='celltype_Monaco')] <- "cell.type"
cell_meta$sample.type <- cell_meta$sample.id
cell_meta[grep("^HNC[0-9]+PBMC", cell_meta$sample.type), "sample.type"] <- "HNC_PBMC"
cell_meta[grep("^HNC[0-9]+TIL", cell_meta$sample.type), "sample.type"] <- "HNC_TIL"
cell_meta[grep("^PBMC", cell_meta$sample.type), "sample.type"] <- "PBMC"
cell_meta[grep("^Tonsil", cell_meta$sample.type), "sample.type"] <- "Tonsil"
[e-mail] <- cell_meta

##提取数据子集
cellsub = rownames(subset(cell_meta, sample.type=="HNC_TIL"|sample.type=="Tonsil"))
scRNA <- scRNA[,cellsub]

##鉴定scRNA中存在的配体受体
ligs <- as.character(unique(ramilowski_pairs$ligand))
recs <- as.character(unique(ramilowski_pairs$receptor))
ligs.present <- rownames(scRNA)[rownames(scRNA) %in% ligs]
recs.present <- rownames(scRNA)[rownames(scRNA) %in% recs]
genes.to.use <- union(ligs.present,recs.present)

##鉴定样本分组之间差异表达的配体受体
Idents(scRNA) <- "sample.type" 
markers <- FindAllMarkers(scRNA,assay="RNA",features=genes.to.use,only.pos=TRUE)
ligs.recs.use <- unique(markers$gene)

##保留差异表达的配体-受体列表
interactions.forward1 <- ramilowski_pairs[as.character(ramilowski_pairs$ligand) %in% ligs.recs.use,]
interactions.forward2 <- ramilowski_pairs[as.character(ramilowski_pairs$receptor) %in% ligs.recs.use,]
interact.for <- rbind(interactions.forward1,interactions.forward2)

##准备celltalker的输入文件
expr.mat <- GetAssayData(scRNA, slot="counts")
defined.clusters <- [e-mail]$cell.type
defined.groups <- [e-mail]$sample.type
defined.replicates <- [e-mail]$sample.id

##重构数据，为每个样本分配相应的表达矩阵
reshaped.matrices <- reshape_matrices(count.matrix=expr.mat, clusters=defined.clusters,
    groups=defined.groups, replicates=defined.replicates, ligands.and.receptors=interact.for)

##分析表达一致性的配体和受体。
# cells.reqd=10：每个分组中至少有10个细胞表达了配体/受体
# freq.pos.reqd=0.5：至少有50％的样本表达的配体/受体
consistent.lig.recs <- create_lig_rec_tib(exp.tib=reshaped.matrices,
                                          clusters=defined.clusters,
                                          groups=defined.groups,
                                          replicates=defined.replicates,
                                          cells.reqd=10,
                                          freq.pos.reqd=0.5,
                                          ligands.and.receptors=interact.for)

##生成Celltalker对象，并绘制circos圈图
# freq.group.in.cluster: 只对细胞数>总细胞数5%的细胞类型分析
put.int<-putative_interactions(ligand.receptor.tibble=consistent.lig.recs,
                               clusters=defined.clusters,
                               groups=defined.groups,
                               freq.group.in.cluster=0.05,
                               ligands.and.receptors=interact.for)


##鉴定分组间特异出现的配体/受体并作图
unique.ints <- unique_interactions(put.int, group1="HNC_TIL", group2="Tonsil", interact.for)
#HNC_TIL组作图
HNC_TIL <- pull(unique.ints[1,2])[[1]]
circos.HNC_TIL <- pull(put.int[1,2])[[1]][HNC_TIL]
par(MAR=c(1,1,1,1))
circos_plot(interactions=circos.HNC_TIL, clusters=defined.clusters)
#Tonsil组作图
Tonsil <- pull(unique.ints[2,2])[[1]]
circos.Tonsil <- pull(put.int[2,2])[[1]][Tonsil]
circos_plot(interactions=circos.Tonsil, clusters=defined.clusters)

实际操作中，celltalker的数据提取和图形调整都比较困难，个人不推荐此包

iTALK的安装使用

 iTALK的特点 

• 集成了多种差异分析和可视化方法；

• 将配体-受体注释为4大类：细胞因子、生长因子、免疫检查点和其他；

• 配体-受体分析的中间数据和最终结果都可以轻松导出。

 iTALK的安装 

install_github("Coolgenome/iTALK", build_vignettes = TRUE)

官方教程：https://github.com/Coolgenome/iTALK/tree/master/example

iTALK用于本例数据 

library(Seurat)
library(tidyverse)
library(iTALK)
rm(list=ls())

##准备iTALK输入文件
# 输入文件格式为数据框，行为细胞名称，列为基因名和细胞注释信息，基因名命名的列值为表达数据。
# 细胞注释信息必选“cell_type”，值为细胞类型注释信息；可选“compare_group”，值为细胞的样本分组信息。
cell_type
compare_group
cell.meta <- subset([e-mail], select=c("orig.ident","celltype_Monaco"))
cell.meta[grep("^HNC[0-9]+PBMC", cell.meta$orig.ident), "orig.ident"] <- "HNC_PBMC"
cell.meta[grep("^HNC[0-9]+TIL", cell.meta$orig.ident), "orig.ident"] <- "HNC_TIL"
cell.meta[grep("^PBMC", cell.meta$orig.ident), "orig.ident"] <- "PBMC"
cell.meta[grep("^Tonsil", cell.meta$orig.ident), "orig.ident"] <- "Tonsil"
names(cell.meta) <- c("compare_group", "cell_type")
cell.expr <- data.frame(t(as.matrix(scRNA@assays$RNA@counts)), check.names=F)
data.italk <- merge(cell.expr, cell.meta, by=0)
rownames(data.italk) <- data.italk$Row.names
data.italk <- data.italk[,-1]

查看iTALK输入文件格式data.italk[1:5,c(1:2,23603:23605)]

##设置绘图颜色和其他变量
mycolor <- c("#92D0E6","#F5949B","#E11E24","#FBB96F","#007AB8","#A2D184","#00A04E","#BC5627","#0080BD","#EBC379","#A74D9D")
cell_type <- unique(data.italk$cell_type)
cell_col <- structure(mycolor[1:length(cell_type)], names=cell_type)
#通讯类型变量
comm_list<-c('growth factor','other','cytokine','checkpoint')
#圈图展示配体-受体对的数量
PN=20

iTALK的可视化结果主要是网络图和弦图，网络图用数字表示所有的配体-受体关系，弦图为了显示效果一般只选前20个，可用通过变量PN修改。

查看样本内配体-受体概况 

##==细胞通讯关系概览==##
dir.create('iTALK/Overview')
setwd("iTALK/Overview")
#实际上一个样本之内的细胞才能通讯，这里提取一组样本分析是为了克服单细胞数据稀疏性造成的误差
data1 <- subset(data.italk, subset=data.italk$compare_group=="Tonsil")
##寻找高表达的配体受体基因,top_genes=50代表提取平均表达量前50%的基因
highly_exprs_genes <- rawParse(data1, top_genes=50, stats="mean")
res<-NULL
for(comm_type in comm_list){
    #comm_type = 'cytokine'   #测试循环代码的临时变量
    res_cat <- FindLR(highly_exprs_genes, datatype='mean count', comm_type=comm_type)
    res_cat <- res_cat[order(res_cat$cell_from_mean_exprs*res_cat$cell_to_mean_exprs, decreasing=T),]
    write.csv(res_cat, paste0('LRpairs_Overview_',comm_type,'.xls'))
  #pdf(paste0('LRpairs_Overview_',comm_type,'.pdf'), width=6, height=7)
  png(paste0('LRpairs_Overview_',comm_type,'_net.png'), width=600, height=650)
  #绘制细胞通讯关系网络图
    NetView(res_cat, col=cell_col, vertex.label.cex=1.2, edge.label.cex=0.9, 
             vertex.size=30, arrow.width=3, edge.max.width=10, margin = 0.2)
  dev.off()
  #绘制topPN的配体-受体圈图
  png(paste0('LRpairs_Overview_',comm_type,'_circ.png'), width=600, height=650)
    LRPlot(res_cat[1:PN,],datatype='mean count',cell_col=cell_col,link.arr.lwd=res_cat$cell_from_mean_exprs[1:PN],link.arr.width=res_cat$cell_to_mean_exprs[1:PN], text.vjust = "0.35cm")
  dev.off()
    res<-rbind(res,res_cat)
}
res <- res[order(res$cell_from_mean_exprs*res$cell_to_mean_exprs, decreasing=T),]
write.csv(res, 'LRpairs_Overview.xls')
res.top <- res[1:PN,]
png('LRpairs_Overview_net.png', width=600, height=650)
NetView(res, col=cell_col, vertex.label.cex=1.2, edge.label.cex=0.9, 
         vertex.size=30, arrow.width=3, edge.max.width=10, margin = 0.2)
dev.off()
png('LRpairs_Overview_circ.png', width=600, height=650)
LRPlot(res.top, datatype='mean count', link.arr.lwd=res.top$cell_from_mean_exprs,
       cell_col=cell_col, link.arr.width=res.top$cell_to_mean_exprs)
dev.off()
setwd("~/project/2020/2007_10xDemo2")

对比分组之间差异的配体-受体 

##==细胞通讯差异分析==##
data2 <- subset(data.italk, subset=data.italk$compare_group=="HNC_TIL"|data.italk$compare_group=="Tonsil")
##各个细胞类型分别执行差异分析，然后合并在一起
deg.genes <- DEG(data2, method='Wilcox', contrast=c("HNC_TIL","Tonsil"))
deg_B <- DEG(data2 %>% filter(cell_type=='B cells'), method='Wilcox',contrast=c("HNC_TIL","Tonsil"))
deg_T <- DEG(data2 %>% filter(cell_type=='T cells'), method='Wilcox',contrast=c("HNC_TIL","Tonsil"))
deg_cd4T <- DEG(data2 %>% filter(cell_type=='CD4+ T cells'), method='Wilcox', contrast=c("HNC_TIL","Tonsil"))
deg_3 <- rbind(deg_B, deg_cd4T, deg_T)

##各个配体-受体分类分别作图
dir.create('iTALK/DEG')
setwd("iTALK/DEG")
res<-NULL
for(comm_type in comm_list){
    #comm_type='other' 
  #多个细胞类型之间显著表达的配体-受体，结果会过滤同一细胞类型内的配体-受体关系
  res_cat <- FindLR(deg_3, datatype='DEG', comm_type=comm_type)
  #单个细胞类型显著表达的配体-受体
  #res_cat <- FindLR(data_1=deg_cd4T, datatype='DEG', comm_type=comm_type)
    res_cat <- res_cat[order(res_cat$cell_from_logFC*res_cat$cell_to_logFC, decreasing=T),]
  write.csv(res_cat, paste0('LRpairs_DEG_',comm_type,'.xls'))
    #plot by ligand category
    png(paste0('LRpairs_DEG_',comm_type,'_circ.png'), width=600, height=650)
  if(nrow(res_cat)==0){
        next
    }else if(nrow(res_cat>=PN)){
        LRPlot(res_cat[1:PN,], datatype='DEG', link.arr.lwd=res_cat$cell_from_logFC[1:PN],
           cell_col=cell_col, link.arr.width=res_cat$cell_to_logFC[1:PN])
    }else{
        LRPlot(res_cat, datatype='DEG', link.arr.lwd=res_cat$cell_from_logFC,
           cell_col=cell_col, link.arr.width=res_cat$cell_to_logFC)
    }
  dev.off()
    png(paste0('LRpairs_DEG_',comm_type,'_net.png'), width=600, height=650)
  NetView(res_cat, col=cell_col, vertex.label.cex=1.2, edge.label.cex=0.9, 
          vertex.size=30, arrow.width=3, edge.max.width=10, margin = 0.2)
    dev.off()
    res<-rbind(res,res_cat)
}

##所有配体-受体分类一起作图
res<-res[order(res$cell_from_logFC*res$cell_to_logFC, decreasing=T),]
write.csv(res, 'LRpairs_DEG_all.xls')
if(is.null(res)){
    print('No significant pairs found')
}else if(nrow(res)>=PN){
    res.top <- res[1:PN,]
    png(paste0('LRpairs_DEG_all_net.png'), width=600, height=650)
  NetView(res, col=cell_col, vertex.label.cex=1.2, edge.label.cex=0.9, 
          vertex.size=30, arrow.width=3, edge.max.width=10, margin = 0.2)
  dev.off()    
    png('LRpairs_DEG_all_circ.png', width=600, height=650)
    LRPlot(res.top, datatype='DEG', link.arr.lwd=res.top$cell_from_logFC,
       cell_col=cell_col, link.arr.width=res.top$cell_to_logFC)
    dev.off()
}else{
    png(paste0('LRpairs_DEG_all_net.png'), width=600, height=650)
  NetView(res, col=cell_col, vertex.label.cex=1.2, edge.label.cex=0.9, 
          vertex.size=30, arrow.width=3, edge.max.width=10, margin = 0.2)
  dev.off()    
    png('LRpairs_DEG_all_circ.png', width=600, height=650)
    LRPlot(res, datatype='DEG', link.arr.lwd=res$cell_from_logFC,
       cell_col=cell_col, link.arr.width=res$cell_to_logFC)
    dev.off()
}
setwd("~/project/2020/2007_10xDemo2")

iTALK图例的问题 

iTALK画的图不带图例，只能从作者原文里截图啦😅

获取帮助

本教程的目的在于把常用的单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以点击 “阅读原文” 找到作者联系方式，获取帮助。

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

• 生信爆款入门-第8期（线上直播4周，马拉松式陪伴，带你入门）你的生物信息入门课
  

• 数据挖掘学习班第6期（线上直播3周，马拉松式陪伴，带你入门） 医学生/医生首选技能提高课

• 生信技能树的2019年终总结 你的生物信息成长宝藏

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注