单细胞测序分析(三)质控和过滤(动态参数)

45 阅读3分钟

大家好,欢迎回到我们的单细胞测序分析系列!在上一篇推文中,我们介绍了如何从不同文件构建Seurat对象,并为进一步的分析做好准备。今天,我们将深入探讨在单细胞数据分析中至关重要的一个环节——质控和过滤。让我们一起来看看为什么质控和过滤对数据分析至关重要,以及如何在Seurat中进行这些操作。

质控和过滤的重要性

单细胞RNA测序能够为我们提供每个细胞的基因表达谱,但同时也会受到技术噪声、细胞质量不佳或其他外部因素的影响。因此,进行质控和过滤是保证数据质量的第一步。通过去除低质量细胞和不可靠的数据点,可以确保我们的分析结果更为精确和可靠。

质控和过滤的常见标准包括:

  • 细胞总RNA数目:极低的RNA数可能意味着细胞已经死亡或没有充分的RNA转录。
  • 基因数目:如果一个细胞只表达极少的基因,它很可能是低质量的细胞。
  • 线粒体基因的比例:线粒体基因的比例过高可能表明细胞正在经历压力或死亡。
  • 核糖体基因的比例:核糖体基因在细胞中占有重要位置,过高的核糖体基因表达可能与细胞的活跃状态和代谢水平有关,反之可能是细胞质量不佳的标志。
  • 红血细胞基因的比例:对于一些类型的样本,红血细胞基因的过高表达可能意味着细胞受到污染或者样本处理不当。

添加线粒体、核糖体、血红细胞信息

如果研究的是人或者小鼠,默认线粒体基因命名前有 MT-,可以用此来识别。 其他物种可以通过 features 参数指定线粒体基因列表。

seob[["percent.mt"]] <- PercentageFeatureSet(
  seob, 
  pattern = "^MT-"
)

seob[["percent.ribo"]] <- PercentageFeatureSet(
  seob, 
  pattern = "^RPS|^RPL"
)


seob[["percent.hb"]] <- PercentageFeatureSet(
  seob, 
  pattern = "^HB"
)

可视化

p0 <- VlnPlot(seob, 
              features = c("nFeature_RNA""nCount_RNA""percent.mt","percent.ribo","percent.hb"), 
              assay"RNA",
              layer"counts",
              group.by  = "samples",
              log = T,
              pt.size = 0)
p0

图片

过滤细胞

固定参数过滤

subset 命令用于对 seurat 对象取子集,具体用法参考 satijalab.org/seurat/arti…

seob_filter1 <- subset(seob,
                       subset = nFeature_RNA > 200 &
                         nFeature_RNA < 1500 &
                         percent.mt < 7 &
                         percent.ribo >3 &
                         percent.hb < 1 )

p1 <- VlnPlot(seob_filter1, 
              features = c("nFeature_RNA""nCount_RNA""percent.mt","percent.ribo","percent.hb"), 
              assay"RNA",
              layer"counts",
              group.by  = "samples",
              log = F,
              pt.size = 0)

图片

动态参数过滤

library(tidyverse)

# 计算中位数和 MAD
cells_filtered <- seob@meta.data %>%
  rownames_to_column(var"cell_id") %>%
  group_by(samples) %>%
  mutate(
    nCount_median = median(nCount_RNA),
    nCount_mad = mad(nCount_RNA),
    nFeature_median = median(nFeature_RNA),
    nFeature_mad = mad(nFeature_RNA),
    percent_mt_median = median(percent.mt),
    percent_mt_mad = mad(percent.mt),
    percent_ribo_median = median(percent.ribo),
    percent_ribo_mad = mad(percent.ribo),
    percent_hb_median = median(percent.hb),
    percent_hb_mad = mad(percent.hb)
  ) %>%
  filter(
    nCount_RNA >= nCount_median - 3*nCount_mad,
    nCount_RNA <= nCount_median + 3*nCount_mad,
    nFeature_RNA >= nFeature_median - 3*nFeature_mad,
    nFeature_RNA <= nFeature_median + 3*nFeature_mad,
    percent.mt >= percent_mt_median - 3*percent_mt_mad,  
    percent.mt <= percent_mt_median + 3*percent_mt_mad,  
    percent.ribo >= percent_ribo_median - 3*percent_ribo_mad,  
    percent.ribo <= percent_ribo_median + 3*percent_ribo_mad,  
    percent.hb >= percent_hb_median - 3*percent_hb_mad,  
    percent.hb <= percent_hb_median + 3*percent_hb_mad  
  ) %>%
  pull(cell_id)


seob_filter2 <- subset(seob, cells = cells_filtered)

p2 <- VlnPlot(seob_filter2, 
              features = c("nFeature_RNA""nCount_RNA""percent.mt","percent.ribo","percent.hb"), 
              assay"RNA",
              layer"counts",
              group.by  = "samples",
              log = F,
              pt.size = 0)

p2