无论怎么做都是错误的单细胞水平相关性分析

生信技能树

发布于 2024-11-21 11:18:35

1.2K0

马拉松授课环节我们的单细胞转录组数据主要是降维聚类分群，以及各个单细胞亚群注释，然后走三个高级分析，拟时序，细胞通讯，转录因子注释。但是很多人有自己的生物学背景，会考虑不同的生物学功能基因集的打分，以及不同的基因的相关性，试图去从单细胞水平提出自己的生物学见解。

不过，我这里可能需要提醒大家一定要注意，如果是真正的单细胞水平相关性分析，无论怎么做都是不合理的。比如大家如果使用代码绘制任意两个基因在单细胞转录组数据里面的相关性散点图，就出现这种多个点排列成一条一条：

多个点排列成一条一条

简单的示例代码是：

library(SeuratData) #加载seurat数据集  
library(Seurat)
getOption('timeout')
options(timeout=10000)
# InstallData("pbmc3k")  
data("pbmc3k")  
sce <- UpdateSeuratObject(pbmc3k) 
table(sce$seurat_annotations)
library(tidyverse)
colnames(sce@meta.data)
dim(sce)
sce <- sce %>% NormalizeData %>% FindVariableFeatures %>% ScaleData

# devtools::install_github("LKremer/ggpointdensity")
library(ggpointdensity) 
selected_genes <- c('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@assays$RNA@scale.data[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:
- expr1是从sce对象的counts槽中提取的原始计数矩阵的转置（t）。这个矩阵包含了每个基因在每个细胞中的原始读段（UMI或读段）计数。这些计数是未经标准化或缩放的，因此它们反映了测序深度和细胞捕获效率的直接测量。
expr2:
- expr2是从sce对象的data槽中提取的表达矩阵的转置。这个矩阵包含了经过预处理的表达值，通常是经过标准化的。在Seurat中，data槽通常存储的是经过预处理步骤（如去除批次效应、标准化和缩放）后的表达数据。这些数据更适合进行下游分析，如差异表达分析和聚类。
expr3:
- expr3是从sce对象的scale.data槽中提取的表达矩阵的转置。这个矩阵包含了经过缩放的表达值，通常用于可视化和聚类。在Seurat中，scale.data槽存储的是经过ScaleData函数处理后的数据，该函数会对数据进行中心化和缩放，使其更适合进行多样本比较和可视化。

总结来说，这三个表达量矩阵的主要区别在于它们在数据预处理流程中的位置和用途：

expr1提供了原始的计数数据，适合进行质量控制和初步探索。
expr2提供了预处理后的数据，适合进行差异表达分析。
expr3提供了缩放后的数据，适合进行可视化和聚类分析。

在使用这些数据时，需要根据分析目的选择合适的表达量矩阵。例如，如果你需要对数据进行标准化处理，可能会选择expr2或expr3。如果你需要评估测序深度或进行质量控制，可能会选择expr1。

但是，实际上，如果大家是做真正的单细胞水平相关性分析，无论取什么样的表达量矩阵怎么做都是不合理的！

跨越细胞亚群的相关性

大家很少会关心这个维度的相关性，因为单细胞转录组最重要的研究单位就是细胞亚群，何苦再把大家混合起来呢！虽然说，从代码角度很容易做出来，比如我有这样的一个单细胞表达量矩阵降维聚类分群结果：

table(phe[,c("orig.ident","celltype"),])
av <-AggregateExpression(sce.all ,
                         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
    P02_PM    288  519  288   668   75      716     12  337   2507
    P02_PT    191  148  439    19   44      125    143  103    815

我们可以简单的看看管家基因（Housekeeping genes）是否有跨越细胞亚群的高度相关性，是一类在细胞中表达水平相对稳定的基因，它们通常参与细胞的基本生理过程，如能量代谢、DNA复制和细胞骨架维持：

# 创建一个包含管家基因名称的向量
housekeeping_genes <- c("ACTB", "GAPDH", "UBC", "HPRT1", "PPIA", "RPLP0", "GUSB", "SDHA", "TFRC", "PGK1", "PEPD", "ALDOA", "ENO1", "G6PD")
housekeeping_genes=housekeeping_genes[housekeeping_genes %in% rownames(av)]
## Top3 genes
load('check-by-celltype/qc-_marker_cosg.Rdata')
top_3 <- unique(as.character(apply(marker_cosg$names,2,head,3)))
df =cor(t(as.matrix(log(av[c(housekeeping_genes,top_3),]+1))))
colnames(df)
pheatmap::pheatmap(df)

请注意，管家基因的概念在近年来受到了一些质疑，因为一些传统上被认为是管家基因的基因在不同条件下的表达水平可能会有显著变化。因此，在使用管家基因进行标准化或数据校正时，需要谨慎，并考虑实验的具体条件和目的。此外，不同物种的管家基因可能不同，上述列表主要基于人类基因。

可以很清晰的看到，多个管家基因在多个细胞亚群的相关性比较好，然后就是每个亚群各自的top3基因：

image-20241115210132544

这个也是合情合理的，上面的热图里面的每个单元格其实都是一个相关性散点图，它们就不会出现上面的多个点排列成一条一条的情况；

df=as.data.frame(t(log(av[c("ACTB", "GAPDH"),]+1)))
head(df)
library(ggpubr)
ggscatter(df,x="ACTB",y="GAPDH",
          color="black",shape=21,size=3,#Pointscolor,shapeandsize
          add="reg.line",#Addregressinline
          add.params=list(color="blue",fill="lightgray"),#Customizereg.line
          conf.int=TRUE,#Addconfidenceinterval
          cor.coef=TRUE,#Addcorrelationcoefficient.see?stat_cor
          cor.coeff.args=list(method="pearson",label.sep="\n")
)

如下所示几乎完美的相关性：