EcoTyper代码实操（五）：在scRNA-seq数据中发现细胞状态和生态型

前言

EcoTyper是一个基于机器学习的工具，能够从Bulk、单细胞、以及空间分辨率的基因表达数据中大规模地识别并验证细胞状态和生态型。我们在前面的推文中介绍了EcoTyper的分析框架和部分实操，感兴趣的小伙伴可以先阅读这一部分哦。

EcoTyper的代码实操主要分为6个部分：

那么，今天让我们一起来学习一下它的第五个部分——De novo Discovery of Cell States and Ecotypes in scRNA-seq Data。

EcoTyper的分析框架：
Extract cell type specific or top variable genes：在分析前，通常需要去除不特定于某一细胞类型的基因（在给定细胞类型和所有其他细胞类型之间差异表达的基因作为特异基因），以减少发现假细胞状态的可能性。然后，通过Wilcox检验和Benjamini-Hochberg校正来筛除Q值大于0.05的基因。
Cell state discovery on correlation matrices：EcoTyper 使用非负矩阵分解（NMF）来从单细胞表达数据中鉴定细胞状态。然而，由于scRNA-seq数据通常是稀疏的，直接应用NMF可能效果不佳。因此，EcoTyper 在给定细胞类型的每对细胞之间的相关矩阵上应用NMF。
Choosing the number of cell states：基于cophenetic系数确定细胞状态。cophenetic系数量化了给定细胞状态数量（或聚类数）的分类稳定性，通常在0到1之间，1表示最稳定。
Extracting cell state information：从步骤2的NMF输出中提取细胞状态信息。
Cell state re-discovery on expression matrices：在步骤2中确定了细胞状态后，EcoTyper进行差异表达分析，以鉴定与每个细胞状态高度相关的基因。然后，选取具有前1000个基因进行新一轮的NMF。
Extracting cell state information：解析第5步NMF结果，提取细胞状态信息用于后续分析。
Cell state QC filter：为了避免假细胞状态的识别，EcoTyper引入了自适应假阳指数（AFI）来过滤具有AFI >= 1的细胞状态。
Ecotype (cellular community) discovery：生态型或细胞社群是通过识别细胞状态在样本之间的共现模式来获得的。EcoTyper使用Jaccard系数来量化每对细胞状态在发现队列中样本之间的重叠程度。然后，根据Jaccard矩阵，通过超几何检验评估是否存在细胞状态之间的重叠。最后，通过层次聚类和轮廓宽度来确定最佳的生态型数量。

代码流程

1.准备环境和数据

下载 EcoTyper

wget https://github.com/digitalcytometry/ecotyper/archive/refs/heads/master.zip
unzip master.zip
cd ecotyper-master
#EcoTyper是一个独立的软件，用R实现，但并不是R包哦。

R环境准备：

这是官网推荐的R和R包版本，但是小编用的最新的R版本，并且安装相应最新版本的R包也是很顺利的。

install.packages(c("RColorBrewer","cluster", "circlize","cowplot","data.table","doParallel","ggplot2","grid", "reshape2", "viridis","config","argparse","colorspace", "plyr"))
BiocManager::install("ComplexHeatmap")
BiocManager::install("Biobase")
BiocManager::install("NMF")

输入数据：

1. 表达谱数据：结直肠癌scRNA-seq数据，数据存储在example_data/scRNA_CRC_data.txt中。

2. 样本注释文件，数据存储在example_data/scRNA_CRC_annotation.txt。

2.discovery scripts

EcoTyper_discovery_scRNA.R脚本用于在单细胞转录组数据中恢复细胞状态和生态型。

Rscript EcoTyper_discovery_scRNA.R -h

此脚本接受 YAML 格式的配置文件作为输入文件（config_discovery_scRNA.yml）。

参数详解：

• Discovery dataset name : "discovery_scRNA_CRC"

  发现数据集名称

• Expression matrix : "example_data/scRNA_CRC_data.txt"

  表达数据的文件路径

• Annotation file : "example_data/scRNA_CRC_annotation.txt"

  注释数据的文件路径

• Annotation file column to scale by : NULL

  指定一个注释文件中的列名，按照这个列的值，对基因进行缩放（均值为0，方差为1）

• Annotation file column(s) to plot : []

  指定注释文件中哪些列将用作输出热图中的颜色条（color bar）

  Annotation file column(s) to plot : ["Histology", "Tissue"]

• Output folder : "DiscoveryOutput_scRNA"

• Filter genes : "cell type specific"

  Accepted values:

  "cell type specific" - select genes overexpressed in a cell type

  "no filter" - use all genes

• Number of threads : 10

  线程数

• Number of NMF restarts : 5

  NMF重新启动次数参数。EcoTyper 使用NMF（非负矩阵分解）来鉴定细胞状态，但NMF的结果可能取决于随机初始化。为了获得稳定的解决方案，NMF通常多次运行，每次使用不同的种子。

• Maximum number of states per cell type : 20

  每种细胞类型可能的状态数量的上限。

• Cophenetic coefficient cutoff : 0.95

  取值范围为[0,1]，较低的值会得到较多的细胞状态数量。

• Jaccard matrix p-value cutoff : 1

  在步骤8中，生态型的识别是通过对衡量每对状态之间的样本重叠的Jaccard矩阵进行聚类来完成的。在执行生态型识别之前，将不显著的样本重叠的Jaccard矩阵值设置为0，以减轻由虚假重叠引入的噪音。值得注意的是，在具有较小样本数量的scRNA-seq数据集中，建议禁用此筛选以避免过度筛选Jaccard矩阵，但如果发现的scRNA-seq数据集包含足够多的样本来可靠评估重叠的显著性，可以将此截断值设置为较低的值（例如0.05）。

• Minimum number of states in ecotypes : 3

  少于3种细胞状态的生态型被过滤掉

运行脚本：

Rscript EcoTyper_discovery_scRNA.R -c config_discovery_scRNA.yml

3.结果解读

输出的结果包括3类，rank_plot图、各种细胞类型的细胞状态以及生态型。我们来逐各学习。

1. rank_plot（Cophenetic系数图）:横向虚线表示配置文件中提供的Cophenetic系数截断值（Cophenetic coefficient cutoff field）。红色的垂直虚线表示基于提供的Cophenetic系数截断值自动选择的状态数量。注意这些图表示在应用步骤6和7中的低质量状态筛选之前获得的状态数量。因此最终结果可能包含更少的状态。

接下来，我们以内皮细胞为例展示ecotyper识别的细胞状态结果

2. state_assignment.txt：每个样本被分配到具有最高推断丰度的细胞状态的信息，也就是每个样本的主要细胞状态。

3. state_assignment_heatmap.png/pdf:展示在每个细胞状态中具有最高 log2 FC的基因的表达情况，这些基因被称为细胞状态特异性标记基因。列代表发现数据集中的样本，而行代表了每个细胞状态的标记基因。

4. "heatmap_top_ann.txt" 是用户提供的注释文件与细胞状态相关的信息的整合。只有被分配到主要细胞状态的样本才会包含在这个文件中。

生态型部分输出的结果

5. "ecotype_abundance.txt" 文件包含了在发现数据集中的每个生态型（ecotype）相对丰度的信息。

6. "ecotype_assignment.txt" 文件包含了将发现数据集中的样本分配到各个生态型的信息。

7. "heatmap_assigned_samples_viridis.pdf/png" 展示了分配了生态型的样本中细胞状态分数的热图。

8. nclusters_jaccard.png/pdf: "initial clusters" 指的是通过对Jaccard指数矩阵进行聚类获得的初始簇的数量。这个数量通常是通过应用平均轮廓（average silhouette）方法来选择的。

9. jaccard_matrix.png/pdf:包含至少3个细胞状态的生态型的Jaccard系数矩阵的热图。

---

小结

在本期推文中，我们介绍了如何使用EcoTyper从单细胞转录组数据中从头发现细胞状态和生态型。小编认为，我们在研究的过程中可以整合较多的单细胞数据集保证更为全面地发现细胞状态和生态型。随后，将模型应用于Bulk转录组数据中恢复细胞状态和生态型，结合Bulk数据样本量大、表型数据多的优势更好的探索相关生物学结果。我们将在下期推文中和大家学习EcoTyper的最后一部分实操--De novo Discovery of Cell States and Ecotypes in Pre-Sorted Data，关于今天的学习内容，感兴趣的小伙伴快来用起来吧~

好啦，本期的分享到这里就结束了，我们下期再会~