学习目标：

• 了解如何确定单个亚群的标记
• 理解聚类和标记识别的迭代过程

单细胞RNA-seq标记鉴定

现在我们已经确定了我们想要的集群，我们可以继续进行标记识别，这将使我们能够验证某些集群的身份，并帮助推测任何未知集群的身份。

image

目标：

• 为了确定每个集群的基因标记
• 使用标记来识别每个集群的细胞类型
• 为了确定是否需要基于细胞类型标记重新聚类，可能需要合并或拆分聚类

挑战：

• 对结果过度解读
• 结合不同类型的标记进行识别

建议：

• 把结果看作是需要验证的假设。夸大的p值可能导致结果的过度解释(基本上每个细胞都被用作一个复制)。得分高的值可信度高。识别每个聚类的所有保存在条件之间的标记
• 识别特定簇之间表达有差异的标记

通过聚类分析产生了以下聚类：

image

请记住，我们在聚类分析中遇到了以下问题：

1. 集群 7 和 20 的细胞类型标识是什么？
2. 对应于相同细胞类型的簇是否具有生物学意义的差异？是否存在这些细胞类型的亚群中？
3. 通过识别这些簇的其他标记基因，我们能否对这些细胞类型身份获得更高的可信度？

我们可以使用 Seurat 探索几种不同类型的识别标记，以获得这些问题的答案。每个都有自己的优点和缺点：

1. 对每个集群的所有标记进行识别：该分析将每个集群与其他所有集群进行比较，并输出差异表达/存在的基因。
   • 用于识别未知簇并提高对假设细胞类型的可信度。
2. 识别每个簇的保守标记：该分析首先寻找在每个条件下差异表达/存在的基因，然后找出所有条件下在簇中保守的基因。这些基因可以帮助确定集群的身份。
   • 用于多个条件以识别不同条件下保守的细胞类型标记。
3. 特定簇之间的标记识别：该分析探讨了特定簇之间差异表达的基因。
   • 用于确定从上述分析中似乎代表相同细胞类型（即具有相似标记）的簇之间的基因表达差异。

识别每个簇的所有标记

这种类型的分析通常推荐用于评估单个样品组/条件。使用 FindAllMarkers() 函数，我们将每个集群与所有其他集群进行比较，以识别潜在的标记基因。每个集群中的细胞被视为重复，本质上通过一些统计检验进行差异表达分析。

注意：默认是 Wilcoxon Rank Sum检验，但也可以选择其他检验方法。

image

该FindAllMarkers()函数有三个重要参数，它们为确定基因是否为标记提供阈值：

• logfc.threshold：相对于所有其他聚类组合中基因平均表达的最小log2倍变化。默认值为 0.25。
  • 缺点：
    • 如果平均logfc不满足阈值，可能会错过那些在感兴趣的集群内的一小部分细胞中表达的细胞标记，而在其他集群中不会
    • 由于不同细胞类型的代谢输出略有差异，可能会返回大量代谢/核糖体基因，这对于区分细胞类型身份没有什么用
• min.diff.pct：表达该基因的细胞百分数与所有其他基因组合中表达基因的细胞百分数之间的最小差异百分数。
  • 缺点：可能会错过那些在所有细胞中表达的细胞标记物，但在这个特定的细胞类型中高度上调
• min.pct：只检测在两个种群中最小部分细胞中检测到的基因。这意味着通过不测试那些很少表达的基因来加快功能。默认是0.1。
  • 缺点：由于不是所有的基因都能在所有的细胞中被检测到(即使是表达出来的)，如果设置一个非常高的值可能会导致许多假阴性

您可以使用这些参数的任意组合，这取决于您想要的严格/宽松程度。此外，默认情况下，这个函数将返回阳性和阴性表达变化的基因，添加参数only.pos来选择只保留阳性表达的基因。为每个集群查找标记的代码如下所示。

## DO NOT RUN THIS CODE ##

# Find markers for every cluster compared to all remaining cells, report only the positive ones
markers <- FindAllMarkers(object = seurat_integrated, 
                          only.pos = TRUE,
                          logfc.threshold = 0.25)                     

注意：此命令可能需要很长时间才能运行，因为它正在针对其他所有细胞处理每个单独的集群。

在所有条件下鉴定保守标记

由于我们的数据集中有不同条件下的样本，因此我们最好的选择是找到保守标记。此函数在内部按样本组/条件分离出细胞，然后针对单个指定簇与所有其他簇（或第二个簇，如果指定）执行差异基因表达测试。计算每个条件下基因的 p 值，然后使用R 包中的 MetaDE 荟萃分析方法跨组组合。

image

在我们开始我们的标记识别之前，我们将明确设置我们的默认检测，我们想要使用标准化数据，而不是整合数据。

DefaultAssay(seurat_integrated) <- "RNA"

注意：默认检测应该已经是RNA，因为我们在之前的聚类质量控制课程中设置了它。但我们鼓励您运行上面的这行代码，以确保万一运行环境在您的分析上游某处发生了更改。请注意，RNA原始计数和归一化计数存储在检测的counts和data插槽。默认情况下，查找标记的函数将使用标准化数据。

函数FindConservedMarkers(), 具有以下结构：

FindConservedMarkers() 句法：

## DO NOT RUN ##

FindConservedMarkers(seurat_integrated,
                     ident.1 = cluster,
                     grouping.var = "sample",
                     only.pos = TRUE,
             min.diff.pct = 0.25,
                     min.pct = 0.25,
             logfc.threshold = 0.25)

您将看到FindAllMarkers()函数描述的一些参数；这是因为它在内部使用该函数首先在每个组中查找标记。在这里，我们列出了一些额外的参数，这些参数在使用FindConservedMarkers()时提供：

• ident.1：这个函数一次只评估一个集群；在这里，您需要指定感兴趣的集群。
• grouping.var：元数据中的变量（列标题），用于指定将细胞分成组

对于我们的分析，只使用大于 0.25 的对数倍数变化阈值。返回每个集群的阳性标记。

让我们在一个集群上测试它，看看它是如何工作的：

cluster0_conserved_markers <- FindConservedMarkers(seurat_integrated,
                              ident.1 = 0,
                              grouping.var = "sample",
                              only.pos = TRUE,
                      logfc.threshold = 0.25)

image

FindConservedMarkers()函数的输出是一个矩阵，其中包含按我们指定的聚类的基因 ID 列出的推定标记的排名列表，以及相关的统计数据。请注意，为每个组（在我们的例子中，Ctrl 和 Stim）计算相同的统计数据集，最后两列对应于两个组的组合 p 值。我们在下面描述了其中的一些列：

• 基因：基因符号
• condition_p_val：未针对条件的多重测试校正调整 p 值
• condition_avg_logFC：条件的平均对数折叠变化。正值表示该基因在簇中的表达更高。
• condition_pct.1：在条件簇中检测到基因的细胞百分比
• condition_pct.2：在其他簇中平均检测到基因的细胞百分比
• condition_p_val_adj：条件的调整 p 值，基于使用数据集中所有基因的 bonferroni 校正，用于确定显著性
• max_pval： 每个组/条件计算的p值的最大p值
• minimump_p_val：组合 p 值

注意：由于每个细胞都被视为一个重复，这将导致每个组内的 p 值膨胀！一个基因可能具有较低 的p 值 < 1e-50，但这并不能转化为高度可靠的标记基因。

在查看输出时，我们建议在pct.1和pct.2之间寻找表达差异很大的标记，以及更大的折叠变化。例如，如果pct.1= 0.90 和pct.2= 0.80，它可能不会像标记一样令人兴奋。但是，如果pct.2改为 = 0.1，则更大的差异会更令人信服。此外，感兴趣的是表达标记的大多数细胞是否在我感兴趣的簇中。如果pct.1低，例如 0.3，则可能不那么有趣。如上所述，这两个也是在运行函数时可能包含的参数。

添加基因注释

添加基因注释列信息会很有帮助。为此，我们将data使用下面提供的代码加载位于您文件夹中的注释文件：

annotations <- read.csv("data/annotation.csv")

注意：如果您有兴趣了解我们如何获得此注释文件，请查看链接的材料。

首先，我们将带有基因标识符的行名变成它自己的列。然后我们将这个注释文件与我们的结果合并FindConservedMarkers()：

# Combine markers with gene descriptions 
cluster0_ann_markers <- cluster0_conserved_markers %>% 
                rownames_to_column(var="gene") %>% 
                left_join(y = unique(annotations[, c("gene_name", "description")]),
                          by = c("gene" = "gene_name"))

View(cluster0_ann_markers)

练习

在上一课中，我们通过检查已知细胞标记物 FCGR3A 和 MS4A7 的表达，将簇 9 确定为 FCGR3A+ 单核细胞。使用FindConservedMarkers()函数找到簇 9 的保守标记。你观察到什么？您是否将 FCGR3A 和 MS4A7 视为簇 9 中的高表达基因？

在多样本上运行

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

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

# 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 家族语法：

## DO NOT RUN ##
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().

评估标记基因

我们想使用这些基因列表来查看我们可以识别这些簇与哪些细胞类型相关。让我们看看每个集群的顶部的基因，看看这是否给我们任何提示。我们可以通过两个组的平均倍数变化查看前 10 个标记，对于每个集群进行快速阅读：

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

# Visualize top 10 markers per cluster
View(top10)

image

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

对于第20 组，富集的基因似乎没有一个让我们感兴趣。我们经常查看pct.1与pct.2良好标记基因相比差异较大的基因。例如，我们可能对基因 TPSB2 感兴趣，该基因显示簇中表达该基因的细胞比例很大，但其他簇中表达该基因的细胞很少。如果我们“谷歌”TPSB2，我们会找到GeneCards 网站。

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

因此，簇20可能代表肥大细胞。肥大细胞是免疫系统的重要细胞，属于造血谱系。研究已经确定肥大细胞特征显着富含丝氨酸蛋白酶，例如TPSAB1 和 TPSB2，这两种蛋白酶都出现在我们的保守标记列表中。另一个不是丝氨酸蛋白酶但是已知肥大细胞特异性基因并出现在我们列表中的基因是 FCER1A（编码 IgE 受体的一个亚基）。此外，我们看到 GATA1 和 GATA2 出现在我们的列表中，它们不是肥大细胞标记基因，但在肥大细胞中大量表达，是已知的转录因子调节各种肥大细胞特异性基因。

可视化标记基因

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

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

image

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

小提琴图类似于箱线图，不同之处在于它们还显示了数据在不同值下的概率密度，通常由核密度评估器平滑。小提琴图比普通箱线图提供更多信息。箱线图仅显示汇总统计数据，例如均值/中位数和四分位距，而小提琴图则显示了数据的完整分布。当数据分布是多峰（多于一个峰）时，这种差异特别有用。在这种情况下，小提琴图显示了不同峰值的存在、它们的位置和相对幅度。

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

image

这些结果和图可以帮助我们确定这些簇的身份，或验证我们之前探索过预期细胞类型的规范标记后假设的身份。

识别每个簇的基因标记

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

image

我们可以尝试所有比较的组合，但我们将从簇 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)

image

在这些顶级基因中，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, 
            order = TRUE,
            min.cutoff = 'q10',
        repel = TRUE
            )

image

由于初始状态和激活状态的标记都显示在标记列表中，因此可视化表达是有帮助的。基于这些图，集群 0 和 2 似乎是可靠的初始 T 细胞。然而，对于激活的 T 细胞，很难说。我们可能会说簇 4 和 18 是激活的 T 细胞，但 CD69 表达不如 CREM 明显。我们将标记初始细胞并将剩余的簇标记为 CD4+ T 细胞。

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

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

# 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)

image

如果我们想删除潜在的压力细胞，我们可以使用以下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)

image

现在我们想要保存我们最终标记的 Seurat 对象和输出sessionInfo()：

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

# Create and save a text file with sessionInfo
sink("sessionInfo_scrnaseq_Oct2020.txt")
sessionInfo()
sink()

您可以在此链接中找到有关该sink()功能的更多信息。

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

• 通过实验验证我们确定的细胞类型的有趣标记。
• 探索细胞类型的子集以发现此处描述的细胞子簇
• 进行条件ctrl和stim条件之间的差异表达分析
  • 进行此分析需要生物学重复，我们还有其他材料来帮助完成此分析。
• 如果试图确定细胞类型或细胞状态之间的进展，则可以进行轨迹分析或谱系追踪。例如，我们可以使用这种类型的分析来探索以下任何一项：
  • 分化过程
  • 表达随时间变化
  • 不同表达下细胞状态的变化