# Comparison with base R

## Stefano Mangiola

### 2023-03-15

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.2.2 (2022-10-31 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.8.1 SummarizedExperiment_1.28.0
    ##  [3] Biobase_2.58.0                 GenomicRanges_1.50.2
    ##  [5] GenomeInfoDb_1.34.9            IRanges_2.32.0
    ##  [7] S4Vectors_0.36.2               BiocGenerics_0.44.0
    ##  [9] MatrixGenerics_1.10.0          matrixStats_0.63.0
    ## [11] tidybulk_1.10.1                ggrepel_0.9.3
    ## [13] ggplot2_3.4.1                  magrittr_2.0.3
    ## [15] tibble_3.2.0                   tidyr_1.3.0
    ## [17] dplyr_1.1.0                    knitr_1.42
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] nlme_3.1-162           bitops_1.0-7           bit64_4.0.5
    ##  [4] httr_1.4.5             SnowballC_0.7.0        tools_4.2.2
    ##  [7] backports_1.4.1        utf8_1.2.3             R6_2.5.1
    ## [10] DBI_1.1.3              lazyeval_0.2.2         mgcv_1.8-42
    ## [13] colorspace_2.1-0       withr_2.5.0            tidyselect_1.2.0
    ## [16] bit_4.0.5              compiler_4.2.2         preprocessCore_1.60.2
    ## [19] cli_3.6.0              DelayedArray_0.24.0    plotly_4.10.1
    ## [22] scales_1.2.1           readr_2.1.4            genefilter_1.80.3
    ## [25] stringr_1.5.0          digest_0.6.31          XVector_0.38.0
    ## [28] pkgconfig_2.0.3        htmltools_0.5.4        fastmap_1.1.1
    ## [31] limma_3.54.2           htmlwidgets_1.6.1      rlang_1.1.0
    ## [34] RSQLite_2.3.0          generics_0.1.3         jsonlite_1.8.4
    ## [37] BiocParallel_1.32.5    tokenizers_0.3.0       RCurl_1.98-1.10
    ## [40] GenomeInfoDbData_1.2.9 Matrix_1.5-3           Rcpp_1.0.10
    ## [43] munsell_0.5.0          fansi_1.0.4            lifecycle_1.0.3
    ## [46] stringi_1.7.12         edgeR_3.40.2           zlibbioc_1.44.0
    ## [49] plyr_1.8.8             Rtsne_0.16             grid_4.2.2
    ## [52] blob_1.2.3             parallel_4.2.2         crayon_1.5.2
    ## [55] lattice_0.20-45        Biostrings_2.66.0      splines_4.2.2
    ## [58] annotate_1.76.0        hms_1.1.2              KEGGREST_1.38.0
    ## [61] locfit_1.5-9.7         pillar_1.8.1           widyr_0.1.5
    ## [64] reshape2_1.4.4         codetools_0.2-19       XML_3.99-0.13
    ## [67] glue_1.6.2             evaluate_0.20          tidytext_0.4.1
    ## [70] data.table_1.14.8      vctrs_0.5.2            png_0.1-8
    ## [73] tzdb_0.3.0             gtable_0.3.1           purrr_1.0.1
    ## [76] cachem_1.0.7           xfun_0.37              xtable_1.8-4
    ## [79] broom_1.0.4            janeaustenr_1.0.0      survival_3.5-5
    ## [82] viridisLite_0.4.1      AnnotationDbi_1.60.2   memoise_2.0.1
    ## [85] sva_3.46.0             ellipsis_0.3.2