Do multiple pairwise comparison with scCODA of pertpy

Ecosystem pertpy
1

Hi, I’m using scCODA model to do composition analysis. As shown in tutorial, if we want to try different reference levels, just modify the treatment level in formula arguement.
Now, if we got A,B,C and D four levels, and want to perform multiple pairwise comparison in between them, how to construct the model?
I’m not sure if my workflow is right. I constucted a model(Mudata) with three ‘identical’ modalities except the fomula which defined three different levels.
Please give your opinion and advice, thanks, everyone!

2

OR LIKE THE FOLLOWING THAT STEPWISELY SUBSETS THE MOD:

#Model setup and inference
#HC as control
sccoda_data_all = sccoda_model.prepare(
sccoda_data_all ,
formula=“disease_cluster”,
reference_cell_type=“automatic”
)
#sampling
sccoda_model.run_nuts(sccoda_data_all )
sccoda_model.set_fdr(sccoda_data_all , est_fdr=0.2)
sccoda_model.summary(sccoda_data_all )

#Compare with AM1----------------------------------------------------------------------
sccoda_data_all .mod[“AM1”] = sccoda_data_all [“coda”][
sccoda_data_all [“coda”].obs[“disease_cluster”].isin([“AM1”,“AM2”,“RM1”,“RM2”])
].copy()

sccoda_data_all = sccoda_model.prepare(
sccoda_data_all,
modality_key= ‘AM1’,
formula=“C(disease_cluster, Treatment(‘AM1’))”,
reference_cell_type=“automatic”
)
print(sccoda_data_all )
#sampling for AM1
sccoda_model.run_nuts(sccoda_data_all , modality_key = “AM1”)
sccoda_model.set_fdr(sccoda_data_all , est_fdr=0.2,modality_key= “AM1”)
sccoda_model.summary(sccoda_data_all ,modality_key = “AM1”)

#Compare with AM2------------------------------------------------------------------------
sccoda_data_all.mod[“AM2”] = sccoda_data_all [“coda”][
sccoda_data_all [“coda”].obs[“disease_cluster”].isin([“AM2”,“RM1”,“RM2”])
].copy()

sccoda_data_all = sccoda_model.prepare(
sccoda_data_all,
modality_key= ‘AM2’,
formula=“C(disease_cluster, Treatment(‘AM2’))”,
reference_cell_type=“automatic”
)
print(sccoda_data_all )
#sampling for AM2
sccoda_model.run_nuts(sccoda_data_all, modality_key = “AM2”)
sccoda_model.set_fdr(sccoda_data_all, est_fdr=0.2,modality_key= “AM2”)
sccoda_model.summary(sccoda_data_all,modality_key = “AM2”)

#Compare with RM1------------------------------------------------------------------------
sccoda_data_all.mod[“RM1”] = sccoda_data_all[“coda”][
sccoda_data_all[“coda”].obs[“disease_cluster”].isin([“RM1”,“RM2”])
].copy()

sccoda_data_all= sccoda_model.prepare(
sccoda_data_all,
modality_key= ‘RM1’,
formula=“C(disease_cluster, Treatment(‘RM1’))”,
reference_cell_type=“automatic”
)
print(sccoda_data_all)
#sampling for RM1
sccoda_model.run_nuts(sccoda_data_all, modality_key = “RM1”)
sccoda_model.set_fdr(sccoda_data_all, est_fdr=0.2,modality_key= “RM1”)
sccoda_model.summary(sccoda_data_all,modality_key = “RM1”)

print(sccoda_model.credible_effects(sccoda_data_all))
print(sccoda_model.credible_effects(sccoda_data_all,modality_key=“AM1”))
print(sccoda_model.credible_effects(sccoda_data_all,modality_key=“AM2”))
print(sccoda_model.credible_effects(sccoda_data_all,modality_key=“RM1”))

3

Dear @shihsama,

thank you for using pertpy!
I assume that you refer to scCODA - Modeling options and result analysis - pertpy.

I guess the easiest answer would be to create subsets where you always have both pairs to run a pairwise comparison. Then you don’t have to mess too much with the formulas. With those subsets you can follow the instructions in the docs that you read already.

4

Thanks for your advice, Sir!