Bioconductor Version

To begin, install Bioconductor and then install BatchQC:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("BatchQC")

Github Version

To install the most up-to-date version of BatchQC, please install directly from github. You will need the devtools package. You can install both of these with the following commands:

if (!require("devtools", quietly = TRUE)) {
    install.packages("devtools")   
}
library(devtools)
install_github("wejlab/BatchQC")

Load BatchQC and Launch Shiny App

You should now be able to load BatchQC and launch the shiny app.

library(BatchQC)

BatchQC()

1. Upload data set

Upon launching the shiny app, you will be on the “Upload Data” screen. From here, you may upload the following data:

1. Your own data w/Count and Metadata file. Count file should have sample ID as column and item of interest as row. Metadata should have sample ID as row and all meta data labeled in the column. Both files should be uploaded as .csv by selecting “Browse” and selecting the appropriate files from your computer.
2. Your own data as a Summarized Experiment object. Upload should be a .RDS file type and can be uploaded by selecting “Browse” and then the appropriate file from your computer
3. An example data set. See the examples vignette for additional information on each example data set

A preview of the selected data will show up under the “Input Summary” and “Full Metadata” tab. Input summary shows the counts file/assay on the Input summary tab and the metadata on the “Full Metadata” tab. If a preview does not load properly, please ensure your dataset is structured properly and a valid file type.

You MUST hit “Upload” for your data to be available to interact with the shiny app.

All shiny app functions can also be run from the command line; all command line execution require that your data is in a Summarized Experiment object. You can create a summarized experiment object from a counts and metadata matrix. The remainder of this documentation will utilize the signature data example data set included in the package. This can be loaded from the command line as a summarized experiment object with the following:

data(signature_data)
data(batch_indicator)
se_object <- BatchQC::summarized_experiment(signature_data, batch_indicator)
SummarizedExperiment::assayNames(se_object) <- "log_intensity"
se_object$batch <- as.factor(se_object$batch)
se_object$condition <- as.factor(se_object$condition)

Likewise, you can upload your own data into a summarized experiment object for use with command line functions by providing a data matrix and sample info matrix to the summarized_experiment function.

2. Apply Normalization methods

After you have successfully uploaded your data set of interest, you can apply one of the following Normalization methods if appropriate:

1. CPM or counts per million. CPM calculates the counts mapped to a feature relative to the total counts mapped to a sample times one million. CPM may be used to adjust expression count biases introduced by sequencing depth. CPM adjusted values are not recommended for differential expression analysis or within sample comparison.
2. DESeq. DESeq calculates the counts mapped to a feature divided by sample- specific size factors. Size factors are determined by the median ratio of gene counts relative to the geometric mean per feature. DESeq may be used to adjust expression count biases introduced by sequencing depth and RNA composition. DESeq adjusted values are not recommended for within sample comparison.

To use one of these two methods, select it from the normalization method drop down, select an assay to perform the normalization on from the drop down, and name the new assay (the prepopulated name can be changed to your preference). You may also choose to log-transform the results of either of these normalization methods or your raw data by selecting the “Log transform” box. Hit “Normalize” to complete the analysis. Your new assay will now be available for further analysis.

Alternatively, this can be called directly from the command line on an SE object and will return the same SE object with an additional assay with the name that you provide as the output_assay_name.

Our signature data set already has a log adjusted assay, so normalization is not needed. Examples are provided below on how you would complete this if needed via the command line:

# CPM Normalization
se_object <- BatchQC::normalize_SE(se = se_object, method = "CPM",
    log_bool = FALSE, assay_to_normalize = "log_intensity",
    output_assay_name = "CPM_normalized_counts")

# CPM Normalization w/log
se_object <- BatchQC::normalize_SE(se = se_object, method = "CPM",
    log_bool = TRUE, assay_to_normalize = "log_intensity",
    output_assay_name = "CPM_log_normalized_counts")

# DESeq Normalization
se_object <- BatchQC::normalize_SE(se = se_object, method = "DESeq",
    log_bool = FALSE, assay_to_normalize = "log_intensity",
    output_assay_name = "DESeq_normalized_counts")

# DESeq Normalization w/log
se_object <- BatchQC::normalize_SE(se = se_object, method = "DESeq",
    log_bool = TRUE, assay_to_normalize = "log_intensity",
    output_assay_name = "DESeq_log_normalized_counts")

# log adjust
se_object <- BatchQC::normalize_SE(se = se_object, method = "none",
    log_bool = TRUE, assay_to_normalize = "log_intensity",
    output_assay_name = "log_normalized_counts")

3. Apply Batch Effect Correction

Select the appropriate Batch Effect correction model to create a batch-corrected assay for your data. You may select one of the following:

1. Combat-Seq. Combat-Seq uses a negative binomial regression to model batch effects. It requires untransformed, raw count data to adjust for batch effect. Only use this model on an untransformed raw counts assay (not normalized assays).
2. Combat. Combat corrects for Batch effect using a parametric empirical Bayes framework and data should be cleaned and normalized. Select a normalized assay to run this on.

To apply your desired batch effect correction model, select the method from the drop down menu, then select the appropriate assay in the second drop down menu, select your batch variable (labeled batch in all built-in examples), and the covariates that you would like to preserve (not including the batch variable). If desired, revise the name for the corrected assay and then hit “Correct.”

A progress bar will pop up in the bottom right hand corner and a message will state “Correction Complete” in the bottom right hand corner when complete.

You are now ready to analyze your raw and batch effect corrected data sets!

The signature data is log normalized, so Combat is the appropriate batch correction tool. Therefore, to run from the command line use the following:

# Combat correction
se_object <- BatchQC::batch_correct(se = se_object, method = "ComBat",
    assay_to_normalize = "log_intensity", batch = "batch",
    covar = c("condition"), output_assay_name = "Combat_corrected")

## Found 53 genes with uniform expression within a single batch (all zeros); these will not be adjusted for batch.

## Found3batches

## Adjusting for1covariate(s) or covariate level(s)

## Standardizing Data across genes

## Fitting L/S model and finding priors

## Finding parametric adjustments

## Adjusting the Data

Combat-Seq can be applied to a data set by providing “Combat-Seq” to the method parameter in the function call.

4. Experimental Design

batch_design_tibble <- BatchQC::batch_design(se = se_object, batch = "batch", 
    covariate = "condition")
batch_design_tibble

## # A tibble: 10 × 4
## # Groups:   condition [10]
##    condition   `1`   `2`   `3`
##    <fct>     <int> <int> <int>
##  1 1             6    12     9
##  2 2             6     0     0
##  3 3             0     6     0
##  4 4             0     6     0
##  5 5             0     5     0
##  6 6             0     6     0
##  7 7             0     6     0
##  8 8             0     0     9
##  9 9             0     0     9
## 10 10            0     0     9

confound_table <- BatchQC::confound_metrics(se = se_object, batch = "batch")
confound_table

##           Pearson Correlation Coefficient Cramer's V
## condition                       0.9179524  0.8005757

pearson_cor_result <- BatchQC::std_pearson_corr_coef(batch_design_tibble)
pearson_cor_result

## [1] 0.9179524

cramers_v_result <- BatchQC::cramers_v(batch_design_tibble)
cramers_v_result

## [1] 0.8005757

5. Variation Analysis

First, select your unadjusted assay, your batch variable and any covariates of interest. Click “Here we go!” to see the variation analysis results. The “Individual Variable” tab shows the individual, or raw variation, explained by each variable. The “Residual” tab displays the type 2, or residual, variation for each variable. A box plot is displayed along with a table with the variation specific to each gene listed. This table is searchable and can be displayed in ascending or descending order. Notice the amount of variation attributed to batch. In our signature data set, most of the variation is attributed to batch rather than the condition. This indicates a need for a batch correction.

The “Individual Variation Variable/Batch Ratio” tab displays the individual, or raw variation, divided by the batch. A ratio greater that 1 indicates that batch has a stronger affect than the variable of interest. And the “Residual Variation Variable/Batch Ratio” displays the residual, or type 2 variation, divided by the batch. A ratio greater that 1 indicates that batch has a stronger affect than the variable of interest.

Next, select the batch corrected assay, with the same batch and condition variable and observe the changes in variation. When batch effect is being properly corrected, you should notice that the variation explained by the batch variable has decreased in the batch corrected assay, indicating that the variation in the data is now more likely to be associated with your condition of interest rather than the batch. The ratio statistics should also be less than 1. This is visible in our example data set and confirms the need for batch adjustment.

To recreate the variation analysis plot for the log_intensity assay, use the following code:

ex_variation <- batchqc_explained_variation(se = se_object,
    batch = "batch", condition = "condition", assay_name = "log_intensity")
EV_boxplot <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[1]])
EV_boxplot

EV_boxplot_type_2 <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[2]])
EV_boxplot_type_2

To recreate the variation analysis table for the log_intensity assay, use the following code:

ex_variation <- batchqc_explained_variation(se = se_object,
    batch = "batch", condition = "condition", assay_name = "log_intensity")
EV_table <- BatchQC::EV_table(batchqc_ev = ex_variation[[1]])
EV_table

##       Gene Name Explained batch condition Unexplained
##          <char>     <num> <num>     <num>       <num>
##    1:     Gene1     36.61 26.52     14.12       63.39
##    2:     Gene2     60.24 57.70     44.63       39.76
##    3:     Gene3     49.31 34.32     36.55       50.69
##    4:     Gene4     17.31  3.01     16.32       82.69
##    5:     Gene5     85.54 81.74     69.30       14.46
##   ---                                                
## 1596:  Gene1596     96.29 90.68     62.02        3.71
## 1597:  Gene1597     92.51 87.16     54.09        7.49
## 1598:  Gene1598     74.53 28.47     56.57       25.47
## 1599:  Gene1599     16.14  5.56     12.94       83.86
## 1600:  Gene1600     27.77 15.21     19.30       72.23

EV_table_type_2 <- BatchQC::EV_table(batchqc_ev = ex_variation[[2]])
EV_table_type_2

##       Gene Name Explained batch condition Unexplained
##          <char>     <num> <num>     <num>       <num>
##    1:     Gene1     36.61 22.49     10.09       63.39
##    2:     Gene2     60.24 15.62      2.54       39.76
##    3:     Gene3     49.31 12.76     14.99       50.69
##    4:     Gene4     17.31  1.00     14.31       82.69
##    5:     Gene5     85.54 16.24      3.80       14.46
##   ---                                                
## 1596:  Gene1596     96.29 34.27      5.61        3.71
## 1597:  Gene1597     92.51 38.41      5.34        7.49
## 1598:  Gene1598     74.53 17.96     46.06       25.47
## 1599:  Gene1599     16.14  3.20     10.58       83.86
## 1600:  Gene1600     27.77  8.47     12.56       72.23

We can then compare our batch corrected assay, using the following code:

ex_variation <- batchqc_explained_variation(se = se_object,
    batch = "batch", condition = "condition", assay_name = "Combat_corrected")
EV_boxplot <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[1]])
EV_boxplot

EV_boxplot_type_2 <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[2]])
EV_boxplot_type_2

EV_table <- BatchQC::EV_table(batchqc_ev = ex_variation[[1]])
EV_table

##       Gene Name Explained batch condition Unexplained
##          <char>     <num> <num>     <num>       <num>
##    1:     Gene1     12.83  0.98      6.99       87.17
##    2:     Gene2      5.37  0.01      4.60       94.63
##    3:     Gene3     23.71  1.65     23.53       76.29
##    4:     Gene4     12.85  0.23     12.72       87.15
##    5:     Gene5     24.47  1.53     22.35       75.53
##   ---                                                
## 1596:  Gene1596     49.51  2.84     43.03       50.49
## 1597:  Gene1597     29.33  0.12     29.29       70.67
## 1598:  Gene1598     70.83 19.44     69.39       29.17
## 1599:  Gene1599     16.14  5.56     12.94       83.86
## 1600:  Gene1600     15.65  0.05     13.21       84.35

EV_table_type_2 <- BatchQC::EV_table(batchqc_ev = ex_variation[[2]])
EV_table_type_2

##       Gene Name Explained batch condition Unexplained
##          <char>     <num> <num>     <num>       <num>
##    1:     Gene1     12.83  5.84     11.85       87.17
##    2:     Gene2      5.37  0.77      5.36       94.63
##    3:     Gene3     23.71  0.18     22.06       76.29
##    4:     Gene4     12.85  0.12     12.61       87.15
##    5:     Gene5     24.47  2.12     22.94       75.53
##   ---                                                
## 1596:  Gene1596     49.51  6.48     46.68       50.49
## 1597:  Gene1597     29.33  0.04     29.21       70.67
## 1598:  Gene1598     70.83  1.44     51.39       29.17
## 1599:  Gene1599     16.14  3.20     10.58       83.86
## 1600:  Gene1600     15.65  2.45     15.60       84.35

6. Heatmaps

Select your assay of interest and the batch and condition(s) of interest. We recommend viewing no more than 500 elements at a time (or the max number in your data set).

heatmaps <- BatchQC::heatmap_plotter(se = se_object, assay = "log_intensity",
    nfeature = 38, annotation_column = c("batch", "condition"),
    log_option = "FALSE")
correlation_heatmap <- heatmaps$correlation_heatmap
correlation_heatmap

heatmaps <- BatchQC::heatmap_plotter(se = se_object, assay = "log_intensity",
    nfeature = 38, annotation_column = c("batch", "condition"),
    log_option = FALSE)
heatmap <- heatmaps$topn_heatmap
heatmap

7. Dendrograms

dendrogram_plot <- BatchQC::dendrogram_plotter(se = se_object,
    assay = "log_intensity",
    batch_var = "batch",
    category_var = "condition")

## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

dendrogram_plot$dendrogram

circular_dendrogram_plot <- BatchQC::dendrogram_plotter(
    se = se_object, assay = "log_intensity", batch_var = "batch",
    category_var = "condition")

## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

circular_dendrogram_plot$circular_dendrogram

8. PCA Analysis

You may select multiple assays to plot PCA analysis side by side. Therefore, select your raw data assay and batch corrected assay, the batch variable as shape, condition as color, and the number of top variable features you are interested in. Review the clustering pattern. In data sets where a strong batch effect is seen, similar shapes (batch) will all cluster together. Whereas if a batch effect is not seen, shapes (batch) should be dispersed throughout and you should see clustering by condition (color).

Use the following code to recreate the PCA plot and associated table that lists the explained variation of each PC for our signature data set:

pca_plot <- BatchQC::PCA_plotter(se = se_object, nfeature = 20, color = "batch",
    shape = "condition", assays = c("log_intensity", "Combat_corrected"),
    xaxisPC = 1, yaxisPC = 2, log_option = FALSE)
pca_plot$plot

## Warning: The shape palette can deal with a maximum of 6 discrete values because more
## than 6 becomes difficult to discriminate
## ℹ you have requested 10 values. Consider specifying shapes manually if you need
##   that many have them.

## Warning: Removed 66 rows containing missing values or values outside the scale range
## (`geom_point()`).

pca_plot$var_explained

##     log_intensity Combat_corrected
## PC1       0.89763          0.86091
## PC2       0.07170          0.08066

9. Differential Expression Analysis

To explore which features are differentially expressed between conditions, use the differential expression tab. DE_Seq2 should be used for count data (non-negative, integer values) while limma can be used for all other data.

Our signature data is logged data, so we will use limma, we will select the assay we are interested in, select the batch variable and select all covariates that we would like to include in the analysis. The “Results Table” return a table for each condition within each variable. The results table displays the log2 fold change, p-value and adjusted p-value for each feature. The “P-Value Analysis” tab displays a box plot of p-values for each analysis. In both the “Results Table” and “P-Value Analysis” tabs, the tables can be sorted by ascending or descending order of a column of your choice. They are also searchable with the “Search” bar (which can be helpful if you are interested in a specific feature or p-value threshold).

To run the limma differential expression analysis from the command line, run the following code (which can be modified to run DESeq2 when appropriate by changing the method argument to “DESeq2”):

differential_expression <- BatchQC::DE_analyze(se = se_object,
    method = "limma", batch = "batch", conditions = c("condition"),
    assay_to_analyze = "log_intensity")

After running the differential expression analysis, you can see all the completed analysis by running this code:

names(differential_expression)

##  [1] "(Intercept)" "condition2"  "condition3"  "condition4"  "condition5" 
##  [6] "condition6"  "condition7"  "condition8"  "condition9"  "condition10"
## [11] "batch2"      "batch3"

You can then view the log2 fold change, p-value and adjusted p-value table for a specific analysis of interest, by running the following:

head(differential_expression[["batch2"]])

##          log2FoldChange        pvalue         padj
## Gene931        5.234550 2.061868e-101 3.298989e-98
## Gene1187       4.726761 4.218134e-101 3.374507e-98
## Gene932        5.710023  1.837234e-95 9.798582e-93
## Gene933        5.665289  2.934568e-94 1.173827e-91
## Gene1361       4.790767  1.559233e-84 4.989545e-82
## Gene330        3.422975  9.181434e-83 2.448382e-80

After running a differential expression analysis, you may replicate the p-value box plot and table from the command line by running the following code:

pval_plotter(differential_expression)

head(pval_summary(differential_expression))

##            (Intercept)   condition2   condition3   condition4   condition5
## Gene1106 2.510236e-145 1.167947e-30 5.161973e-73 2.772024e-80 2.221457e-20
## Gene689  1.541911e-136 5.513534e-17 6.200256e-33 3.196720e-30 2.407150e-19
## Gene806  1.067757e-133 1.413820e-11 2.936216e-27 5.936389e-24 2.701257e-19
## Gene189  2.020368e-133 2.292188e-11 9.935786e-26 5.962531e-24 6.143495e-18
## Gene715  5.591466e-133 4.793547e-11 1.641991e-22 5.175309e-20 1.708665e-16
## Gene284  8.463674e-133 6.122471e-11 1.688931e-20 4.465234e-18 3.569192e-16
##            condition6   condition7   condition8   condition9  condition10
## Gene1106 4.621614e-45 6.362206e-38 1.241689e-59 1.269818e-58 2.469975e-58
## Gene689  1.539260e-41 1.554160e-37 4.520514e-38 7.800657e-42 2.278908e-50
## Gene806  3.234945e-41 4.388781e-35 3.260273e-23 1.141698e-40 3.503531e-39
## Gene189  1.488640e-37 1.719055e-34 4.907893e-23 6.171553e-40 3.201013e-37
## Gene715  1.786228e-37 1.829450e-31 5.039421e-22 5.356080e-35 1.429255e-36
## Gene284  8.294929e-37 5.381066e-29 6.166698e-22 5.775666e-31 3.746781e-30
##                 batch2        batch3
## Gene1106 2.061868e-101 1.446943e-100
## Gene689  4.218134e-101  1.637307e-96
## Gene806   1.837234e-95  1.406129e-94
## Gene189   2.934568e-94  1.492472e-92
## Gene715   1.559233e-84  4.770485e-84
## Gene284   9.181434e-83  2.736390e-81

value <- round((max(abs(
    differential_expression[[length(differential_expression)]][, 1]))
    + min(
        abs(differential_expression[[length(differential_expression)]][, 1])))
    / 2)
volcano_plot(DE_results = differential_expression[["batch2"]],
    pslider = 0.05,
    fcslider = value)

10. Data Download

You may download the data set that you have been working with int he shiny app, including the assays you added in the “Normalization” or “Batch Effect Correction” steps when you uploaded your data. A preview is provided that shows the Summarized Experiment (SE) object, including the dimensions, metadata table, assays, rownames and colNames.

Click the “Download” button and you can save this SE object to your computer for additional analysis outside of the BatchQC shiny app. The object is saved as a .RDS Summarized Experiment object in the folder of your choice. You should give it a descriptive name (default is “se.RDS”).

If you run code from the command line according to the above instructions, all assays will be saved to your original SE object and available for your use in other applications. You may save it as an SE object by running the following:

file_location <- "location/to/save/object"

saveRDS(object = se_object, file = file_location)