单细胞数据挖掘2-QC

1. 写在前面

• 参考教程：https://hbctraining.github.io/scRNA-seq/

2. 搭建环境

• 在 Rproj 的位置创建 data 文件夹
• 在 data 文件夹中为每个样本创建一个链接到 outs/raw_feature_bc_matrix 文件夹的软链接

3. 读入数据

• 可选 readMM() 或 Read10X()，后者比较方便
• 使用 for 循环读入
  • 代码

  for (file in c("ctrl_raw_feature_bc_matrix", "stim_raw_feature_bc_matrix")){
          seurat_data <- Read10X(data.dir = paste0("data/", file))
          seurat_obj <- CreateSeuratObject(counts = seurat_data,
                                           min.features = 100,
                                           project = file)
          assign(file, seurat_obj)
  }
  

  • 最后的 assign 这个技巧很实用，使得每次循环创造的 Seurat 对象不会重复写入同一个变量
• 然后把这些对象 merge 起来方便下一步操作

# Create a merged Seurat object
merged_seurat <- merge(x = sigaf1,
                       y = c(sigag1,sigah1),
                       add.cell.id = c("sigaf1","sigag1","sigah1"))

4. 质控

• 在质控之前，首先要对自己的细胞有一定的了解，主要是是否属于低复杂度细胞以及是否属于高线粒体表达的细胞，比如我挖掘的数据集的 LK 细胞就不属于这两类细胞

4.1. 添加 metadata

• 首先需要添加这两项
  • number of genes detected per UMI: this metric with give us an idea of the complexity of our dataset (more genes detected per UMI, more complex our data)
  • mitochondrial ratio: this metric will give us a percentage of cell reads originating from the mitochondrial genes
  • 分别用于评估细胞的复杂度和线粒体表达情况
• 代码

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

# Compute percent mito ratio
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^mt-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100

• 人是 MT-开头，鼠是 mt- 开头
• 接下来添加细胞名和样本名信息，并将 Seurat 对象储存为 Rdata

# 为了不改动 Seurat 对象，我们把 metadata 单独取出来添加信息
# Create metadata dataframe
metadata <- merged_seurat@meta.data
# Add cell IDs to metadata
metadata$cells <- rownames(metadata)

# Rename columns
metadata <- metadata %>%
  dplyr::rename(seq_folder = orig.ident,
                nUMI = nCount_RNA,
                nGene = nFeature_RNA)

# Create sample column
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^sigaf1_"))] <- "sigaf1"
metadata$sample[which(str_detect(metadata$cells, "^sigag1_"))] <- "sigag1"
metadata$sample[which(str_detect(metadata$cells, "^sigah1_"))] <- "sigah1"
metadata = metadata[,c(1,2,3,5,6,7,4)]

# Add metadata back to Seurat object
merged_seurat@meta.data <- metadata

# Create .RData object to load at any time
save(merged_seurat, file="data/merged_filtered_seurat.RData")

4.2. 针对双细胞

• Currently, we recommend not including any thresholds at this point in time. When we have identified markers for each of the clusters, we suggest exploring the markers to determine whether the markers apply to more than one cell type.

4.3. 细胞计数

# Visualize the number of cell counts per sample
metadata %>%
    ggplot(aes(x=sample, fill=sample)) +
    geom_bar() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    theme(plot.title = element_text(hjust=0.5, face="bold")) +
    ggtitle("NCells")

4.4. UMI counts (transcripts) per cell

# Visualize the number UMIs/transcripts per cell
metadata %>%
    ggplot(aes(color=sample, x=nUMI, fill= sample)) +
    geom_density(alpha = 0.2) +
    scale_x_log10() +
    theme_classic() +
    ylab("Cell density") +
    geom_vline(xintercept = 500)

• 最少 500，500-1000 勉强能用，最好大于 1000

  
  

4.5. Genes detected per cell

# Visualize the distribution of genes detected per cell via histogram
metadata %>%
    ggplot(aes(color=sample, x=nGene, fill= sample)) +
    geom_density(alpha = 0.2) +
    theme_classic() +
    scale_x_log10() +
    geom_vline(xintercept = 300)

# Visualize the distribution of genes detected per cell via boxplot
metadata %>%
    ggplot(aes(x=sample, y=log10(nGene), fill=sample)) +
    geom_boxplot() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    theme(plot.title = element_text(hjust=0.5, face="bold")) +
    ggtitle("NCells vs NGenes")

• 高质量的数据应是单峰，且峰的位置与细胞类型相匹配

• 如果出现大峰边的小峰或双峰现象，则既可能是技术上有部分细胞检测失败，也可能是真实的生物学差异

  
  
  

4.6. UMIs vs. genes detected

# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
metadata %>%
  ggplot(aes(x=nUMI, y=nGene, color=mitoRatio)) +
  geom_point() +
  scale_colour_gradient(low = "gray90", high = "black") +
  stat_smooth(method=lm) +
  scale_x_log10() +
  scale_y_log10() +
  theme_classic() +
  geom_vline(xintercept = 500) +
  geom_hline(yintercept = 300) +
  facet_wrap(~sample)

• 线粒体基因比例高，且基因数和 UMI count 都低的细胞为低质量细胞。

  
  

4.7. Mitochondrial counts ratio

# Visualize the distribution of mitochondrial gene expression detected per cell
metadata %>%
  ggplot(aes(color=sample, x=mitoRatio, fill=sample)) +
  geom_density(alpha = 0.2) +
  scale_x_log10() +
  theme_classic() +
  geom_vline(xintercept = 0.2)

4.8. Complexity

# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI
metadata %>%
  ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
  geom_density(alpha = 0.2) +
  theme_classic() +
  geom_vline(xintercept = 0.8)

4.9. 细胞水平过滤

# Filter out low quality reads using selected thresholds - these will change with experiment
filtered_seurat <- subset(x = merged_seurat,
                          subset= (nUMI >= 500) &
                            (nGene >= 300) &
                            (log10GenesPerUMI > 0.80) &
                            (mitoRatio < 0.20))

4.10. 基因水平过滤

• 仅保留在10个及以上细胞中表达的基因

# Output a logical vector for every gene on whether the more than zero counts per cell
# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")

# Output a logical vector for every gene on whether the more than zero counts per cell
nonzero <- counts > 0

# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10

# Only keeping those genes expressed in more than 10 cells
filtered_counts <- counts[keep_genes, ]

# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)

4.11. Re-assess QC metrics

• 重跑一遍前面的可视化以评估质控效果

友情宣传

• 全国巡讲全球听（买一得五），第二期 ，你的生物信息学入门课
• 生信技能树的2019年终总结 ，你的生物信息学成长宝藏
• 2020学习主旋律，B站74小时免费教学视频为你领路