马拉松授课环节我们的单细胞转录组数据主要是降维聚类分群,以及各个单细胞亚群注释,然后走三个高级分析,拟时序,细胞通讯,转录因子注释。但是很多人有自己的生物学背景,会考虑不同的生物学功能基因集的打分,以及不同的基因的相关性,试图去从单细胞水平提出自己的生物学见解。
不过,我这里可能需要提醒大家一定要注意, 如果是真正的单细胞水平相关性分析,无论怎么做都是不合理的。比如大家如果使用代码绘制任意两个基因在单细胞转录组数据里面的相关性散点图,就出现这种多个点排列成一条一条:
多个点排列成一条一条
简单的示例代码是:
library(SeuratData) #加载seurat数据集
library(Seurat)
getOption('timeout')
options(timeout=10000)
# InstallData("pbmc3k")
data("pbmc3k")
sce table(sce$seurat_annotations)
library(tidyverse)
colnames([email protected])
dim(sce)
sce % NormalizeData %>% FindVariableFeatures %>% ScaleData
# devtools::install_github("LKremer/ggpointdensity")
library(ggpointdensity)
selected_genes 'GJA4', 'GJA5', 'EFNB2', # Art_marker
'DLL4', 'NOTCH1', 'JAG1', 'JAG2', 'RBPJ', # Notch_marker
'BCL2', 'NR2F2') # 治病基因
selected_genes = selected_genes[selected_genes %in% rownames(sce)]
expr1 = t(as.data.frame(sce@assays$RNA@counts[c('CD3D','CD4'),]))
expr2 = t(as.data.frame(sce@assays$RNA@data[c('CD3D','CD4'),]))
expr3 = t(as.data.frame(sce@[email protected][c('CD3D','CD4'),]))
draw function(expr){
expr %>%
ggplot(aes(x = CD3D, y = CD4)) +
ggrastr::rasterise(ggpointdensity::geom_pointdensity(size = 0.5), dpi = 600) +
geom_smooth(method = "lm", formula = y ~ x,
color = "#6b76ae", fill = "#e5a323",
size = 0.5, alpha = 0.2) +
theme_bw() +
theme(aspect.ratio = 1, legend.margin = margin(l= -8),
axis.title = element_text(face = "italic", size = 12)) +
guides(color = guide_colorbar(
frame.colour = "black", frame.linewidth = 0.5,
ticks.colour = "black", title = "Density")) +
scale_color_viridis_c() +
ggpubr::stat_cor(method = "pearson") # 添加相关性
}
p1=draw(expr1)
p2=draw(expr2)
p3=draw(expr3)
p1+p2+p3
初学者往往是纠结于到底使用哪个表达量矩阵,比如上面的代码里面的 ,
expr1
、
expr2
和
expr3
代表了不同类型的表达量矩阵,它们各自有不同的特点和用途:
-
-
expr1
是从
sce
对象的
counts
槽中提取的原始计数矩阵的转置(
t
)。这个矩阵包含了每个基因在每个细胞中的原始读段(UMI或读段)计数。这些计数是未经标准化或缩放的,因此它们反映了测序深度和细胞捕获效率的直接测量。
-
expr2
是从
sce
对象的
data
槽中提取的表达矩阵的转置。这个矩阵包含了经过预处理的表达值,通常是经过标准化的。在Seurat中,
data
槽通常存储的是经过预处理步骤(如去除批次效应、标准化和缩放)后的表达数据。这些数据更适合进行下游分析,如差异表达分析和聚类。
-
expr3
是从
sce
对象的
scale.data
槽中提取的表达矩阵的转置。这个矩阵包含了经过缩放的表达值,通常用于可视化和聚类。在Seurat中,
scale.data
槽存储的是经过
ScaleData
函数处理后的数据,该函数会对数据进行中心化和缩放,使其更适合进行多样本比较和可视化。
总结来说,这三个表达量矩阵的主要区别在于它们在数据预处理流程中的位置和用途:
-
expr1
提供了原始的计数数据,适合进行质量控制和初步探索。
-
expr2
提供了预处理后的数据,适合进行差异表达分析。
-
expr3
提供了缩放后的数据,适合进行可视化和聚类分析。
在使用这些数据时,需要根据分析目的选择合适的表达量矩阵。例如,如果你需要对数据进行标准化处理,可能会选择
expr2
或
expr3
。如果你需要评估测序深度或进行质量控制,可能会选择
expr1
。
但是,实际上,如果大家是做真正的单细胞水平相关性分析,无论取什么样的表达量矩阵怎么做都是不合理的!
跨越细胞亚群的相关性
大家很少会关心这个维度的相关性, 因为单细胞转录组最重要的研究单位就是细胞亚群,何苦再把大家混合起来呢!虽然说,从代码角度很容易做出来,比如我有这样的一个单细胞表达量矩阵降维聚类分群结果:
table(phe[,c("orig.ident","celltype"),])
av group.by = c("orig.ident","celltype"),
assays = "RNA")
av=as.data.frame(av[[1]])
> table(phe[,c("orig.ident","celltype"),])
celltype
orig.ident Bcells endo epi fibro mast myeloids plasma SMC Tcells
P01_OM 14 349 678 3786 1 240 2 409 204
P01_PM 2349 80 7 252 19 419 29 56 2791
P01_PT 8 20 1005 39 22 52 60 14 343
P02_AS 191 0 701 1214 1 2291 0 10 496
P02_OM 9 237 472 1667 7 492 9 277 330