scRNA-seq marker identification(二)

分享是一种态度

回顾

运行多个样本

函数FindConservedMarkers()一次接受一个群集，我们可以运行这个函数，有多少群集就运行多少次。然而，这并不是很有效。相反，我们将首先创建一个函数来查找保守的标记，包括我们想要包含的所有参数。我们还将添加几行代码来修改输出。我们的函数将会：

1. 运行 FindConservedMarkers() 函数
2. 使用 rowames_to_column() 函数将行名传输到列
3. 合并注释
4. 使用 cbind() 函数创建群集IDs列

# Create function to get conserved markers for any given cluster
get_conserved <- function(cluster){
  FindConservedMarkers(seurat_integrated,
                       ident.1 = cluster,
                       grouping.var = "sample",
                       only.pos = TRUE) %>%
    rownames_to_column(var = "gene") %>%
    left_join(y = unique(annotations[, c("gene_name", "description")]),
               by = c("gene" = "gene_name")) %>%
    cbind(cluster_id = cluster, .)
  }

现在我们已经创建了这个函数，我们可以将其用作适当的 map 函数的参数。我们希望map系列函数的输出是一个数据框，我们将使用map_dfr()函数将每个集群输出通过行合并在一起。

map family syntax:

map_dfr(inputs_to_function, name_of_function)

现在，让我们尝试用此函数来查找未识别细胞类型的群集（群集7和群集20）的保守标记。

# Iterate function across desired clusters
conserved_markers <- map_dfr(c(7,20), get_conserved)

查找所有群集的标记

对于您的数据，您可能希望在所有群集上运行此函数，在这种情况下，您可以输入 0：20 而不是 c(7，20) ；但是，运行该函数需要相当长的时间。此外，当您在所有群集上运行此函数时，在某些情况下，您的群集可能没有足够的细胞用于特定的组-因此您的函数将失败。对于这些群集，您可以使用 FindAllMarkers() 。

评估标记基因

我们想要使用这些基因列表来看看我们能不能识别出这些细胞群与哪些细胞类型一致。我们可以通过查看每个群集的top基因，看看这是否能我们一些提示。我们可以按两组的平均fold change来查看前10个标记，对于每个群集，我们可以快速浏览：

# Extract top 10 markers per cluster
top10 <- conserved_markers %>% 
  mutate(avg_fc = (ctrl_avg_logFC + stim_avg_logFC) /2) %>% 
  group_by(cluster_id) %>% 
  top_n(n = 10, 
        wt = avg_fc)

# Visualize top 10 markers per cluster
View(top10)

我们看到第7簇出现了许多热休克和DNA损伤相关基因。根据这些标记，这些很可能是应激或濒临死亡的细胞。然而，如果我们更详细地研究这些细胞的质量指标(即群集的mitoRatio和nUMI)，我们并没有真正看到支持这一论点的数据。如果我们仔细看一看标记基因列表，我们还会发现一些与T细胞相关的基因和标记的激活。这些可能是激活的(细胞毒性)T细胞。有广泛的研究支持热休克蛋白与反应性T细胞在慢性炎症中诱导抗炎细胞因子的联系。在这个群集中，我们需要对免疫细胞有更深入的了解，才能真正梳理结果并得出最终结论。

对于群集20，富集基因似乎没有一个突出的共同主题。为了寻找好的标记基因，我们经常查看 pct.1 与 pct.2 差异较大的基因。例如，我们可能对TPSB2基因感兴趣，它显示该群集中有很大比例的细胞表达该基因，但其他群集中表达该基因的细胞很少。如果我们用谷歌搜索TPSB2，就会找到GeneCard网站（https://www.genecards.org/cgi-bin/carddisp.pl?gene=TPSB2&keywords=TPSB2）。

“β-类胰蛋白酶似乎是肥大细胞表达的主要同工酶，而在嗜碱性细胞中，α-类胰蛋白酶占主导地位。类胰蛋白酶被认为是哮喘和其他过敏性和炎症性疾病的致病介质。“

因此，群集20可能表示肥大细胞。肥大细胞是免疫系统的重要细胞，属于造血系统。研究（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5045264/）已经发现肥大细胞的特征明显富含丝氨酸蛋白酶，如TPSAB1和TPSB2，这两种蛋白酶都出现在我们保守的标志物列表中。另一个基因不是丝氨酸蛋白酶，而是已知的肥大细胞特异性基因，出现在我们的列表中是FCER1A(编码IgE受体的一个亚单位)。此外，我们发现GATA1和GATA2出现在我们的列表中，它们不是肥大细胞标记基因，但在肥大细胞中大量表达，是已知的调节各种肥大细胞特异性基因的转录因子。（https://mcb.asm.org/content/34/10/1812）

可视化标记基因

为了更好地了解群集20的细胞类型标识，我们可以使用 FeaturePlot() 函数按群集探索不同标识标记的表达。

# Plot interesting marker gene expression for cluster 20
FeaturePlot(object = seurat_integrated, 
                        features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"),
                         sort.cell = TRUE,
                         min.cutoff = 'q10', 
                         label = TRUE,
    repel = TRUE)

我们还可以通过使用小提琴图来探索特定标记的表达范围：

小提琴图类似于箱线图，不同之处在于它们还显示不同值的数据的概率密度，通常由kernel density estimator进行平滑处理。小提琴图比普通的箱线图更能提供信息。箱线图仅显示平均值/中位数和四分位数范围等汇总统计数据，而小提琴图则显示数据的完整分布。当数据分布是多模式(多个峰值)时，这种差异特别有用。在这种情况下，小提琴曲线图显示了不同峰值的存在、它们的位置和相对振幅。

# Vln plot - cluster 20
VlnPlot(object = seurat_integrated, 
        features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"))

这些结果和曲线图可以帮助我们确定这些群的身份，或者在之前探索了预期细胞类型的规范标记后，验证是否是我们假设的身份。

识别每个群集的基因标记

关于分析的最后一组问题涉及到同一细胞类型相对应的群集是否具有生物学意义上的差异。有时返回的标记列表不能充分分隔某些群集。例如，我们之前已经确定第0、2、4、10和18簇为CD4+T细胞，但是这些细胞群集之间是否存在生物学上的差异？我们可以使用 FindMarkers() 函数来确定两个特定群集之间差异表达的基因。

我们可以尝试所有组合的比较，但我们将从群集2与所有其他CD4+T细胞群集开始：

# Determine differentiating markers for CD4+ T cell
cd4_tcells <- FindMarkers(seurat_integrated,
                          ident.1 = 2,
                          ident.2 = c(0,4,10,18))                  

# Add gene symbols to the DE table
cd4_tcells <- cd4_tcells %>%
  rownames_to_column(var = "gene") %>%
  left_join(y = unique(annotations[, c("gene_name", "description")]),
             by = c("gene" = "gene_name"))

# Reorder columns and sort by padj      
cd4_tcells <- cd4_tcells[, c(1, 3:5,2,6:7)]

cd4_tcells <- cd4_tcells %>%
  dplyr::arrange(p_val_adj) 

# View data
View(cd4_tcells)

在这些top基因中，CREM基因作为激活的标志脱颖而出。我们知道另一个激活的标志是CD69，而原始细胞或记忆细胞的标志包括SELL和CCR7基因。有趣的是，SELL基因也位于列表顶部。让我们使用这些新的细胞状态标记直观地了解一下激活状态：

# Plot gene markers of activated and naive/memory T cells
FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("CREM", "CD69", "CCR7", "SELL"),
            label = TRUE, 
            sort.cell = TRUE,
            min.cutoff = 'q10',
     repel = TRUE
            )

由于原始状态和激活状态的标记都出现在标记列表中，因此有助于可视化表达。根据这些图，似乎群集0和2确实是原始T细胞。然而，对于被激活的T细胞来说，就很难说了。我们可以说群集4和群集18是激活的T细胞，但CD69的表达不像CREM那样明显。我们将标记原始细胞，并将剩余的群集标记为CD4+T细胞。

现在利用所有这些信息，我们可以推测不同簇的细胞类型，并用细胞类型标签绘制细胞图。

Cluster ID	Cell Type
0	Naive or memory CD4+ T cells
1	CD14+ monocytes
2	Naive or memory CD4+ T cells
3	CD14+ monocytes
4	CD4+ T cells
5	CD8+ T cells
6	B cells
7	Stressed cells / Activated T cells
8	NK cells
9	FCGR3A+ monocytes
10	CD4+ T cells
11	B cells
12	NK cells
13	CD8+ T cells
14	CD14+ monocytes
15	Conventional dendritic cells
16	Megakaryocytes
17	B cells
18	CD4+ T cells
19	Plasmacytoid dendritic cells
20	Mast cells

然后，我们可以将群集的标识重新分配给这些细胞类型：

# Rename all identities
seurat_integrated <- RenameIdents(object = seurat_integrated, 
                               "0" = "Naive or memory CD4+ T cells",
                               "1" = "CD14+ monocytes",
                               "2" = "Naive or memory CD4+ T cells",
                               "3" = "CD14+ monocytes",
                               "4" = "CD4+ T cells",
                               "5" = "CD8+ T cells",
                               "6" = "B cells",
                               "7" = "Stressed cells / Activated T cells",
                               "8" = "NK cells",
                               "9" = "FCGR3A+ monocytes",
                               "10" = "CD4+ T cells",
                               "11" = "B cells",
                               "12" = "NK cells",
                               "13" = "CD8+ T cells",
                               "14" = "CD14+ monocytes",
                               "15" = "Conventional dendritic cells",
          "16" = "Megakaryocytes",
          "17" = "B cells", 
          "18" = "CD4+ T cells", 
          "19" = "Plasmacytoid dendritic cells", 
          "20" = "Mast cells")


# Plot the UMAP
DimPlot(object = seurat_integrated, 
        reduction = "umap", 
        label = TRUE,
        label.size = 3,
        repel = TRUE)

如果我们想要移除应激或死亡细胞，可以使用 subset() 函数：

# Remove the stressed or dying cells
seurat_subset_labeled <- subset(seurat_integrated,
                               idents = "Stressed cells / Activated T cells", invert = TRUE)

# Re-visualize the clusters
DimPlot(object = seurat_subset_labeled, 
        reduction = "umap", 
        label = TRUE,
        label.size = 3,
 repel = TRUE)

现在，我们保存最终标记的Seurat对象：

# Save final R object
write_rds(seurat_integrated,
          path = "results/seurat_labelled.rds")       

---

现在我们已经定义了群集，并且每个群集的标记都已定义，我们有以下几个不同的选择：

• 通过实验验证我们识别的细胞类型的有趣标记。
• 进行 ctrl 和 stim 条件之间的差异表达分析
• 生物重复是进行这项分析所必需的，我们有额外的资料（https://hbctraining.github.io/scRNA-seq/lessons/pseudobulk_DESeq2_scrnaseq.html）来帮助我们完成这项分析。
• 如果试图确定细胞类型或细胞状态之间的发育，可以进行轨迹分析或谱系追踪。例如，我们可以使用此类型的分析来探索以下任何一项：
• 分化过程
• 随时间的表达变化
• 表达中的细胞状态的变化

---

注：以上内容来自哈佛大学生物信息中心(HBC)_的教学团队的生物信息学培训课程。原文链接：https://hbctraining.github.io/scRNA-seq/schedule/  点击 “阅读原文” 可直达

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注