15.1 前言
除了基因表达模式的变化之外，细胞组成（例如细胞类型的比例）也会在不同条件下发生变化。例如，特定药物可以诱导细胞类型的转分化，这将反映在细胞身份组合物中。需要足够的细胞和样本数量才能准确确定细胞同一性簇比例和背景变异。可以在已知细胞类型或对应于例如最近受扰动影响的细胞的细胞状态的形式的细胞身份簇的水平上进行组成分析。
差异丰度分析比较两种条件下细胞类型的组成。两种模式的样品含有不同比例的细胞类型，可以测试其丰度的显着变化。
本篇将介绍这两种方法并将其应用于Haber数据集[Haber et al., 2017]。该数据集包含来自小鼠小肠和类器官的53193个单个上皮细胞。一些细胞还受到细菌或蠕虫感染，例如分别通过沙门氏菌和Heligmosomoides polygyrus感染。在本教程中，我们使用完整Haber数据集的子集，其中仅包括专门为此目的收集的对照细胞和受感染细胞。值得注意的是，我们排除了仅收集大量细胞的附加数据集，以加快计算速度并降低复杂性。
第一步，我们加载数据集。
15.2 加载数据
import warnings
import pandas as pd
warnings.filterwarnings("ignore")
warnings.simplefilter("ignore")
import scanpy as sc
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import altair as alt
import pertpy as pt
adata = pt.dt.haber_2017_regions()
adata
AnnData object with n_obs × n_vars = 9842 × 15215
obs: 'batch', 'barcode', 'condition', 'cell_label'
adata.obs
batch barcode condition cell_label
index
B1_AAACATACCACAAC_Control_Enterocyte.Progenitor B1 AAACATACCACAAC Control Enterocyte.Progenitor
B1_AAACGCACGAGGAC_Control_Stem B1 AAACGCACGAGGAC Control Stem
B1_AAACGCACTAGCCA_Control_Stem B1 AAACGCACTAGCCA Control Stem
B1_AAACGCACTGTCCC_Control_Stem B1 AAACGCACTGTCCC Control Stem
B1_AAACTTGACCACCT_Control_Enterocyte.Progenitor B1 AAACTTGACCACCT Control Enterocyte.Progenitor
... ... ... ... ...
B10_TTTCACGACAAGCT_Salmonella_TA B10 TTTCACGACAAGCT Salmonella TA
B10_TTTCAGTGAGGCGA_Salmonella_Enterocyte B10 TTTCAGTGAGGCGA Salmonella Enterocyte
B10_TTTCAGTGCGACAT_Salmonella_Stem B10 TTTCAGTGCGACAT Salmonella Stem
B10_TTTCAGTGTGACCA_Salmonella_Endocrine B10 TTTCAGTGTGACCA Salmonella Endocrine
B10_TTTCAGTGTTCTCA_Salmonella_Enterocyte.Progenitor B10 TTTCAGTGTTCTCA Salmonella Enterocyte.Progenitor
9842 rows × 4 columns
数据分10批收集。独特的条件是对照、沙门氏菌、Hpoly.Day3和Hpoly.Day10，分别对应于健康对照状态、沙门氏菌感染、Heligmosomoides polygyrus感染3天后的细胞和Heligmosomoides polygyrus感染10天后的细胞。cell_label 对应于细胞类型。
15.3 为什么细胞类型计数数据是组合的
在分析细胞计数数据的成分变化时，需要考虑多种技术和方法学限制。一项挑战是实验重复次数较少，这会导致在使用频率统计检验进行差异丰度分析时出现较大的置信区间。更重要的是，单细胞测序自然受到每个样本细胞数量的限制——我们无法对组织或器官中的每个细胞进行测序，而是使用一个小的、有代表性的快照来代替。然而，这迫使我们将细胞类型计数视为纯粹的比例，即样本中的细胞总数只是一个比例因子。在统计文献中，此类数据被称为成分数据，其特征是一个样本中所有特征（在我们的例子中为细胞类型）的相对丰度总和为1。
由于这种总和为一的限制，导致细胞类型丰度之间存在负相关。为了说明这一点，让我们考虑以下示例：
在病例对照研究中，我们想要比较健康器官和患病器官的细胞类型组成。在这两种情况下，我们都有三种细胞类型（A、B 和 C），但它们的丰度不同：
健康器官由每种类型2000个细胞组成（总共6000个细胞）。
这种疾病导致A型细胞倍增，而B型和C型细胞不受影响，因此患病器官有8000个细胞。
healthy_tissue = [2000, 2000, 2000]
diseased_tissue = [4000, 2000, 2000]
example_data_global = pd.DataFrame(
data=np.array([healthy_tissue, diseased_tissue]),
index=[1, 2],
columns=["A", "B", "C"],
)
example_data_global["Disease status"] = ["Healthy", "Diseased"]
example_data_global
A B C Disease status
1 2000 2000 2000 Healthy
2 4000 2000 2000 Diseased
plot_data_global = example_data_global.melt(
"Disease status", ["A", "B", "C"], "Cell type", "count"
)
fig, ax = plt.subplots(1, 2, figsize=(12, 6))
sns.barplot(
data=plot_data_global, x="Disease status", y="count", hue="Cell type", ax=ax[0]
)
ax[0].set_title("Global abundances, by status")
sns.barplot(
data=plot_data_global, x="Cell type", y="count", hue="Disease status", ax=ax[1]
)
ax[1].set_title("Global abundances, by cell type")
plt.show()
我们想找出患病器官中哪些细胞类型的丰度增加或减少。如果我们能够确定两个器官中每个细胞的类型，情况就会很清楚，如上右图所示。不幸的是，这是不可能的。由于我们的测序过程容量有限，我们只能从两个群体中抽取600个细胞的代表性样本。为了模拟此步骤，我们可以使用numpy的random.multinomial函数从群体中采样600个细胞，无需放回：
np.random.seed(1234)
healthy_sample = np.random.multinomial(
pvals=healthy_tissue / np.sum(healthy_tissue), n=600
)
diseased_sample = np.random.multinomial(
pvals=diseased_tissue / np.sum(diseased_tissue), n=600
)
example_data_sample = pd.DataFrame(
data=np.array([healthy_sample, diseased_sample]),
index=[1, 2],
columns=["A", "B", "C"],
)
example_data_sample["Disease status"] = ["Healthy", "Diseased"]
example_data_sample
A B C Disease status
1 193 201 206 Healthy
2 296 146 158 Diseased
plot_data_sample = example_data_sample.melt(
"Disease status", ["A", "B", "C"], "Cell type", "count"
)
fig, ax = plt.subplots(1, 2, figsize=(12, 6))
sns.barplot(
data=plot_data_sample, x="Disease status", y="count", hue="Cell type", ax=ax[0]
)
ax[0].set_title("Sampled abundances, by status")
sns.barplot(
data=plot_data_sample, x="Cell type", y="count", hue="Disease status", ax=ax[1]
)
ax[1].set_title("Sampled abundances, by cell type")
plt.show()
现在图片已经不太清楚了。虽然A型细胞的计数仍然增加（大约从200个到300个），但其他两种细胞类型似乎从大约200个减少到150个这种明显的减少是由于我们对600个细胞的限制造成的-如果样品被A型细胞占据，B型和C型细胞的比例必须较低。因此，如果不考虑其他细胞类型，就不可能确定一种细胞类型的丰度变化。
如果我们忽略数据的组合性，并使用像 Wilcoxon 秩和检验或 scDC 这样的单变量方法，我们可能会错误地将细胞类型群体变化视为统计上的声音效应，尽管它们是由细胞类型比例的固有负相关性引起的。
此外，二次采样数据不仅为我们的问题提供了一种有效的解决方案。如果患病病例中B型和C型细胞均减少1000个细胞，我们将获得与上述相同的600个细胞的代表性样本。为了获得唯一的结果，我们可以固定数据的参考点，假设该参考点在所有样本中都保持不变。这可以是单一细胞类型、多种细胞类型的聚合（例如几何平均值）。
scCODA已经解决了这个问题，我们将在下一节中介绍它并将其应用于我们的数据集。
15.4 带有标记的簇
scCODA属于需要预定义簇（最常见的细胞类型）来统计得出成分变化的工具系列。受微生物组数据成分分析方法的启发，scCODA提出了一种贝叶斯方法来解决单细胞分析中常见的低重复问题。它使用分层狄利克雷多项式模型对细胞类型计数进行建模，该模型通过对所有测量的细胞类型比例进行联合建模来解释细胞类型比例的不确定性和负相关偏差。为了确保唯一可识别的解决方案和易于解释，scCODA中的参考被选择为特定的细胞类型。因此，scCODA检测到的任何成分变化始终必须相对于所选参考进行查看。
然而，scCODA假设协变量和细胞丰度之间存在对数线性关系，这在使用连续协变量时可能并不总是反映潜在的生物过程。 scCODA的另一个限制是无法推断除成分效应之外的细胞成分之间的相关结构。此外，scCODA仅模拟平均丰度的变化，但无法检测响应变异性的变化。
第一步，我们实例化scCODA模型。
然后，我们使用load函数准备MuData对象以供后续处理，并根据输入数据创建成分分析数据集。在我们的例子中，我们指定cell_type_identifier为cell_label，sample_identifier为batch，covariate_obs为condition。
sccoda_model = pt.tl.Sccoda()
sccoda_data = sccoda_model.load(
adata,
type="cell_level",
generate_sample_level=True,
cell_type_identifier="cell_label",
sample_identifier="batch",
covariate_obs=["condition"],
)
sccoda_data
MuData object with n_obs × n_vars = 9852 × 15223
2 modalities
rna: 9842 x 15215
obs: 'batch', 'barcode', 'condition', 'cell_label', 'scCODA_sample_id'
coda: 10 x 8
obs: 'condition', 'batch'
var: 'n_cells'
为了概述不同条件下的细胞类型分布，我们可以使用scCODA的箱线图。为了更好地了解数据的分布方式，红点显示实际的数据点。
pt.pl.coda.boxplots(
sccoda_data,
modality_key="coda",
feature_name="condition",
figsize=(12, 5),
add_dots=True,
args_swarmplot={"palette": ["red"]},
)
plt.show()
箱线图突出显示了细胞类型分布的一些差异。明显值得注意的是沙门氏菌病症的肠上皮细胞比例很高。但其他细胞类型，例如转运扩增（TA）细胞，与对照相比，沙门氏菌条件下的丰度也表现出明显差异。必须正确评估这些差异是否具有统计显着性。
另一种可视化是scCODA提供的堆叠条形图。该可视化很好地显示了成分数据的特征：如果我们比较对照组和沙门氏菌组，我们可以看到受感染小鼠中肠上皮细胞的比例大大增加。由于数据是成比例的，这会导致满足总和为一约束的所有其他细胞类型的份额减少。
pt.pl.coda.stacked_barplot(
sccoda_data, modality_key="coda", feature_name="condition", figsize=(4, 2)
)
plt.show()
除了细胞计数AnnData对象之外，scCODA还需要两个主要参数：公式和引用细胞类型。该公式描述了使用R-style指定的协变量。在我们的例子中，我们将条件指定为唯一的协变量。由于它是具有四个水平（对照和三种疾病状态）的离散协变量，因此该模型对每个状态与其他样本的比较进行了建模。如果我们想一次对多个协变量进行建模，只需将它们添加到公式中（即
formula = "covariate_1 + covariate_2"）就足够了。如上所述，scCODA需要参考细胞类型进行比较，据信协变量不会改变该参考细胞类型。scCODA可以自动选择适当的细胞类型作为参考，该细胞类型在所有样品中具有几乎恒定的相对丰度，也可以使用用户指定的参考细胞类型运行。在这里，我们将内分泌细胞设置为参考，因为从视觉上看它们的丰度似乎相当恒定。手动设置参考单元类型的另一种方法是将reference_cell_type设置为“automatic”，这将强制scCODA本身选择合适的参考单元类型。如果参考细胞类型的选择不清楚，我们建议使用此选项来获取指标甚至最终选择。
sccoda_data = sccoda_model.prepare(
sccoda_data,
modality_key="coda",
formula="condition",
reference_cell_type="Endocrine",
)
sccoda_model.run_nuts(sccoda_data, modality_key="coda", rng_key=1234)
sample: 100%|██████████| 11000/11000 [01:08<00:00, 161.54it/s, 255 steps of size 1.72e-02. acc. prob=0.85]
sccoda_data["coda"].varm["effect_df_condition[T.Salmonella]"]
Final Parameter HDI 3% HDI 97% SD Inclusion probability Expected Sample log2-fold change
Cell Type
Endocrine 0.0000 0.000 0.000 0.000 0.0000 32.598994 -0.526812
Enterocyte 1.5458 0.985 2.071 0.283 0.9996 382.634978 1.703306
Enterocyte.Progenitor 0.0000 -0.475 0.566 0.143 0.2817 126.126003 -0.526812
Goblet 0.0000 -0.345 1.013 0.290 0.4354 52.735108 -0.526812
Stem 0.0000 -0.742 0.297 0.173 0.3092 135.406509 -0.526812
TA 0.0000 -0.876 0.331 0.211 0.3358 78.986854 -0.526812
TA.Early 0.0000 -0.338 0.615 0.151 0.3033 152.670412 -0.526812
Tuft 0.0000 -1.221 0.548 0.342 0.4098 23.041143 -0.526812
接受率描述了在初始老化阶段后接受的建议样本的比例，并且可以作为不良优化运行的临时指标。对于scCODA，所需的接受率在0.4到0.9之间。接受率过高或过低表明抽样过程存在问题。
sccoda_data
MuData object with n_obs × n_vars = 9852 × 15223
2 modalities
rna: 9842 x 15215
obs: 'batch', 'barcode', 'condition', 'cell_label', 'scCODA_sample_id'
coda: 10 x 8
obs: 'condition', 'batch'
var: 'n_cells'
uns: 'scCODA_params'
obsm: 'covariate_matrix', 'sample_counts'
varm: 'intercept_df', 'effect_df_condition[T.Hpoly.Day3]', 'effect_df_condition[T.Hpoly.Day10]', 'e
scCODA根据包含概率选择可信效应。可信效应和不可信效应之间的界限取决于所需的错误发现率 (FDR)。较小的FDR值会产生更保守的结果，但可能会遗漏一些效应，而较大的FDR值会选择更多的效应，但会产生更多的错误发现。推理后可以通过 sim_results.set_fdr() 轻松设置所需的FDR级别。默认情况下，该值为0.05。由于根据数据集的不同，FDR可能会对结果产生重大影响，因此我们建议尝试高达0.2的不同 FDR，以获得最显着的效果。
在我们的例子中，我们使用不太严格的FDR 0.2。
sccoda_model.set_fdr(sccoda_data, 0.2)
为了获得每种细胞类型的成分变化的二元分类，我们在结果对象上使用scCODA的credible_effects函数。标记为“True”的每种细胞类型或多或少都存在。倍数变化描述了细胞类型是否更多或更少存在。因此，我们将把它们与下面的二元分类一起绘制。
sccoda_model.credible_effects(sccoda_data, modality_key="coda")
Covariate Cell Type
condition[T.Hpoly.Day3] Endocrine False
Enterocyte False
Enterocyte.Progenitor False
Goblet False
Stem False
TA False
TA.Early False
Tuft False
condition[T.Hpoly.Day10] Endocrine False
Enterocyte True
Enterocyte.Progenitor False
Goblet False
Stem False
TA False
TA.Early False
Tuft True
condition[T.Salmonella] Endocrine False
Enterocyte True
Enterocyte.Progenitor False
Goblet False
Stem False
TA False
TA.Early False
Tuft False
Name: Final Parameter, dtype: bool
为了将倍数变化与二元分类一起绘制出来，我们可以轻松地使用effects_bar_plot函数。
pt.pl.coda.effects_barplot(sccoda_data, "coda", "condition")
plt.show()
这些图很好地显示了条件对细胞类型的显着且可信的影响。这些效应在很大程度上与 Haber论文中的发现一致，该论文使用非组合泊松回归模型得出了他们的发现：
“沙门氏菌感染后，成熟肠上皮细胞的频率大幅增加。”[Haber et al., 2017]
“Heligmosomoides polygyrus 导致杯状细胞和簇状细胞丰度增加。”[Haber 等人，2017]
15.5 带有标记的簇和层次结构
除了每种细胞类型的丰度之外，典型的单细胞数据集还以基于树的分层排序的形式包含有关不同细胞相似性的信息。这些层次结构可以通过基因表达的聚类（通常用于发现属于同一细胞类型的细胞簇）自动确定，也可以通过生物信息层次结构（如细胞谱系）自动确定。 tascCODA是scCODA的扩展，它将分层信息和实验协变量数据集成到成分计数数据的生成模型中。这对于提高分辨率的细胞图谱工作特别有益。
从本质上讲，它使用与scCODA几乎相同的狄利克雷多项式设置，但扩展了模型，从而对细胞类型集产生影响，这些细胞类型被定义为树结构中的内部节点。
import schist
warnings.filterwarnings("ignore")
warnings.simplefilter("ignore")
要使用tascCODA，我们首先必须定义细胞类型的分层排序。一种可能的层次聚类使用八种细胞类型，并根据它们在sc.tl.dendrogram的PCA表示中的相似性（皮尔逊相关性）对它们进行排序。由于这种结构在我们的数据中非常简单，因此不会给我们带来很多新的见解，因此我们希望有一个更复杂的聚类。获得此类簇的最新方法是schist包，它使用嵌套随机块模型，以不同的分辨率级别对细胞群进行聚类。使用标准设置运行该方法需要一些时间（在我们的数据上大约需要 15 分钟），并为我们提供了将每个单元格分配给adata.obs中的层次聚类的信息。首先，我们需要通过PCA嵌入定义细胞之间的距离度量：
# use logcounts to calculate PCA and neighbors
adata.layers["counts"] = adata.X.copy()
adata.layers["logcounts"] = sc.pp.log1p(adata.layers["counts"]).copy()
adata.X = adata.layers["logcounts"].copy()
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=30, random_state=1234)
# Calculate UMAP for visualization purposes
sc.tl.umap(adata)
WARNING: You’re trying to run this on 15215 dimensions of `.X`, if you really want this, set `use_rep='X'`.
Falling back to preprocessing with `sc.pp.pca` and default params.
然后，我们可以在AnnData对象上运行schist，这会产生通过adata.obs中的一组列“nsbm_level_{i}”定义的聚类：
schist.inference.nested_model(adata, samples=100, random_seed=5678)
adata.obs
objc[13409]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x12f0f1c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x1448ce6b0). One of the two will be used. Which one is undefined.
objc[13410]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x1265f3c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x13bdb76b0). One of the two will be used. Which one is undefined.
objc[13408]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x125a9ec30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x13b1576b0). One of the two will be used. Which one is undefined.
objc[13411]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x129969c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x13f0d36b0). One of the two will be used. Which one is undefined.
objc[13407]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x127cb9c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x13d4106b0). One of the two will be used. Which one is undefined.
objc[13414]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x12ee9ac30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x1446806b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[13490]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x124cf9c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x13a4136b0). One of the two will be used. Which one is undefined.
objc[13492]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x131710c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x146e806b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[13660]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x13455ec30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x149c2f6b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[13699]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x131764c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x146ee86b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[13757]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x12ad09c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x1404af6b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[14239]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x1278c5c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x13d0326b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[14327]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x132a2ec30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x1481d76b0). One of the two will be used. Which one is undefined.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
Numba: Attempted to fork from a non-main thread, the TBB library may be in an invalid state in the child process.
objc[14348]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x1343f7c30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x149ba86b0). One of the two will be used. Which one is undefined.
objc[14356]: Class GNotificationCenterDelegate is implemented in both /opt/anaconda3/lib/libgio-2.0.0.dylib (0x12f11cc30) and /usr/local/Cellar/glib/2.74.4/lib/libgio-2.0.0.dylib (0x1448ce6b0). One of the two will be used. Which one is undefined.
batch barcode condition cell_label scCODA_sample_id nsbm_level_0 nsbm_level_1 nsbm_level_2 nsbm_level_3 nsbm_level_4 nsbm_level_5
index
B1_AAACATACCACAAC_Control_Enterocyte.Progenitor B1 AAACATACCACAAC Control Enterocyte.Progenitor B1 0 0 0 0 0 0
B1_AAACGCACGAGGAC_Control_Stem B1 AAACGCACGAGGAC Control Stem B1 1 5 3 1 0 0
B1_AAACGCACTAGCCA_Control_Stem B1 AAACGCACTAGCCA Control Stem B1 10 2 2 1 0 0
B1_AAACGCACTGTCCC_Control_Stem B1 AAACGCACTGTCCC Control Stem B1 34 3 3 1 0 0
B1_AAACTTGACCACCT_Control_Enterocyte.Progenitor B1 AAACTTGACCACCT Control Enterocyte.Progenitor B1 91 35 0 0 0 0
... ... ... ... ... ... ... ... ... ... ... ...
B10_TTTCACGACAAGCT_Salmonella_TA B10 TTTCACGACAAGCT Salmonella TA B10 6 5 3 1 0 0
B10_TTTCAGTGAGGCGA_Salmonella_Enterocyte B10 TTTCAGTGAGGCGA Salmonella Enterocyte B10 142 36 4 1 0 0
B10_TTTCAGTGCGACAT_Salmonella_Stem B10 TTTCAGTGCGACAT Salmonella Stem B10 112 1 1 1 0 0
B10_TTTCAGTGTGACCA_Salmonella_Endocrine B10 TTTCAGTGTGACCA Salmonella Endocrine B10 146 36 4 1 0 0
B10_TTTCAGTGTTCTCA_Salmonella_Enterocyte.Progenitor B10 TTTCAGTGTTCTCA Salmonella Enterocyte.Progenitor B10 77 14 6 3 0 0
9842 rows × 11 columns
UMAP图很好地显示了schist的聚类如何与细胞类型分配相关联。层次结构第1层的表示是对上面级别的严格细化，即第2层的每个簇都分为多个较小的簇：
sc.pl.umap(
adata, color=["nsbm_level_1", "nsbm_level_2", "cell_label"], ncols=3, wspace=0.5
)
现在，我们将细胞级数据转换为样本级数据并创建树。我们以与scCODA相同的方式创建
tasccoda_model对象，但使用由schist和树级别定义的聚类。
Tasccoda的加载函数将准备一个MuData对象，并将我们的树表示转换为ete树结构并将其保存为tasccoda_data['coda'].uns["tree"]。为了获得一些不太小的簇，我们在最后一层之前砍掉了树，省略了“nsbm_level_0”。
tasccoda_model = pt.tl.Tasccoda()
tasccoda_data = tasccoda_model.load(
adata,
type="cell_level",
cell_type_identifier="nsbm_level_1",
sample_identifier="batch",
covariate_obs=["condition"],
levels_orig=["nsbm_level_4", "nsbm_level_3", "nsbm_level_2", "nsbm_level_1"],
add_level_name=True,
)
tasccoda_data
MuData object with n_obs × n_vars = 9852 × 15256
2 modalities
rna: 9842 x 15215
obs: 'batch', 'barcode', 'condition', 'cell_label', 'scCODA_sample_id', 'nsbm_level_0', 'nsbm_level_1', 'nsbm_level_2', 'nsbm_level_3', 'nsbm_level_4', 'nsbm_level_5'
uns: 'neighbors', 'umap', 'schist', 'nsbm_level_1_colors', 'nsbm_level_2_colors', 'cell_label_colors'
obsm: 'X_pca', 'X_umap', 'CM_nsbm_level_0', 'CM_nsbm_level_1', 'CM_nsbm_level_2', 'CM_nsbm_level_3', 'CM_nsbm_level_4', 'CM_nsbm_level_5'
layers: 'counts', 'logcounts'
obsp: 'distances', 'connectivities'
coda: 10 x 41
obs: 'condition', 'batch'
var: 'n_cells'
uns: 'tree'
pt.pl.coda.draw_tree(tasccoda_data)
tascCODA中的模型设置和执行与scCODA类似，参考的自由参数和公式也相同。此外，我们可以通过
pen_args参数中的参数phi和lambda_1调整树聚合和模型选择。在这里，我们使用无偏设置phi=0和模型选择，该模型选择比默认值lambda_1=1.7稍微宽松一些。我们使用簇 18 作为参考，因为它与内分泌细胞组几乎相同。
tasccoda_model.prepare(
tasccoda_data,
modality_key="coda",
reference_cell_type="18",
formula="condition",
pen_args={"phi": 0, "lambda_1": 3.5},
tree_key="tree",
)
Zero counts encountered in data! Added a pseudocount of 0.5.
MuData object with n_obs × n_vars = 9852 × 15256
2 modalities
rna: 9842 x 15215
obs: 'batch', 'barcode', 'condition', 'cell_label', 'scCODA_sample_id', 'nsbm_level_0', 'nsbm_level_1', 'nsbm_level_2', 'nsbm_level_3', 'nsbm_level_4', 'nsbm_level_5'
uns: 'neighbors', 'umap', 'schist', 'nsbm_level_1_colors', 'nsbm_level_2_colors', 'cell_label_colors'
obsm: 'X_pca', 'X_umap', 'CM_nsbm_level_0', 'CM_nsbm_level_1', 'CM_nsbm_level_2', 'CM_nsbm_level_3', 'CM_nsbm_level_4', 'CM_nsbm_level_5'
layers: 'counts', 'logcounts'
obsp: 'distances', 'connectivities'
coda: 10 x 41
obs: 'condition', 'batch'
var: 'n_cells'
uns: 'tree', 'scCODA_params'
obsm: 'covariate_matrix', 'sample_counts'
tasccoda_model.run_nuts(
tasccoda_data, modality_key="coda", rng_key=1234, num_samples=10000, num_warmup=1000
)
sample: 100%|██████████| 11000/11000 [04:50<00:00, 37.83it/s, 127 steps of size 3.18e-02. acc. prob=0.97]
tasccoda_model.summary(tasccoda_data, modality_key="coda")
Compositional Analysis summary
┌────────────────────────────────────────────┬────────────────────────────────────────────────────────────────────┐
│ Name │ Value │
├────────────────────────────────────────────┼────────────────────────────────────────────────────────────────────┤
│ Data │ Data: 10 samples, 41 cell types │
│ Reference cell type │ 18 │
│ Formula │ condition │
└────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Intercepts │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Final Parameter Expected Sample │
│ Cell Type │
│ 0 1.313 53.195 │
│ 1 1.098 42.904 │
│ 2 1.205 47.749 │
│ 3 0.526 24.215 │
│ 4 -0.707 7.057 │
│ 5 0.634 26.976 │
│ 6 -0.432 9.290 │
│ 7 1.038 40.405 │
│ 8 1.276 51.263 │
│ 9 1.345 54.925 │
│ 10 0.625 26.735 │
│ 11 0.817 32.394 │
│ 12 -0.359 9.994 │
│ 13 0.260 18.559 │
│ 14 0.851 33.514 │
│ 15 0.524 24.166 │
│ 16 0.934 36.414 │
│ 17 -0.142 12.416 │
│ 18 0.684 28.360 │
│ 19 0.857 33.716 │
│ 20 0.198 17.443 │
│ 21 0.209 17.636 │
│ 22 -0.159 12.206 │
│ 23 0.913 35.658 │
│ 24 1.190 47.038 │
│ 25 0.057 15.149 │
│ 26 -0.086 13.131 │
│ 27 -0.002 14.281 │
│ 28 0.786 31.405 │
│ 29 -0.589 7.940 │
│ 30 -0.713 7.014 │
│ 31 0.210 17.654 │
│ 32 -0.797 6.449 │
│ 33 -0.806 6.391 │
│ 34 -0.839 6.184 │
│ 35 -0.104 12.897 │
│ 36 1.443 60.580 │
│ 37 0.215 17.742 │
│ 38 -1.062 4.948 │
│ 39 -0.879 5.942 │
│ 40 0.084 15.564 │
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Effects │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Effect Expected Sample log2-fold change │
│ Covariate Cell Type │
│ conditionT.Hpoly.Day3 0 0.000 51.027 -0.060 │
│ 1 0.000 41.155 -0.060 │
│ 2 -0.257 35.423 -0.431 │
│ 3 0.439 36.030 0.573 │
│ 4 0.000 6.769 -0.060 │
│ 5 0.439 40.139 0.573 │
│ 6 0.000 8.912 -0.060 │
│ 7 0.000 38.759 -0.060 │
│ 8 0.439 76.276 0.573 │
│ 9 -0.257 40.746 -0.431 │
│ 10 0.000 25.645 -0.060 │
│ 11 0.000 31.073 -0.060 │
│ 12 0.000 9.586 -0.060 │
│ 13 0.000 17.803 -0.060 │
│ 14 0.000 32.148 -0.060 │
│ 15 0.000 23.181 -0.060 │
│ 16 0.000 34.930 -0.060 │
│ 17 0.000 11.910 -0.060 │
│ 18 0.000 27.204 -0.060 │
│ 19 0.000 32.342 -0.060 │
│ 20 0.000 16.733 -0.060 │
│ 21 0.439 26.242 0.573 │
│ 22 0.000 11.709 -0.060 │
│ 23 -0.257 26.453 -0.431 │
│ 24 0.000 45.121 -0.060 │
│ 25 0.000 14.532 -0.060 │
│ 26 0.000 12.596 -0.060 │
│ 27 0.000 13.699 -0.060 │
│ 28 0.000 30.125 -0.060 │
│ 29 0.000 7.617 -0.060 │
│ 30 0.000 6.729 -0.060 │
│ 31 0.000 16.935 -0.060 │
│ 32 -0.257 4.784 -0.431 │
│ 33 0.000 6.131 -0.060 │
│ 34 0.000 5.932 -0.060 │
│ 35 0.000 12.371 -0.060 │
│ 36 0.000 58.111 -0.060 │
│ 37 0.000 17.019 -0.060 │
│ 38 0.000 4.746 -0.060 │
│ 39 0.000 5.699 -0.060 │
│ 40 0.439 23.158 0.573 │
│ conditionT.Hpoly.Day10 0 -1.759 12.539 -2.085 │
│ 1 -0.786 26.759 -0.681 │
│ 2 -1.637 12.716 -1.909 │
│ 3 0.000 33.144 0.453 │
│ 4 0.373 14.025 0.991 │
│ 5 0.000 36.924 0.453 │
│ 6 0.000 12.716 0.453 │
│ 7 0.000 55.305 0.453 │
│ 8 0.000 70.166 0.453 │
│ 9 -1.637 14.627 -1.909 │
│ 10 0.000 36.593 0.453 │
│ 11 -0.242 34.808 0.104 │
│ 12 -0.242 10.739 0.104 │
│ 13 0.000 25.403 0.453 │
│ 14 -0.242 36.012 0.104 │
│ 15 -0.242 25.968 0.104 │
│ 16 0.000 49.842 0.453 │
│ 17 0.000 16.994 0.453 │
│ 18 0.000 38.817 0.453 │
│ 19 0.000 46.148 0.453 │
│ 20 -0.242 18.744 0.104 │
│ 21 0.000 24.140 0.453 │
│ 22 -0.242 13.116 0.104 │
│ 23 -1.637 9.496 -1.909 │
│ 24 -1.597 13.038 -1.851 │
│ 25 0.000 20.736 0.453 │
│ 26 -0.242 14.110 0.104 │
│ 27 0.000 19.548 0.453 │
│ 28 -0.242 33.746 0.104 │
│ 29 0.000 10.868 0.453 │
│ 30 0.000 9.601 0.453 │
│ 31 0.000 24.164 0.453 │
│ 32 1.217 29.810 2.209 │
│ 33 0.564 15.377 1.267 │
│ 34 1.186 27.712 2.164 │
│ 35 0.000 17.652 0.453 │
│ 36 -1.716 14.907 -2.023 │
│ 37 0.000 24.285 0.453 │
│ 38 0.000 6.772 0.453 │
│ 39 0.000 8.132 0.453 │
│ 40 0.000 21.303 0.453 │
│ conditionT.Salmonella 0 0.000 34.663 -0.618 │
│ 1 0.000 27.957 -0.618 │
│ 2 0.000 31.114 -0.618 │
│ 3 0.000 15.779 -0.618 │
│ 4 0.000 4.598 -0.618 │
│ 5 0.000 17.578 -0.618 │
│ 6 0.000 6.054 -0.618 │
│ 7 0.000 26.329 -0.618 │
│ 8 0.000 33.404 -0.618 │
│ 9 0.213 44.286 -0.311 │
│ 10 0.000 17.421 -0.618 │
│ 11 0.000 21.108 -0.618 │
│ 12 0.000 6.512 -0.618 │
│ 13 0.000 12.094 -0.618 │
│ 14 2.173 191.842 2.517 │
│ 15 1.547 73.971 1.614 │
│ 16 0.000 23.728 -0.618 │
│ 17 0.000 8.090 -0.618 │
│ 18 0.000 18.480 -0.618 │
│ 19 0.000 21.970 -0.618 │
│ 20 0.000 11.367 -0.618 │
│ 21 0.000 11.492 -0.618 │
│ 22 0.000 7.954 -0.618 │
│ 23 0.000 23.235 -0.618 │
│ 24 0.000 30.651 -0.618 │
│ 25 0.000 9.872 -0.618 │
│ 26 1.547 40.192 1.614 │
│ 27 0.000 9.306 -0.618 │
│ 28 1.547 96.127 1.614 │
│ 29 0.000 5.174 -0.618 │
│ 30 0.000 4.571 -0.618 │
│ 31 0.000 11.504 -0.618 │
│ 32 0.000 4.202 -0.618 │
│ 33 0.000 4.165 -0.618 │
│ 34 0.000 4.030 -0.618 │
│ 35 0.000 8.404 -0.618 │
│ 36 0.000 39.475 -0.618 │
│ 37 0.000 11.561 -0.618 │
│ 38 0.000 3.224 -0.618 │
│ 39 0.000 3.872 -0.618 │
│ 40 0.000 10.142 -0.618 │
└─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Nodes │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Covariate=condition[T.Hpoly.Day10]_node │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Final Parameter Is credible │
│ Node │
│ nsbm_level_4_0 0.00 False │
│ nsbm_level_3_2 0.00 False │
│ nsbm_level_3_0 0.00 False │
│ nsbm_level_3_1 0.00 False │
│ nsbm_level_3_3 -0.24 True │
│ nsbm_level_2_8 0.00 False │
│ nsbm_level_2_10 0.00 False │
│ 10 0.00 False │
│ 31 0.00 False │
│ nsbm_level_2_0 0.00 False │
│ nsbm_level_2_7 0.00 False │
│ nsbm_level_2_11 0.00 False │
│ nsbm_level_2_3 0.00 False │
│ nsbm_level_2_2 -1.64 True │
│ nsbm_level_2_13 0.00 False │
│ nsbm_level_2_1 0.00 False │
│ nsbm_level_2_4 0.00 False │
│ nsbm_level_2_6 0.00 False │
│ nsbm_level_2_14 0.00 False │
│ 11 0.00 False │
│ 16 0.00 False │
│ 37 0.00 False │
│ 19 0.00 False │
│ 27 0.00 False │
│ 30 0.00 False │
│ 0 -1.76 True │
│ 35 0.00 False │
│ 17 0.00 False │
│ 4 0.37 True │
│ 25 0.00 False │
│ 13 0.00 False │
│ 29 0.00 False │
│ 38 0.00 False │
│ 5 0.00 False │
│ 3 0.00 False │
│ 8 0.00 False │
│ 40 0.00 False │
│ 21 0.00 False │
│ 2 0.00 False │
│ 23 0.00 False │
│ 9 0.00 False │
│ 32 2.85 True │
│ 6 0.00 False │
│ 34 1.19 True │
│ 7 0.00 False │
│ 1 -0.79 True │
│ 24 -1.60 True │
│ 18 0.00 False │
│ 36 -1.72 True │
│ 33 0.56 True │
│ 39 0.00 False │
│ 26 0.00 False │
│ 14 0.00 False │
│ 28 0.00 False │
│ 15 0.00 False │
│ 12 0.00 False │
│ 20 0.00 False │
│ 22 0.00 False │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Covariate=condition[T.Hpoly.Day3]_node │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Final Parameter Is credible │
│ Node │
│ nsbm_level_4_0 0.00 False │
│ nsbm_level_3_2 0.00 False │
│ nsbm_level_3_0 0.00 False │
│ nsbm_level_3_1 0.00 False │
│ nsbm_level_3_3 0.00 False │
│ nsbm_level_2_8 0.00 False │
│ nsbm_level_2_10 0.00 False │
│ 10 0.00 False │
│ 31 0.00 False │
│ nsbm_level_2_0 0.00 False │
│ nsbm_level_2_7 0.00 False │
│ nsbm_level_2_11 0.00 False │
│ nsbm_level_2_3 0.44 True │
│ nsbm_level_2_2 -0.26 True │
│ nsbm_level_2_13 0.00 False │
│ nsbm_level_2_1 0.00 False │
│ nsbm_level_2_4 0.00 False │
│ nsbm_level_2_6 0.00 False │
│ nsbm_level_2_14 0.00 False │
│ 11 0.00 False │
│ 16 0.00 False │
│ 37 0.00 False │
│ 19 0.00 False │
│ 27 0.00 False │
│ 30 0.00 False │
│ 0 0.00 False │
│ 35 0.00 False │
│ 17 0.00 False │
│ 4 0.00 False │
│ 25 0.00 False │
│ 13 0.00 False │
│ 29 0.00 False │
│ 38 0.00 False │
│ 5 0.00 False │
│ 3 0.00 False │
│ 8 0.00 False │
│ 40 0.00 False │
│ 21 0.00 False │
│ 2 0.00 False │
│ 23 0.00 False │
│ 9 0.00 False │
│ 32 0.00 False │
│ 6 0.00 False │
│ 34 0.00 False │
│ 7 0.00 False │
│ 1 0.00 False │
│ 24 0.00 False │
│ 18 0.00 False │
│ 36 0.00 False │
│ 33 0.00 False │
│ 39 0.00 False │
│ 26 0.00 False │
│ 14 0.00 False │
│ 28 0.00 False │
│ 15 0.00 False │
│ 12 0.00 False │
│ 20 0.00 False │
│ 22 0.00 False │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Covariate=condition[T.Salmonella]_node │
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Final Parameter Is credible │
│ Node │
│ nsbm_level_4_0 0.00 False │
│ nsbm_level_3_2 0.00 False │
│ nsbm_level_3_0 0.00 False │
│ nsbm_level_3_1 0.00 False │
│ nsbm_level_3_3 0.00 False │
│ nsbm_level_2_8 0.00 False │
│ nsbm_level_2_10 0.00 False │
│ 10 0.00 False │
│ 31 0.00 False │
│ nsbm_level_2_0 0.00 False │
│ nsbm_level_2_7 0.00 False │
│ nsbm_level_2_11 0.00 False │
│ nsbm_level_2_3 0.00 False │
│ nsbm_level_2_2 0.00 False │
│ nsbm_level_2_13 0.00 False │
│ nsbm_level_2_1 0.00 False │
│ nsbm_level_2_4 0.00 False │
│ nsbm_level_2_6 1.55 True │
│ nsbm_level_2_14 0.00 False │
│ 11 0.00 False │
│ 16 0.00 False │
│ 37 0.00 False │
│ 19 0.00 False │
│ 27 0.00 False │
│ 30 0.00 False │
│ 0 0.00 False │
│ 35 0.00 False │
│ 17 0.00 False │
│ 4 0.00 False │
│ 25 0.00 False │
│ 13 0.00 False │
│ 29 0.00 False │
│ 38 0.00 False │
│ 5 0.00 False │
│ 3 0.00 False │
│ 8 0.00 False │
│ 40 0.00 False │
│ 21 0.00 False │
│ 2 0.00 False │
│ 23 0.00 False │
│ 9 0.21 True │
│ 32 0.00 False │
│ 6 0.00 False │
│ 34 0.00 False │
│ 7 0.00 False │
│ 1 0.00 False │
│ 24 0.00 False │
│ 18 0.00 False │
│ 36 0.00 False │
│ 33 0.00 False │
│ 39 0.00 False │
│ 26 0.00 False │
│ 14 0.63 True │
│ 28 0.00 False │
│ 15 0.00 False │
│ 12 0.00 False │
│ 20 0.00 False │
│ 22 0.00 False │
└────────────────────────────────────────────────────────────────────────────────────────────────
同样，接受概率恰好在tascCODA的期望值0.85左右，表明优化没有明显问题。
tascCODA的结果首先应该被解释为对树节点的影响。节点上的非零参数意味着该节点下所有细胞类型的聚合计数显着变化。我们可以轻松地将这三种疾病状态的树状图可视化。蓝色圆圈表示增加，红色圆圈表示减少：
pt.pl.coda.draw_effects(
tasccoda_data,
modality_key="coda",
tree="tree",
covariate="condition[T.Salmonella]",
show_leaf_effects=False,
show_legend=False,
)
pt.pl.coda.draw_effects(
tasccoda_data,
modality_key="coda",
tree="tree",
covariate="condition[T.Hpoly.Day3]",
show_leaf_effects=False,
show_legend=False,
)
pt.pl.coda.draw_effects(
tasccoda_data,
modality_key="coda",
tree="tree",
covariate="condition[T.Hpoly.Day10]",
show_leaf_effects=False,
show_legend=False,
)
或者，对内部节点的影响也可以通过树转换到细胞类型级别，从而允许像scCODA 一样计算对数倍数变化。为了可视化细胞类型的对数倍数变化，我们做了与scCODA相同的绘图，灵感来自“高分辨率单细胞图谱揭示了非小细胞肺癌中组织驻留中性粒细胞的多样性和可塑性”。
pt.pl.coda.effects_barplot(tasccoda_data, modality_key="coda", covariates="condition")
<seaborn.axisgrid.FacetGrid at 0x19ad2cd90>
通过绘制UMAP嵌入上每个条件的效应大小，并将其与细胞类型分配进行比较，可以获得另一种富有洞察力的表示：
kwargs = {"ncols": 3, "wspace": 0.25, "vcenter": 0, "vmax": 1.5, "vmin": -1.5}
pt.pl.coda.effects_umap(
tasccoda_data,
effect_name=[
"effect_df_condition[T.Salmonella]",
"effect_df_condition[T.Hpoly.Day3]",
"effect_df_condition[T.Hpoly.Day10]",
],
cluster_key="nsbm_level_1",
**kwargs
)
sc.pl.umap(
tasccoda_data["rna"], color=["cell_label", "nsbm_level_1"], ncols=2, wspace=0.5
)
结果与scCODA的发现非常相似：
对于沙门氏菌感染，我们得到簇的聚集增加，大约代表细胞类型簇中的肠细胞。对于簇12，这种增加甚至更强，如对叶水平的额外积极影响所示
对于Heligmosomoides polygyrus感染，3 天后我们没有得到可信的变化。10 天后，我们发现包含干细胞和转运扩增细胞的细胞簇减少，肠上皮细胞和肠上皮细胞祖细胞的减少也不太明显，scCODA也发现了这一点。
15.6 最后
由于本人最近要毕业了，要写毕业论文找gz等，很多事情堆到了一起，没有办法抽出足够的时间来更新，接下来一段时间，更新的频率会比较随缘，希望大家谅解。但是，我也保证，等忙完这段时间，我会回归的，谢谢大家。
小鱼“美照”奉上。
