# Comparison with base R

## Stefano Mangiola

### 2023-04-26

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

    rowData(tt)$gene_name = rownames(tt)
    tt.aggr = tt %>% aggregate_duplicates(.transcript = gene_name)

Base 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

    tt.norm = tt.aggr %>% identify_abundant(factor_of_interest = condition) %>% scale_abundance()

Base 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

    tt.norm.variable = tt.norm %>% keep_variable()

Base 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

    tt.norm.MDS =
      tt.norm %>%
      reduce_dimensions(method="MDS", .dims = 2)

Base 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

    tt.norm.PCA =
      tt.norm %>%
      reduce_dimensions(method="PCA", .dims = 2)

Base 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

    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

    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

    tt.norm.MDS.rotated =
      tt.norm.MDS %>%
    	rotate_dimensions(`Dim1`, `Dim2`, rotation_degrees = 45, action="get")

Base R

    rotation = function(m, d) {
    	r = d * pi / 180
    	((bind_rows(
    		c(`1` = cos(r), `2` = -sin(r)),
    		c(`1` = sin(r), `2` = cos(r))
    	) %>% 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

    tt.de =
    	tt %>%
    	test_differential_abundance( ~ condition, action="get")
    tt.de

Base 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

    tt.norm.adj =
    	tt.norm %>% adjust_abundance(	~ condition + time)

Base 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

    tt.cibersort =
    	tt %>%
    	deconvolve_cellularity(action="get", cores=1)

Base 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

    tt.norm.cluster = tt.norm.MDS %>%
      cluster_elements(method="kmeans",	centers = 2, action="get" )

Base 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

    tt.norm.SNN =
    	tt.norm.tSNE %>%
    	cluster_elements(method = "SNN")

Base 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

    tt.norm.non_redundant =
    	tt.norm.MDS %>%
      remove_redundancy(	method = "correlation" )

Base 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

    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

    # 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

    # 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

    # 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.3.0 RC (2023-04-13 r84269 ucrt)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2022 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
    ##
    ## time zone: America/New_York
    ## tzcode source: internal
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] tidySummarizedExperiment_1.10.0 SummarizedExperiment_1.30.0
    ##  [3] Biobase_2.60.0                  GenomicRanges_1.52.0
    ##  [5] GenomeInfoDb_1.36.0             IRanges_2.34.0
    ##  [7] S4Vectors_0.38.0                BiocGenerics_0.46.0
    ##  [9] MatrixGenerics_1.12.0           matrixStats_0.63.0
    ## [11] tidybulk_1.12.0                 ggrepel_0.9.3
    ## [13] ggplot2_3.4.2                   magrittr_2.0.3
    ## [15] tibble_3.2.1                    tidyr_1.3.0
    ## [17] dplyr_1.1.2                     knitr_1.42
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] DBI_1.1.3               bitops_1.0-7            widyr_0.1.5
    ##  [4] rlang_1.1.0             tidytext_0.4.1          compiler_4.3.0
    ##  [7] RSQLite_2.3.1           mgcv_1.8-42             reshape2_1.4.4
    ## [10] png_0.1-8               vctrs_0.6.2             sva_3.48.0
    ## [13] stringr_1.5.0           pkgconfig_2.0.3         crayon_1.5.2
    ## [16] fastmap_1.1.1           backports_1.4.1         XVector_0.40.0
    ## [19] ellipsis_0.3.2          utf8_1.2.3              tzdb_0.3.0
    ## [22] preprocessCore_1.62.0   purrr_1.0.1             bit_4.0.5
    ## [25] xfun_0.39               zlibbioc_1.46.0         cachem_1.0.7
    ## [28] jsonlite_1.8.4          blob_1.2.4              SnowballC_0.7.1
    ## [31] DelayedArray_0.26.0     BiocParallel_1.34.0     broom_1.0.4
    ## [34] parallel_4.3.0          R6_2.5.1                stringi_1.7.12
    ## [37] limma_3.56.0            genefilter_1.82.0       Rcpp_1.0.10
    ## [40] readr_2.1.4             Matrix_1.5-4            splines_4.3.0
    ## [43] tidyselect_1.2.0        codetools_0.2-19        plyr_1.8.8
    ## [46] lattice_0.21-8          withr_2.5.0             KEGGREST_1.40.0
    ## [49] evaluate_0.20           Rtsne_0.16              survival_3.5-5
    ## [52] Biostrings_2.68.0       pillar_1.9.0            janeaustenr_1.0.0
    ## [55] plotly_4.10.1           generics_0.1.3          RCurl_1.98-1.12
    ## [58] hms_1.1.3               munsell_0.5.0           scales_1.2.1
    ## [61] xtable_1.8-4            glue_1.6.2              lazyeval_0.2.2
    ## [64] tools_4.3.0             tokenizers_0.3.0        data.table_1.14.8
    ## [67] annotate_1.78.0         locfit_1.5-9.7          XML_3.99-0.14
    ## [70] grid_4.3.0              AnnotationDbi_1.62.0    edgeR_3.42.0
    ## [73] colorspace_2.1-0        nlme_3.1-162            GenomeInfoDbData_1.2.10
    ## [76] cli_3.6.1               fansi_1.0.4             viridisLite_0.4.1
    ## [79] gtable_0.3.3            digest_0.6.31           htmlwidgets_1.6.2
    ## [82] memoise_2.0.1           htmltools_0.5.5         lifecycle_1.0.3
    ## [85] httr_1.4.5              bit64_4.0.5