scRNA-seq marker identification(一)

分享是一种态度

回顾

学习目标

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

Single-cell RNA-seq marker identification

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

目标

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

挑战

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

建议

• 将结果视为需要验证的假设。过高的p值可能会导致对结果的过度解释(基本上每个细胞都有重复)。排名靠前的标记是值得信赖的。确定每个聚类条件之间所有的保守标记。
• 识别在特定群集之间差异表达的标记

我们的聚类分析产生了以下群集：

我们在聚类分析中有以下问题：

1. 群集7和20的细胞类型是什么？
2. 与相同细胞类型相对应的群集是否有生物学意义的差异？
3. 这些细胞类型有亚群吗？我们是否可以通过识别这些群集的其他标记基因来获得对这些细胞类型特性的更准确的识别？

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

1. 识别每个群集的所有标记：该分析将每个群集与所有其他群集进行比较，并输出差异表达/存在的基因。

• 对于识别未知群集和提高假设细胞类型的置信度非常有用。

识别每个群集的保守标记：该分析首先寻找在每个条件下差异表达/存在的基因，然后报告那些在所有条件下都在群集中保守的基因。这些基因可以帮助鉴定群集的身份。

• 与多个条件配合使用时，可用于标识跨条件保留的细胞类型标记。

特定聚类之间的标记识别：该项分析探索了特定簇之间差异表达的基因。

• 对于从上述分析中确定的似乎代表相同细胞类型(即具有相似标记)的群集之间的基因表达差异很有用。

识别每个群集的所有标记

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

注：默认值为Wilcoxon Rank Sum检验，但也有其他可用选项。

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

• logfc.threshold ：相对于所有其他群集组合中的平均表达，群集中基因的平均表达的最小log2倍数变化。默认值为0.25。
• 如果平均log2FC不满足阈值，则可能会错过那些在目标簇中的一小部分细胞中表达但不在其他簇中表达的细胞标记
• 由于不同细胞类型的代谢输出存在细微差异，可能会返回许多代谢/核糖体基因，这对于区分细胞类型身份没有帮助
• 缺点：
• 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值，然后使用MetaDE R软件包中的meta分析方法进行跨组组合。

在开始标记鉴定之前，我们将明确设置默认测定，我们希望使用原始计数，而不是集成数据。

DefaultAssay(seurat_integrated) <- "RNA"

注意：虽然此函数的默认设置是从“RNA”插槽获取数据，但我们建议您运行上面的代码行，以绝对确保万一活动插槽在您的分析中的上游某处发生更改。原始计数和归一化计数存储在此槽中，用于查找标记的函数将自动提取原始计数。

函数 FindConservedMarkers() 具有以下结构：FindConservedMarkers() 语法:

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的log2倍数变化的阈值。我们还将指定仅返回每个群集的正标记。

在一个群集上进行测试，看看它是如何工作的：

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

FindConservedMarkers() 函数的输出是一个矩阵，其中包含按我们指定的群集的基因ID列出的假定标记的排序列表，以及相关的统计数据。请注意，为每个组（在我们的 Case,Ctrl和Stim）计算相同的统计信息集，最后两列对应于这两个组中的组合p值。

• gene:  gene symbol
• condition_p_val:  未针对条件进行多次测试校正调整的P值
• condition_avg_logFC: 条件的平均log2FC。正值表示该基因在簇中的表达量更高。
• condition_pct.1: 在条件群集中检测到该基因的细胞百分比
• condition_pct.2: 在条件的其他群集中平均检测到该基因的细胞百分比
• condition_p_val_adj:  基于使用数据集中所有基因的Bonferroni校正。调整后的p值，用于确定显著性
• max_pval:  按group/condition计算的p值的最大p值
• minimump_p_val:   combined p value

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

在查看输出时，我们建议寻找在 pct.1 和 pct.2 之间表达差异较大且logFC较大的标记。例如，如果 pct.1 =0.90和 pct.2 =0.80，则可能没有那么令人兴奋的标记。然而，如果 pct.2 =0.1，那么更大的差异将会很有说服力。同样，我们感兴趣的是表达该标记的大多数细胞是否在我感兴趣的群集中。如果 pct.1 较低，例如0.3，则可能没有那么有意义。如上所述，这两个参数也是运行函数时可能包括的参数。如上所述，这两个参数也是运行函数时可能包括的参数。

添加基因注释

添加带有基因注释信息的列可能会很有帮助。为此，将此文件（https://github.com/hbctraining/scRNA-seq/raw/master/data/annotation.csv）下载到您的数据文件夹。然后将其加载到R环境中：

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

注意：如果您有兴趣了解我们是如何获得此注释文件的，请查看链接(https://hbctraining.github.io/scRNA-seq/lessons/fetching_annotations.html)。

首先，我们将带有基因标识符的行名转换为自己的列。然后，我们将此注释文件与来自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)

---

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

• 生信爆款入门-第9期（线上直播4周，马拉松式陪伴，带你入门）你的生物信息入门课
  

• 数据挖掘第7期（线上直播3周，马拉松式陪伴，带你入门） 医学生/医生首选技能提高课

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

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

长按扫码可关注