Foundations R Fundamentals

R Fundamentals for Bioinformatics

When to Use

• Reading or adapting R code from a published pipeline, DESeq2/edgeR/Seurat vignette, or supplementary script
• Translating a Python data-analysis snippet into R (or vice versa) and need the syntax mapping
• Debugging an off-by-one or "wrong subset" bug caused by R's 1-based, exclude-on-negative indexing
• Setting up a minimal Bioconductor workflow (install, DESeqDataSet, filtering results)
• Choosing/running a basic statistical test in R (t-test, Wilcoxon, chi-squared, correlation)

Version Compatibility

R ≥ 4.2, Bioconductor ≥ 3.17 (BiocManager ≥ 1.30), DESeq2 ≥ 1.40, dplyr ≥ 1.1, ggplot2 ≥ 3.4. Base-R syntax here (<-, data.frame, apply) is stable across R versions back to R 3.x.

Prerequisites

• R installed (R --version), or an R kernel (IRkernel) for Jupyter
• install.packages("BiocManager") for any Bioconductor package (DESeq2, GenomicRanges, ...)
• Comfort with Python's data types helps — this skill is written as a Python→R translation

Python vs R Syntax Traps

Goal: avoid the bugs that come from assuming R behaves like Python. Approach: memorize this table before writing/reading any R code — every row is a real, common bug source when porting between the two languages.

Vectors and Data Frames

Goal: manipulate vectors and data frames the R way. Approach: use vectorized operations and logical indexing instead of loops; use data.frame for tabular data and filter with boolean masks on columns.

# --- Vectors ---
gene_expression <- c(2.5, 3.1, 4.2, 1.8, 5.6)
gene_names <- c("BRCA1", "TP53", "EGFR", "MYC", "KRAS")
positions <- 1:10
coverage <- seq(from = 0, to = 100, by = 10)
groups <- rep(c("control", "treatment"), each = 3)

# Named vector + logical indexing
gc_content <- c(BRCA1 = 0.42, TP53 = 0.38, EGFR = 0.55, MYC = 0.61)
gc_rich <- gc_content[gc_content > 0.5]

# Vectorized ops -- no explicit loops needed
log2fc <- log2(gene_expression + 1)          # pseudocount avoids log(0)
pvalue <- c(0.01, 0.15, 0.001, 0.02, 0.4)
sig_up <- log2fc > 1 & pvalue < 0.05          # vectorized filter
cat("Upregulated:", gene_names[sig_up], "\n")
which(sig_up)                                  # indices of TRUE values

# --- Data frames ---
gene_data <- data.frame(
    gene = c("BRCA1", "TP53", "EGFR"),
    log2fc = c(1.2, -0.5, 3.8),
    pvalue = c(0.01, 0.15, 0.001),
    stringsAsFactors = FALSE   # keep strings as character, not factor
)

gene_data$gene                  # column access
gene_data[1, ]                  # first row
gene_data[2, "log2fc"]          # specific cell
gene_data[gene_data$log2fc > 1 & gene_data$pvalue < 0.05, ]  # filter rows

gene_data$padj <- p.adjust(gene_data$pvalue, method = "BH")   # FDR (standard for genomics)
gene_data[order(gene_data$pvalue), ]                            # sort by p-value

Matrices and Statistical Distributions

Goal: work with expression matrices and distribution functions. Approach: use matrix() with dimnames for gene-by-sample data and apply() for row/column stats; use the d/p/q/r prefix family for any distribution.

# --- Matrices (expression matrices) ---
expr_matrix <- matrix(
    c(5.2, 3.1, 8.5, 6.2, 4.8, 2.9, 7.1, 5.8),
    nrow = 2, byrow = TRUE,
    dimnames = list(c("Sample1", "Sample2"), c("BRCA1", "TP53", "EGFR", "MYC"))
)

expr_matrix[2, "EGFR"]          # specific cell
expr_matrix[1, ]                # full row
expr_matrix[, "BRCA1"]          # full column
apply(expr_matrix, 1, mean)     # row means (margin = 1)
apply(expr_matrix, 2, sd)       # column SDs (margin = 2)

# --- Statistical distributions: d(ensity) p(CDF) q(uantile) r(andom) ---
set.seed(42)                     # reproducibility
rnorm(n = 100, mean = 10, sd = 2)    # random samples
pnorm(q = 1.96, lower.tail = TRUE)   # CDF -> 0.975
qnorm(p = 0.975)                     # quantile -> 1.96
rnbinom(n = 100, mu = 50, size = 5)  # negative binomial -- what DESeq2 models RNA-seq counts with

# t-test / Wilcoxon / chi-squared -- the basic hypothesis-test toolkit
control <- rnorm(30, mean = 10, sd = 2)
treatment <- rnorm(30, mean = 12, sd = 2)
t.test(control, treatment)$p.value
wilcox.test(control, treatment)$p.value      # non-parametric alternative
chisq.test(table(c("Mut","WT","Mut","Mut"), c("Ctrl","Ctrl","Drug","Drug")))

Minimal DESeq2 Workflow

Goal: run a minimal DESeq2 differential-expression analysis. Approach: build a DESeqDataSet from a raw count matrix + sample metadata, run DESeq(), extract a contrast, then filter by FDR and effect size.

# Install once
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("DESeq2")

library(DESeq2)

#' Run DESeq2 and return significant genes for a two-group contrast.
#' @param counts   integer matrix, genes x samples (raw counts, NOT normalized)
#' @param metadata data.frame, one row per sample, must include `condition`
#' @param alpha    FDR threshold (default 0.05)
#' @param lfc      minimum absolute log2 fold-change to call significant
#' @return data.frame of significant genes sorted by adjusted p-value
run_deseq2 <- function(counts, metadata, alpha = 0.05, lfc = 1) {
    dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~condition)
    dds <- DESeq(dds)
    res <- results(dds, contrast = c("condition", "treatment", "control"))
    res_df <- as.data.frame(res)
    sig <- res_df[!is.na(res_df$padj) & res_df$padj < alpha & abs(res_df$log2FoldChange) > lfc, ]
    sig[order(sig$padj), ]
}

# Reading/writing tabular data
df <- read.csv("data.csv", stringsAsFactors = FALSE)
df <- read.delim("data.tsv", sep = "\t")
write.csv(df, "output.csv", row.names = FALSE)

# Self-check: run_deseq2() end to end on synthetic count data
demo_deseq2 <- function() {
    set.seed(1)
    counts <- matrix(rnbinom(400, mu = 100, size = 5), nrow = 100, ncol = 4,
                      dimnames = list(paste0("gene", 1:100), paste0("s", 1:4)))
    metadata <- data.frame(condition = factor(c("control", "control", "treatment", "treatment")),
                            row.names = colnames(counts))
    sig <- run_deseq2(counts, metadata)
    stopifnot(is.data.frame(sig))
    stopifnot(all(c("log2FoldChange", "padj") %in% colnames(sig)))
    cat("OK:", nrow(sig), "significant genes out of", nrow(counts), "\n")
}
demo_deseq2()

Pitfalls

• R indexing starts at 1: x[1] is first; x[-1] means "all except first" (NOT last element like Python)
• Factors look like strings but are not: a factor stores integer codes + labels; functions expecting strings will misbehave — use stringsAsFactors = FALSE in data.frame() or as.character() to convert
• library() vs require(): use library() in scripts — it errors loudly if a package is missing; require() returns FALSE silently and can hide broken pipelines
• Data frames vs matrices: Bioconductor functions often require numeric matrices; a data frame with character columns cannot be used directly — convert with as.matrix(df[, numeric_cols])
• <- inside function arguments: f(x <- 1) assigns x in the enclosing scope AND passes the value as the argument; use = for keyword arguments to avoid the side effect
• Multiple testing: p.adjust(..., method = "BH") (Benjamini-Hochberg FDR) is the genomics standard; Bonferroni is overly conservative for large gene sets
• Scalar vs vectorized AND/OR: && and || only evaluate the first element — always use & and | when filtering vectors/data frames

See Also

• foundations-r-hypothesis-testing-and-nonparametrics — Shapiro-Wilk, Wilcoxon, ANOVA, correlation tests in depth
• foundations-r-regression-correlation-and-diagnostics — linear models, lm(), diagnostic plots
• foundations-biostatistics-fundamentals — statistical concepts underlying these tests
• bio-differential-expression-deseq2-basics — full DESeq2 workflow beyond this quick-start