_adjusted`. At the moment just an unwanted covariates is allowed at a time.
> Batch effect correction is important for removing technical variation that could confound biological signals.
TidyTranscriptomics
```{r adjust, message=FALSE, warning=FALSE, results='hide'}
airway.norm.adj =
airway.norm |> adjust_abundance( .factor_unwanted = cell, .factor_of_interest = dex, method="combat")
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
library(sva)
count_m_log = log(count_m + 1)
design =
model.matrix(
object = ~ factor_of_interest + batch,
data = annotation
)
count_m_log.sva =
ComBat(
batch = design[,2],
mod = design,
...
)
count_m_log.sva = ceiling(exp(count_m_log.sva) -1)
count_m_log.sva$cell_type = counts[
match(counts$sample, rownames(count_m_log.sva)),
"cell"
]
```
## Deconvolve `Cell type composition`
We may want to infer the cell type composition of our samples (with the algorithm Cibersort; Newman et al., 10.1038/nmeth.3337). `deconvolve_cellularity` takes as arguments a tibble, column names (as symbols; for `sample`, `transcript` and `count`) and returns a tibble with additional columns for the adjusted cell type proportions.
*Note: Cellularity deconvolution requires gene symbols that match the reference data. The airway dataset uses Ensembl IDs, so this example is commented out.*
TidyTranscriptomics
```{r cibersort, eval=FALSE}
# Requires gene symbols that match reference data
# airway.cibersort =
# airway |>
# deconvolve_cellularity( cores=1, prefix = "cibersort__") |>
# pivot_sample()
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
source("CIBERSORT.R")
count_m |> write.table("mixture_file.txt")
results <- CIBERSORT(
"sig_matrix_file.txt",
"mixture_file.txt",
perm=100, QN=TRUE
)
results$cell_type = tibble_counts[
match(tibble_counts$sample, rownames(results)),
"cell"
]
```
*Note: The plotting code is commented out as the deconvolution step is not executed.*
```{r plot_cibersort, eval=FALSE}
# airway.cibersort |>
# pivot_longer(
# names_to= "Cell_type_inferred",
# values_to = "proportion",
# names_prefix ="cibersort__",
# cols=contains("cibersort__")
# ) |>
# ggplot(aes(x=Cell_type_inferred, y=proportion, fill=cell)) +
# geom_boxplot() +
# facet_wrap(~cell) +
# my_theme +
# theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), aspect.ratio=1/5)
```
## Step 7: Cluster Samples
We may want to cluster our samples based on the transcriptomic profiles. `cluster_elements` takes as arguments a tibble, column names (as symbols; for `sample`, `transcript` and `count`) and returns a tibble with additional columns for the cluster labels.
> Clustering helps identify groups of samples with similar expression profiles.
TidyTranscriptomics
```{r cluster}
airway.norm.cluster =
airway.norm |>
cluster_elements(method="kmeans", centers = 2)
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
library(stats)
count_m_log = log(count_m + 1)
count_m_log_centered = count_m_log - rowMeans(count_m_log)
kmeans_result = kmeans(t(count_m_log_centered), centers = 2)
cluster_labels = kmeans_result$cluster
```
## Step 9: Test Differential Abundance (Alternative Method)
We may want to test for differential abundance between conditions. `test_differential_abundance` takes as arguments a tibble, column names (as symbols; for `sample`, `transcript` and `count`) and a formula representing the desired linear model, and returns a tibble with additional columns for the statistics from the hypothesis test (e.g., log fold change, p-value and false discovery rate).
> This demonstrates an alternative approach to differential expression analysis.
TidyTranscriptomics
```{r differential}
airway.norm.de =
airway.norm |>
test_differential_abundance(~ dex, method="edgeR_quasi_likelihood")
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
library(edgeR)
count_m_log = log(count_m + 1)
design = model.matrix(~ dex, data = annotation)
dge = DGEList(counts = count_m)
dge = calcNormFactors(dge)
dge = estimateDisp(dge, design)
fit = glmQLFit(dge, design)
qlf = glmQLFTest(fit, coef=2)
results = topTags(qlf, n = Inf)
```
## Conclusion
Tidybulk provides a unified interface for comprehensive transcriptomic data analysis with seamless integration of SummarizedExperiment objects and tidyverse principles. It streamlines the entire workflow from raw data to biological insights while maintaining compatibility with standard Bioconductor practices.
```{r session-info}
sessionInfo()
```