专栏名称: 生信技能树
生物信息学学习资料分析,常见数据格式及公共数据库资料分享。常见分析软件及流程,基因检测及癌症相关动态。
目录
相关文章推荐
中国基金报  ·  重大重组,复牌! ·  15 小时前  
中国基金报  ·  华为,AI大消息! ·  昨天  
中国基金报  ·  “元宇宙第一股”,暴跌! ·  2 天前  
哔哩哔哩  ·  用300天从零制作游戏,最后上架Steam ·  2 天前  
中国基金报  ·  这家“独角兽”与OpenAI“分手”了! ·  3 天前  
51好读  ›  专栏  ›  生信技能树

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

生信技能树  · 公众号  ·  · 2024-11-15 21:21

正文

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

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

多个点排列成一条一条

简单的示例代码是:

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 代表了不同类型的表达量矩阵,它们各自有不同的特点和用途:

  1. 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                          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






    请到「今天看啥」查看全文


    推荐文章
    中国基金报  ·  重大重组,复牌!
    15 小时前
    中国基金报  ·  华为,AI大消息!
    昨天
    中国基金报  ·  “元宇宙第一股”,暴跌!
    2 天前
    中国基金报  ·  这家“独角兽”与OpenAI“分手”了!
    3 天前
    Python开发者  ·  Git 版本控制与工作流
    7 年前