MarioniLab / miloR

R package implementation of Milo for testing for differential abundance in KNN graphs
https://bioconductor.org/packages/release/bioc/html/miloR.html
GNU General Public License v3.0
316 stars 20 forks source link

Multiple comparisons gives identical results by using model.contrasts #333

Closed mainharryHR closed 1 month ago

mainharryHR commented 1 month ago

I am trying to study how different cancer types and sex affects the cell abundances in different cell types.

design_matrix <- data.frame(colData(Bcells_milo))[, c("patientNumber", "Cancer", "Sex")]

contrast.all <- c("CancerControl - CancerNo", "CancerControl - CancerCRC", "CancerNo - CancerCRC") model <- model.matrix(~ 0 + Cancer, data=design_matrix) mod.contrasts <- makeContrasts(contrasts=contrast.all, levels=model)

ControlVSNo <- testNhoods(Bcells_milo, design=~ Cancer, design.df=design_matrix, fdr.weighting="graph-overlap", model.contrasts = mod.contrasts[c("CancerControl"), c("CancerControl - **CancerNo**")])

ControlVSNo %>% arrange(SpatialFDR) %>% head() ControlVSNo <- annotateNhoods(Bcells_milo, ControlVSNo, coldata_col = "hairuV2") plotDAbeeswarm(ControlVSNo, group.by = "Cell types", alpha = 0.1)

I tried different combination of following: model.contrasts = mod.contrasts[c("CancerControl"), c("CancerControl - CancerCRC")] , they give me the same results. I feel the model.contrasts is not working properly.

Any comments are welcome.

Screenshot 2024-05-29 at 13 43 33

Many thanks. Harry

MikeDMorgan commented 1 month ago

A good place to learn is this guide from the developers of edgeR: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7873980/

Milo uses largely the same syntax and format for model contrasts.

mainharryHR commented 1 month ago

A good place to learn is this guide from the developers of edgeR: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7873980/

Milo uses largely the same syntax and format for model contrasts.

Thank you for reply. I went through the paper and multiple comparison vignette many times. https://bioconductor.org/packages/release/bioc/vignettes/miloR/inst/doc/milo_contrasts.html
My above pipeline always gives me the identical results. I am almost sure testNhoods(Bcells_milo, design=~ Cancer, design.df=design_matrix, fdr.weighting="graph-overlap", model.contrasts = mod.contrasts) , model.contrasts = mod.contrasts might be not working. In addition, the published papers (2 decent papers) only use the design, but not model.contrasts, for their multiple comparisons by subsetting.

I have been struggling with this part for several days. It will be great if someone could help me to figure out.

I really appreciate your kind help and great packages

mainharryHR commented 1 month ago

Dear Mike, I am wondering if you have good strategies to compare 3 conditions: A, B, C?

I tried the subsetting strategies, the results make sense. But subsetting will not consider the multiple comparisons, which might be not a perfect way. Here comes results from the subsetting strategies.

Screenshot 2024-05-31 at 13 49 08

Here comes the mod.contrast

contrast.all <- c("CancerControl - CancerNo", "CancerControl - CancerCRC", "CancerNo - CancerCRC") model <- model.matrix(~ 0 + Cancer, data=design_matrix) mod.contrasts <- makeContrasts(contrasts=contrast.all, levels = model)

Screenshot 2024-05-31 at 13 25 25

But after testNhoods, I got the identical results no matter how I change the mod.contrasts[c("CancerControl"), c("CancerControl - CancerCRC"), Plus the results are all negative LogFC , which does not make senses as well.

ControlVSNo <- testNhoods(Bcells_milo, design=~ Cancer, design.df=design_matrix, fdr.weighting="graph-overlap", model.contrasts = mod.contrasts[c("CancerControl"), c("CancerControl - CancerNo")])

Screenshot 2024-05-31 at 13 48 17

Any comments are welcome! Thank you very much.

MikeDMorgan commented 1 month ago

GitHub issues are for reporting issues with software, such as bugs and errors. I suggest you post an issue on the Bioconductor forum to explain what you seek to do and ask for advice there.

mainharryHR commented 1 month ago

Thanks Mike,

I think this is a bug in the MiloR package. Even in the official multiple comparison vignette, the results are identical after multiple comparisons . I am experiencing the exact same problems.

Thank you very much for looking into that with your team member.

Have a great day! Best, Harry

Screenshot 2024-06-03 at 10 02 43
MikeDMorgan commented 1 month ago

If you try running each contrast separately rather than a vector do you get a different result? Fundamentally, Milo doesn't do anything to the model contrasts, it just hands them straight over to edgeR for the QL F-test.

MikeDMorgan commented 1 month ago

OK - the issue is in the vignette, not Milo - the vignette passes mod.contrasts to testNhoods when this should be all.contrasts. Milo converts the vector of contrasts into the correct format internally using the limma function makeContrasts. Therefore, you should pass contrast.all to testNhoods, not mod.contrasts. The vignettes have been updated on the devel branch.

mainharryHR commented 1 month ago

If you try running each contrast separately rather than a vector do you get a different result? Fundamentally, Milo doesn't do anything to the model contrasts, it just hands them straight over to edgeR for the QL F-test.

Thank you for information. When subsetting, I have different results.

I am curious how to compute the spacialFDR for multiple comparison?

For multiply comparision, P value can be adjusted by "BH", I guess each row will be assumed as one condition? or cell type?

Many thanks

mainharryHR commented 1 month ago

OK - the issue is in the vignette, not Milo - the vignette passes mod.contrasts to testNhoods when this should be all.contrasts. Milo converts the vector of contrasts into the correct format internally using the limma function makeContrasts. Therefore, you should pass contrast.all to testNhoods, not mod.contrasts. The vignettes have been updated on the devel branch.

contrast.all <- c("CancerControl - CancerNo", "CancerControl - CancerCRC", "CancerNo - CancerCRC") model <- model.matrix(~ 0 + Cancer, data=design_matrix) mod.contrasts <- makeContrasts(contrasts=contrast.all, levels = model)

contrast.res <- testNhoods(Bcells_milo, design=~ Cancer, design.df=design_matrix, fdr.weighting="graph-overlap", model.contrasts = contrast.all),

I made the change accordingly, but it gives me errors. Could you pinpoint which code should I modify?

MikeDMorgan commented 1 month ago

You need to provide the code AND the errors - it's impossible to help and debug any problems without the appropriate information.

Next, you need to make sure that what you are doing is consistent. The model contrasts uses the convention of ~ 0 + Variable but you're passing ~Variable into testNhoods <- these are not equivalent behaviours.

I strongly recommend you re-read the guide to making good design matrices, as well as the limma/edgeR guides on contrast, before proceeding further.

dravichandar commented 3 weeks ago

Hi Mike

I am in a similar boat. And I think your suggestion to past contrast.all instead of mod.contrast fixed the spurious log fold change calculations. However I want to make sure its correctly set up. I am pasting a mock example


contrast_list <- c('Group_varGroup1 - Group_varGroup2','Group_varGroup4 - Group_varGroup3')
contrast.res_Group1vsGroup2 <- testNhoods(milo_obj, design=~ 0 + Group_var, design.df=design_matrix, fdr.weighting="graph-overlap", model.contrasts = contrast_list[1])#first comparison
contrast.res_Group4vsGroup2 <- testNhoods(milo_obj, design=~ 0 + Group_var, design.df=design_matrix, fdr.weighting="graph-overlap", model.contrasts = contrast_list[2])#first comparison```
MikeDMorgan commented 2 weeks ago

Hi @dravichandar - yes that looks correct. You can also pass contrast_list in as a single vector to test both contrasts - but this will return an omnibus test p-value/SpatialFDR; LFCs for each test are returned on each nhood.