_adjusted`. At the moment just an unwanted covariated is allowed at a time.
TidyTranscriptomics
```{r adjust, message=FALSE, warning=FALSE, results='hide'}
counts_SE.norm.adj =
counts_SE.norm |> adjust_abundance( ~ factor_of_interest + batch)
```
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.type"
]
```
## 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.
TidyTranscriptomics
```{r cibersort}
counts_SE.cibersort =
counts_SE |>
deconvolve_cellularity(action="get", cores=1, prefix = "cibersort__")
```
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.type"
]
```
With the new annotated data frame, we can plot the distributions of cell types across samples, and compare them with the nominal cell type labels to check for the purity of isolation. On the x axis we have the cell types inferred by Cibersort, on the y axis we have the inferred proportions. The data is facetted and coloured by nominal cell types (annotation given by the researcher after FACS sorting).
```{r plot_cibersort, eval=FALSE}
counts_SE.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.type`)) +
geom_boxplot() +
facet_wrap(~`Cell.type`) +
my_theme +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), aspect.ratio=1/5)
```
## Test differential cell-type abundance
We can also perform a statistical test on the differential cell-type abundance across conditions
```{r DC, eval=FALSE}
counts_SE |>
test_differential_cellularity(. ~ condition )
```
We can also perform regression analysis with censored data (coxph).
```{r DC_censored, eval=FALSE}
counts_SE |>
test_differential_cellularity(survival::Surv(time, dead) ~ .)
```
## Cluster `samples`
We may want to cluster our data (e.g., using k-means sample-wise). `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 annotation. At the moment only k-means clustering is supported, the plan is to introduce more clustering methods.
**k-means**
TidyTranscriptomics
```{r cluster}
counts_SE.norm.cluster = counts_SE.norm.MDS |>
cluster_elements(method="kmeans", centers = 2, action="get" )
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
count_m_log = log(count_m + 1)
k = kmeans(count_m_log, iter.max = 1000, ...)
cluster = k$cluster
cluster$cell_type = tibble_counts[
match(tibble_counts$sample, rownames(cluster)),
c("Cell.type", "Dim1", "Dim2")
]
```
We can add cluster annotation to the MDS dimension reduced data set and plot.
```{r plot_cluster}
counts_SE.norm.cluster |>
ggplot(aes(x=`Dim1`, y=`Dim2`, color=`cluster_kmeans`)) +
geom_point() +
my_theme
```
**SNN**
Matrix package (v1.3-3) causes an error with Seurat::FindNeighbors used in this method. We are trying to solve this issue. At the moment this option in unaviable.
TidyTranscriptomics
```{r SNN, eval=FALSE, cache=TRUE, message=FALSE, warning=FALSE, results='hide'}
counts_SE.norm.SNN =
counts_SE.norm.tSNE |>
cluster_elements(method = "SNN")
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
library(Seurat)
snn = CreateSeuratObject(count_m)
snn = ScaleData(
snn, display.progress = TRUE,
num.cores=4, do.par = TRUE
)
snn = FindVariableFeatures(snn, selection.method = "vst")
snn = FindVariableFeatures(snn, selection.method = "vst")
snn = RunPCA(snn, npcs = 30)
snn = FindNeighbors(snn)
snn = FindClusters(snn, method = "igraph", ...)
snn = snn[["seurat_clusters"]]
snn$cell_type = tibble_counts[
match(tibble_counts$sample, rownames(snn)),
c("Cell.type", "Dim1", "Dim2")
]
```
```{r SNN_plot, eval=FALSE}
counts_SE.norm.SNN |>
pivot_sample() |>
select(contains("tSNE"), everything())
counts_SE.norm.SNN |>
pivot_sample() |>
gather(source, Call, c("cluster_SNN", "Call")) |>
distinct() |>
ggplot(aes(x = `tSNE1`, y = `tSNE2`, color=Call)) + geom_point() + facet_grid(~source) + my_theme
# Do differential transcription between clusters
counts_SE.norm.SNN |>
mutate(factor_of_interest = `cluster_SNN` == 3) |>
test_differential_abundance(
~ factor_of_interest,
action="get"
)
```
## Drop `redundant` transcripts
We may want to remove redundant elements from the original data set (e.g., samples or transcripts), for example if we want to define cell-type specific signatures with low sample redundancy. `remove_redundancy` takes as arguments a tibble, column names (as symbols; for `sample`, `transcript` and `count`) and returns a tibble with redundant elements removed (e.g., samples). Two redundancy estimation approaches are supported:
+ removal of highly correlated clusters of elements (keeping a representative) with method="correlation"
+ removal of most proximal element pairs in a reduced dimensional space.
**Approach 1**
TidyTranscriptomics
```{r drop}
counts_SE.norm.non_redundant =
counts_SE.norm.MDS |>
remove_redundancy( method = "correlation" )
```
Standard procedure (comparative purpose)
```{r, eval=FALSE}
library(widyr)
.data.correlated =
pairwise_cor(
counts,
sample,
transcript,
rc,
sort = TRUE,
diag = FALSE,
upper = FALSE
) |>
filter(correlation > correlation_threshold) |>
distinct(item1) |>
rename(!!.element := item1)
# Return non redudant data frame
counts |> anti_join(.data.correlated) |>
spread(sample, rc, - transcript) |>
left_join(annotation)
```
We can visualise how the reduced redundancy with the reduced dimentions look like
```{r plot_drop}
counts_SE.norm.non_redundant |>
pivot_sample() |>
ggplot(aes(x=`Dim1`, y=`Dim2`, color=`Cell.type`)) +
geom_point() +
my_theme
```
**Approach 2**
```{r drop2}
counts_SE.norm.non_redundant =
counts_SE.norm.MDS |>
remove_redundancy(
method = "reduced_dimensions",
Dim_a_column = `Dim1`,
Dim_b_column = `Dim2`
)
```
We can visualise MDS reduced dimensions of the samples with the closest pair removed.
```{r plot_drop2}
counts_SE.norm.non_redundant |>
pivot_sample() |>
ggplot(aes(x=`Dim1`, y=`Dim2`, color=`Cell.type`)) +
geom_point() +
my_theme
```
## Other useful wrappers
The above wrapper streamline the most common processing of bulk RNA sequencing data. Other useful wrappers are listed above.
## From BAM/SAM to tibble of gene counts
We can calculate gene counts (using FeatureCounts; Liao Y et al., 10.1093/nar/gkz114) from a list of BAM/SAM files and format them into a tidy structure (similar to counts).
```{r eval=FALSE}
counts = tidybulk_SAM_BAM(
file_names,
genome = "hg38",
isPairedEnd = TRUE,
requireBothEndsMapped = TRUE,
checkFragLength = FALSE,
useMetaFeatures = TRUE
)
```
## From ensembl IDs to gene symbol IDs
We can add gene symbols from ensembl identifiers. This is useful since different resources use ensembl IDs while others use gene symbol IDs. This currently works for human and mouse.
```{r ensembl}
counts_ensembl |> ensembl_to_symbol(ens)
```
## From gene symbol to gene description (gene name in full)
We can add gene full name (and in future description) from symbol identifiers. This currently works for human and mouse.
```{r description}
counts_SE |>
describe_transcript() |>
select(feature, description, everything())
```
## Appendix
```{r}
sessionInfo()
```