单细胞转录组+空间转录组联合分析系列教程
大家好！今日给大家分享新出炉的软件NeuronChat，看名字大概知道这是一类预测神经细胞间的细胞通讯软件。针对单细胞转录组分析，大家较为广泛使用CellChat、CellPhoneDB、NicheNet等来预测细胞通讯信号，然而大家在用于特定的组织或细胞类别，可能会觉得分析得到的配受体对不够或不足以表达想说明的生物学意义等困惑，笔者近期分析也是遇到如此情况，尤其是在具有神经元的组织中。而这篇新出炉的比较特异性预测神经元细胞间通讯信号的软件能部分解决此类涉及到神经组织细胞内通讯的问题，所以分享给大家。该文章于2023年2月28日发表在Nature Communications上标题为《Inferring neuron-neuron communications from single-cell transcriptomics through NeuronChat》。作者将该软件功能归结于3块，笔者逐篇分享，供大家参考。
本软件实质性还是作者构建了一个特异性由神经细胞间信号组成的数据库，比如小分子神经递质、神经肽、缝隙连接蛋白、气体递质和突触黏附分子，而靶标通常但不限于受体。例如，神经递质的靶蛋白也可以是摄取转运蛋白或失活酶；缝隙连接蛋白的靶蛋白是其他相容的缝隙连接蛋白；对于非肽类神经递质，条目中包括相应的合成酶和/或囊泡转运蛋白；对于含有多个不同亚基的异构体受体，对应的亚基被归类到具有相同配体的不同条目中。在373个配体-靶标相互作用对条目中，分别有221个、73个、39个、16个和24个条目涉及小分子神经递质、神经肽、缝隙连接蛋白、气体递质和突触黏附分子。
接下来分享第一篇代码示例 ：
NeuronChat（一）：特异性神经细胞间通讯推断和分析
devtools::install_github("Wei-BioMath/NeuronChat")
library(NeuronChat)
library(CellChat)
library(ggalluvial)
###Part I: 创建NeuronChat对象
# 作者使用的是皮质层数据来自Yao 2021;也可用自己的rds数据
data(list='cortex_data') #or data <- readRDS("xx.rds")
# 作者取了数据中VISp 区域的子集
region_name <- 'VISp'; cell_class <- names(table(meta$class_label))
cell_idx <- which(meta$region_label %in% region_name & meta$class_label %in% cell_class & !(meta$subclass_label %in%c('Car3','CR','DG','L2/3 IT PPP','L5/6 IT TPE-ENT')))
target_df_single <- target_df[cell_idx,] # a data frame: row cell, column gene (the last column is the cell subclass)
meta_tmp <- meta[cell_idx,];rownames(meta_tmp) <- meta_tmp$sample_name # subset meta
df_group <- meta_tmp[!duplicated(meta_tmp$subclass_label),c('class_label','subclass_label')]
group <- structure(df_group$class_label,names=df_group$subclass_label) # create the cell class for the subclasses, used for plot
# creat NeuronChat object; choose the database 'mouse' for mouse data; 'human' for human data
# note that the first parameter should be a numeric matrix with row gene and column cell
x <- createNeuronChat(t(as.matrix(target_df_single[,1:(dim(target_df_single)[2]-1)])),DB='mouse',group.by = target_df_single$cell_subclass);
#> Create a NeuronChat object from a data matrix
###Part2 II：执行NeuronChat来推断神经元特异性的细胞间通讯网络
# M is for the permutation test; typically ~4 mins when M=100, 依据数据集大小和细胞数目而定
#作者文章中测试过M=100和M=1000准确性类似，而且计算速度还快
x <- run_NeuronChat(x,M=100)
#> Time difference of 5.181635 mins
# the the communication networks for individual interaction pairs are stored in slot 'net'
# aggregate the communication networks over all interaction pairs, method can be 'weight', 'count' and so on
net_aggregated_x <- net_aggregation(x@net,method = 'weight')
###Part III：可视化神经特异性细胞间的通讯网络
#for aggregated network
par(mfrow=c(1,2))
# Visualization, circle plot, for the aggregated network
netVisual_circle_neuron(net_aggregated_x,group=group,vertex.label.cex = 1)
# Visualization, chordDiagram, for the aggregated network; also using cellchat function netVisual_chord_cell_internal(net_aggregated_x, group = group,lab.cex=1)
netVisual_chord_neuron(x,method = 'weight',group=group,lab.cex = 1)
图1.png
# Visualization, heatmap for the aggregated network
heatmap_aggregated(x, method='weight',group=group)
图2.png
#for individual network
par(mfrow=c(1,2))
# Visualization for the single interaction pair, circle plot
netVisual_circle_neuron(x@net$Glu_Gria2,group=group,vertex.label.cex = 1)
# Visualization for the single interaction pair, chord diagram
netVisual_chord_neuron(x,interaction_use='Glu_Gria2',group=group,lab.cex = 1)
图3.png
# Visualization for the single interaction pair, heatmap
heatmap_single(x,interaction_name='Glu_Gria2',group=group)
图4.png
# Visualization for the single interaction pair, heatmap with violin plots showing expression of genes realted to ligand and target
lig_tar_heatmap(x,interaction_name='Glu_Gria2',width.vector=c(0.38,0.35,0.27))
图5.png
###Part IV：细胞通讯配受体对分析
#barplot show the count of links or information flow for all interaction pairs
g1 <- rankNet_Neuron(x,slot.name = "net",measure = c("weight"),mode='single',font.size = 5)
g2 <- rankNet_Neuron(x,slot.name = "net",measure = c("count"),mode='single',font.size = 5)
g1+g2
图6.png
#outgoing/incoming pattern
selectK_Neuron(x,pattern = "outgoing")
selectK_Neuron(x,pattern = "incoming")
x<- identifyCommunicationPatterns_Neuron(x, slot.name = "net", pattern = c("outgoing"), k=4,height = 18)
x<- identifyCommunicationPatterns_Neuron(x, slot.name = "net", pattern = c("incoming"), k=4,height = 18)
图7.png
library(ggalluvial)
g3 < netAnalysis_river_Neuron(x,slot.name = "net", pattern = c("outgoing"),font.size = 2.5,cutoff.1 = 0.5,cutoff.2=0.5)
g4 < netAnalysis_river_Neuron(x,slot.name = "net", pattern = c("incoming"),font.size = 2.5,cutoff.1 = 0.5,cutoff.2=0.5)
g3+g4
#> Please make sure you have load `library(ggalluvial)` when running this function
图合并.png
#manifold learning and classification
x <- computeNetSimilarity_Neuron(x,type='functional')
x <- CellChat::netEmbedding(x, slot.name = "net_analysis", type = "functional")
#> Manifold learning of the signaling networks for a single dataset
x <- CellChat::netClustering(x, type='functional',slot.name = "net_analysis",k=5)
#> Classification learning of the signaling networks for a single dataset
netVisual_embedding_Neuron(x, type = "functional", label.size = 5,pathway.remove.show = F)
图10.png
netVisual_embeddingZoomIn_Neuron(x, type = "functional", nCol = 2,label.size = 3)
图11.png
###R语言依赖包版本：
R >= 4.1.0
dplyr >= 1.0.9
data.table >= 1.14.2
CellChat >= 1.1.3
Seurat >= 4.1.0
SeuratObject >= 4.1.0
NMF >= 0.23.0
igraph >= 1.3.4
ggplot2 >= 3.3.6
ComplexHeatmap >= 2.8.0
circlize >= 0.4.14
ggalluvial >= 0.12.3
