repository·master·Indexed 41 days ago
https://github.com/tidyomics/tidybulkAn R framework for modular transcriptomic data analysis designed for the tidyverse ecosystem. It provides a unified interface for SummarizedExperiment objects, supporting workflows from abundance normalization and filtering to dimensionality reduction, clustering, and gene enrichment. It features multiplexed differential expression analysis supporting methods such as edgeR, DESeq2, and limma-voom.
tidybulk is an R package designed for modular transcriptomic data analysis. It aims to bring transcriptomics workflows into the tidyverse ecosystem by providing a unified interface that integrates seamlessly with SummarizedExperiment objects.
It is designed to streamline the entire analysis workflow, from raw data processing to biological insight generation, using tidy principles.
When using tidybulk, the fit objects produced by underlying methods (like edgeR or DESeq2) are stored as attributes in the metadata. You can access these to perform standard quality control plots.
For edgeR:
plotBCV() to plot the biological coefficient of variation (BCV) trend.plotMD() to plot the log-fold change vs mean plot.For DESeq2:
plotDispEsts() to plot the mean-variance trend.plotMA() to plot the log-fold change vs the mean plot.# edgeR example
metadata(airway)$tidybulk$edgeR_quasi_likelihood_object |> plotBCV()
# DESeq2 example
metadata(airway)$tidybulk$DESeq2_object |> plotMA()You can use test_differential_expression() to run multiple differential expression methods on a SummarizedExperiment object. By providing a prefix argument, you can multiplex the results, allowing for easy comparison of different statistical approaches (e.g., edgeR, DESeq2, limma-voom) within the same object.
Supported methods include:
edgeR_quasi_likelihood (edgeR QL method)edger_robust_likelihood_ratio (edgeR robust likelihood ratio)DESeq2 (DESeq2 method)limma_voom (limma-voom)limma_voom_sample_weights (limma-voom with sample weights)# Standard differential expression analysis
airway = airway |>
# Use QL method
test_differential_expression(~ dex, method = "edgeR_quasi_likelihood", prefix = "ql__") |>
# Use edger_robust_likelihood_ratio
test_differential_expression(~ dex, method = "edger_robust_likelihood_ratio", prefix = "lr_robust__") |>
# Use DESeq2 method
test_differential_expression(~ dex, method = "DESeq2", prefix = "deseq2__") |>
# Use limma_voom
test_differential_expression(~ dex, method = "limma_voom", prefix = "voom__") |>
# Use limma_voom_sample_weights
test_differential_expression(~ dex, method = "limma_voom_sample_weights", prefix = "voom_weights__") After multiplexing differential expression methods, you can visualize the distribution of p-values using histograms. A flat histogram indicates no differential expression, while a spike near zero indicates genuine signal. Use pivot_transcript() to convert the SummarizedExperiment into a tidy format for plotting with ggplot2.
airway |>
pivot_transcript() |>
select(
ql__PValue,
lr_robust__PValue,
voom__P.Value,
voom_weights__P.Value,
deseq2__pvalue
) |>
pivot_longer(everything(), names_to = "method", values_to = "pvalue") |>
ggplot(aes(x = pvalue, fill = method)) +
geom_histogram(binwidth = 0.01) +
facet_wrap(~method) +
my_theme +
labs(title = "Histogram of p-values across methods")To compare significance and effect size across different methods, you can generate volcano plots. This involves pivoting the transcript data, harmonizing column names (e.g., mapping log2FoldChange and logFC to a single logFC name), and using tidybulk::log10_reverse_trans() for the y-axis scale.
# Create volcano plots
airway |
# Select the columns we want to plot
pivot_transcript() |>
select(
.feature,
ql__logFC, ql__PValue,
lr_robust__logFC, lr_robust__PValue,
voom__logFC, voom__P.Value,
voom_weights__logFC, voom_weights__P.Value,
deseq2__log2FoldChange, deseq2__pvalue
) |>
# Pivot longer to get a tidy data frame
pivot_longer(
- .feature,
names_to = c("method", "stat"),
values_to = "value", names_sep = "__"
) |>
# Harmonize column names
mutate(stat = case_when(
stat %in% c("logFC", "log2FoldChange") ~ "logFC",
stat %in% c("PValue", "pvalue", "P.Value", "p.value") ~ "PValue"
)) |>
pivot_wider(names_from = "stat", values_from = "value") |>
unnest(c(logFC, PValue)) |>
# Plot
ggplot(aes(x = logFC, y = PValue)) +
geom_point(aes(color = PValue < 0.05, size = PValue < 0.05)) +
scale_y_continuous(trans = tidybulk::log10_reverse_trans()) +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "black")) +
scale_size_manual(values = c("TRUE" = 0.5, "FALSE" = 0.1)) +
facet_wrap(~method) +
my_theme +
labs(title = "Volcano Plots by Method")tidybulk provides a way to automatically generate a bibliography of all methods and packages used in your analysis workflow. This is useful for ensuring proper citation of the statistical methods and software employed during your transcriptomic data analysis.
Use the get_bibliography() function on your data object (e.g., a tidybulk object) to retrieve the BibTeX entries.
# Get bibliography of all methods used in our workflow
airway |> get_bibliography()The following utility and validation functions are available:
Utility Functions:
describe_transcript(): Describe transcript characteristics.get_bibliography(): Get bibliography.resolve_complete_confounders_of_non_interest(): Resolve confounders.Validation Functions:
check_if_counts_is_na(): Check if counts contain NA values.check_if_duplicated_genes(): Check for duplicated genes.check_if_wrong_input(): Validate input data.log10_reverse_trans(): Log10 reverse transformation.logit_trans(): Logit transformation.The following function is available for testing differential expression:
test_differential_expression(): Test differential expression using various methods.The following functions are available for normalizing abundance data within the tidybulk framework:
scale_abundance(): Scale abundance data.quantile_normalise_abundance(): Perform quantile normalization.adjust_abundance(): Adjust abundance for unwanted variation.fill_missing_abundance(): Fill missing abundance values.impute_missing_abundance(): Impute missing abundance values.The following functions are available for testing gene enrichment:
test_gene_enrichment(): Test gene enrichment.test_gene_overrepresentation(): Test gene overrepresentation.test_gene_rank(): Test gene rank.The following functions are available for filtering and selecting transcripts:
identify_abundant(): Identify abundant transcripts without removing them.keep_abundant(): Keep abundant transcripts.identify_abundant_per_category(): Identify abundant transcripts per experimental category.keep_abundant_per_category(): Keep transcripts abundant in at least one experimental category.keep_variable(): Keep variable transcripts.filterByExpr(): Filter by expression.The following functions are available for cellularity analysis:
deconvolve_cellularity(): Deconvolve cellularity using various methods.cibersort(): CIBERSORT analysis.