Practice Session & Exercise

Introduction to R and Basic Programming Concepts

Author

Affiliation

Bioinformatics Core Facility CECAD

1 Practice by Sessions

1.1 Session 1

1.1.1 Exercise 1:

• Use R to calculate the following:
  • multiply 31 and 78
  • divide 697 with 41
• Assign the value of 39 to x
• Assign the value of 22 to y

Solution 1.

31 * 78

[1] 2418

697 / 41

[1] 17

x <- 39
x

[1] 39

y <- 22
y

[1] 22

1.2 Session 2

Before You Start : Some Raw Material

Please copy the below code:

Assume you observed the expression of six genes under two conditions:

C1 <- c(g1=0,  g2=1,  g3=10,  g4=6,  g5=4,  g6=23)
C2 <- c(g1=20, g2=11, g3=1,   g4=7,  g5=3,  g6=5)

In addition you have a grouping variable for the genes:

G  <- c(g1="A",g2="A",g3="B", g4="C",g5="C",g6="B")

1.2.1 Exercise 1: Create a Matrix

Use the cbind() function (“column-wise bind”) to bind C1 and C2 into a matrix (if you are unsure, check ?cbind).

Assign the result to a variable named “M” and print it.

Can you get the number of rows and columns of M using R code?  

Solution 2.

M <- cbind(C1,C2) # note that the variable names are now column names!
M

   C1 C2
g1  0 20
g2  1 11
g3 10  1
g4  6  7
g5  4  3
g6 23  5

nrow(M) # the number of rows

[1] 6

ncol(M) # the number of columns

[1] 2

1.2.2 Exercise 2: Data Type of a Matrix

What is the data type of M? Do you have an idea why actually two types are listed?  

Solution 3.

class(M)

[1] "matrix" "array" 

# M is both a matrix and also an array with a single slice!

1.2.3 Exercise 3: Inspect the Matrix

Extract the first column by index and by name. What is its data type? 
Use the is.vector() function to confirm that it is a vector. 

Extract the second row by index and by name. What is its data type? 
Use the is.vector() function to confirm that it is a vector, too. 

Solution 4.

class(M[,1])

[1] "numeric"

is.vector(M[,1])

[1] TRUE

class(M[2,])

[1] "numeric"

is.vector(M[2,])

[1] TRUE

1.2.4 Exercise 4a: Row Means by Hand

The operators “+” and “/” are vectorizing, that is, they can operate on entire vectors. 

Get the mean expression of the 6 genes in the two conditions by adding C1 and C2, and then dividing the result by 2. 

Assign the result to a variable named Cmean. 

Hint: Put parentheses where needed, to ensure the correct order of computations!  

Solution 5.

C1C2_mean <- (C1+C2) / 2 
C1C2_mean # note: just putting the variable name is often equivalent to using print()

  g1   g2   g3   g4   g5   g6 
10.0  6.0  5.5  6.5  3.5 14.0 

1.2.5 Exercise 4b: Row (and Column) Means Computationally

R has a function rowMeans(), which does exactly the computation of Exercise 4a. Try it and compare to C1C2_mean of Exercise 4a, using function identical().

What do you think the pendant colMeans() does? Try it.

Solution 6.

# Observe that each index in the vectors C1 and C2 correspond to a ROW in M
C1C2_mean_rowMeans <- rowMeans(M)
identical(C1C2_mean, C1C2_mean_rowMeans)

[1] TRUE

columnMeans <- colMeans(M)
columnMeans 

      C1       C2 
7.333333 7.833333 

1.2.6 Exercise 5: Type Conversion (Not Always Intended)

Try to cbind() all 3 vectors C1, C2, and G. 
It will work, but something goes wrong. Do you observe something unexpected?  

Solution 7.

cbind(C1,C2,G)

   C1   C2   G  
g1 "0"  "20" "A"
g2 "1"  "11" "A"
g3 "10" "1"  "B"
g4 "6"  "7"  "C"
g5 "4"  "3"  "C"
g6 "23" "5"  "B"

# A matrix is a 2D generalization of a vector, and like in a vector, all its
# elements must be of the same data type. If we try to mix truly 
# incompatible data types, we get an error. 
# But by an internal conversion rule numbers can be converted to text.
# Hence R, trying to be obedient and to make our matrix, converts all
# numbers to character values! 
# Such things can be hard to spot if they happen unintentionally.

1.2.7 Exercise 6: Make a List

The list() function combines heterogeneous R objects into a common list object. 

The names of the resulting list elements can be specified like so: list(my_name=my_object).  

Use the function to create a list of C1, C2, and G,  
with names “condition1”, “condition2”, and “group”.`  

Assign the list object to a variable named “L” and print it.

Solution 8.

L <- list(condition1 = C1,
          condition2 = C2,
          group = G) 
# note that you can break a long command line anywhere,
# except in the middle of a word

L

$condition1
g1 g2 g3 g4 g5 g6 
 0  1 10  6  4 23 

$condition2
g1 g2 g3 g4 g5 g6 
20 11  1  7  3  5 

$group
 g1  g2  g3  g4  g5  g6 
"A" "A" "B" "C" "C" "B" 

# All 3 objects have been "hooked onto" L, keeping their 
# entire, original form

1.2.8 Exercise 7: Make a Dataframe

A list is convenient because of its permissiveness: you can store (nearly) anything there. But it can also be cumbersome to handle. 

List L has a special property: all of its elements are vectors of a fixed length (although of different data types).  
Such lists can be transformed such that they look very similar to matrices, yet internally stay lists with all their benefits.

This matrix-like data structure is the dataframe (data type data.frame).  
Convert L to a dataframe using as.data.frame(), assign the result to a variable named df, and print it.  

Solution 9.

df <- data.frame(L)
df

   condition1 condition2 group
g1          0         20     A
g2          1         11     A
g3         10          1     B
g4          6          7     C
g5          4          3     C
g6         23          5     B

1.2.9 Exercise 8a: Inspect the Dataframe

A dataframe can be indexed like a matrix.  
Extract the first column by index or by name, and use class() and is.vector() on it. 
Then do the same with the second row. Do both dimensions behave in the same way?

Solution 10.

df[,1]

[1]  0  1 10  6  4 23

class(df[,1])

[1] "numeric"

is.vector(df[,1])

[1] TRUE

# equivalent formulations:
# is.vector(df[,"condition1"])
# is.vector(df$condition1)

df[2,]

   condition1 condition2 group
g2          1         11     A

class(df[2,])

[1] "data.frame"

is.vector(df[2,])

[1] FALSE

1.2.10 Exercise/ Question 8b: Why Can Dataframe Rows Not Be Vectors?

Solution 11. The dataframe structure is meant for cases where a tabular format is convenient for analysis and readability, 
yet the columns are not all of the same data type, and hence a matrix representation is not possible. 
Therefore the definition of a dataframe must not make assumptions which require all columns (= all entries of each row) to have the same type.
But exactly this assumption is made by the vector data type.

1.2.11 Exercise 9: Make a Factor

Vector G (which is also the “group” column of df) is a grouping variable for the ratio of the count vectors C1 and C2:  
C2 is higher than C1 in group “A” of G, lower than C1 in group “B”, and about the same in group “C”. 

We can make it explicit that it is a grouping variable by coding it as a factor. 

In a much larger, real dataset, such a grouping factor for genes could be used to statistically relate gene induction or repression to some other attribute of the genes, for example to pathway membership.  
Here we simply practice to set up a factor variable:

Create an unordered factor with base level “C” from vector G and assign it to variable f1. 
Print f1.  
Create an ordered factor with levels c(“C”, “B”, “A”) from G and assign it to variable f2. 
Print f2.  
Create an unordered factor with base level “C” from vector df$group, and assign it back to df$group. 
Column “group” should now be a factor. Check it using class().

Solution 12.

f1 <- factor(G, levels = c("C","A","B"))
f1

g1 g2 g3 g4 g5 g6 
 A  A  B  C  C  B 
Levels: C A B

f2 <- factor(G, levels = c("C","B","A"))
f2

g1 g2 g3 g4 g5 g6 
 A  A  B  C  C  B 
Levels: C B A

df$group <- factor(G, levels = c("C","A","B"))
class(df$group)

[1] "factor"

1.3 Session 3

1.3.1 Exercise 1: R Control Flow (if, if...else)

You are analyzing the expression level of a target cancer gene ((p53)). You need to write a script that checks a scalar expression value and prints a diagnostic classification.

Task: Create a variable named p53_expr. Write an if...else if...else structure that:

1. Prints “Critical Overexpression” if p53_expr is strictly greater than 10.
2. Prints “Normal Expression” if p53_expr is between 3 and 10 (inclusive).
3. Prints “Underexpressed” if p53_expr is less than 3.

Bonus Challenge: Test your code by converting the structure into a vectorized function using ifelse() to handle an entire vector of expressions: c(12, 5, 1, 8, 15).

Solution 13.

p53_expr <- 8
if (p53_expr > 10) {
  print("Critical Overexpression")
} else if (p53_expr >= 3) {
  print("Normal Expression")
} else {
  print("Underexpressed")
}

[1] "Normal Expression"

# Bonus Vectorized Check
p53_vec <- c(12, 5, 1, 8, 15)
ifelse(p53_vec > 10, "Critical Overexpression", 
       ifelse(p53_vec >= 3, "Normal Expression", "Underexpressed"))

[1] "Critical Overexpression" "Normal Expression"      
[3] "Underexpressed"          "Normal Expression"      
[5] "Critical Overexpression"

1.3.2 Exercise 2: R Loops (for, next, break)

You have a list of DNA primer sequences. You need to screen them sequentially. However, your laboratory pipeline must skip short sequences and completely halt if an invalid or contaminated sequence is detected.

Task: Write a for loop to iterate through primers and evaluate the text length using nchar():

primers <- c("ATGGC", "AT", "GTCGATC", "UNKNOWN", "CTGA","AYTGC")

1. If a sequence matches unknown character, print “Contamination detected! Halting pipeline.” and exit the loop immediately.
2. If a sequence has fewer than 4 characters, skip it , print the message: "Skipping short primer: [sequence_here]" and continue to the next iteration.
3. For all valid sequences, print the message: "Processing primer: [sequence_here]".

Solution 14.

primers <- c("ATGGC", "AT", "GTCGATC", "UNKNOWN", "CTGA","AYTGC")
for (p in primers) {
  if (p == "UNKNOWN") {
    print("Contamination detected! Halting pipeline.")
    break
  }
  if (nchar(p) < 4) {
    print(paste("Skipping short primer:", p))
    next
  }
  print(paste("Processing primer:", p))
}

[1] "Processing primer: ATGGC"
[1] "Skipping short primer: AT"
[1] "Processing primer: GTCGATC"
[1] "Contamination detected! Halting pipeline."

1.3.3 Exercise 3: R Loops (while or repeat)

A bacterial culture starts with a population of 150 cells and doubles every hour. You need to determine exactly how many hours it takes for the colony to breach a critical biohazard threshold of 20,000 cells.

Task: Write a for loop to iterate through primers and evaluate the text length using nchar():

1. Initialize two variables: population <- 150 and hours <- 0.
2. Write a while loop that continues running as long as population is less than 20000.
3. Inside the loop, increment the hours counter by 1 and double the population.
4. After the loop terminates, print a final summary sentence showing the total hours and ending population.

Bonus Challenge: Modify your loop to include a safety check that halts the process if it exceeds 24 hours, printing “Process halted: Time limit exceeded.” and breaking the loop. Plot the population growth over time using a simple line graph.

Solution 15.

population <- 150
hours <- 0
while (population < 20000) {
  hours <- hours + 1
  population <- population * 2
  
  if (hours > 24) {
    print("Process halted: Time limit exceeded.")
    break
  }
}

print(paste("Threshold breached at hour", hours, "with", population, "cells."))

[1] "Threshold breached at hour 8 with 38400 cells."

# Bonus: Plotting population growth
population <- 150
hours <- 0

pop_growth <- data.frame(Hours = hours, Population = population)
while (population < 20000) {
  hours <- hours + 1
  population <- population * 2
  pop_growth <- rbind(pop_growth, data.frame(Hours = hours, Population = population))
  if (hours > 24) {
    print("Process halted: Time limit exceeded.")
    break
  }
}

plot(pop_growth$Hours, pop_growth$Population, type = "o", col = "blue",
     xlab = "Hours", ylab = "Population Size",
     main = "Bacterial Population Growth Over Time")

1.3.4 Exercise 4: apply() family function

You receive a 2D dataset tracking the concentrations of 4 different metabolites (columns) across 6 laboratory tissue samples (rows). You must extract the baseline variability for each metabolite.

Task: Copy this code to generate your matrix:

set.seed(42)
metabolite_matrix <- matrix(rnorm(24, mean = 5, sd = 2), nrow = 6)
colnames(metabolite_matrix) <- c("Glucose", "Lactate", "Pyruvate", "Alanine")

1. Use a single apply() function call to compute the median concentration for each row (per tissue sample).
2. Use a second apply() function call to compute the standard deviation (sd) for each column (per metabolite type).

Solution 16.

set.seed(42)
metabolite_matrix <- matrix(rnorm(24, mean = 5, sd = 2), nrow = 6)
colnames(metabolite_matrix) <- c("Glucose", "Lactate", "Pyruvate", "Alanine")

# Row medians (Margin = 1)
apply(metabolite_matrix, 1, median)

[1] 4.982098 4.626552 5.229807 5.570149 5.232351 6.108550

# Column standard deviations (Margin = 2)
apply(metabolite_matrix, 2, sd)

 Glucose  Lactate Pyruvate  Alanine 
1.317636 2.043521 2.324981 3.057981 

1.4 Session 4

1.4.1 Exercise 1: CSV I/O & the factor gotcha

Estimated time: 20 minutes

A collaborator sends you a small gene expression table. The decimals are written with a comma and columns are separated by semicolons. You want to round-trip the file (write it, read it back) and make sure nothing changes silently — in particular, that character columns do not get converted to factors.

Task:

expr <- data.frame(
  gene      = c("TP53", "BRCA1", "MYC", "EGFR", "KRAS"),
  condition = c("ctrl", "ctrl", "treat", "treat", "ctrl"),
  log2FC    = c(-0.42, 1.87, 2.31, -1.12, 0.05)
)

1. Write expr to a file expr_de.csv using write.csv2() (German-style: ; separator, , decimal). Inspect the raw file content with readLines("expr_de.csv").
2. Read the file back with read.csv2() into an object expr_back. Check class(expr_back$gene) — is it character or factor?
3. Re-read the file passing stringsAsFactors = FALSE (or convert with as.character()) so that gene and condition stay as character vectors.
4. Compare the two data frames with identical(expr, expr_back) and all.equal(expr, expr_back). If they differ, explain why.

Bonus Challenge: Write the same data with the default write.csv() and read it back with read.csv(). Which separator and decimal symbol does each pair use? When would the German variant matter in practice?

Solution 17.

expr <- data.frame(
  gene      = c("TP53", "BRCA1", "MYC", "EGFR", "KRAS"),
  condition = c("ctrl", "ctrl", "treat", "treat", "ctrl"),
  log2FC    = c(-0.42, 1.87, 2.31, -1.12, 0.05)
)

# 1. Write German-style CSV
write.csv2(expr, "expr_de.csv", row.names = FALSE)
readLines("expr_de.csv")

[1] "\"gene\";\"condition\";\"log2FC\"" "\"TP53\";\"ctrl\";-0,42"          
[3] "\"BRCA1\";\"ctrl\";1,87"           "\"MYC\";\"treat\";2,31"           
[5] "\"EGFR\";\"treat\";-1,12"          "\"KRAS\";\"ctrl\";0,05"           

# 2. Read back and check class
expr_back <- read.csv2("expr_de.csv")
class(expr_back$gene)

[1] "character"

# 3. Force character columns
expr_back <- read.csv2("expr_de.csv", stringsAsFactors = FALSE)
class(expr_back$gene)

[1] "character"

# 4. Compare
identical(expr, expr_back)

[1] TRUE

all.equal(expr, expr_back)

[1] TRUE

# Bonus: default CSV
write.csv(expr, "expr_us.csv", row.names = FALSE)
readLines("expr_us.csv")

[1] "\"gene\",\"condition\",\"log2FC\"" "\"TP53\",\"ctrl\",-0.42"          
[3] "\"BRCA1\",\"ctrl\",1.87"           "\"MYC\",\"treat\",2.31"           
[5] "\"EGFR\",\"treat\",-1.12"          "\"KRAS\",\"ctrl\",0.05"           

expr_us <- read.csv("expr_us.csv", stringsAsFactors = FALSE)
all.equal(expr, expr_us)

[1] TRUE

1.4.2 Exercise 2: Binary formats — saveRDS vs save

Estimated time: 15 minutes

After a long overnight analysis you want to checkpoint your results so a colleague can reload them tomorrow without rerunning anything. R offers two binary formats: saveRDS() / readRDS() (one object, assigned on read) and save() / load() (multiple objects, restored under their original names).

Task:

deg_results <- data.frame(
  gene   = paste0("g", 1:6),
  log2FC = c(1.2, -0.8, 2.5, -1.7, 0.3, 0.0),
  padj   = c(0.001, 0.04, 0.0001, 0.02, 0.6, 0.9)
)
sample_info <- data.frame(
  sample_id = paste0("s", 1:4),
  group     = c("ctrl", "ctrl", "treat", "treat")
)

1. Save deg_results to deg.rds with saveRDS(). Reload it under a different name, deg2 <- readRDS("deg.rds"). Verify with identical(deg_results, deg2).
2. Save both deg_results and sample_info together to session.RData with save(). Remove them from your workspace with rm(deg_results, sample_info) and confirm with ls().
3. Reload them with load("session.RData"). Notice they reappear under their original names — you do not assign the result of load().
4. Explain in one sentence: when would you prefer .rds over .RData?

Solution 18.

deg_results <- data.frame(
  gene   = paste0("g", 1:6),
  log2FC = c(1.2, -0.8, 2.5, -1.7, 0.3, 0.0),
  padj   = c(0.001, 0.04, 0.0001, 0.02, 0.6, 0.9)
)
sample_info <- data.frame(
  sample_id = paste0("s", 1:4),
  group     = c("ctrl", "ctrl", "treat", "treat")
)

# 1. Single object
saveRDS(deg_results, "deg.rds")
deg2 <- readRDS("deg.rds")
identical(deg_results, deg2)

[1] TRUE

# 2. Multiple objects
save(deg_results, sample_info, file = "session.RData")
rm(deg_results, sample_info)
ls()

 [1] "C1"                 "C1C2_mean"          "C1C2_mean_rowMeans"
 [4] "C2"                 "colorize"           "columnMeans"       
 [7] "deg2"               "df"                 "expr"              
[10] "expr_back"          "expr_us"            "f1"                
[13] "f2"                 "G"                  "hours"             
[16] "L"                  "M"                  "metabolite_matrix" 
[19] "p"                  "p53_expr"           "p53_vec"           
[22] "pop_growth"         "population"         "primers"           
[25] "x"                  "y"                 

# 3. Reload — original names restored
load("session.RData")
ls()

 [1] "C1"                 "C1C2_mean"          "C1C2_mean_rowMeans"
 [4] "C2"                 "colorize"           "columnMeans"       
 [7] "deg_results"        "deg2"               "df"                
[10] "expr"               "expr_back"          "expr_us"           
[13] "f1"                 "f2"                 "G"                 
[16] "hours"              "L"                  "M"                 
[19] "metabolite_matrix"  "p"                  "p53_expr"          
[22] "p53_vec"            "pop_growth"         "population"        
[25] "primers"            "sample_info"        "x"                 
[28] "y"                 

# 4. Prefer .rds when you want to control the variable name on read
#    (e.g. inside a function or a pipeline). Use .RData to checkpoint
#    several related objects at once.

1.4.3 Exercise 3: Reshaping a time-course experiment

Estimated time: 25 minutes

You have measured the expression of three genes at four time points. Your raw table is in wide format (one column per time point), but the plotting function you want to use later expects long format (one row per measurement).

Task:

wide <- data.frame(
  gene = c("TP53", "BRCA1", "MYC"),
  t0   = c(1.0, 0.8, 1.2),
  t1   = c(1.4, 0.9, 2.1),
  t2   = c(2.1, 1.1, 3.5),
  t3   = c(3.0, 1.3, 5.0)
)

1. Reshape wide to long format with reshape() using gene as idvar, the four time columns as varying, v.names = "expression", timevar = "time", and direction = "long". Store it in long.
2. Inspect long — how many rows does it have? Why?
3. Reshape long back to wide format and store it in wide_back.
4. Compare with all.equal(wide, wide_back). If attributes differ, try all.equal(wide, wide_back, check.attributes = FALSE).

Bonus Challenge: Use stack() on a subset of the wide table (e.g. wide[, -1]) and compare the result with what reshape() produced. When is stack() simpler? When does it fall short?

Solution 19.

wide <- data.frame(
  gene = c("TP53", "BRCA1", "MYC"),
  t0   = c(1.0, 0.8, 1.2),
  t1   = c(1.4, 0.9, 2.1),
  t2   = c(2.1, 1.1, 3.5),
  t3   = c(3.0, 1.3, 5.0)
)

# 1. Wide -> long
long <- reshape(wide,
                varying  = c("t0", "t1", "t2", "t3"),
                v.names  = "expression",
                timevar  = "time",
                times    = c("t0", "t1", "t2", "t3"),
                idvar    = "gene",
                direction = "long")
long

          gene time expression
TP53.t0   TP53   t0        1.0
BRCA1.t0 BRCA1   t0        0.8
MYC.t0     MYC   t0        1.2
TP53.t1   TP53   t1        1.4
BRCA1.t1 BRCA1   t1        0.9
MYC.t1     MYC   t1        2.1
TP53.t2   TP53   t2        2.1
BRCA1.t2 BRCA1   t2        1.1
MYC.t2     MYC   t2        3.5
TP53.t3   TP53   t3        3.0
BRCA1.t3 BRCA1   t3        1.3
MYC.t3     MYC   t3        5.0

# 2. 3 genes x 4 time points = 12 rows
nrow(long)

[1] 12

# 3. Long -> wide
wide_back <- reshape(long,
                     idvar    = "gene",
                     timevar  = "time",
                     direction = "wide")
wide_back

          gene expression.t0 expression.t1 expression.t2 expression.t3
TP53.t0   TP53           1.0           1.4           2.1           3.0
BRCA1.t0 BRCA1           0.8           0.9           1.1           1.3
MYC.t0     MYC           1.2           2.1           3.5           5.0

# 4. Compare (column names will differ: t0 vs expression.t0)
all.equal(wide, wide_back, check.attributes = FALSE)

[1] TRUE

# Bonus
stack(wide[, -1])

   values ind
1     1.0  t0
2     0.8  t0
3     1.2  t0
4     1.4  t1
5     0.9  t1
6     2.1  t1
7     2.1  t2
8     1.1  t2
9     3.5  t2
10    3.0  t3
11    1.3  t3
12    5.0  t3

# stack() is simpler when you only need a single value/group pair,
# but it loses the gene identifier — reshape() preserves it via idvar.

1.4.4 Exercise 4: Combining tables — merge, cbind, rbind

Estimated time: 25 minutes

You have two tables from the same experiment: gene-level counts and sample metadata. They share a sample_id column. You want to attach metadata to counts so that each row of counts knows which group it belongs to.

Task:

counts <- data.frame(
  sample_id = c("s1", "s2", "s3", "s4"),
  TP53      = c(120, 95, 300, 280),
  MYC       = c(45,  60, 410, 390)
)
meta <- data.frame(
  sample_id = c("s1", "s2", "s3", "s4"),
  group     = c("ctrl", "ctrl", "treat", "treat"),
  batch     = c(1, 1, 2, 2)
)

1. Merge counts and meta by sample_id using merge(). Store the result in full. How many rows does it have?
2. Now duplicate the first row of meta (meta_bad <- rbind(meta, meta[1, ])) and merge again. How many rows do you get now? Explain the Cartesian product.
3. Use cbind() to add a new column total <- counts$TP53 + counts$MYC to full. Why does cbind() work safely here but would fail if you tried to bind counts and meta directly?
4. Build a second batch of two samples (s5, s6) with the same columns as counts, then rbind() it onto counts. Verify the result has 6 rows.

Solution 20.

counts <- data.frame(
  sample_id = c("s1", "s2", "s3", "s4"),
  TP53      = c(120, 95, 300, 280),
  MYC       = c(45,  60, 410, 390)
)
meta <- data.frame(
  sample_id = c("s1", "s2", "s3", "s4"),
  group     = c("ctrl", "ctrl", "treat", "treat"),
  batch     = c(1, 1, 2, 2)
)

# 1. Clean merge
full <- merge(counts, meta, by = "sample_id")
nrow(full)   # 4

[1] 4

# 2. Cartesian blow-up from a duplicated key
meta_bad <- rbind(meta, meta[1, ])
full_bad <- merge(counts, meta_bad, by = "sample_id")
nrow(full_bad)   # 5 — s1 appears twice on the right side

[1] 5

# 3. cbind a derived column
full$total <- full$TP53 + full$MYC
full

  sample_id TP53 MYC group batch total
1        s1  120  45  ctrl     1   165
2        s2   95  60  ctrl     1   155
3        s3  300 410 treat     2   710
4        s4  280 390 treat     2   670

# cbind on counts + meta would just glue columns positionally,
# duplicating sample_id and *assuming* the row order matches.
# merge() uses the key, which is safer.

# 4. rbind a new batch
new_batch <- data.frame(
  sample_id = c("s5", "s6"),
  TP53      = c(150, 200),
  MYC       = c(80,  110)
)
counts2 <- rbind(counts, new_batch)
nrow(counts2)   # 6

[1] 6

1.5 Session 5

Before You Start

Please re-create species_tables, species_means, species_colors, and trait_colors as given at the beginning of Session5.R. Try to understand each step of the computations!

1.5.1 Exercise 1:

We want to study the distances between members of the versicolor and the setosa population, and compare it to the distances between members of the versicolor and the virginia population. Given the assumed genotypes (virginia 4n, setosa 2n, and the hybrid versicolor 4n+2n = 6n), versicolor and virginica should be closer to each other than versicolor and setosa, because they share more alleles. Here we do it for one flower trait at a time (assigned to variable “feature”).

Take the difference between the vectors species_tables$versicolor[[feature]] species_tables$setosa[[feature]] and assign it to variable “d_vs”.

Take the difference between the vectors species_tables$versicolor[[feature]] species_tables$virginica[[feature]] and assign it to variable “d_vv”.

We want to produce a plot which shows histograms of both difference vectors on the same canvas. See lines 347-388 in Session5.R for hints.

Solution 21.

species_tables <- 
    split(iris[,-5],  ## what to split
          iris[, 5])  ## split by what


feature <- "Sepal.Length"
d_vs <- 
  abs(species_tables$versicolor[[feature]] -
  species_tables$setosa[[feature]])
d_vv <- 
  abs(species_tables$versicolor[[feature]] -
  species_tables$virginica[[feature]])

# total range of both vectors:
use_range <- c(min(c(d_vs,d_vv)),
               max(c(d_vs,d_vv)))

hist(d_vs,
     col=rgb(1,0,0,1/4),
     
     add=FALSE, ## this is the default
     
     ## initialize to full x range !
     xlim=use_range,
     
     ## full y range you usually 
     ## only know after some trials ..
     ylim=c(0,22),
     
     ## x-axis label
     xlab="abs. trait diff. between random pairs of plants from the two species",
     
     ## larger axis labels:
     cex.lab = 2,
     
         main = 
       paste("Distance histograms for ",
           feature, 
             "\nVersicolor/Setosa (red) and ",                    "Versicolor/Virginica (blue)"),
     
     ## larger title:
     cex.main = 2
)


## ----Histogram_added2---------------------------------------------------------

hist(d_vv,
     add=TRUE,
     col=rgb(0,0,1,1/4)
)

1.5.2 Exercise 2: Label Extreme Points

Copy the code of Slide 17/Session5.html / Session5.R/lines 216-234.

Try the following: (1) Identify the 3 highest y values. (2) find all points which reach these values (these can be more than the 3 unique highest values, if multiple points have the same value). (3) mark all points found with a fat red point.

Hints: (1) We need function sort() for sorting all values. The sort should be decreasing, such that the highest values are on top. Check ?sort to learn how a decreasing sort is selected. Then run unique on the sorted values, to remove repeated values. Finally assign the top 3 values of the sorted list to variable “top3”. (2) We saw the operator %in%, which is used as "x %in% y" and finds those points in vector x which are in vector y. Can you transfer this pattern to the task at hand? (3) Use the points() function to color the points red. Check ?points to learn about parameters for choosing symbol, color and expansion factor. Try different values of the cex parameter.

Solution 22.

df <- subset(iris, Species=="setosa")

plot(
# x is now the row number in df
x=1:nrow(df),
xlab="individual plant",
y=df$Petal.Width,
ylab="Petal.Width",
## show both points and
## connecting lines:
type="b",
## line width:
lwd = 2,
## line style = dashed:
lty=2,
main="Iris setosa"
)
sort(df$Petal.Width, decreasing = TRUE)[1:3]

[1] 0.6 0.5 0.4

highest <- sort(df$Petal.Width, decreasing = TRUE)[1:3]
i <- which(df$Petal.Width %in% highest)
i

[1]  6 16 17 22 24 27 32 44 45

points(i, df$Petal.Width[i], pch=21, col="red", bg="red", cex=3)

1.6 Session 6

1.6.1 Exercise 1: Measures of Central Tendency (mean, median, summary)

You are measuring the concentration of a signaling protein (pg/mL) in 10 wild-type mice. However, due to a pipetting error or an extreme biological anomaly, one sample has an massive outlier value.

Task

protein_levels <- c(22, 24, 21, 25, 23, 22, 26, 24, 22, 450)

1. Calculate the arithmetic mean and the median of the raw data.
2. Recalculate the mean and median excluding the outlier value of 450 (Hint: you can index it out or use protein_levels[protein_levels < 100]).
3. Look at how much each metric shifted. Which measure of central tendency proved robust against the extreme outlier

Solution 23.

protein_levels <- c(22, 24, 21, 25, 23, 22, 26, 24, 22, 450)
# 1. With Outlier
mean(protein_levels)   # Output: 65.7

[1] 65.9

median(protein_levels) # Output: 23.5

[1] 23.5

# 2. Without Outlier
clean_protein <- protein_levels[protein_levels < 100]
mean(clean_protein)   # Output: 23.11

[1] 23.22222

median(clean_protein) # Output: 23

[1] 23

Interpretation

• The median barely moved (23.5 to 23), making it highly robust.

• The mean collapsed from 65.7 down to 23.11 because it was heavily skewed by 450.

1.6.2 Exercise 2: Measures of Variability (sd, var, range, IQR)

An ecologist measures the root depth of a desert plant species in meters. For a publication layout, they need to convert all measurements to centimeters (multiply by 100).

Task

root_m <- c(1.2, 1.5, 0.8, 2.1, 1.4, 1.7)

1. Calculate the standard deviation (sd()) and the variance (var()) of root_m.
2. Create a new vector called root_cm by multiplying root_m by 100.
3. Calculate the standard deviation and variance of root_cm.
4. Observe the mathematical changes: By what scaling factor did the standard deviation increase? By what factor did the variance increase?

Solution 24.

root_m <- c(1.2, 1.5, 0.8, 2.1, 1.4, 1.7)

# 1. Original in meters
sd_m <- sd(root_m)
var_m <- var(root_m) 
sd_m

[1] 0.441588

var_m

[1] 0.195

# 2. Convert to centimeters
root_cm <- root_m * 100

# 3. Standard deviation and variance in centimeters
sd_cm <- sd(root_cm)   
var_cm <- var(root_cm) 
sd_cm

[1] 44.1588

var_cm

[1] 1950

# 4. Scaling factors
sd_factor <- sd_cm / sd_m   # Output: 100
var_factor <- var_cm / var_m # Output: 10000
sd_factor

[1] 100

var_factor

[1] 10000

Interpretation

• Standard deviation scaled perfectly by 100 (same as the data).
• Variance scaled by 10,000 (which is 100^2), because variance uses squared units.

1.6.3 Exercise 3: Contingency Tables (table, prop.table, addmargins)

You are breeding fruit flies (Drosophila) to study classical genetics. You cross two flies and count the offspring phenotypes, categorizing them by eye color (Red vs. White) and wing type (Normal vs. Vestigial).

Task

flies <- data.frame(
  eyes = c("Red", "Red", "White", "Red", "Red", "White", "Red", "White"),
  wings = c("Normal", "Normal", "Vestigial", "Vestigial", "Normal", "Normal", "Normal", "Vestigial")
)

1. Generate a 2D contingency table of counts using fly_table <- table(flies$eyes, flies$wings). Print it out.
2. Convert this count table into conditional relative frequencies based on the row totals (Hint: look at the margin argument inside prop.table()). What percentage of Red-eyed flies have Normal wings?
3. Use addmargins(fly_table) to quickly add row and column sums to your original count table.

Solution 25.

flies <- data.frame(
  eyes = c("Red", "Red", "White", "Red", "Red", "White", "Red", "White"),
  wings = c("Normal", "Normal", "Vestigial", "Vestigial", "Normal", "Normal", "Normal", "Vestigial")
)

# 1. Count Table
fly_table <- table(flies$eyes, flies$wings)
print(fly_table)

       
        Normal Vestigial
  Red        4         1
  White      1         2

# 2. Row proportions (margin = 1 targets the rows: Red vs White)
prop.table(fly_table, margin = 1) 

       
           Normal Vestigial
  Red   0.8000000 0.2000000
  White 0.3333333 0.6666667

# Interpretation: Output shows 0.80 for Red-Normal, meaning 80% of Red flies have Normal wings.

# 3. Add Margin Totals
addmargins(fly_table)

       
        Normal Vestigial Sum
  Red        4         1   5
  White      1         2   3
  Sum        5         3   8

1.6.4 Exercise 4: Correlation Coefficients (cor, cor.test)

A limnologist tracks water dissolved oxygen levels (mg/L) and water temperature (°C) in an alpine lake over several weeks to find if there is a stable relationship.

Task

temp <- c(4.2, 5.5, 7.1, 9.3, 12.1, 14.5, 16.2)
oxygen <- c(12.5, 11.8, 11.2, 10.1, 9.4, 8.7, 8.1)

1. Use cor(temp, oxygen) to calculate the basic Pearson correlation coefficient. Is the linear relationship positive or negative?
2. Use cor.test(temp, oxygen) to extract deeper inferential parameters. 3 Read the console output from step 2: What are the exact lower and upper bounds of the 95 percent confidence interval for this correlation coefficient?

Solution 26.

temp <- c(4.2, 5.5, 7.1, 9.3, 12.1, 14.5, 16.2)
oxygen <- c(12.5, 11.8, 11.2, 10.1, 9.4, 8.7, 8.1)

# 1. Simple vector correlation
cor(temp, oxygen) # Output: -0.9935147 (Strong negative linear correlation)

[1] -0.993692

# 2 & 3. Full testing
cor.test(temp, oxygen)

    Pearson's product-moment correlation

data:  temp and oxygen
t = -19.813, df = 5, p-value = 6.049e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.999109 -0.956064
sample estimates:
      cor 
-0.993692 

Focus

• Look at the bottom of the output text block under ‘95 percent confidence interval:’ Output bounds: -0.9991272 to -0.9529061

2 Practice across Sessions

2.1 Exercise 1: Environmental Monitoring Pipeline

You are an environmental data scientist processing seasonal weather anomalies. You need to read raw files, reshape and clean data, apply programmatic thresholds, isolate metrics using the apply family, generate analytical base R plots, and automatically compile everything into a professional HTML report.

For this project, you will use the airquality dataset, which is built directly into R. It tracks daily air quality measurements in New York, including solar radiation, wind, temperature, and ozone levels—factors that heavily influence biological ecosystems.

Step 1: Data IO & Reshaping

• Save the built-in R dataset to a temporary local CSV file on your hard drive using write.csv().
• Read that exact file back into a new R variable named raw_env using read.csv().
• Extract only rows matching the summer months of July and August into a filtered data frame called summer_env.

Step 2: Loop Execution & Custom Control Flow

• Write a for loop that iterates sequentially through every row of your summer_env data frame to check the Temp (Temperature) value for that day using an if...else if...else structure:
  • If Temp is strictly greater than 90°F, print: “Extreme Heat Detected ([Temp] degrees)”.
  • If Temp is between 80°F and 90°F (inclusive), print: “Moderate Warning”.
  • For all other temperatures, skip processing entirely using the next statement without printing anything.

Step 3: High-Throughput Metrics via the apply Family

• Use apply family function to find the min, maximum and standard deviation for Temp from the summer_env

Step 4: Multi-Panel Visualizations with Base R Plots

• Configure your R graphics engine using par(mfrow = c(1, 2)) to arrange your upcoming plots side-by-side in a 1-row, 2-column grid.
• Plot 1: Draw a shaded hist() histogram of summer_env$Temp. Color it "coral", add clear axis text, and label the title "Summer Temperature Profile".
• Plot 2: Draw a plot() scatter plot showing Ozone on the y-axis and Solar.R (Solar Radiation) on the x-axis. Use solid circular indicators (pch = 16), color the dots "darkblue", and add a line of best fit through the plot utilizing abline(lm(Ozone ~ Solar.R, data = summer_env), col = "red", lwd = 2).
• Reset your plotting layout back to standard single-panel parameters using par(mfrow = c(1, 1)).

Step 5: Wrap it into an Automated HTML Report

To practice modern literate programming as covered in our very first discussion, compile your steps into an automatically generated standalone HTML report document using Quarto (or R Markdown).

1. In RStudio, create a new document: File > New File > Quarto Document… (Select HTML output).
2. ipe out the default placeholder text and paste the complete solution script into clean execution code chunks.
3. Add professional Markdown headers (## Step 1…) between chunks to describe your analysis steps.
4. Click the Render (or Knit) button at the top of your screen to run the entire pipeline and generate your final report.html file automatically.

Solution 27.

## --- STEP 1: DATA IO & RESHAPING ---
# Write built-in data out to file
write.csv(airquality, "temp_airquality.csv", row.names = FALSE)

# Read it back in
raw_env <- read.csv("temp_airquality.csv")


# Reshape/Filter for summer months (July = 7, August = 8)
summer_env <- raw_env[raw_env$Month == 7 | raw_env$Month == 8, ]



## --- STEP 2: LOOPS & CONTROL FLOW ---
for (i in 1:nrow(summer_env)) {
  current_temp <- summer_env$Temp[i]
  
  if (current_temp > 90) {
    print(paste("Row", i, ": Extreme Heat Detected (", current_temp, "degrees)"))
  } else if (current_temp >= 80 & current_temp <= 90) {
    print(paste("Row", i, ": Moderate Warning"))
  } else {
    next # Skip everything else safely
  }
}

[1] "Row 1 : Moderate Warning"
[1] "Row 2 : Moderate Warning"
[1] "Row 3 : Moderate Warning"
[1] "Row 4 : Moderate Warning"
[1] "Row 5 : Moderate Warning"
[1] "Row 6 : Moderate Warning"
[1] "Row 7 : Moderate Warning"
[1] "Row 8 : Extreme Heat Detected ( 92 degrees)"
[1] "Row 9 : Extreme Heat Detected ( 92 degrees)"
[1] "Row 10 : Moderate Warning"
[1] "Row 11 : Moderate Warning"
[1] "Row 13 : Moderate Warning"
[1] "Row 14 : Extreme Heat Detected ( 91 degrees)"
[1] "Row 15 : Moderate Warning"
[1] "Row 16 : Moderate Warning"
[1] "Row 17 : Moderate Warning"
[1] "Row 18 : Moderate Warning"
[1] "Row 19 : Moderate Warning"
[1] "Row 20 : Moderate Warning"
[1] "Row 22 : Moderate Warning"
[1] "Row 23 : Moderate Warning"
[1] "Row 24 : Moderate Warning"
[1] "Row 25 : Moderate Warning"
[1] "Row 26 : Moderate Warning"
[1] "Row 27 : Moderate Warning"
[1] "Row 28 : Moderate Warning"
[1] "Row 29 : Moderate Warning"
[1] "Row 30 : Moderate Warning"
[1] "Row 31 : Moderate Warning"
[1] "Row 32 : Moderate Warning"
[1] "Row 33 : Moderate Warning"
[1] "Row 34 : Moderate Warning"
[1] "Row 35 : Moderate Warning"
[1] "Row 36 : Moderate Warning"
[1] "Row 37 : Moderate Warning"
[1] "Row 38 : Moderate Warning"
[1] "Row 39 : Moderate Warning"
[1] "Row 40 : Moderate Warning"
[1] "Row 41 : Extreme Heat Detected ( 92 degrees)"
[1] "Row 42 : Moderate Warning"
[1] "Row 43 : Moderate Warning"
[1] "Row 44 : Moderate Warning"
[1] "Row 45 : Moderate Warning"
[1] "Row 56 : Moderate Warning"
[1] "Row 57 : Moderate Warning"
[1] "Row 58 : Moderate Warning"
[1] "Row 59 : Extreme Heat Detected ( 97 degrees)"
[1] "Row 60 : Extreme Heat Detected ( 94 degrees)"
[1] "Row 61 : Extreme Heat Detected ( 96 degrees)"
[1] "Row 62 : Extreme Heat Detected ( 94 degrees)"

# --- STEP 3: THE APPLY FAMILY ---
tapply(summer_env$Temp, summer_env$Month, min)

 7  8 
73 72 

tapply(summer_env$Temp, summer_env$Month, max)

 7  8 
92 97 

tapply(summer_env$Temp, summer_env$Month, mean)

       7        8 
83.90323 83.96774 

tapply(summer_env$Temp, summer_env$Month, sd)

       7        8 
4.315513 6.585256 

# --- STEP 4: BASE R PLOTS ---
# Setup side-by-side canvas panels
par(mfrow = c(1, 2))

# Plot 1: Temperature Frequency Distribution
hist(summer_env$Temp, 
     col = "coral", 
     border = "white",
     main = "Summer Temperature Profile", 
     xlab = "Temperature (°F)")

# Plot 2: Scatter plot tracking Solar influence on Ozone levels
plot(summer_env$Solar.R, summer_env$Ozone,
     pch = 16, 
     col = "darkblue",
     main = "Ozone vs. Solar Radiation",
     xlab = "Solar Radiation (langleys)", 
     ylab = "Ozone (ppb)")

# Overlay best fit linear regression path
abline(lm(Ozone ~ Solar.R, data = summer_env), col = "red", lwd = 2)

# Reset graphic canvas parameters
par(mfrow = c(1, 1))

2.2 Exercise 2: Mini differential-expression pipeline

Estimated time: 75 minutes (Phases A–C: ~45 min · Phase D + Bonus: ~30 min)

This exercise ties together everything from Sessions 1–6. You will save and reload a small gene-expression dataset, classify genes by their fold change (control flow), summarise the groups two different ways (loop vs functional style), reshape the per-sample measurements, merge in sample metadata, order by biological importance, round-trip the result through CSV and RDS to see how the two formats differ, and finally peek ahead at visualisation.

Setup — generate the two data frames in memory:

set.seed(2026)
genes <- paste0("g", 1:12)

# Per-gene baseline expression and treatment effect.
# Effects are a mix of up-regulation (+), down-regulation (-) and flat (~0)
# so the classification step later produces all three categories.
baseline <- round(rnorm(12, mean = 5, sd = 1), 2)
effect   <- c( 3.0,  2.5, -2.2, -1.6,  0.1, -0.2,
               2.8, -1.9,  0.3,  3.5, -2.0,  0.0)
noise <- function() round(rnorm(12, 0, 0.3), 2)

expr_wide <- data.frame(
  gene = genes,
  s1   = baseline           + noise(),  # ctrl
  s2   = baseline           + noise(),  # ctrl
  s3   = baseline + effect  + noise(),  # treat
  s4   = baseline + effect  + noise()   # treat
)

meta <- data.frame(
  sample_id = c("s1", "s2", "s3", "s4"),
  group     = c("ctrl", "ctrl", "treat", "treat")
)

The two objects expr_wide and meta now live in your R session. From here on, you are responsible for getting them onto disk and back.

2.2.1 Phase A — Read & Write

1. Save expr_wide and meta to disk as expr_wide.csv and meta.csv. Read them back into new objects expr and meta_in. (Since R 4.0, read.csv() keeps strings as character by default — no extra argument needed, but you may pass stringsAsFactors = FALSE for clarity.)
2. Compute a per-gene log2 fold change between the treat samples (s3, s4) and the ctrl samples (s1, s2). Hint: log2(mean_treat / mean_ctrl). Store it as a new column log2FC in expr. Then run hist(expr$log2FC) — does the distribution suggest where the ±0.5 threshold should sit?

2.2.2 Phase B — Classify & Summarise

3. Classify each gene using ifelse() into one of three categories, stored in expr$direction:
   • "up" if log2FC > 0.5
   • "down" if log2FC < -0.5
   • "flat" otherwise
   Then run barplot(table(expr$direction)) to verify you have genes in all three groups.
4. Summarise the groups two different ways so you can compare imperative vs functional R:
   1. Write a for loop that, for each direction in c("up", "down", "flat"), prints the number of genes and the mean log2FC.
   2. Reproduce the same summary in a single call using tapply() (or sapply() + subsetting). Confirm the numbers match.

2.2.3 Phase C — Reshape, Merge, Order

5. Reshape the four sample columns of expr from wide to long with reshape(), producing a data frame long with columns gene, sample_id, expression.
6. Merge long with meta_in by sample_id so each row also knows its group. Then merge in log2FC and direction from expr so every row carries its gene-level annotation. Run boxplot(expression ~ group, data = long) as a sanity check — can you already see the treatment effect?
7. Order the merged table by abs(log2FC) descending (most strongly affected genes first), then by group. This is the canonical ordering in differential-expression output.

2.2.4 Phase D — Save, Round-trip, Visualise

8. Save the final ordered table both as results.csv (text) and results.rds (binary). Read both back and compare each reload to the original using identical() and all.equal(). Which format round-trips exactly? Where does the other one differ (row order? row names? numeric precision? column types?).
9. Peek at Session 5: use stripchart(expression ~ factor(direction, levels = c("down", "flat", "up")), data = long[long$group == "treat", ], vertical = TRUE, method = "jitter", pch = 19, col = c("tomato", "grey60", "steelblue"), main = "Treated samples: expression by direction") to see the individual data points for each gene category. Do the up and down genes separate cleanly from the flat ones?

Bonus Challenge: Use aggregate() to compute the standard deviation of expression per (gene, group) pair on the long table. Which gene is the most variable inside the treat group? Does its direction (up/down/flat) explain that variability, or is something else going on?

Solution 28.

# ---- Setup (same as in the task) ----
set.seed(2026)
genes <- paste0("g", 1:12)
baseline <- round(rnorm(12, mean = 5, sd = 1), 2)
effect   <- c( 3.0,  2.5, -2.2, -1.6,  0.1, -0.2,
               2.8, -1.9,  0.3,  3.5, -2.0,  0.0)
noise <- function() round(rnorm(12, 0, 0.3), 2)
expr_wide <- data.frame(
  gene = genes,
  s1 = baseline + noise(), s2 = baseline + noise(),
  s3 = baseline + effect + noise(), s4 = baseline + effect + noise()
)
meta <- data.frame(
  sample_id = c("s1", "s2", "s3", "s4"),
  group     = c("ctrl", "ctrl", "treat", "treat")
)

# ---- Phase A: Read & Write ----
# 1. Save and reload
write.csv(expr_wide, "expr_wide.csv", row.names = FALSE)
write.csv(meta,      "meta.csv",      row.names = FALSE)
expr    <- read.csv("expr_wide.csv")
meta_in <- read.csv("meta.csv")

# 2. log2 fold change (treat vs ctrl)
ctrl_mean  <- rowMeans(expr[, c("s1", "s2")])
treat_mean <- rowMeans(expr[, c("s3", "s4")])
expr$log2FC <- round(log2(treat_mean / ctrl_mean), 3)
hist(expr$log2FC, breaks = 8, main = "Distribution of log2FC",
     xlab = "log2 Fold Change", col = "lightblue")
abline(v = c(-0.5, 0.5), col = "red", lty = 2)  # show the thresholds

# ---- Phase B: Classify & Summarise ----
# 3. Classify
expr$direction <- ifelse(expr$log2FC >  0.5, "up",
                  ifelse(expr$log2FC < -0.5, "down", "flat"))
barplot(table(expr$direction),
        main = "Genes per direction", col = c("tomato", "grey70", "steelblue"))

# 4a. for-loop summary (imperative)
for (d in c("up", "down", "flat")) {
  sub <- expr[expr$direction == d, ]
  cat(d, ":", nrow(sub), "genes, mean log2FC =",
      round(mean(sub$log2FC), 3), "\n")
}

up : 4 genes, mean log2FC = 0.712 
down : 4 genes, mean log2FC = -0.797 
flat : 4 genes, mean log2FC = -0.121 

# 4b. Same summary in one shot (functional)
tapply(expr$log2FC, expr$direction, length)   # n genes

down flat   up 
   4    4    4 

tapply(expr$log2FC, expr$direction, mean)     # mean log2FC

    down     flat       up 
-0.79700 -0.12125  0.71175 

# ---- Phase C: Reshape, Merge, Order ----
# 5. Wide -> long
long <- reshape(expr[, c("gene", "s1", "s2", "s3", "s4")],
                varying  = c("s1", "s2", "s3", "s4"),
                v.names  = "expression",
                timevar  = "sample_id",
                times    = c("s1", "s2", "s3", "s4"),
                idvar    = "gene",
                direction = "long")
rownames(long) <- NULL

# 6. Merge with sample metadata and gene-level annotation
long <- merge(long, meta_in, by = "sample_id")
long <- merge(long, expr[, c("gene", "log2FC", "direction")], by = "gene")
boxplot(expression ~ group, data = long,
        main = "Sanity check: expression by group",
        col = c("steelblue", "tomato"))  # treat group should be higher

# 7. Order by |log2FC| desc, then group — most affected genes first
long <- long[order(-abs(long$log2FC), long$group), ]
head(long)

   gene sample_id expression group log2FC direction
41   g8        s1       4.19  ctrl -0.976      down
42   g8        s2       4.03  ctrl -0.976      down
43   g8        s3       2.17 treat -0.976      down
44   g8        s4       2.01 treat -0.976      down
5   g10        s1       4.48  ctrl  0.830        up
6   g10        s2       4.63  ctrl  0.830        up

# ---- Phase D: Save, Round-trip, Visualise ----
# 8. Compare CSV vs RDS round-trip
write.csv(long, "results.csv", row.names = FALSE)
saveRDS(long,   "results.rds")

long_csv <- read.csv("results.csv")
long_rds <- readRDS("results.rds")

identical(long, long_rds)   # TRUE — binary preserves everything

[1] TRUE

identical(long, long_csv)   # FALSE — text loses row names / row order info

[1] FALSE

all.equal(long, long_csv)   # describes the differences

[1] "Attributes: < Component \"row.names\": Mean relative difference: 0.6639004 >"

# Typical CSV-vs-original differences:
#  - row.names are character "1","2",... instead of integer
#  - row order is preserved here, but numeric precision can shift on
#    values with many decimals
# Lesson: use RDS when you need an exact reload; use CSV for portability.

# 9. Session 5 preview: individual points per direction in treated samples
stripchart(expression ~ factor(direction, levels = c("down", "flat", "up")),
           data     = long[long$group == "treat", ],
           vertical = TRUE,        # put direction on the x-axis
           method   = "jitter", pch = 19,
           col      = c("tomato", "grey60", "steelblue"),
           main     = "Treated samples: expression by direction",
           xlab     = "Direction", ylab = "Expression")

# up and down genes should sit clearly apart from flat ones

# ---- Bonus: per-(gene, group) variability ----
var_tab <- aggregate(expression ~ gene + group, data = long, FUN = sd)
var_tab[order(-var_tab$expression), ]

   gene group  expression
6    g3  ctrl 0.905096680
14  g10 treat 0.509116882
5    g2  ctrl 0.452548340
18   g3 treat 0.367695526
3   g11  ctrl 0.332340187
9    g6  ctrl 0.318198052
17   g2 treat 0.282842712
7    g4  ctrl 0.268700577
19   g4 treat 0.134350288
16  g12 treat 0.127279221
15  g11 treat 0.120208153
22   g7 treat 0.120208153
11   g8  ctrl 0.113137085
23   g8 treat 0.113137085
2   g10  ctrl 0.106066017
24   g9 treat 0.084852814
21   g6 treat 0.070710678
4   g12  ctrl 0.063639610
1    g1  ctrl 0.056568542
10   g7  ctrl 0.056568542
20   g5 treat 0.049497475
13   g1 treat 0.042426407
12   g9  ctrl 0.014142136
8    g5  ctrl 0.007071068

# The most variable gene in `treat` is usually one with a strong effect
# (large |log2FC|), because the +effect amplifies sample-to-sample noise.