本文来源于哈佛大学的单细胞课程系列，在此做一些学习，不当之处请指正。

scRNA-seq/09_merged_SC_marker_identification.md at master · hbctraining/scRNA-seq · GitHubhttps://github.com/hbctraining/scRNA-seq/blob/master/lessons/09_merged_SC_marker_identification.md

学习目标：

• 了解如何确定单个细胞群的标记

• 了解聚类和标记识别的迭代过程

单细胞RNA-seq标记鉴定

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

目标：

• 确定用于每个集群的基因标记

• 使用标记每个集群的鉴定细胞类型

• 为了确定是否需要基于细胞类型标记进行重新聚类，可能需要合并或拆分聚类

挑战：

• 对结果的过度解释

• 结合不同类型的标记识别

建议：

• 将结果视为需要验证的假设。夸大的p值可能导致对结果的过度解释（基本上每个细胞都用作重复项）。最佳标记是最值得信赖, 识别每个簇的条件之间保守的所有标记。

• 识别在特定簇之间差异表达的标记

我们的聚类分析得出以下聚类：

请记住，从聚类分析中我们有以下问题：

1. 集群7和20的类型标识是什么？

2. 对应于相同细胞类型的簇是否具有生物学上有意义的差异？这些细胞类型有亚群吗？

3. 通过为这些簇鉴定其他标记基因，我们能否对这些细胞类型的身份获得更高的信心？

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

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

• 对于识别未知簇和提高对假定细胞类型的置信度很有用。

2. 鉴定每个簇的保守标记：该分析首先寻找在每种条件下差异表达/存在的基因，然后报告在所有条件下在簇中保守的那些基因。这些基因可以帮助弄清楚簇的身份。

• 在不止一种条件下有用，可识别在各种条件下保守的细胞类型标记。

3. 特定簇之间的标记鉴定：该分析探索了特定簇之间差异表达的基因。

• 根据上述分析，可用于确定似乎代表相同细胞类型（即标记相似）的簇之间基因表达的差异。

识别每个聚类的所有标记

在评估单个样本组/条件时，通常建议使用这种类型的分析。通过该FindAllMarkers()功能，我们将每个聚类与所有其他聚类进行比较，以识别潜在的标记基因。将每个簇中的细胞视为重复样品，并通过一些统计学检验进行差异表达分析。

注意：默认是Wilcoxon秩和检验，但是还有其他选项可用。

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

• logfc.threshold：相对于所有其他集群组合中的平均表达，集群中基因的平均表达的最小log2倍变化。默认值为0.25。

• 如果平均log2FC不满足阈值，则可能会错过那些在目标簇中的一小部分细胞中表达但不在其他簇中表达的细胞标记

• 由于不同细胞类型的代谢输出存在细微差异，可能会返回许多代谢/核糖体基因，这对于区分细胞类型身份没有帮助

• 缺点：

• min.diff.pct：簇中表达基因的细胞百分比与所有其他簇中表达基因的细胞百分比之和之间的最小百分比差异。

• 缺点：可能会错过在所有细胞中表达但在该特定细胞类型中高度上调的那些细胞标志物

• min.pct：仅在两个种群中的任何一个中至少有一部分细胞中检测到的测试基因。旨在通过不测试很少表达的基因来加快功能。默认值为0.1。

• 缺点：如果设置为很高的值，可能会导致许多假阴性，原因是并非在所有细胞中都检测到了所有基因（即使它被表达了）

您可以根据需要的严格程度使用这些参数的任意组合。同样，默认情况下，此功能将返回显示正向和负向表达变化的基因。通常，我们添加一个参数only.pos以选择仅保留积极的更改。查找每个集群标记的代码如下所示。

# 查找每个群集的显著基因
markers <- FindAllMarkers(object = seurat_integrated, 
                          only.pos = TRUE,
                          logfc.threshold = 0.25)  

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

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

由于我们的数据集中有代表不同条件的样本，因此我们最好的选择是找到保守标记。此功能在内部按样本组/条件分离出细胞，然后针对所有其他聚类（或第二个聚类（如果指定））针对单个指定聚类执行差异基因表达测试。计算每种条件的基因水平p值，然后使用MetaDE R软件包中的荟萃分析方法跨组进行组合。在开始标记识别之前，我们将明确设置默认测定，我们要使用原始计数而不是积分数据。

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列出的推定标记的排名列表，以及相关的统计信息。请注意，为每个组（在我们的示例中为Ctrl和Stim）计算一组相同的统计信息，最后两列对应于两个组中的组合p值。我们在下面介绍其中的一些列：

• 基因：基因符号

• condition_p_val： p值未针对条件的多次测试校正进行调整

• condition_avg_logFC：条件的平均log2倍变化, 正值表明该基因在簇中表达更高。

• condition_pct.1：在条件中的簇中检测到该基因的细胞百分比

• condition_pct.2：在其他集群中平均检测到该基因的细胞所占的百分比

• condition_p_val_adj：基于使用数据集中所有基因的邦费罗尼校正，条件调整后p值，用于确定显著性

• 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文件夹中。然后将其加载到您的R环境中：

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

注意：如果您想知道我们如何获得此批注文件，请查看链接的资料。

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

# 将标记与基因描述结合起来
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)

在多个样本上运行

该函数一次FindConservedMarkers() 只接受一个集群，我们可以在有集群的情况下运行该函数多次。但是，这不是很有效。相反，我们将首先创建一个函数来查找保守标记，包括我们要包括的所有参数。我们还将添加几行代码来修改输出。我们的职能将：

1. 运行FindConservedMarkers()功能

2. 使用rownames_to_column()功能将行名称传输到列

3. 合并注释

4. 使用cbind()函数创建集群ID列

# 创建功能，可以为任何给定的集群保守标记
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**** 家庭语法：**

map_dfr(inputs_to_function, name_of_function)

现在，让我们尝试使用此功能为未识别的细胞类型的集群找到保守标记：cluster7和cluster 20。

#簇迭代函数
conserved_markers <- map_dfr(c(7,20), get_conserved)

查找所有集群的标记

对于您的数据，您可能想在所有群集上运行此功能，在这种情况下，您可以输入0:20而不是c(7,20); 但是，要花很长时间才能运行。同样，有可能在所有群集上运行此功能时，在某些情况下，群集中的细胞不足以容纳特定组-功能将失败。对于这些群集，您将需要使用FindAllMarkers()。

评估标记基因

我们希望使用这些基因列表来确定这些簇可以识别哪些细胞类型。让我们看一下每个簇的顶级基因，看看是否能给我们任何提示。我们可以按两组的平均倍数变化查看前10个标记，以便快速浏览每个集群：

# 提取每个簇的前10个标记
top10 <- conserved_markers %>% 
         mutate(avg_fc = (ctrl_avg_logFC + stim_avg_logFC) /2) %>% 
         group_by(cluster_id) %>% 
         top_n(n = 10, 
         wt = avg_fc)

# 可视化每个群集的前10个标记
View(top10)

我们看到簇7出现了许多热休克和DNA损伤基因。基于这些标记，可能是这些细胞处于应激或垂死状态。但是，如果我们更详细地研究这些细胞的质量指标（即，覆盖在集群上的mitoRatio和nUMI），我们实际上看不到支持该参数的数据。如果我们仔细观察标记基因列表，我们还会发现一些与T细胞相关的基因和激活标记。这些可能是活化的（细胞毒性）T细胞。有广泛的研究支持热休克蛋白与反应性T细胞在慢性炎症中诱导抗炎细胞因子的结合。在这个簇中，我们需要对免疫细胞有更深入的了解，才能真正弄清结果并得出最终结论。对于聚类20，富集基因似乎没有一个共同的细胞类型，我们经常会查看pct.1与pct.2好标记基因相比，差异较大的基因。例如，我们可能对基因TPSB2感兴趣，该基因显示了表达该基因的簇中有很大比例的细胞，而其他簇中却很少表达该基因的细胞。如果我们使用Google 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"),
                         sort.cell = TRUE,
                         min.cutoff = 'q10', 
                         label = TRUE,
			 repel = TRUE)

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

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

＃ VLN情节-集群20 

# cluster 20韦恩图
VlnPlot(object = seurat_integrated, 
        features = c("TPSAB1", "TPSB2", "FCER1A", "GATA1", "GATA2"))

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

识别每个簇的基因标记

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

# 确定分化为CD4 + T细胞标记
cd4_tcells <- FindMarkers(seurat_integrated,
                          ident.1 = 2,
                          ident.2 = c(0,4,10,18))                  

# 将基因符号添加到DE表
cd4_tcells <- cd4_tcells %>%
  rownames_to_column(var = "gene") %>%
  left_join(y = unique(annotations[, c("gene_name", "description")]),
             by = c("gene" = "gene_name"))

# 重新排序列和排序padj   
cd4_tcells <- cd4_tcells[, c(1, 3:5,2,6:7)]

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

#查看数据
View(cd4_tcells)

在这些最重要的基因中，CREM基因是激活标记。我们知道激活的另一个标记是CD69，而幼稚或记忆细胞的标记包括SELL和CCR7基因。有趣的是，SELL基因也位于列表的顶部。让我们使用这些新的细胞状态标记在可视化上稍微探索一下激活状态：

# 绘制激活的和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
            )

由于Naive和Activated 状态的标记都显示在标记列表中，因此可视化表达很有帮助。根据这些图，似乎簇0和2确实是Naive 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)

如果我们要删除Stressed cells / Activated T cells细胞，可以使用以下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

• 要进行此分析，必须进行生物学复制，并且我们还有其他材料可以帮助您完成此分析。

• 如果试图确定细胞类型或细胞状态之间的进程，可以执行轨迹分析或谱系追踪。例如，使用这种类型的分析，我们可以探索以下任何一种：

• 分化过程

• 表达随时间变化

• 表达中的细胞状态变化

本课程由哈佛大学生物信息学核心（HBC）的教学团队成员开发。这些是根据知识共享署名许可（CC BY 4.0）的条款分发的开放获取材料，只要注明原始作者和出处，就可以在任何介质中进行不受限制的使用，分发和复制。

• 这些材料和动手操作的一部分改编自Satija Lab的Seurat-_ [引导聚类教程_](https://satijalab.org/seurat/pbmc3k_tutorial.html)

##############################################
####09_merged_SC_marker_identification.md
## 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)  
####新版本wt = avg_logFC，改为wt = avg_log2FC
markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

top5 <- markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC)

DefaultAssay(seurat_integrated) <- "RNA"

DoHeatmap(seurat_integrated, pt.size=1, features = top5$gene) + NoLegend()


DotPlot(seurat_integrated, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
top3<- markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_log2FC)
DotPlot(seurat_integrated, features =unique(top3$gene),cols=c('#deebf7', "#6baed6"))+
  theme_bw()+
theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5))


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)
??FindConservedMarkers()

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

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

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

# 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 family syntax:

map_dfr(inputs_to_function, name_of_function)


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

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


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

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


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


# 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
)
# 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)+scale_color_npg()

# 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)+
        theme_bw()+
        scale_color_npg()


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