单细胞Seurat的umi矩阵-与feature、counts(用于质控)

时间:2024-12-10 21:39:43

目录

关于umi矩阵学习

用umi计算feature、counts值

①meta数据查看

②Count和Feature计算(生成Seurat时自动计算)

1)提取UMI矩阵

2)计算

其他指标

评估质量指标(重点)

1)UMI计数

2)基因计数

3)UMIs vs. genes detected

4)线粒体计数比率

5)综合过滤

过滤提取子集

关于umi矩阵学习

10X数据因为有UMI,不需要考虑基因长度的影响,但仍然需要考虑测序深度在不同细胞之间的差异,所以需要用函数LogNormalize进行处理,具体方法是用该基因的UMI/该细胞的全部UMI,再乘以10000,在按列LogNormalize后,又可以按行进行scale​​,以去除极大值极小值基因对数据的影响。关于单细胞TPM、Count数据的处理_rds文件的umicount和readcount是什么-****博客

用umi计算feature、counts值

scRNA-seq质量控制流程scRNA-seq—质量控制 ()

①meta数据查看

③单细胞学习-pbmc的Seurat 流程_seurat 删除离散细胞-****博客

Seurat会自动为每个细胞创建一些元数据 data <- pbmc@

  1. rm(list=ls())
  2. library(dplyr)
  3. library(Seurat)
  4. library(patchwork)
  5. ## 读取数据
  6. pbmc.data <- Read10X(data.dir = "F:/##24年单细胞处理##/pbmc3k_filtered_gene_bc_matrices")
  7. ## 创建Seruat对象
  8. pbmc <- CreateSeuratObject(counts = pbmc.data,
  9. project = "pbmc3k",
  10. = 3, # :每个feature至少在多少个细胞中表达(feature=gene)
  11. = 200) # :每个细胞中至少有多少个feature被检测到
  12. pbmc
  13. data <- pbmc@meta.data
> head(data)
                  nCount_RNA nFeature_RNA
AAACATACAACCAC-1     pbmc3k       2419          779
AAACATTGAGCTAC-1     pbmc3k       4903         1352
AAACATTGATCAGC-1     pbmc3k       3147         1129
AAACCGTGCTTCCG-1     pbmc3k       2639          960
AAACCGTGTATGCG-1     pbmc3k        980          521
AAACGCACTGGTAC-1     pbmc3k       2163          781

  •  :通常包含样本标识,通常默认project为我们为其分配的身份

  • nCount_RNA :每个细胞的UMI数量

  • nFeature_RNA :每个细胞检测到的基因数量(非0数量)


②Count和Feature计算(生成Seurat时自动计算
1)提取UMI矩阵
  1. #提取UMI矩阵
  2. exp <- GetAssayData(pbmc, slot = "counts", assay = "RNA")
  3. umi_df <- data.frame(exp)
  4. umi_df[1:5,1:3]
umi_df[1:5,1:3]
              AAACATACAACCAC.1 AAACATTGAGCTAC.1 AAACATTGATCAGC.1
AL627309.1                   0                0                0
AP006222.2                   0                0                0
RP11-206L10.2                0                0                0
RP11-206L10.9                0                0                0
LINC00115                    0                0                0
2)计算
  1. dat <- cbind(colnames(umi_df),
  2. colSums(umi_df), ##每列UMI的和
  3. colSums(umi_df != 0))#每列检测到的基因数量feature
  4. dat1 <- apply(dat[,c(2,3)],2,as.numeric) #转化为数值型
  5. rownames(dat1) <- (dat)
  6. colnames(dat1) <- c("UMI1","UMI2")
  7. dat1 <- as.data.frame(dat1)
> head(dat1)
                 UMI1 UMI2
AAACATACAACCAC.1 2419  779
AAACATTGAGCTAC.1 4903 1352
AAACATTGATCAGC.1 3147 1129
AAACCGTGCTTCCG.1 2639  960
AAACCGTGTATGCG.1  980  521
AAACGCACTGGTAC.1 2163  781

其他指标
  • number of genes detected per UMI: 这个度量让我们对数据集的复杂性有了一个概念(每个UMI检测到的基因越多,我们的数据就越复杂)

  • mitochondrial ratio: 这个度量将给我们提供一个源自线粒体基因的细胞读数的百分比。

每个细胞的每个UMI的基因数量非常容易计算,我们将对结果进行log10转换,以便更好地在样本之间进行比较。

  1. # Add number of genes per UMI for each cell to metadata
  2. pbmc$log10GenesPerUMI <- log10(pbmc$nFeature_RNA) / log10(pbmc$nCount_RNA)

PercentageFeatureSet()将采用一定模型并搜索基因标识符。此功能可以轻松计算属于每个细胞的可能功能子集的所有计数的百分比。

  1. # Compute percent mito ratio
  2. pbmc$mitoRatio <- PercentageFeatureSet(object = pbmc, pattern = "^MT-")
  3. pbmc$mitoRatio <- pbmc@meta.data$mitoRatio / 100

注意:提供的模式(“ ^ MT-”)适用于人类基因名称。

新的原数据查看

  1. data1 <- pbmc@meta.data
  2. head(data1)
> head(data1)
                  nCount_RNA nFeature_RNA log10GenesPerUMI
AAACATACAACCAC-1     pbmc3k       2419          779        0.8545652
AAACATTGAGCTAC-1     pbmc3k       4903         1352        0.8483970
AAACATTGATCAGC-1     pbmc3k       3147         1129        0.8727227
AAACCGTGCTTCCG-1     pbmc3k       2639          960        0.8716423
AAACCGTGTATGCG-1     pbmc3k        980          521        0.9082689
AAACGCACTGGTAC-1     pbmc3k       2163          781        0.8673469

评估质量指标(重点)
  • 细胞计数

  • 每个细胞的UMI计数

  • 每个细胞检测到的基因

  • UMI与检测到的基因

  • 线粒体比率

  • Novelty

scRNA-seq—质量控制 ()

③单细胞学习-pbmc的Seurat 流程_seurat 删除离散细胞-****博客

1)UMI计数

nCount_RNA :每个细胞的UMI数量:这里没有做相关质控

可视化

  1. # Visualize the number UMIs/transcripts per cell
  2. library(ggplot2)
  3. data1 %>%
  4. ggplot(aes(color='', x= nCount_RNA, fill= '')) +
  5. geom_density(alpha = 0.2) +
  6. scale_x_log10() +
  7. theme_classic() +
  8. ylab("Cell density") +
  9. geom_vline(xintercept = 500)

2)基因计数

nFeature_RNA :每个细胞检测到的基因数量

  1. data1 %>%
  2. ggplot(aes(color='', x= nFeature_RNA, fill= '')) +
  3. geom_density(alpha = 0.2) +
  4. scale_x_log10() +
  5. theme_classic() +
  6. ylab("Cell density") +
  7. geom_vline(xintercept = 500)

3)UMIs vs. genes detected

通常一起评估的两个指标是UMI的数量和每个细胞检测到的基因数量。在这里,我们绘制了基因数量与线粒体reads所占比例UMI数量的关系图。线粒体reads部分只在检测到很少基因的特别低计数的细胞中才高(浅蓝色)。这可能是损伤/死亡的细胞,其细胞质的mRNA已经通过破裂的膜泄漏出来,因此,只有位于线粒体的mRNA仍然是保守的。这些细胞被我们的计数和基因数量阈值过滤掉。联合可视化计数和基因阈值可显示联合过滤效果

质量差的细胞很可能每个细胞的基因和UMI都很低,并且与图左下象限的数据点相对应。好的细胞通常会表现为每个细胞有更多的基因和更高数量的UMI。

通过此图,我们还评估了线的斜率,以及图的右下角象限中数据点的任何散布情况。这些细胞有大量的UMI,但只有几个基因。这些可能是濒临死亡的细胞,但也可能代表一个低复杂性细胞类型的群体(即红细胞)。

  1. # Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
  2. p <- data1 %>%
  3. ggplot(aes(x=nCount_RNA, y=nFeature_RNA, color=mitoRatio)) +
  4. geom_point() +
  5. scale_colour_gradient(low = "gray90", high = "black") +
  6. stat_smooth(method=lm) +
  7. scale_x_log10() +
  8. scale_y_log10() +
  9. theme_classic() +
  10. geom_vline(xintercept = 500) +
  11. geom_hline(yintercept = 250) +
  12. facet_wrap(~)
  13. p

4)线粒体计数比率

这一指标可以确定死亡或濒临死亡的细胞是否存在大量的线粒体污染。我们将线粒体计数质量差的样品定义为超过0.2线粒体比率标记的细胞,除非您希望样品中有这种情况。

  1. # Visualize the distribution of mitochondrial gene expression detected per cell
  2. p1 <- data1 %>%
  3. ggplot(aes(color=, x=mitoRatio, fill=)) +
  4. geom_density(alpha = 0.2) +
  5. scale_x_log10() +
  6. theme_classic() +
  7. geom_vline(xintercept = 0.2)
  8. p1

5)综合过滤

我们可以看到,我们对每个细胞测序较少的样本具有更高的整体复杂性,这是因为我们还没有开始对这些样本的任何给定基因进行饱和测序。这些样本中的异常值细胞可能是RNA种类比其他细胞简单的细胞。有时,我们可以通过这个指标检测到低复杂性细胞类型的污染,比如红细胞。一般来说,我们预计novelty得分在0.80以上。

  1. # Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI
  2. p3 <- data1 %>%
  3. ggplot(aes(x=log10GenesPerUMI, color = , fill=)) +
  4. geom_density(alpha = 0.2) +
  5. theme_classic() +
  6. geom_vline(xintercept = 0.8)
  7. p3

过滤提取子集
  1. #过滤
  2. pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & < 5)
  3. pbmc

scRNA-seq—质量控制 ()

[1] 代码可用于自行计算该指标: /hbctraining/scRNA-seq/blob/master/lessons/
[2] Scrublet: /AllonKleinLab/scrublet