Experimental Design - exercises

Scientific Inquiry 2024

1 Introduction

This vignette contains examples generated with ConfoundingExplorer.

The package can be installed from GitHub:

install.packages("remotes")
remotes::install_github("csoneson/ConfoundingExplorer")

After installing the package, an app can be launched by typing

library(ConfoundingExplorer)
ConfoundingExplorer()

in the R console.

2 Overview of the interface

The ConfoundingExplorer interface consists of two parts - a sidebar to the left where the user can change the settings and the experimental design, and a collection of panels displaying various aspects of the generated data and results from statistical tests. The main purpose of the app is to evaluate the ability to correctly detect differentially expressed features between two conditions, when the samples may be collected in different batches.

A detailed description of how the data is generated can be found by clicking on the Data description button in the bottom of the sidebar, but briefly, a base matrix with 1,000 features (rows) and the number of samples (columns) equal to the sum of the numbers in the table in the top of the sidebar is first generated by random draws from a normal distribution with mean 10 and standard deviation 2. Next, condition effects and batch effects are added or subtracted from the samples from a given condition or batch. The strength of the condition and batch signals can be modified using the sliders in the sidebar, and the number of samples in each condition and batch can be changed by modifying the table in the top-left corner.

The Analysis approach section of the sidebar determines how the generated data will be analysed. In each case, a linear model will be used to compare the values for each feature between the two conditions, but depending on the choice, the batch effect will be accounted for in different ways. Again, more details can be obtained by clicking on the Data description button.

The first row of panels in the main window illustrate the results of the statistical test. The left panel shows the number of features in each category (affected by condition effect and/or by batch effect), and the features that are detected as significantly differentially expressed between the conditions are colored in red. The right panel shows a histogram of all the 1,000 p-values. If there is no signal in the data (no truly differentially expressed features) we expect this histogram to show a uniform distribution. If there are some truly differential features, we expect additionally a peak close to p=0. A peak close to p=1 often indicates a problem with the data, e.g., a factor that affects the data but that is not accounted for in the analysis.

The next panel (Summary) provides a summary of the statistical test results. More precisely, it shows the TPR (true positive rate, or power - the fraction of the features that are truly affected by the condition that are detected as significantly differential), the FDP (false discovery proportion - the fraction of the features detected as significant that are false positives), and the FPR (false positive rate - the fraction of the truly non-differential features that are nevertheless found to be significant).

The Data heatmap shows the simulated data in the form of a heatmap. It also indicates the condition and batch assignment of the samples (here shown in the rows), and the true designation of the features (batch- and/or condition-affected) as well as whether they are detected as significantly differential between the two conditions.

Finally, the last row of panels show examples from each of the four categories of features (same categories as in the top-right panel). A new set of examples can be obtained by clicking on the New selection button.

3 Non-confounded setup

First, we consider the setup where there is no confounding between the condition (group) and the batch. In other words, each condition is equally represented in each of the batches. We will assume that 25% of the features are truly different between the two conditions, and that 50% of the features are affected by the batch. The two sets of genes overlap partly. Thus, some features that are truly different between conditions are also affected by the batch.

3.1 Don’t adjust for the batch effect

We first consider the results if we don’t adjust for the batch effect at all, and just compare the two conditions. Since each condition contains samples from both batches, the within-condition variance estimate will also incorporate the differences induced by the batch effect, and thus be artificially inflated compared to the ‘true’ (within-batch) within-condition variance. The effect on the results of the test comparing the two conditions is that we will have lower power to detect true differences (due to the overestimated variance).

app <- ConfoundingExplorer(
    sampleSizes = matrix(rep(5, 4), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "dontAdjust"
)

3.2 Increase the size of the batch effect, don’t adjust

To see the effect described above even clearer, let’s increase the size of the batch effect. This simulates a situation where, for example, five replicates of each condition were obtained by each of two technologies, or using completely different protocols. Now we see the loss of power much more strikingly, and we notice a strangely shaped p-value histogram.

app <- ConfoundingExplorer(
    sampleSizes = matrix(rep(5, 4), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 10,
    analysisApproach = "dontAdjust"
)

3.3 Adjust for the batch effect by including it as a covariate

The recommended way of adjusting for batch effects in statistical modeling is to include them as a covariate (an ‘additional predictor variable’) in the model. Here we follow this approach, and thus model the observed values as a linear combination of the estimated batch effect and the estimated condition effect. As before, the p-values correspond to the test of the condition effect (but now we’re accounting for the batch effect). We note that by including the extra predictor corresponding to the batch effect, we can explain some of the between-sample variance and thus increase the power to detect differences.

app <- ConfoundingExplorer(
    sampleSizes = matrix(rep(5, 4), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "inclBatch"
)

4 Fully confounded setup

Next, we consider a fully confounded setup. Here, the 10 replicates of each condition were all acquired in the same batch, and only samples from one condition were included in each batch.

4.1 Don’t adjust for the batch effect

If we do not adjust for the batch effect in this case, both the batch and condition effects are interpreted as the condition effect, which leads to a high false positive rate. In other words, many features are interpreted as being different between the two conditions, whereas the difference in reality is caused by the batch effect, which is completely confounded with the condition. In addition, some features that are truly differential between the two conditions are no longer detected, since the difference is counteracted by the batch effect.

app <- ConfoundingExplorer(
    sampleSizes = matrix(c(10, 0, 0, 10), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "dontAdjust"
)

4.2 Attempt to adjust for the batch effect

In the non-confounded situation above, we could adjust for the batch effect by including it as a covariate in the model. When the batch and condition factors are fully confounded, however, this is not possible, since it would mean including two perfectly correlated predictors in the model.

app <- ConfoundingExplorer(
    sampleSizes = matrix(c(10, 0, 0, 10), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "inclBatch"
)

5 Partly confounded setup

Finally, we consider a partly confounded setup. Here, each batch contains 9 samples of one condition and 1 sample of the other condition.

5.1 Adjust for the batch effect by including it as a covariate

As in the non-confounded case, we attempt to adjust for the batch effect by including it as a covariate in the model. While it is now statistically feasible (as opposed to the fully confounded case), we notice that the power to detect true differences between the conditions is lower than in the non-confounded case.

app <- ConfoundingExplorer(
    sampleSizes = matrix(c(9, 1, 1, 9), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "inclBatch"
)

5.2 Remove the batch effect in advance

Another approach that is sometimes used, especially when the data is going to be used for exploratory analysis rather than statistical modeling, is to estimate the batch effect in advance, and subtract it from the observed values. Here we illustrate this approach for the partly confounded setup. We notice that almost all the true signal is removed by this procedure, since we do not inform the batch effect correction step of the condition effect, and thus any signal that could potentially be explained by the batch variable will be removed.

app <- ConfoundingExplorer(
    sampleSizes = matrix(c(9, 1, 1, 9), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "removeBatch"
)

5.3 Remove batch effect in advance, accounting for the condition

As an alternative to the approach outlined in the previous section, we could inform the batch effect adjustment step of the condition factor, and instruct it to not remove any effect that could be explained by the latter. Since the batch and condition factors are highly similar, the effect of this is that very little of the batch effect is in fact eliminated.

app <- ConfoundingExplorer(
    sampleSizes = matrix(c(9, 1, 1, 9), nrow = 2,
                         dimnames = list(c("group1", "group2"),
                                         c("batch1", "batch2"))),
    fracVarCond = 0.25,
    fracVarBatch = 0.5,
    condEffectSize = 4,
    batchEffectSize = 4,
    analysisApproach = "removeBatchAccCond"
)

6 Conclusion

The examples above illustrate the complications that arise when the condition of interest is (partly or fully) confounded with a batch effect or another source of unwanted variation affecting the data. It is important to note that in a fully confounded experimental setup, statistical modeling can never ‘rescue’ the experiment, and there is no reliable way of distinguishing the effect of the batch and the condition, respectively. In a non-confounded experiment, even if there is a batch effect, it can typically be accounted for in the statistical modeling.

7 Session info

sessionInfo()
# R version 4.3.2 (2023-10-31)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS Sonoma 14.3
# 
# Matrix products: default
# BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
# LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# time zone: Europe/Zurich
# tzcode source: internal
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] ConfoundingExplorer_0.5.0 BiocStyle_2.30.0         
# 
# loaded via a namespace (and not attached):
#   [1] bitops_1.0-7                gridExtra_2.3              
#   [3] rlang_1.1.3                 magrittr_2.0.3             
#   [5] shinydashboard_0.7.2        clue_0.3-65                
#   [7] GetoptLong_1.0.5            matrixStats_1.2.0          
#   [9] compiler_4.3.2              mgcv_1.9-1                 
#  [11] png_0.1-8                   vctrs_0.6.5                
#  [13] reshape2_1.4.4              stringr_1.5.1              
#  [15] pkgconfig_2.0.3             shape_1.4.6                
#  [17] crayon_1.5.2                fastmap_1.1.1              
#  [19] XVector_0.42.0              ellipsis_0.3.2             
#  [21] fontawesome_0.5.2           utf8_1.2.4                 
#  [23] promises_1.2.1              rmarkdown_2.25             
#  [25] shinyAce_0.4.2              UpSetR_1.4.0               
#  [27] purrr_1.0.2                 xfun_0.41                  
#  [29] zlibbioc_1.48.0             cachem_1.0.8               
#  [31] GenomeInfoDb_1.38.5         jsonlite_1.8.8             
#  [33] highr_0.10                  later_1.3.2                
#  [35] DelayedArray_0.28.0         parallel_4.3.2             
#  [37] cluster_2.1.6               R6_2.5.1                   
#  [39] bslib_0.6.1                 stringi_1.8.3              
#  [41] RColorBrewer_1.1-3          limma_3.58.1               
#  [43] GenomicRanges_1.54.1        jquerylib_0.1.4            
#  [45] shinyMatrix_0.6.0           Rcpp_1.0.12                
#  [47] bookdown_0.37               SummarizedExperiment_1.32.0
#  [49] iterators_1.0.14            knitr_1.45                 
#  [51] IRanges_2.36.0              splines_4.3.2              
#  [53] igraph_1.6.0                httpuv_1.6.14              
#  [55] Matrix_1.6-5                tidyselect_1.2.0           
#  [57] rstudioapi_0.15.0           abind_1.4-5                
#  [59] yaml_2.3.8                  miniUI_0.1.1.1             
#  [61] doParallel_1.0.17           codetools_0.2-19           
#  [63] lattice_0.22-5              tibble_3.2.1               
#  [65] plyr_1.8.9                  shiny_1.8.0                
#  [67] Biobase_2.62.0              ROCR_1.0-11                
#  [69] evaluate_0.23               circlize_0.4.15            
#  [71] pillar_1.9.0                BiocManager_1.30.22        
#  [73] MatrixGenerics_1.14.0       DT_0.31                    
#  [75] foreach_1.5.2               stats4_4.3.2               
#  [77] shinyjs_2.1.0               generics_0.1.3             
#  [79] iSEE_2.14.0                 RCurl_1.98-1.14            
#  [81] S4Vectors_0.40.2            ggplot2_3.4.4              
#  [83] munsell_0.5.0               scales_1.3.0               
#  [85] xtable_1.8-4                glue_1.7.0                 
#  [87] tools_4.3.2                 colourpicker_1.3.0         
#  [89] cowplot_1.1.3               grid_4.3.2                 
#  [91] tidyr_1.3.1                 iCOBRA_1.30.0              
#  [93] shinyBS_0.61.1              colorspace_2.1-0           
#  [95] SingleCellExperiment_1.24.0 nlme_3.1-164               
#  [97] GenomeInfoDbData_1.2.11     vipor_0.4.7                
#  [99] cli_3.6.2                   fansi_1.0.6                
# [101] viridisLite_0.4.2           S4Arrays_1.2.0             
# [103] ComplexHeatmap_2.18.0       dplyr_1.1.4                
# [105] gtable_0.3.4                rintrojs_0.3.4             
# [107] sass_0.4.8                  digest_0.6.34              
# [109] BiocGenerics_0.48.1         ggrepel_0.9.5              
# [111] SparseArray_1.2.3           rjson_0.2.21               
# [113] htmlwidgets_1.6.4           memoise_2.0.1              
# [115] htmltools_0.5.7             lifecycle_1.0.4            
# [117] shinyWidgets_0.8.1          GlobalOptions_0.1.2        
# [119] statmod_1.5.0               mime_0.12