专栏首页生信菜鸟团如果你觉得相关性热图不好看,或者太简陋

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

前面的教程:混合到同一个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$RNA@var.features
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

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

本文分享自微信公众号 - 生信菜鸟团(bio_123456789),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-07-23

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 使用SingleR构建自定义细胞亚群数据库

    如何很多朋友留言问,为什么不使用现成的工具呢,比如SingleR就构建自定义细胞亚群数据库。我们当然知道这样的工具很好用,但是我们要分享的是技术细节,如果一切都...

    生信菜鸟团
  • 曾经我是一个只会Excel的数据分析师,直到我遇到了……

    算法与数据结构 我是一个数据分析师。 准确来说我是一个当年只会excel数据透视表,就天不怕地不怕地来当数据分析师的人。当年的某一天,我的老板Q我: 小刘啊,我...

    前朝楚水
  • 月末总结与推书

    潘成涛
  • 国外程序员陋习,写在农历狗年前国产程序员陋习,写在农历猴年前

    为了呼应下面以前写的文章,今天来写写国外程序员的陋习(当然不是所有人都是这样,但是也应该是我碰到比较多的) “国产程序员陋习,写在农历猴年前” http://w...

    麦克-堂
  • 数据可视化三步理论,好的图表会说话

    数据可视化,是目前非常热门的方向,所谓“文不如表,表不如图”,“一图胜千言”,在工作中,数据可视化变得越来越重要。

    挖数
  • 你为什么需要 Kotlin

    导语 在当今的互联网时代,新技术犹如雨后春笋般层出不穷。精神哥之前也和开发同学一起讨论过程序员的成长离不开哪些软技能?当时很多人都有提到探究新技术对程序员的重要...

    腾讯Bugly
  • PyTorch 到底好用在哪里?

    提问内容如下: 之前非常熟悉 Tensorflow,后来都说 PyTorch 简单易上手,自己就去试了试。 PyTorch 连最基本的 maximum, min...

    AI研习社
  • #夏日编程团#天这么热,没法出去浪,不如趁着宅在空调房里的这两个月,跟我们组团点亮你的编程技能点

    最近这天热得不像样子,一非洲游客在天安门广场中暑晕倒……还是选择老老实实窝在房间里吹空调吧。 不过也有些人,学习的热情比气温还高,不写代码就不舒服!(因为会被助...

    Crossin先生
  • 如何学习计算OpenCV

    如何学习OpenCV 一:学习OpenCV三个阶段 人工智能带火了计算机视觉的人才需求,作为计算机视觉应用开发框架OpenCV也越来越受到欢迎,市场需求大增,...

    OpenCV学堂
  • 业界 | 人工智能看走眼的图像都长什么样?

    选自the Verge 作者:James Vincent 机器之心编译 参与:Ellen Han、黄小天 威廉·吉布森(William Gibson)写于 2...

    机器之心
  • 【实战案例】如何利用大数据思维在北京租到好房子?

    PPV课大数据 第一步:精准定位。 确定找房地点,精确到小区。每个小区在任意时间,至少有三五间空房待租。大的小区,有几十间。完全不要担心没房。如果没有,基本是因...

    小莹莹
  • 「推荐」从openresty谈到rust

    本文转载自知乎:https://zhuanlan.zhihu.com/p/87922545?utm_source=wechat_session&utm_medi...

    MikeLoveRust
  • 代码命名:僧敲月下门

    忽一日於驴上吟得:‘鸟宿池中树,僧敲月下门。’初欲著‘推’字,或欲著‘敲’字,炼之未定,遂于驴上作‘推’字手势,又作‘敲’字手势。 - 《鉴戒录·贾忤旨》 两句...

    tyrchen
  • 加强版!如果编程语言是车,那么你开的是……

    世界上,总是充满活力的人,热衷于创造新语言,并不遗余力地推介,开大会,开专栏,立项目,开论坛,只求开发者能注意到:“嗨,这儿有一玩杂耍的,看着飞刀嗖嗖嗖~~~喷...

    程序员小助手
  • 拼凑6个网页工具图表还不够那就再加上6个组学

    TCGA数据挖掘真的是绵绵不绝,这里就不再赘述了,从基因集到ceRNA,到可变剪切,肿瘤免疫, 再到现在的m6A和自噬基因, 马上缺氧,代谢应该是也要出来了,每...

    生信技能树
  • 敢挑战吗?这30个以太坊开发示例,让你成为80万都挖不走的区块链人才!

    我曾经买过加密货币,曾试图使用一些丑陋矿机挖矿,看过一些稀稀拉拉的Solidity教程。但不得不承认,在当时,我更偏爱前者,我切身体会到了加密货币的狂热,急切需...

    区块链大本营
  • iOS开发工具篇-AppStore统计工具 (转载)

    随着iOS开发的流行,针对iOS开发涉及的方方面面,早有一些公司提供了专门的解决方案或工具。这些解决方案或工具包括:用户行为统计工具(友盟,Flurry,Goo...

    tandaxia
  • 什么样的代码是好代码?

    即Single Responsibility (单一职责),Open Close(开闭),Liskov Substitution(里氏替换),Interface...

    NaughtyCat
  • 免费创建个人静态网站最佳实践:hugo+github+netlify

    关于搭建一个博客或个人网站的好处不用我多说,但创建网站的难度可能会让人望而却步。本人从网络上获得过很多帮助,学到很多。很早就萌发了自己建网站并分享知识的想法,但...

    风不止

扫码关注云+社区

领取腾讯云代金券