1
$\begingroup$

I've posted this question on the Bioconductor support site but haven't received any responses so I figured I'd try here. The post was quite detailed so I'm copying it over mainly as-written. Please let me know if there's anything I can add or clarify that would help!

I'm working on analyzing a proteomic dataset in which I am interested in evaluating the difference in protein expression between case and control and the effect and interaction of gender on expression differences by case status.

I also have a factor (site) and a covariate (age) that I need to adjust for, neither of which I am interested in knowing the effects of, and that is where I'm having trouble. In a typical regression, I would just add them to the model when estimating my coefficients of interest, but, following this extremely helpful guide, I think I can/should include site as a random effects variable (if it's possible to do that with age, a continuous covariate, I would not know how to do that, which is why I don't list it as well even though the effect of age is not of interest to me). One of the issues I'm having with that is that it results in site being included as a blocking variable, and from my experience with blocking in clinical trials, this would not be a blocking variable; it's not that strict, it's just the sites across which observations are distributed.

So far, I've set up how to look at: differential expression by case status; differential expression by gender; differential expression by case status within genders; differential expression by gender within case status; and interaction between case status and gender. I believe that's more or less been coded correctly, but I will still include it and please do let me know if there are any improvements that could be made. I've also set up those same comparisons, but with two different attempts adjusting for age and site. All of that is included below. I am asking if either of those methods for adjusting is appropriate in this case, and if not, how should I adjust for age and site?

Any assistance is greatly appreciated!

This is the rough structure of the data I'm working with, except rather than "expression" being a single column, it's more like 8000 columns, and I have 153 rather than 10 observations, though the distribution of variables is similar.

# Create mock dataset
pheno_dat <- data.frame(matrix(NA, nrow = 10, ncol = 5))
pheno_dat[1,] <- c("case", "female", "A", 44, "case_female")
pheno_dat[2,] <- c("control", "female", "B", 36, "control_female")
pheno_dat[3,] <- c("case", "female", "A", 45, "case_female")
pheno_dat[4,] <- c("control", "female", "C", 38, "control_female")
pheno_dat[5,] <- c("case", "female", "B", 46, "case_female")
pheno_dat[6,] <- c("control", "female", "A", 38, "control_female")
pheno_dat[7,] <- c("case", "male",  "A", 48, "case_male")
pheno_dat[8,] <- c("control", "male", "B", 39, "control_male")
pheno_dat[9,] <- c("case", "male",  "C", 47, "case_male")
pheno_dat[10,] <- c("control",  "male", "C", 41, "control_male")
colnames(pheno_dat) <- c("status", "gender", "site", "age", "stat_gender")
pheno_dat$status <- as.factor(pheno_dat$status)
pheno_dat$status <- relevel(pheno_dat$status, ref="control")
pheno_dat$gender <- as.factor(pheno_dat$gender)
pheno_dat$gender <- relevel(pheno_dat$gender, ref="male")
pheno_dat$age <- as.numeric(pheno_dat$age)
pheno_dat$site <- as.factor(pheno_dat$site)
pheno_dat$stat_gender <- as.factor(pheno_dat$stat_gender)
pheno_dat$stat_gender <- relevel(pheno_dat$stat_gender, ref="control_male")

> pheno_dat
>          status gender site age    stat_gender
> 1         case female    A  44    case_female
> 2      control female    B  36 control_female
> 3         case female    A  45    case_female
> 4      control female    C  38 control_female
> 5         case female    B  46    case_female
> 6      control female    A  38 control_female
> 7         case   male    A  48      case_male
> 8      control   male    B  39   control_male
> 9         case   male    C  47      case_male
> 10     control   male    C  41   control_male

exprs_dat <- as.data.frame(c(9.7424779,
                             8.541871,
                             7.6822924,
                             9.6101021,
                             9.2693605,
                             8.8050989,
                             8.8638765,
                             7.9937876,
                             7.3689429,
                             8.9164766))
colnames(exprs_dat) <- "expression"

# Convert to expression set
exp_set <- as.ExpressionSet(exprs_dat)
exp_set@phenoData@data <- pheno_dat

The modeling strategy I've used is similar between the following two methods of covariate adjustment. I believe this part is done correctly (but again, please correct me if I'm wrong!).

This is my attempt at covariate adjustment by adding age to the design matrix and treating site as a random effects variable.

# Design matrix adjusting for age
means_stat_gender_age <- model.matrix(~ 0 + exp_set$stat_gender + exp_set$age)
colnames(means_stat_gender_age) <- c("control_male", "case_female", "case_male", "control_female", "age")

# Setting up site as a blocking (random effects) variable
cor <- duplicateCorrelation(exp_set, means_stat_gender_age, block=exp_set$site)

> cor$consensus.correlation
> [1] -0.3233333

# Edited to add: the actual consensus correlation is 0.20, so for this example, 
# I'm ignoring that a negative consensus correlation would be removed.

# Model blocking on site 
fit_stat_gender_age_block <- lmFit(object=exp_set, design=means_stat_gender_age, block=exp_set$site, correlation=cor$consensus.correlation)
fit_stat_gender_age_block <- eBayes(fit_stat_gender_age_block)

> topTable(fit_stat_gender_age_block)
>   control_male case_female case_male control_female       age  AveExpr        F      P.Value    adj.P.Val
> 1     -3.97873   -4.984572 -6.532959      -2.566777 0.3094411 8.679429 3894.656 5.733608e-09 5.733608e-09

# Contrasts for comparisons of interest, adjusting for age and blocking by site
contrasts_comp_int_stat_gender_age_block <- makeContrasts(
  caseVScontrol=(case_male+case_female)/2-(control_male+control_female)/2, # Differential expression by case status
  femaleVSmale=(case_female+control_female)/2-(case_male+control_male)/2, # Differential expression by gender
  caseVScontrol_female=case_female-control_female, caseVScontrol_male=case_male-control_male, # Differential expression by case status within gender
  femaleVSmale_case=case_female-case_male, femaleVSmale_control=control_female-control_male, # Differential expression by gender within case status
  statusXgender=(case_female-control_male)-(control_female-control_male)-(case_male-control_male), # Interaction between case status and gender
  levels=colnames(means_stat_gender_age))
fit_contrasts_comp_int_stat_gender_age_block <- contrasts.fit(fit_stat_gender_age_block, contrasts_comp_int_stat_gender_age_block)
fit_contrasts_comp_int_stat_gender_age_block <- eBayes(fit_contrasts_comp_int_stat_gender_age_block)

> topTable(fit_contrasts_comp_int_stat_gender_age_block)
>   caseVScontrol femaleVSmale caseVScontrol_female caseVScontrol_male femaleVSmale_case femaleVSmale_control statusXgender
> 1     -2.486012      1.48017            -2.417794          -2.554229          1.548387             1.411952     0.1364348
>    AveExpr         F   P.Value adj.P.Val
> 1 8.679429 0.9186268 0.4954753 0.4954753

This is my attempt at covariate adjustment by adding age and site to the design matrix.

# Design matrix adjusting for age and site
means_stat_gender_cor <- model.matrix(~ 0 + exp_set$stat_gender + exp_set$age + exp_set$site)
colnames(means_stat_gender_cor) <- c("control_male", "case_female", "case_male", "control_female", "age", "B", "C")
fit_stat_gender_cor <- lmFit(exp_set, means_stat_gender_cor)
fit_stat_gender_cor <- eBayes(fit_stat_gender_cor)

> topTable(fit_stat_gender_cor)
>   control_male case_female case_male control_female       age          B           C  AveExpr        F     P.Value   adj.P.Val
> 1     -2.58555   -3.574291 -5.064257      -1.339882 0.2779798 -0.1102714 -0.04674809 8.679429 83.64332 0.001962154 0.001962154

# Contrasts for comparisons of interest, adjusting for age and site
contrasts_comp_int_stat_gender_cor <- makeContrasts(
  caseVScontrol=(case_male+case_female)/2-(control_male+control_female)/2, # Differential expression by case status
  femaleVSmale=(case_female+control_female)/2-(case_male+control_male)/2, # Differential expression by gender
  caseVScontrol_female=case_female-control_female, caseVScontrol_male=case_male-control_male, # Differential expression by case status within gender
  femaleVSmale_case=case_female-case_male, femaleVSmale_control=control_female-control_male, # Differential expression by gender within case status
  statusXgender=(case_female-control_male)-(control_female-control_male)-(case_male-control_male), # Interaction between case status and gender
  levels=colnames(means_stat_gender_cor))
fit_contrasts_comp_int_stat_gender_cor <- contrasts.fit(fit_stat_gender_cor, contrasts_comp_int_stat_gender_cor)
fit_contrasts_comp_int_stat_gender_cor <- eBayes(fit_contrasts_comp_int_stat_gender_cor)

> topTable(fit_contrasts_comp_int_stat_gender_cor)
>   caseVScontrol femaleVSmale caseVScontrol_female caseVScontrol_male femaleVSmale_case femaleVSmale_control statusXgender
> 1     -2.356557     1.367817            -2.234408          -2.478707          1.489967             1.245668     0.2442984
>    AveExpr         F   P.Value adj.P.Val
> 1 8.679429 0.3691084 0.7826167 0.7826167

If anything is missing or unclear, please let me know, and I'll address it as best as possible!

Cheers!


# please also include the results of running the following in an R session 

sessionInfo( )

R version 4.1.1 (2021-08-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Monterey 12.5

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] statmod_1.4.36       BiocManager_1.30.18  leukemiasEset_1.30.0 destiny_3.8.1        limma_3.50.3         ROCR_1.0-11         
 [7] glmnet_4.1-4         Matrix_1.4-1         gridExtra_2.3        janitor_2.1.0        UniprotR_2.2.0       magrittr_2.0.3      
[13] glue_1.6.2           Gmisc_3.0.0          htmlTable_2.4.1      Rcpp_1.0.9           robust_0.7-1         fit.models_0.64     
[19] sva_3.42.0           BiocParallel_1.28.3  genefilter_1.76.0    mgcv_1.8-40          nlme_3.1-158         report_0.5.1        
[25] broom_1.0.0          flextable_0.7.2      rempsyc_0.0.4.5      devtools_2.4.4       usethis_2.1.6        forcats_0.5.1       
[31] table1_1.4.2         psych_2.2.5          dplyr_1.0.9          tidyr_1.2.0          Hmisc_4.7-0          Formula_1.2-4       
[37] survival_3.3-1       lattice_0.20-45      ggplot2_3.3.6        lmSupport_2.9.13     readxl_1.4.0         Biobase_2.54.0      
[43] BiocGenerics_0.40.0 

loaded via a namespace (and not attached):
#Removed because of character limit
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

The potential advantage of treating a predictor as a random rather than a fixed effect is that you are often able to use up fewer degrees of freedom with a random effect. For a random effect you estimate the distribution of the effects among the multiple categorical sources of the random effects (e.g., among mice in the guide to which you link), rather than one coefficient for each source beyond the first. The potential downside is that a random effect is generally forced to take a Gaussian distribution, which might not be appropriate particularly when there are only a few sources of the random effect.

For random effects in limma, don't worry about the "block" terminology. I think that carries over from old agricultural statistics terminology and simply represents a random intercept term that might be expressed as (1|site) in a model from the lme4 package.

With respect to age, modeling it linearly as you are already only uses up one degree of freedom. It also makes no sense to treat age as a "random effect" as it's a continuous predictor. If anything, you might consider modeling it more flexibly, for example with regression splines, to better capture its association with gene expression for which you want to account.

For site it's perhaps more a question of how many are involved. If you only have 3, as in the example you present, there isn't much advantage to treating site as a random effect. If you have more than half a dozen or so, it might make sense. That said, in your example you seem to have lost some power by including site as a fixed effect in the model matrix versus as a random effect. Note that in your example the consensus.correlation associated with site is fairly negative, in which case the guide to which you link seems to recommend removing the "block" argument entirely and ignoring the random effect.

Your contrasts seem fine, given the way that you have set up the model matrix with 4 categorical predictors for the combinations of sex and treatment.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you @EdM , this is helpful. Regarding the negative consensus correlation, I've edited the text to reflect that while it is negative in this example, in my dataset, it's 0.20, so I'm ignoring its negativity here since it's just for illustrative purposes. Overall, it sounds like you're saying that modeling as a random effect for site doesn't make much difference when there are only three sites (which is true in the full data as well), so it would be acceptable to model the covariates as fixed effects, though this may lead to some power loss? Would that be about right? $\endgroup$
    – Sam
    Commented Aug 24, 2022 at 17:53
  • $\begingroup$ @Sam It's important to account for potential differences among sites. With 3 site values you only use up 2 degrees of freedom treating site as a fixed effect, not a big problem with power loss when you have 153 observations. Add site into the model to correct for potential site/batch effects. You might want to read up further on batch-effect correction with limma as I haven't done that in a while and have no experience using limma with proteomic (as opposed to RNA expression) data. $\endgroup$
    – EdM
    Commented Aug 24, 2022 at 21:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.