statOmics / msqrob2

Implementation of the MSqRob analysis of differentially expressed proteins using the Features infrastructure
9 stars 10 forks source link

formula for all pairwise comparisons #48

Open Aintzane99 opened 1 year ago

Aintzane99 commented 1 year ago

Hello,

I am using msqrob2 package to perform differential expression analysis on a protein dataset that has 8 conditions. I would like to know if the makeContrast() fuction from MSqRob2 allows to specify testing for all pairwise comparisons (as in the makeContrast() function from MSqRob previous package), without setting a particular condition as the reference or 'intercept'.

Thanks in advance!

Aintzane

leonardblaschek commented 1 year ago

This is far from ideal and can most definitely be done more elegantly, but the following code works to calculate all pairwise comparisons. I think the only thing you'd need to change in the 'Compare genotypes' section is the name of your Qfeatures object. The rest should be universal. I haven't done any tests though.

    ## Estimate protein levels
    pe <- msqrob(object = pe, i = "protein", formula = ~ genotype, overwrite = TRUE)

    ## Compare genotypes
    # Get coef names
    coefs <- names(getCoef(rowData(pe[["protein"]])$msqrobModels[[1]]))[-1] 
    # Generate all combinations of coefs
    coef_comb_data <- expand.grid(coefs, coefs) |>
                      filter(Var1 != Var2)

    # List duplicate combinations (same combination with switched positions)
    indx <- !duplicated(t(apply(coef_comb_data, 1, sort)))

    # Remove duplicates and format combinations
    coef_comb <- coef_comb_data[indx, ] |>
                unite(contrast, Var1, Var2, sep = " - ") |>
                pull(var = contrast, name = NULL)

    # Combine single coefs and coef combinations and add the '=0' for the null hypothesis
    contrast_names <- c(coefs, coef_comb)
    contrast_list <- paste0(c(coefs, coef_comb), "=0")

    # Make contrasts
    contrasts <- makeContrast(
      contrast_list,
      parameterNames = coefs
    )

    # Test all pairwise comparisons
    pe <- hypothesisTest(object = pe, i = "protein", contrast = contrasts, overwrite = TRUE)

Let me know if it doesn't work, maybe we can figure it out together ;)