R Programming Guide

Getting Started with R

# we start by assigning the number 10 to a variable called 'a'
# why do we use '<-' instead of '='? honestly it's just old-school R habit, 
# and it means we're stuffing the value 10 into 'a' so we can reuse it later without retyping numbers
a <- 10

# now we use print() to actually show the value on the screen
# why? because in R scripts, if you just write 'a', it won't always show up in the console, 
# so printing forces it to display
print(a)

Essential Data Structures

# here we are creating a vector using the combine function c()
# why use c()? because R needs a way to glue multiple individual numbers together into a single list object
numbers <- c(1, 2, 3, 4, 5)

# print out the vector to check it worked correctly
print(numbers)

Data Wrangling with dplyr

# load the dplyr package so we can actually use data manipulation tools
# why? because base R can be super clunky for filtering, so dplyr makes it way easier
library(dplyr)

# we are using the built-in 'starwars' dataset for practice
# what's with the '%>%' symbol? that's a pipe, it just takes whatever data is on the left and shoves it into the next step
starwars %>% 
  filter(species == "Droid") %>%     # keep only rows where the species is specifically a Droid
  select(name, height, mass) %>%     # pick just these 3 columns because we don't care about the rest
  arrange(desc(height))              # sort them from tallest to shortest so the biggest is at the top

Writing Custom Functions

# write our own custom function to calculate Body Mass Index (BMI)
# why make a function? if we have to calculate this for 50 different patients, writing the math out every single time is a nightmare
calculate_bmi <- function(weight_kg, height_meters) {
  bmi <- weight_kg / (height_meters ^ 2) # the actual math formula
  return(round(bmi, 2))                  # round it nicely to 2 decimal places so it looks clean
}

# test out our new function with a weight of 75kg and height of 1.8m
patient_bmi <- calculate_bmi(weight_kg = 75, height_meters = 1.8)
print(patient_bmi)

Hypothesis Testing (T-tests)

# perform an independent t-test using R's built-in sleep dataset
# why 'extra ~ group'? the tilde (~) basically means "test extra sleep hours grouped by the drug groups"
t.test(extra ~ group, data = sleep)

Advanced Visualisation: Volcano Plots

# before we start we need to load the packages we need for this, as they dont come with R default.
# why ggplot2? because it's basically the industry standard for making clean graphs in R
# why readxl? because typing out raw data by hand is way too much effort, so we just read the excel file directly
library(ggplot2)
library(readxl)

# grab the excel file and save it as 'df' (dataframe)
# why save as df? because writing 'OX1.xlsx' every single time is super annoying
df <- read_excel("OX1.xlsx")

# force p-value and log2 fold change to be treated strictly as numbers
# why? because if R accidentally reads them as text strings, math and plotting will completely break
df$pvalue <- as.numeric(df$pvalue)
df$log2FoldChange <- as.numeric(df$log2FoldChange)

# multiply the log2 fold change by -1 to flip the direction
# why? because the original comparison was backwards (WT vs ERFOX1 instead of ERFOX1 vs WT), 
# so multiplying by -1 fixes it so positive values actually mean higher expression in our gene of interest
df$log2FoldChange <- df$log2FoldChange * -1

# convert the p-value into -log10 scale
# why do this? because raw p-values are tiny decimals (like 0.0001), 
# and taking -log10 turns them into big positive numbers so the most significant genes shoot right up to the top
df$neg_log10_p <- -log10(df$pvalue)

# make a new column called expression status and set everything to 'Not Significant' by default
df$Expression_Status <- "Not Significant"

# now apply our cutoffs to change the status of specific genes
# why > 1 and < 0.05? anything above a fold change of 1 (a 2-fold change) 
# and a p-value under 0.05 is what we actually care about as upregulated or downregulated
df$Expression_Status[df$log2FoldChange > 1 & df$pvalue < 0.05] <- "Upregulated"
df$Expression_Status[df$log2FoldChange < -1 & df$pvalue < 0.05] <- "Downregulated"

# finally, build the volcano plot using ggplot
ggplot(df, aes(x = log2FoldChange, y = neg_log10_p, color = Expression_Status)) +
  
  # add the genes as little points on the plot
  # alpha makes them slightly see-through so overlapping dots are visible
  geom_point(alpha = 0.6, size = 1.5) +

  # set custom colors so the significant groups pop out clearly against the grey background
  scale_color_manual(values = c(
    "Upregulated" = "#D55E00",
    "Downregulated" = "#0072B2",
    "Not Significant" = "gray75"
  )) +

  # add vertical dashed lines at -1 and 1 to clearly show our fold change cutoffs
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray40", linewidth = 0.4) +

  # add a horizontal dashed line for our p-value cutoff (0.05 on a -log10 scale)
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray40", linewidth = 0.4) +

  # add clean titles and axis labels to explain what we are looking at
  labs(
    title = "ERFOX1 vs WT",
    x = "Log2 Fold Change",
    y = "-log10(p-value)",
    color = "Gene Status"
  ) +

  # use a classic theme to keep the background clean and simple
  theme_classic(base_size = 14) +

  # fine-tune layout details like making the plot square and bolding the title
  theme(
    aspect.ratio = 1,
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    axis.line = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

The R Ecosystem & Helpful Packages

Base R is cool and all, but the real magic happens when you start pulling in community packages. Here are the main tools that will save you a ton of headache: