<!– badges: start –> <!– badges: end –>

In this article we show some examples of the differences in coding between tidybulk/tidyverse and base R. We noted a decrease > 10x of assignments and a decrease of > 2x of line numbers.

## Create `tidybulk` tibble.

    tt = se_mini

## Aggregate duplicated `transcripts`

Tidy transcriptomics “{.r .yellow} rowData(tt)$gene_name = rownames(tt) tt.aggr = tt %>% aggregate_duplicates(.transcript = gene_name) ”
Base R “r temp = data.frame( symbol = dge_list$genes$symbol, dge_list$counts ) dge_list.nr <- by(temp, temp$symbol, function(df) if(length(df[1,1])>0) matrixStats:::colSums(as.matrix(df[,-1])) ) dge_list.nr <- do.call("rbind", dge_list.nr) colnames(dge_list.nr) <- colnames(dge_list) ```

## Scale `counts`

Tidy transcriptomics ”r tt.norm = tt.aggr %>% identify_abundant(factor_of_interest = condition) %>% scale_abundance() “
Base R ”r library(edgeR) dgList <- DGEList(count_m=x,group=group) keep <- filterByExpr(dgList) dgList <- dgList[keep,,keep.lib.sizes=FALSE] [...] dgList <- calcNormFactors(dgList, method="TMM") norm_counts.table <- cpm(dgList) ```

## Filter `variable transcripts`

We may want to identify and filter variable transcripts.

Tidy transcriptomics “r tt.norm.variable = tt.norm %>% keep_variable() ”
Base R “r library(edgeR) x = norm_counts.table s <- rowMeans((x-rowMeans(x))^2) o <- order(s,decreasing=TRUE) x <- x[o[1L:top],,drop=FALSE] norm_counts.table = norm_counts.table[rownames(x)] norm_counts.table$cell_type = tibble_counts[ match( tibble_counts$sample, rownames(norm_counts.table) ), "Cell type" ] ```

## Reduce `dimensions`

Tidy transcriptomics ”r tt.norm.MDS = tt.norm %>% reduce_dimensions(method=“MDS”, .dims = 2) “
Base R ”r library(limma) count_m_log = log(count_m + 1) cmds = limma::plotMDS(ndim = .dims, plot = FALSE) cmds = cmds %$% cmdscale.out %>% setNames(sprintf(“Dim%s”, 1:6)) cmds$cell_type = tibble_counts[ match(tibble_counts$sample, rownames(cmds)), “Cell type” ] “

### PCA

Tidy transcriptomics ”r tt.norm.PCA = tt.norm %>% reduce_dimensions(method=“PCA”, .dims = 2) “
Base R ”r count_m_log = log(count_m + 1) pc = count_m_log %>% prcomp(scale = TRUE) variance = pc$sdev^2 variance = (variance / sum(variance))[1:6] pc$cell_type = counts[ match(counts$sample, rownames(pc)), “Cell type” ] “

### tSNE

Tidy transcriptomics ”r tt.norm.tSNE = breast_tcga_mini_SE %>% tidybulk( sample, ens, count_scaled) %>% identify_abundant() %>% reduce_dimensions( method = “tSNE”, perplexity=10, pca_scale =TRUE ) “
Base R ”r count_m_log = log(count_m + 1) tsne = Rtsne::Rtsne( t(count_m_log), perplexity=10, pca_scale =TRUE )$Y tsne$cell_type = tibble_counts[ match(tibble_counts$sample, rownames(tsne)), “Cell type” ] “

## Rotate `dimensions`

Tidy transcriptomics ”r tt.norm.MDS.rotated = tt.norm.MDS %>% rotate_dimensions(Dim1, Dim2, rotation_degrees = 45, action=“get”) “
Base R ”r rotation = function(m, d) { r = d * pi / 180 ((bind_rows( c(1 = cos®, 2 = -sin®), c(1 = sin®, 2 = cos®) ) %>% as_matrix) %*% m) } mds_r = pca %>% rotation(rotation_degrees) mds_r$cell_type = counts[ match(counts$sample, rownames(mds_r)), “Cell type” ] “

## Test `differential abundance`

Tidy transcriptomics ”r tt.de = tt %>% test_differential_abundance( ~ condition, action=“get”) tt.de “
Base R ”r library(edgeR) dgList <- DGEList(counts=counts_m,group=group) keep <- filterByExpr(dgList) dgList <- dgList[keep,,keep.lib.sizes=FALSE] dgList <- calcNormFactors(dgList) design <- model.matrix(~group) dgList <- estimateDisp(dgList,design) fit <- glmQLFit(dgList,design) qlf <- glmQLFTest(fit,coef=2) topTags(qlf, n=Inf) ```

## Adjust `counts`

Tidy transcriptomics “r tt.norm.adj = tt.norm %>% adjust_abundance( ~ condition + time) ”
Base R “r library(sva) count_m_log = log(count_m + 1) design = model.matrix( object = ~ condition + time, 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`

Tidy transcriptomics ”r tt.cibersort = tt %>% deconvolve_cellularity(action=“get”, cores=1) “
Base R ”r 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" ] ```

## Cluster `samples`

### k-means

Tidy transcriptomics “r tt.norm.cluster = tt.norm.MDS %>% cluster_elements(method="kmeans”, centers = 2, action=“get” ) “
Base R ”r 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”) ] “

### 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.

Tidy transcriptomics ”r tt.norm.SNN = tt.norm.tSNE %>% cluster_elements(method = “SNN”) “
Base R ”r 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”) ] “

## Drop `redundant` transcripts

Tidy transcriptomics ”r tt.norm.non_redundant = tt.norm.MDS %>% remove_redundancy( method = “correlation” ) “
Base R ”r 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 redundant data frame counts %>% anti_join(.data.correlated) %>% spread(sample, rc, - transcript) %>% left_join(annotation) “

## Draw `heatmap`

tidytranscriptomics ”r tt.norm.MDS %>% # filter lowly abundant keep_abundant() %>% # extract 500 most variable genes keep_variable( .abundance = count_scaled, top = 500) %>% # create heatmap heatmap(sample, transcript, count_scaled, transform = log1p) %>% add_tile(Cell type) “
Base R ”r # Example taken from airway dataset from BioC2020 workshop. dgList <- SE2DGEList(airway) group <- factor(dgList$samples$`Cell type`) keep.exprs <- filterByExpr(dgList, group=group) dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE] dgList <- calcNormFactors(dgList) logcounts <- cpm(dgList, log=TRUE) var_genes <- apply(logcounts, 1, var) select_var <- names(sort(var_genes, decreasing=TRUE))[1:500] highly_variable_lcpm <- logcounts[select_var,] colours <- c("#440154FF", "#21908CFF", "#fefada" ) col.group <- c("red","grey")[group] gplots::heatmap.2(highly_variable_lcpm, col=colours, trace="none", ColSideColors=col.group, scale="row") ```

## Draw `density plot`

tidytranscriptomics “r # Example taken from airway dataset from BioC2020 workshop. airway %>% tidybulk() %>% identify_abundant() %>% scale_abundance() %>% pivot_longer(cols = starts_with("counts”), names_to = “source”, values_to = “abundance”) %>% filter(!lowly_abundant) %>% ggplot(aes(x=abundance + 1, color=sample)) + geom_density() + facet_wrap(~source) + scale_x_log10() “
Base R ”r # Example taken from airway dataset from BioC2020 workshop. dgList <- SE2DGEList(airway) group <- factor(dgList$samples$dex) keep.exprs <- filterByExpr(dgList, group=group) dgList <- dgList[keep.exprs,, keep.lib.sizes=FALSE] dgList <- calcNormFactors(dgList) logcounts <- cpm(dgList, log=TRUE) var_genes <- apply(logcounts, 1, var) select_var <- names(sort(var_genes, decreasing=TRUE))[1:500] highly_variable_lcpm <- logcounts[select_var,] colours <- c("#440154FF", "#21908CFF", "#fefada" ) col.group <- c("red","grey")[group] gplots::heatmap.2(highly_variable_lcpm, col=colours, trace="none", ColSideColors=col.group, scale="row") ```

## Appendix

    sessionInfo()

    ## R version 4.2.1 (2022-06-23 ucrt)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] tidySummarizedExperiment_1.6.1 SummarizedExperiment_1.26.1
    ##  [3] Biobase_2.56.0                 GenomicRanges_1.48.0
    ##  [5] GenomeInfoDb_1.32.4            IRanges_2.30.1
    ##  [7] S4Vectors_0.34.0               BiocGenerics_0.42.0
    ##  [9] MatrixGenerics_1.8.1           matrixStats_0.62.0
    ## [11] tidybulk_1.8.2                 ggrepel_0.9.1
    ## [13] ggplot2_3.3.6                  magrittr_2.0.3
    ## [15] tibble_3.1.8                   tidyr_1.2.1
    ## [17] dplyr_1.0.10                   knitr_1.40
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] nlme_3.1-159           bitops_1.0-7           bit64_4.0.5
    ##  [4] httr_1.4.4             SnowballC_0.7.0        backports_1.4.1
    ##  [7] tools_4.2.1            utf8_1.2.2             R6_2.5.1
    ## [10] DBI_1.1.3              lazyeval_0.2.2         mgcv_1.8-40
    ## [13] colorspace_2.0-3       withr_2.5.0            tidyselect_1.1.2
    ## [16] bit_4.0.4              compiler_4.2.1         preprocessCore_1.58.0
    ## [19] cli_3.4.1              DelayedArray_0.22.0    plotly_4.10.0
    ## [22] scales_1.2.1           readr_2.1.2            genefilter_1.78.0
    ## [25] stringr_1.4.1          digest_0.6.29          XVector_0.36.0
    ## [28] pkgconfig_2.0.3        htmltools_0.5.3        fastmap_1.1.0
    ## [31] limma_3.52.4           htmlwidgets_1.5.4      rlang_1.0.6
    ## [34] RSQLite_2.2.17         generics_0.1.3         jsonlite_1.8.0
    ## [37] BiocParallel_1.30.3    tokenizers_0.2.3       RCurl_1.98-1.8
    ## [40] GenomeInfoDbData_1.2.8 Matrix_1.5-1           Rcpp_1.0.9
    ## [43] munsell_0.5.0          fansi_1.0.3            lifecycle_1.0.2
    ## [46] stringi_1.7.8          edgeR_3.38.4           zlibbioc_1.42.0
    ## [49] plyr_1.8.7             Rtsne_0.16             grid_4.2.1
    ## [52] blob_1.2.3             parallel_4.2.1         crayon_1.5.2
    ## [55] lattice_0.20-45        Biostrings_2.64.1      splines_4.2.1
    ## [58] annotate_1.74.0        hms_1.1.2              KEGGREST_1.36.3
    ## [61] locfit_1.5-9.6         pillar_1.8.1           widyr_0.1.5
    ## [64] reshape2_1.4.4         codetools_0.2-18       XML_3.99-0.10
    ## [67] glue_1.6.2             evaluate_0.16          tidytext_0.3.4
    ## [70] data.table_1.14.2      vctrs_0.4.2            png_0.1-7
    ## [73] tzdb_0.3.0             gtable_0.3.1           purrr_0.3.4
    ## [76] assertthat_0.2.1       cachem_1.0.6           xfun_0.33
    ## [79] broom_1.0.1            xtable_1.8-4           janeaustenr_1.0.0
    ## [82] survival_3.4-0         viridisLite_0.4.1      AnnotationDbi_1.58.0
    ## [85] memoise_2.0.1          sva_3.44.0             ellipsis_0.3.2