# To ensure Rstudio looks up BioConductor packages run:
setRepositories(ind = c(1:6, 8))
# Then install package with
devtools::install_github("biosurf/cyCombine")Thanks to gartician for the conda implementation. Source.
conda install -c gartician r-cycombineThe article introducing cyCombine is published in Nature Communications.
Vignettes are available at biosurf.
There is a total of eight vignettes, which cover a range of different topics:
- A reference manual showing general commands for using cyCombine.
- A performance benchmarking vignette that shows an example of batch correction using five different tools, including calculation of EMD reductions and MAD scores.
- A vignette covering how to detect batch effects.
- A panel merging vignette including an extended discussion of how to use and evaluate merging.
- An example vignette covering how to process a 1-panel mass cytometry (CyTOF) dataset.
- An example vignette covering how to process a 2-panel mass cytometry (CyTOF) dataset.
- A vignette showing an example with integration of datasets from spectral flow cytometry (SFC) and mass cytometry (CyTOF).
- A vignette covering an example with integration of datasets from spectral flow cytometry (SFC), CITE-seq, and mass cytometry (CyTOF).
If you have any issues or questions regarding the use of cyCombine, please do not hesitate to raise an issue on GitHub. In this way, others may also benefit from the answers and discussions.
library("cyCombine")
library("magrittr")
# Directory containing .fcs files
data_dir <- "data/raw"
# Markers of interest
markers <- c("CD20", "CD3", "CD27", "CD45RA", "CD279", "CD5", "CD19", "CD14", "CD45RO", "GranzymeA", "GranzymeK", "FCRL6", "CD355", "CD152", "CD69", "CD33", "CD4", "CD337", "CD8", "CD197", "LAG3", "CD56", "CD137", "CD161", "FoxP3", "CD80", "CD270", "CD275", "CD134", "CD278", "CD127", "KLRG1", "CD25", "HLADR", "TBet", "XCL1")
# The list of markers can also be imported from a panel file (See the reference manual for an example)
# Compile fcs files, down-sample, and preprocess
uncorrected <- prepare_data(data_dir = data_dir,
markers = markers,
metadata = file.path(data_dir, "metadata.xlsx"), # Can also be .csv file or data.frame object
sample_ids = NULL,
batch_ids = "Batch",
filename_col = "FCS_name",
condition = "Set",
down_sample = TRUE,
sample_size = 500000,
seed = 473,
cofactor = 5)
saveRDS(uncorrected, file = "_data/cycombine_raw_uncorrected.RDS")
# Run batch correction
corrected <- uncorrected %>%
batch_correct(markers = markers,
norm_method = "scale", # "rank" is recommended when combining data with heavy batch effects
rlen = 10, # Consider a larger value, if results are not convincing (e.g. 100)
covar = "condition")
saveRDS(corrected, file = "_data/cycombine_raw_corrected.RDS")If your data is in another format than FCS files or a flowset, please convert your data to a tibble, add the relevant columns (sample, batch, covar/condition/anchor), and begin from transform_asinh() (if your data is not yet transformed; otherwise, skip that step as well).
library(cyCombine)
library(magrittr)
# Directory containing .fcs files
data_dir <- "data/raw"
# Markers of interest
markers <- c("CD20", "CD3", "CD27", "CD45RA", "CD279", "CD5", "CD19", "CD14", "CD45RO", "GranzymeA", "GranzymeK", "FCRL6", "CD355", "CD152", "CD69", "CD33", "CD4", "CD337", "CD8", "CD197", "LAG3", "CD56", "CD137", "CD161", "FoxP3", "CD80", "CD270", "CD275", "CD134", "CD278", "CD127", "KLRG1", "CD25", "HLADR", "TBet", "XCL1")
# Compile fcs files, down-sample, and preprocess
flowset <- compile_fcs(data_dir = data_dir,
pattern = "\\.fcs")
# Convert the generated flowset into a tibble
df <- convert_flowset(metadata = file.path(data_dir, "metadata.xlsx"),
sample_ids = NULL,
batch_ids = "Batch",
filename_col = "FCS_name",
condition = "Set",
down_sample = TRUE,
sample_size = 500000,
seed = 473)
# Transform data
uncorrected <- df %>%
transform_asinh(markers = markers)
saveRDS(uncorrected, file = "_data/cycombine_raw_uncorrected.RDS")
# Run batch correction
labels <- uncorrected %>%
normalize(markers = markers,
norm_method = "scale") %>%
create_som(markers = markers,
rlen = 10)
corrected <- uncorrected %>%
correct_data(label = labels,
covar = "condition")
saveRDS(corrected, file = "_data/cycombine_raw_corrected.RDS")# Full analysis - type ?run_analysis to see how you can modify the analysis
run_analysis(tool = "cycombine", data = "raw", data_dir = "_data", markers = markers)
# Otherwise, plots can be made like so:
plot_density(uncorrected = uncorrected,
corrected = corrected,
markers = markers,
filename = 'figs/densities_withcovar.pdf',
markers_per_page = NULL, # Markers per pdf page - thanks to asongggg for the suggestion
ncol = 6 # Number of columns of plots
)
# PCA plot uncorrected
pca1 <- uncorrected %>%
plot_dimred('uncorrected', type = 'pca')
# PCA plot corrected
pca2 <- corrected %>%
plot_dimred('corrected', type = 'pca')
plot_save_two(pca1, pca2, filename = 'figs/pca.png')
# UMAP
# UMAP plot uncorrected
set.seed(473)
sample <- sample(1:nrow(uncorrected), 20000)
plot1 <- plot_dimred(uncorrected[sample,], type = 'umap', name = 'Uncorrected')
plot2 <- plot_dimred(corrected[sample,], type = 'umap', name = 'Corrected')
plot_save_two(plot1, plot2, filename = 'figs/umap.png')After user requests, we have looked into a method by which anchor-based
(i.e. using replicates) correction can be performed with cyCombine. We
have implemented an experimental option for this in the function
correct_data(). Anchor sample information will be used in the same way
as any other covariate for ComBat. Be aware that anchor status is very
likely to be confounded with the experimental condition, and as such,
these two parameters typically cannot be used together. If you want to
give anchors a try with cyCombine, you need to specify which samples are
replicates and which are not. Therefore, one needs to give
“non-replicates” a unique value for each sample. Otherwise, it will
assume “non-replicates” are replicates of the other “non-replicates”,
which of course they are not.
Example code to specify the anchor variable is shown below. The example is a scenario with three batches containing data from 10 healthy donors. The first sample (HD01) is replicated in all batches, and the others are unique to a batch, making the sample IDs: HD01_batch1, HD01_batch2, HD01_batch3, HD02_batch1, HD03_batch1, HD04_batch1, HD05_batch2, HD06_batch2, HD07_batch2, HD08_batch3, HD09_batch3, and HD10_batch3.
uncorrected$anchor <- ifelse(grepl("^HD01", uncorrected$sample), "HD01", uncorrected$sample)If there are 10000 cells in each sample, the table of the batch vs. anchor parameters will look like the following:
table(uncorrected$anchor, uncorrected$batch)
# 1 2 3
# HD01 10000 10000 10000
# HD02_batch1 10000 0 0
# HD03_batch1 10000 0 0
# HD04_batch1 10000 0 0
# HD05_batch2 0 10000 0
# HD06_batch2 0 10000 0
# HD07_batch2 0 10000 0
# HD08_batch3 0 0 10000
# HD09_batch3 0 0 10000
# HD10_batch3 0 0 10000