Seurat对象的QC和filter

brief

• seurat提供的教学里面包含了Standard pre-processing workflow,workflow包括QC，normalization，scale data ，detection of highly variable features。
• 其中 normalization就有蛮多方法的，seurat自己就提供了两种，一种是"LogNormalization"，另一种是目前正在大力推荐的"SCTransform",而且每一种方法参数众多，我觉得有必要对其进行仔细的探究，所以需要分开学习和记录。
• 还有 detection of highly variable features以及scale data，也是方法和参数众多，我觉得有必要对其进行仔细的探究，所以需要分开学习和记录
• 概要图以及系列博文可以参见链接。
  

这里主要记录了使用 UMI，gene feature，percentage of mitochondrial expression进行细胞过滤和提取的操作，以及一些seurat提供的可视化方法。

• 其中seurat提供的思路如下：
  

实例

• 构建seurat对象

library(dplyr)
library(Seurat)
library(patchwork)

rm(list=ls())

# 使用read10X读取output of the cellranger pipeline from 10X，包括barcodes，genes，matrix.mtx三个文件
pbmc.data <- Read10X(data.dir = "D:/djs/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19")
# 使用 CreateSeuratObject函数构造seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200,
                           names.delim = "-",names.field = 1)

• QC信息的构建和过滤

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
# 计算 a percentage of cell reads originating from the mitochondrial genes
# 这里利用的是正则表达式去map 人的线粒体基因，然后进行计算比例，其他物种的话线粒体基因是不是以MT开头的我不太清楚，所以可能需要手动计算了

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

# 手动计算部分percentage of cell reads originating from the mitochondrial genes可以参考下面的链接
# https://github.com/hbctraining/scRNA-seq/blob/master/lessons/mitoRatio.md


# 计算 complexity of the RNA species
pbmc@meta.data$log10GenesPerUMI <- log10(pbmc$nFeature_RNA) / log10(pbmc$nCount_RNA)

# 下面一步把pre-processing workflow中的QC阶段走完了，
# 其中线粒体含量最好根据组织特性设定，比如心肌细胞线粒体基因含量超过20%很正常啊
# 这里的过滤参数设置规定每个细胞包含的gene应该在200-2500个，线粒体基因含量占比小于5%
# 在第一步构建seurat 对象的时候，就把那些少于3个细胞表达的gene以及少于200个gene的cell丢弃了
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

• 可视化探索部分

# Visualize QC metrics as a violin plot
# 这里 VlnPlot()函数的参数很多，可视化的内容绝不仅仅是下面这个图形
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

# Visualize the number UMIs/transcripts per cell
pbmc@meta.data %>% 
  ggplot(aes(color=orig.ident, x=nCount_RNA, fill= orig.ident)) + 
  geom_density(alpha = 0.2) + 
  scale_x_log10() + 
  theme_classic() +
  ylab("Cell density") +
  geom_vline(xintercept = 500)


# Visualize the distribution of genes detected per cell via histogram
pbmc@meta.data %>% 
  ggplot(aes(color=orig.ident, x=nFeature_RNA, fill= orig.ident)) + 
  geom_density(alpha = 0.2) + 
  theme_classic() +
  scale_x_log10() + 
  geom_vline(xintercept = 300)

# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI (novelty score)
pbmc@meta.data %>%
  	ggplot(aes(x=log10GenesPerUMI, color=orig.ident, fill= orig.ident)) +
  	geom_density(alpha = 0.2) +
  	theme_classic() +
  	geom_vline(xintercept = 0.8)

# Visualize the distribution of mitochondrial gene expression detected per cell
pbmc@meta.data %>% 
  	ggplot(aes(x=percent.mt, color=orig.ident, fill= orig.ident)) + 
  	geom_density(alpha = 0.2) + 
  	scale_x_log10() + 
  	theme_classic() +
  	geom_vline(xintercept = 0.2)