15.1 前言
除了基因表達模式的變化之外，細胞組成（例如細胞類型的比例）也會在不同條件下發(fā)生變化。例如奸例，特定藥物可以誘導(dǎo)細胞類型的轉(zhuǎn)分化，這將反映在細胞身份組合物中向楼。需要足夠的細胞和樣本數(shù)量才能準確確定細胞同一性簇比例和背景變異查吊。可以在已知細胞類型或?qū)?yīng)于例如最近受擾動影響的細胞的細胞狀態(tài)的形式的細胞身份簇的水平上進行組成分析蜜自。
差異豐度分析比較兩種條件下細胞類型的組成菩貌。兩種模式的樣品含有不同比例的細胞類型，可以測試其豐度的顯著變化重荠。
本篇將介紹這兩種方法并將其應(yīng)用于Haber數(shù)據(jù)集[Haber et al., 2017]箭阶。該數(shù)據(jù)集包含來自小鼠小腸和類器官的53193個單個上皮細胞虚茶。一些細胞還受到細菌或蠕蟲感染，例如分別通過沙門氏菌和Heligmosomoides polygyrus感染仇参。在本教程中嘹叫，我們使用完整Haber數(shù)據(jù)集的子集，其中僅包括專門為此目的收集的對照細胞和受感染細胞诈乒。值得注意的是罩扇，我們排除了僅收集大量細胞的附加數(shù)據(jù)集，以加快計算速度并降低復(fù)雜性怕磨。
第一步喂饥，我們加載數(shù)據(jù)集。
15.2 加載數(shù)據(jù)
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
數(shù)據(jù)分10批收集肠鲫。獨特的條件是對照员帮、沙門氏菌、Hpoly.Day3和Hpoly.Day10导饲，分別對應(yīng)于健康對照狀態(tài)捞高、沙門氏菌感染、Heligmosomoides polygyrus感染3天后的細胞和Heligmosomoides polygyrus感染10天后的細胞渣锦。cell_label 對應(yīng)于細胞類型硝岗。
15.3 為什么細胞類型計數(shù)數(shù)據(jù)是組合的
在分析細胞計數(shù)數(shù)據(jù)的成分變化時，需要考慮多種技術(shù)和方法學限制袋毙。一項挑戰(zhàn)是實驗重復(fù)次數(shù)較少型檀，這會導(dǎo)致在使用頻率統(tǒng)計檢驗進行差異豐度分析時出現(xiàn)較大的置信區(qū)間。更重要的是听盖，單細胞測序自然受到每個樣本細胞數(shù)量的限制——我們無法對組織或器官中的每個細胞進行測序贱除，而是使用一個小的月幌、有代表性的快照來代替。然而悬蔽，這迫使我們將細胞類型計數(shù)視為純粹的比例扯躺，即樣本中的細胞總數(shù)只是一個比例因子。在統(tǒng)計文獻中蝎困，此類數(shù)據(jù)被稱為成分數(shù)據(jù)录语，其特征是一個樣本中所有特征（在我們的例子中為細胞類型）的相對豐度總和為1。
由于這種總和為一的限制禾乘，導(dǎo)致細胞類型豐度之間存在負相關(guān)澎埠。為了說明這一點，讓我們考慮以下示例：
在病例對照研究中始藕，我們想要比較健康器官和患病器官的細胞類型組成蒲稳。在這兩種情況下氮趋，我們都有三種細胞類型（A、B 和 C）江耀，但它們的豐度不同：
健康器官由每種類型2000個細胞組成（總共6000個細胞）昵观。
這種疾病導(dǎo)致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函數(shù)從群體中采樣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()
現(xiàn)在圖片已經(jīng)不太清楚了。雖然A型細胞的計數(shù)仍然增加（大約從200個到300個），但其他兩種細胞類型似乎從大約200個減少到150個這種明顯的減少是由于我們對600個細胞的限制造成的-如果樣品被A型細胞占據(jù)了嚎，B型和C型細胞的比例必須較低露筒。因此，如果不考慮其他細胞類型鞍历，就不可能確定一種細胞類型的豐度變化扇救。
如果我們忽略數(shù)據(jù)的組合性装畅，并使用像 Wilcoxon 秩和檢驗或 scDC 這樣的單變量方法，我們可能會錯誤地將細胞類型群體變化視為統(tǒng)計上的聲音效應(yīng)迅诬，盡管它們是由細胞類型比例的固有負相關(guān)性引起的。
此外，二次采樣數(shù)據(jù)不僅為我們的問題提供了一種有效的解決方案搏屑。如果患病病例中B型和C型細胞均減少1000個細胞解幼，我們將獲得與上述相同的600個細胞的代表性樣本害晦。為了獲得唯一的結(jié)果鲫剿，我們可以固定數(shù)據(jù)的參考點，假設(shè)該參考點在所有樣本中都保持不變枚抵。這可以是單一細胞類型、多種細胞類型的聚合（例如幾何平均值）逼泣。
scCODA已經(jīng)解決了這個問題魏蔗，我們將在下一節(jié)中介紹它并將其應(yīng)用于我們的數(shù)據(jù)集廓鞠。
15.4 帶有標記的簇
scCODA屬于需要預(yù)定義簇（最常見的細胞類型）來統(tǒng)計得出成分變化的工具系列砌们。受微生物組數(shù)據(jù)成分分析方法的啟發(fā)，scCODA提出了一種貝葉斯方法來解決單細胞分析中常見的低重復(fù)問題。它使用分層狄利克雷多項式模型對細胞類型計數(shù)進行建模峻堰，該模型通過對所有測量的細胞類型比例進行聯(lián)合建模來解釋細胞類型比例的不確定性和負相關(guān)偏差闹击。為了確保唯一可識別的解決方案和易于解釋狰腌，scCODA中的參考被選擇為特定的細胞類型。因此，scCODA檢測到的任何成分變化始終必須相對于所選參考進行查看。
然而洲赵，scCODA假設(shè)協(xié)變量和細胞豐度之間存在對數(shù)線性關(guān)系绪商，這在使用連續(xù)協(xié)變量時可能并不總是反映潛在的生物過程。 scCODA的另一個限制是無法推斷除成分效應(yīng)之外的細胞成分之間的相關(guān)結(jié)構(gòu)锣尉。此外落蝙，scCODA僅模擬平均豐度的變化，但無法檢測響應(yīng)變異性的變化。
第一步荡陷，我們實例化scCODA模型。
然后，我們使用load函數(shù)準備MuData對象以供后續(xù)處理，并根據(jù)輸入數(shù)據(jù)創(chuàng)建成分分析數(shù)據(jù)集。在我們的例子中浸赫，我們指定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的箱線圖浇坐。為了更好地了解數(shù)據(jù)的分布方式觉渴，紅點顯示實際的數(shù)據(jù)點哎迄。
pt.pl.coda.boxplots(
sccoda_data,
modality_key="coda",
feature_name="condition",
figsize=(12, 5),
add_dots=True,
args_swarmplot={"palette": ["red"]},
)
plt.show()
箱線圖突出顯示了細胞類型分布的一些差異蹬屹。明顯值得注意的是沙門氏菌病癥的腸上皮細胞比例很高弧腥。但其他細胞類型更鲁，例如轉(zhuǎn)運擴增（TA）細胞，與對照相比，沙門氏菌條件下的豐度也表現(xiàn)出明顯差異塘慕。必須正確評估這些差異是否具有統(tǒng)計顯著性筋夏。
另一種可視化是scCODA提供的堆疊條形圖蒂胞。該可視化很好地顯示了成分數(shù)據(jù)的特征：如果我們比較對照組和沙門氏菌組图呢，我們可以看到受感染小鼠中腸上皮細胞的比例大大增加赴叹。由于數(shù)據(jù)是成比例的是辕，這會導(dǎo)致滿足總和為一約束的所有其他細胞類型的份額減少贞谓。
pt.pl.coda.stacked_barplot(
sccoda_data, modality_key="coda", feature_name="condition", figsize=(4, 2)
)
plt.show()
除了細胞計數(shù)AnnData對象之外暇矫，scCODA還需要兩個主要參數(shù)：公式和引用細胞類型夯接。該公式描述了使用R-style指定的協(xié)變量芯丧。在我們的例子中驼唱，我們將條件指定為唯一的協(xié)變量叶雹。由于它是具有四個水平（對照和三種疾病狀態(tài)）的離散協(xié)變量编饺，因此該模型對每個狀態(tài)與其他樣本的比較進行了建模。如果我們想一次對多個協(xié)變量進行建模视译，只需將它們添加到公式中（即
formula = "covariate_1 + covariate_2"）就足夠了。如上所述省古，scCODA需要參考細胞類型進行比較呼股，據(jù)信協(xié)變量不會改變該參考細胞類型。scCODA可以自動選擇適當?shù)募毎愋妥鳛閰⒖迹摷毎愋驮谒袠悠分芯哂袔缀鹾愣ǖ南鄬ωS度破托，也可以使用用戶指定的參考細胞類型運行诫龙。在這里卑吭，我們將內(nèi)分泌細胞設(shè)置為參考笋额，因為從視覺上看它們的豐度似乎相當恒定。手動設(shè)置參考單元類型的另一種方法是將reference_cell_type設(shè)置為“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
接受率描述了在初始老化階段后接受的建議樣本的比例夕玩，并且可以作為不良優(yōu)化運行的臨時指標你弦。對于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根據(jù)包含概率選擇可信效應(yīng)】常可信效應(yīng)和不可信效應(yīng)之間的界限取決于所需的錯誤發(fā)現(xiàn)率 (FDR)旷偿。較小的FDR值會產(chǎn)生更保守的結(jié)果，但可能會遺漏一些效應(yīng)碍沐，而較大的FDR值會選擇更多的效應(yīng)狸捅，但會產(chǎn)生更多的錯誤發(fā)現(xiàn)衷蜓。推理后可以通過 sim_results.set_fdr() 輕松設(shè)置所需的FDR級別累提。默認情況下，該值為0.05磁浇。由于根據(jù)數(shù)據(jù)集的不同斋陪，F(xiàn)DR可能會對結(jié)果產(chǎn)生重大影響，因此我們建議嘗試高達0.2的不同 FDR置吓，以獲得最顯著的效果无虚。
在我們的例子中，我們使用不太嚴格的FDR 0.2衍锚。
sccoda_model.set_fdr(sccoda_data, 0.2)
為了獲得每種細胞類型的成分變化的二元分類友题，我們在結(jié)果對象上使用scCODA的credible_effects函數(shù)。標記為“True”的每種細胞類型或多或少都存在戴质。倍數(shù)變化描述了細胞類型是否更多或更少存在度宦。因此，我們將把它們與下面的二元分類一起繪制告匠。
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
為了將倍數(shù)變化與二元分類一起繪制出來戈抄，我們可以輕松地使用effects_bar_plot函數(shù)。
pt.pl.coda.effects_barplot(sccoda_data, "coda", "condition")
plt.show()
這些圖很好地顯示了條件對細胞類型的顯著且可信的影響后专。這些效應(yīng)在很大程度上與 Haber論文中的發(fā)現(xiàn)一致划鸽，該論文使用非組合泊松回歸模型得出了他們的發(fā)現(xiàn)：
“沙門氏菌感染后，成熟腸上皮細胞的頻率大幅增加÷惴蹋”[Haber et al., 2017]
“Heligmosomoides polygyrus 導(dǎo)致杯狀細胞和簇狀細胞豐度增加嫂用。”[Haber 等人丈冬，2017]
15.5 帶有標記的簇和層次結(jié)構(gòu)
除了每種細胞類型的豐度之外尸折，典型的單細胞數(shù)據(jù)集還以基于樹的分層排序的形式包含有關(guān)不同細胞相似性的信息。這些層次結(jié)構(gòu)可以通過基因表達的聚類（通常用于發(fā)現(xiàn)屬于同一細胞類型的細胞簇）自動確定殷蛇，也可以通過生物信息層次結(jié)構(gòu)（如細胞譜系）自動確定实夹。 tascCODA是scCODA的擴展，它將分層信息和實驗協(xié)變量數(shù)據(jù)集成到成分計數(shù)數(shù)據(jù)的生成模型中粒梦。這對于提高分辨率的細胞圖譜工作特別有益亮航。
從本質(zhì)上講，它使用與scCODA幾乎相同的狄利克雷多項式設(shè)置匀们，但擴展了模型缴淋，從而對細胞類型集產(chǎn)生影響，這些細胞類型被定義為樹結(jié)構(gòu)中的內(nèi)部節(jié)點泄朴。
import schist
warnings.filterwarnings("ignore")
warnings.simplefilter("ignore")
要使用tascCODA重抖，我們首先必須定義細胞類型的分層排序。一種可能的層次聚類使用八種細胞類型祖灰，并根據(jù)它們在sc.tl.dendrogram的PCA表示中的相似性（皮爾遜相關(guān)性）對它們進行排序钟沛。由于這種結(jié)構(gòu)在我們的數(shù)據(jù)中非常簡單，因此不會給我們帶來很多新的見解局扶，因此我們希望有一個更復(fù)雜的聚類恨统。獲得此類簇的最新方法是schist包，它使用嵌套隨機塊模型三妈，以不同的分辨率級別對細胞群進行聚類畜埋。使用標準設(shè)置運行該方法需要一些時間（在我們的數(shù)據(jù)上大約需要 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模燥，這會產(chǎn)生通過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的聚類如何與細胞類型分配相關(guān)聯(lián)咖祭。層次結(jié)構(gòu)第1層的表示是對上面級別的嚴格細化，即第2層的每個簇都分為多個較小的簇：
sc.pl.umap(
adata, color=["nsbm_level_1", "nsbm_level_2", "cell_label"], ncols=3, wspace=0.5
)
現(xiàn)在涧窒，我們將細胞級數(shù)據(jù)轉(zhuǎn)換為樣本級數(shù)據(jù)并創(chuàng)建樹心肪。我們以與scCODA相同的方式創(chuàng)建
tasccoda_model對象，但使用由schist和樹級別定義的聚類纠吴。
Tasccoda的加載函數(shù)將準備一個MuData對象硬鞍，并將我們的樹表示轉(zhuǎn)換為ete樹結(jié)構(gòu)并將其保存為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中的模型設(shè)置和執(zhí)行與scCODA類似，參考的自由參數(shù)和公式也相同伐坏。此外怔匣，我們可以通過
pen_args參數(shù)中的參數(shù)phi和lambda_1調(diào)整樹聚合和模型選擇。在這里桦沉，我們使用無偏設(shè)置phi=0和模型選擇每瞒，該模型選擇比默認值lambda_1=1.7稍微寬松一些。我們使用簇 18 作為參考纯露，因為它與內(nèi)分泌細胞組幾乎相同剿骨。
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左右埠褪，表明優(yōu)化沒有明顯問題浓利。
tascCODA的結(jié)果首先應(yīng)該被解釋為對樹節(jié)點的影響。節(jié)點上的非零參數(shù)意味著該節(jié)點下所有細胞類型的聚合計數(shù)顯著變化钞速。我們可以輕松地將這三種疾病狀態(tài)的樹狀圖可視化贷掖。藍色圓圈表示增加，紅色圓圈表示減少：
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,
)
或者渴语，對內(nèi)部節(jié)點的影響也可以通過樹轉(zhuǎn)換到細胞類型級別苹威，從而允許像scCODA 一樣計算對數(shù)倍數(shù)變化。為了可視化細胞類型的對數(shù)倍數(shù)變化，我們做了與scCODA相同的繪圖丈钙，靈感來自“高分辨率單細胞圖譜揭示了非小細胞肺癌中組織駐留中性粒細胞的多樣性和可塑性”。
pt.pl.coda.effects_barplot(tasccoda_data, modality_key="coda", covariates="condition")
<seaborn.axisgrid.FacetGrid at 0x19ad2cd90>
通過繪制UMAP嵌入上每個條件的效應(yīng)大小，并將其與細胞類型分配進行比較腰鬼，可以獲得另一種富有洞察力的表示：
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
)
結(jié)果與scCODA的發(fā)現(xiàn)非常相似：
對于沙門氏菌感染，我們得到簇的聚集增加浓恳，大約代表細胞類型簇中的腸細胞琢感。對于簇12，這種增加甚至更強糕殉，如對葉水平的額外積極影響所示
對于Heligmosomoides polygyrus感染亩鬼，3 天后我們沒有得到可信的變化。10 天后阿蝶，我們發(fā)現(xiàn)包含干細胞和轉(zhuǎn)運擴增細胞的細胞簇減少雳锋，腸上皮細胞和腸上皮細胞祖細胞的減少也不太明顯，scCODA也發(fā)現(xiàn)了這一點羡洁。
15.6 最后
由于本人最近要畢業(yè)了玷过，要寫畢業(yè)論文找gz等，很多事情堆到了一起，沒有辦法抽出足夠的時間來更新辛蚊，接下來一段時間粤蝎，更新的頻率會比較隨緣，希望大家諒解袋马。但是初澎，我也保證，等忙完這段時間虑凛，我會回歸的碑宴，謝謝大家。
小魚“美照”奉上桑谍。
