如果你觉得相关性热图不好看，或者太简陋

前面的教程：混合到同一个10X样品里面的多个细胞系如何注释，我们提到了可以使用细胞系的表达量矩阵去跟细胞亚群表达量矩阵进行相关性计算。

就有粉丝提问，把单细胞亚群使用   AverageExpression 函数做成为了亚群矩阵，是不是忽略了单细胞亚群的异质性呢？毕竟每个单细胞亚群背后都是成百上千个具体的细胞啊。代码如下所示：

rm(list = ls())
library(SingleR)
library(Seurat)
library(ggplot2)

# 读入scRNA数据 -------
scRNA <- readRDS("../step1_聚类/sce_all.Rds")
table(Idents(scRNA) )
Idents(scRNA) <- "RNA_snn_res.0.2"
table(Idents(scRNA) )


可以看到，我们的8个细胞亚群，确实是都是有几百个细胞的，使用   AverageExpression 函数会被融合成为了一个样品：

> table(Idents(scRNA) )

0   1   2   3   4   5   6   7 
976 881 511 510 493 487 440 316 


而且原文很明显使用的是具体的单细胞亚群里面的全部细胞去和各个细胞系进行计算相关性，然后绘制密度图：

相关性密度图

其实也很容易实现啊！

首先拿到全部的细胞系和全部的具体的每个单细胞的表达相关性矩阵（Pearson correlation coefficient），代码如下：

sce <- readRDS("../step1_聚类/sce_all.Rds")
load(file = 'for_cor.Rdata')
head(reference)
head(query)
cg=sce@assays$[e-mail]
ids=intersect(rownames(reference),cg)
ids=rownames(reference)
ct=sce@assays$RNA@data[ids,]
ct[1:4,1:4]
# 这个时候不能使用   AverageExpression 函数把表达量矩阵融合

table(sce$RNA_snn_res.0.2)
data_merged=cbind(ct,reference)
str(data_merged)
tmp=apply(data_merged, 2,as.numeric)
colnames(tmp)=colnames(data_merged)
rownames(tmp)=rownames(data_merged)
cor_df=cor(tmp)
cor_df[1:4,1:4]
dim(cor_df)
dim(ct)
tt_df = cor_df[1:4614,4615:ncol(cor_df)]
tt_df=as.data.frame(tt_df)
tt_df$clusterID = sce$RNA_snn_res.0.2
head(tt_df)
colnames(tt_df)
colnames(tt_df)[c(1:6 )]

data=tt_df


如下所示，每个具体的单细胞，都有各自对应的多个细胞系的表达相关性矩阵（Pearson correlation coefficient），：

 

然后就可以对这个 Pearson correlation coefficient 矩阵，进行密度图可视化，构造一个函数：

plot_density <- function(clusterID){
  samplenames <- colnames(data)[c(1:6)]
  col <- c('#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02') 
  plot(density(data[which(data$clusterID == clusterID),1]), 
       xlim = c(0.2,0.8),ylim = c(0, 20), lwd=5,
       las=2, main="", xlab="", col=col[1])
  title(main=paste0("cluster", clusterID), xlab="Pearson correlation coefficient")
  den <- density(data[which(data$clusterID == clusterID),2])
  lines(den$x, den$y, col=col[2], lwd=5)
  den <- density(data[which(data$clusterID == clusterID),3])
  lines(den$x, den$y, col=col[3], lwd=5)
  den <- density(data[which(data$clusterID == clusterID),4])
  lines(den$x, den$y, col=col[4], lwd=5)
  den <- density(data[which(data$clusterID == clusterID),5])
  lines(den$x, den$y, col=col[5], lwd=5)
  den <- density(data[which(data$clusterID == clusterID),6])
  lines(den$x, den$y, col=col[6], lwd=5)
  # legend("right", samplenames, text.col=col, bty="n")
}

使用上面的 plot_density 函数 先看一个细胞亚群里面的全部的细胞的相关性系数的密度图：

#画图例并与上图手动合并
col <- c('#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02')
samplenames <- colnames(data)[ 1:6]
plot_density(clusterID = 0)
legend("right", samplenames, text.col=col, bty="n")

如下所示：

单个细胞亚群的全部细胞的相关性密度

可以看到，很清晰的显示了这个cluster0的细胞亚群就是MCF7（我截图的时候，眼神有点飘，大家忽略哈）

我们前面步骤的 混合到同一个10X样品里面的多个细胞系如何注释，手动注释如下：

# st ep4. 手动注释细胞类型------
# cluster0  MCF7
# cluster1  HEK293
# cluster2  T47D
# cluster3  BCK4 # 排除法
# cluster4  T47D
# cluster5  MCF10A
# cluster6  MM134
# cluster7  SUM44


基本上是一致的！

然后画全部的细胞亚群：

pdf(file = "density_cluster.pdf", width = 16, height = 6)
par(mfrow=c(2,4))
plot_density(clusterID = 0)
plot_density(clusterID = 1)
plot_density(clusterID = 2)
plot_density(clusterID = 3)
plot_density(clusterID = 4)
plot_density(clusterID = 5)
plot_density(clusterID = 6)
plot_density(clusterID = 7)
dev.off()


如下所示，基本上跟前面的步骤是一样的结果啦；

细胞亚群

手动注释如下：

# st ep4. 手动注释细胞类型------
# cluster0  MCF7
# cluster1  HEK293
# cluster2  T47D
# cluster3  BCK4 # 排除法
# cluster4  T47D
# cluster5  MCF10A
# cluster6  MM134
# cluster7  SUM44


基本上是一致的！

也可以使用小提琴图来展示：


# 二、可视化之——小提琴图 -------
# 1.数据
data_used <- data[,c(1:6,13)]
# 宽变长
data_used <- tidyr::gather(data_used, key = "celltype", value = "avg_expr", -clusterID)
head(data_used)
# 2.绘图
library(ggplot2)
p <- ggplot(data_used, aes(x = clusterID, y = avg_expr, fill = celltype)) +
  geom_violin() +
  theme(plot.margin=margin(t = 0, b = 0, unit="pt"),
        plot.subtitle = element_text(family = "serif",  colour = "gray0"),
        plot.background = element_rect(fill = "aliceblue"),
        plot.title = element_text(family = "serif", hjust = 1, vjust = 0, colour = "forestgreen"),
        #  背景板
        panel.grid.minor = element_line(linetype = "blank"),
        panel.grid.major.y = element_line(colour = "coral", linetype = "dashed"),
        panel.background = element_rect(fill = "white"),
        # 坐标轴标题
        axis.title.x = element_text(family = "serif",face = "bold",colour = "chocolate4", hjust = 0.5, vjust = 0),
        axis.title.y = element_text(family = "serif",face = "bold",colour = "chocolate4", hjust = 0.5, vjust = 0.5),
        # 图例
        legend.title = element_text(hjust = 0.55,face = "bold", colour = "chocolate4",family = "serif"),
        legend.direction = "vertical",
        legend.text = element_text(family = "serif"),
        legend.key = element_rect(fill = "aliceblue"),
        legend.background = element_rect(fill = "aliceblue")) +
  facet_wrap(vars(clusterID), nrow = 2, scales = "free")
p
ggsave(p, filename = "violin.pdf", width = 14, height = 6)

如下所示：

小提琴图

手动注释如下：

# st ep4. 手动注释细胞类型------
# cluster0  MCF7
# cluster1  HEK293
# cluster2  T47D
# cluster3  BCK4 # 排除法
# cluster4  T47D
# cluster5  MCF10A
# cluster6  MM134
# cluster7  SUM44


基本上是一致的！（看到这里，你有没有觉得这个流程跟某个软件的算法类似呢？）

这个文献及其数据处理也会纳入我们的两个b站单细胞栏目，持续更新半年，基本上涵盖了大家需要的技能：

• https://www.bilibili.com/video/BV1DK4y1X7bb/  更新至第8篇，「生信技能树」100个单细胞文献解读(8/100)
• https://www.bilibili.com/video/BV19Q4y1R7cu/  section 3已更新，「生信技能树」单细胞公开课2021

如果是简单的降维聚类分群，可以参考前面的例子：人人都能学会的单细胞聚类分群注释  ，我们演示了第一层次的分群。

如果你对单细胞数据分析还没有基础认知，可以看基础10讲：

• 01. 上游分析流程
• 02.课题多少个样品，测序数据量如何
• 03. 过滤不合格细胞和基因（数据质控很重要）
• 04. 过滤线粒体核糖体基因
• 05. 去除细胞效应和基因效应
• 06.单细胞转录组数据的降维聚类分群
• 07.单细胞转录组数据处理之细胞亚群注释
• 08.把拿到的亚群进行更细致的分群
• 09.单细胞转录组数据处理之细胞亚群比例比较