📦 R Data Manipulation & Management

🔁 Review: Lecture 5

Before new structures, ten minutes on last week’s ground: syntax, variables, types, math, strings, operators, control flow.

⚡ Rapid Fire 1 — Syntax & Variables

1. You type print("hello) and press Enter. What happens?
2. How do you store the value 150 in a variable called x?
3. Two of these names crash R. One is legal but dangerous. One is perfectly safe. Which is which?

• 1gene
• cell count
• gene_A
• Gene_A

✅ Rapid Fire 1 — Syntax & Variables

1. Nothing crashes and nothing runs. The console sits on a + prompt, waiting for the rest of the string — it stays stuck until every quote and bracket matches
2. x <- 150
3. 1gene and cell count → syntax errors, R stops. Gene_A → legal but dangerous: a different variable from gene_A. gene_A → safe

⚡ Rapid Fire 2 — Types, Math & Strings

1. What type is 150? What type is "EGFR"?
2. 10 %% 3 versus 10 %/% 3?
3. read.csv("data\visium\counts.csv") crashes. Why, and what is the rule?

✅ Rapid Fire 2 — Types, Math & Strings

1. 150 is numeric. "EGFR" is character — the quotes make it text, so it can never be counted or averaged
2. 1 and 3: %% returns the remainder, %/% floor-divides
3. \ is the escape character, so \v and \c are invalid codes. The rule: always use forward slashes / in R paths

⚡ Rapid Fire 3 — Operators & Control Flow

1. Your if tested one condition and it came back FALSE. How do you check a second condition?
2. Write the R test for 500 <= x <= 2500.
3. In a for loop, how do next and break differ?

✅ Rapid Fire 3 — Operators & Control Flow

1. Add an else if — it runs only when the condition above it was FALSE, and the first TRUE branch wins
2. Two comparisons joined: (x >= 500) & (x <= 2500). R has no chained comparisons
3. next skips the rest of the current iteration and keeps going; break aborts the loop entirely

🎯 Lecture 5, Closed

• The Basics — comments, <-, naming
• Values — types, NA / NaN / Inf
• Numbers — %%, %/%, log2, CPM
• Text — quotes, cat(), escapes
• Logic — ==, &, |, !, no chaining
• Flow — if / else, for, next, break, while

Everything from here builds on top of this.

🗺️ Where We Are in the Course

📂 Case Study: GSE129221

📂 The Design

📂 The Question

Which genes change when apatinib hits resistant cells — and do the same genes change in the sensitive ones?

• To analyze, we need: a count matrix (39,376 genes × 12 samples)

📂 Before We Download

Let’s make a folder for the data.

mkdir -p counts      # create it once
cd counts

📂 curl: Fetch the Counts Matrix

curl -o counts.tsv.gz \
  "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE129221&format=file&file=GSE129221_raw_counts_GRCh38.p13_NCBI.tsv.gz"

gunzip counts.tsv.gz     # .tsv.gz -> .tsv

📂 Why a Script?

That was two lines — and one was a 136-character URL.

📂 download.sh

#!/usr/bin/env bash
set -euo pipefail

OUT_DIR="./counts"
mkdir -p "$OUT_DIR"

curl -o "$OUT_DIR/counts.tsv.gz" \
  "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE129221&format=file&file=GSE129221_raw_counts_GRCh38.p13_NCBI.tsv.gz"

gunzip -f "$OUT_DIR/counts.tsv.gz"

echo "File downloaded at $OUT_DIR/counts.tsv:"
ls -lh "$OUT_DIR"

echo "Head of the file:"
head -n 5 "$OUT_DIR/counts.tsv"

📂 Making the script executable

chmod +x download.sh    # RUN ONLY ONCE!
./download.sh           # Script is executable now!

📂 counts/README.md

# counts/

Raw RNA-seq data for **[GSE129221](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129221)** (Song et al. 2019, Cancer Cell Int 19:117).

| File | Source | Notes |
|---|---|---|
| `counts.tsv` | GEO (via download.sh) | 39,376 genes x 12 samples, raw counts |

Run `./download.sh` to fetch it

🧬 So… What Now?

You have the count matrix. It is not yet useful information:

• Rows are Entrez IDs — 100287102 tells you nothing
• Columns are GSM accessions — PC9 and PC9GR look identical
• Many genes are zero in every sample

📋 1. R Data Structures

The Family Tree

Vectors — Same Type, One Dimension

A vector is a 1D chain of values that all share the same type.

Vectors — Same Type, One Dimension

A plain sequence:

and a custom one

Or, repetition:

Subsetting Vectors

Counting Starts at 1

Matrices — Same Type, Two Dimensions

A matrix is a 2D grid where every cell holds the same type. Spot metadata matrices are exactly this shape.

Naming Rows and Columns

Arrays — Matrices Stacked

Lists — A Box of Anything

Elements can have completely different types and lengths

Data Frames — The Workhorse

A data frame is a 2D table where each column keeps its own type. This is your spot-metadata table.

Getting Oriented in a Data Frame

Subsetting Data Frames

Factors — Categories, With an Order

A factor stores category labels plus the list of allowed levels — the R way to hold annotation columns.

Tibbles — Modern Data Frames

A tibble (tbl_df) is an enhanced modern data frame from the tidyverse.

• Better printing: Previews only the first 10 rows and columns that fit your screen.
• Safer subsetting: Never silently drops down to a vector when you select one column with [ , 1].
• Strict typing: Never converts strings to factors automatically.

Creating a Tibble

✍️ Exercise 6.1: What Kind of Object Should This Table Be?

Four genes, three samples:

Two data types live here: gene IDs (text) and read counts (numbers).

• 1. Vector
• 2. Matrix
• 3. Array
• 4. List
• 5. Data frame

✅ Exercise 6.1: What Kind of Object Should This Table Be?

e) Data frame — the only structure that is both rectangular and lets each column keep its own type.

counts <- data.frame(
  gene_id = c("gene_A", "gene_B", "gene_C", "gene_D"),
  S1 = c(2, 0, 310, 89),
  S2 = c(3, 45, 85, 12),
  S3 = c(1, 12, 0, 640),
  row.names = 1
)

✍️ Exercise 6.2: What Shape Is counts?

counts <- data.frame(
  gene_id = c("gene_A", "gene_B", "gene_C", "gene_D"),
  S1 = c(2, 0, 310, 89),
  S2 = c(3, 45, 85, 12),
  S3 = c(1, 12, 0, 640),
  row.names = 1
)

What does dim(counts) give you?

• 1. 3 4
• 2. 4 3
• 3. 1 3
• 4. 4 4

✅ Exercise 6.2: What Shape Is counts?

b) 4 3

dim() is always rows, then columns.

• Rows = genes, running down the page
• Columns = samples, running across the top

dim(counts)              # 4 3
head(rownames(counts))   # "gene_A" "gene_B" "gene_C" ...

✍️ Exercise 6.3: Pull One Sample

You want the read counts for S1 only, as a vector.

• 1. counts[1, ]
• 2. counts$S1
• 3. counts[S1]
• 4. counts["S1"]

✅ Exercise 6.3: Pull One Sample

b) counts$S1 — $ pulls the column out by name and hands you the bare values.

counts$S1        # 2 0 310 89 — all 4 genes, that one sample

📋 2. Reading, Writing & Merging Files

Where Am I?

getwd()                      # print the current working directory
dir()                        # list what is inside it
setwd("~/projects/spatial")  # move to a project folder

Paths Are Relative to Where You Stand

One path, two folders. Walk into a subdirectory and the same line flips its answer:

f <- "FCWT_2026/lecture-06/counts/counts.tsv"

From the repo root:

getwd()          # Spatial_Transcriptomics_Course
file.exists(f)   # TRUE

Now step one level deeper — nothing about f changes:

setwd("FCWT_2026/lecture-06")
getwd()          # Spatial_Transcriptomics_Course/FCWT_2026/lecture-06
file.exists(f)   # FALSE

The Fix: Anchor to the Project Root

here() builds the path from the project root, so the answer stops depending on where you stand:

library(here)
f <- here("FCWT_2026", "lecture-06", "counts", "counts.tsv")

file.exists(f)   # TRUE — from the repo root
file.exists(f)   # TRUE — from analysis/code/ too

Read: the Real Count Matrix

The file you downloaded in the first half of this lecture. It is tab-separated, so read.delim() is the right call:

library(here)
f <- here("FCWT_2026", "lecture-06", "counts", "counts.tsv")

counts <- read.delim(f, row.names = 1)

dim(counts)     # 39376 12
head(counts)    # first 6 genes, all 12 samples

Base R Table Readers Compared

All read.* functions wrap read.table(), but differ in default parameters:

Write: Save a Slice to Disk

Export a subset to disk with write.csv():

small <- head(counts, 5)

# Export with row labels (gene IDs) preserved:
write.csv(small, "counts_head.csv", row.names = TRUE)

The Problem: Opaque Sample IDs

Look at the column headers of your count matrix:

colnames(counts)
# [1] "GSM3702795" "GSM3702796" "GSM3702797" "GSM3702798" ...

• GEO accessions (GSM...) tell you nothing about the biology.
• Which samples are untreated controls? Which received Apatinib? Which are drug-resistant?
• We need sample metadata mapped to these accession IDs.

Create & Save a Sample Metadata Table

Build the lookup mapping and save it to disk with write.csv():

sample_info <- data.frame(
  sample_id = colnames(counts),
  cell_line = rep(c("PC9", "PC9GR"), each = 6),
  treatment = rep(c("Control", "Apatinib", "Control", "Apatinib"), each = 3)
)

# Save sample metadata table for your project:
write.csv(sample_info, "sample_metadata.csv", row.names = FALSE)
head(sample_info, 4)

Merge Metadata with Expression Summaries

Calculate total reads per sample, then link them with merge():

# Total sequenced reads per sample:
lib_sizes <- data.frame(
  sample_id   = colnames(counts),
  total_reads = colSums(counts)
)

# Merge metadata with library sizes by matching sample_id:
sample_summary <- merge(sample_info, lib_sizes, by = "sample_id")
head(sample_summary, 4)

Large Files: Skip read.csv()

# Base R is correct but slow on big matrices
slow <- read.csv("visium_counts.csv")

# data.table::fread() is 10-100x faster on the same file
library(data.table)
fast <- fread("visium_counts.csv")

# Spot coordinates
coords <- fread("tissue_positions_list.csv")

📋 3. Vectorization & the apply Family

The Problem

How many genes read zero in every sample?

counts <- read.delim("counts/counts.tsv", row.names = 1)

silent_gene <- 0
for (i in 1:nrow(counts)) {
  s <- sum(counts[i, ])
  if (s == 0) {
    silent_gene <- silent_gene + 1
  }
}
silent_gene   # 10261

The Solution: One Operation, Whole Table

The same 39,376 numbers, in one call:

gene_means_fast <- rowMeans(mat)

# Identical values — the loop's vector had no names, this one does
identical(unname(gene_means_loop), unname(gene_means_fast))   # TRUE

head(gene_means_fast, 3)
head(names(gene_means_fast), 3)   # "100287102" "653635" "102466751" — free gene labels

Don’t take my word for the speedup — run it. This uses the same shape as the real matrix (39,376 genes × 12 samples):

Your absolute numbers will be larger than mine — this is R compiled to WebAssembly, not native R. The ratio is the point: one trip, versus 39,376.

The rest of the family was already there:

rowSums(mat)     # total reads per gene
colMeans(mat)    # mean expression per sample
range(mat)       # smallest and largest count in the whole table

apply(): The Margin Decides Everything

🎭 apply() Is a Loop Wearing a Mask

apply(X, 1, mean) slices out each row, calls mean() once per row, and pastes the answers back together. It is the for loop from two slides ago — with better manners.

At three rows you cannot see the cost. Scale up:

lapply() vs sapply()

tapply(): Summarise by Group

Choosing the Right Tool

📋 4. Strings & Regular Expressions

The Problem: Messy Gene Labels

Your gene list looks like this:

You need the suffix stripped and the cases normalised — across 20,000 rows.

grep(): Find the Matches

sub() vs gsub(): Replace Text

Case, Length & Slicing

Pattern Cheat Sheet

📋 5. Cleaning Missing Data (NA)

Spot the Gap: is.na()

Why mean() Returns NA

Three Ways to Handle NA

Replace or Remove?

🧬 Practice: Clean a QC Table

1. Add a pass column that is TRUE when umi > 500.
2. Replace the missing mito values with 0 in a copy of the table.
3. Report how many spots remain after both filters.

✅ Practice: Clean a QC Table — Solution

📋 6. Custom Functions

Anatomy of a Function

A function is a named recipe: inputs, a body, and a returned value.

return() or the Last Line?

Validate or Stop

🧬 Practice: Wrap Yesterday’s CPM

1. Write calc_cpm(reads, total) that returns counts per million.
2. Add a mito_pct column for spot_id, mito_reads = c(12, 4).
3. Use sapply() to turn the CPM column into a numeric vector of values.

✅ Practice: Wrap Yesterday’s CPM — Solution

📋 7. Practice Exercises

✍️ Exercise 1: Lists & Double Brackets

A QC report came back from the pipeline as a list:

1. Extract umi_counts as a vector (not wrapped in a list).
2. Compute the mean UMI count.
3. How many spots fall below 100 reads?

✅ Exercise 1: Lists & Double Brackets — Solution

✍️ Exercise 2: apply() & Missing Data

A gene expression matrix arrived with holes in it:

1. Replace every NA with 0.
2. Compute the median of each column with apply().
3. Print the column with the highest median.

✅ Exercise 2: apply() & Missing Data — Solution

✍️ Exercise 3: Clean the Gene List

1. Trim the whitespace and force everything to upper case.
2. Drop the -mouse suffix.
3. Report how many genes remain and print them.

✅ Exercise 3: Clean the Gene List — Solution

🧠 Recap in Seven Lines

• Vectors hold one type; data frames hold one type per column — count everything else you meet inside a list
• [[ ]] unwraps a list element, [ ] wraps it again
• Vectors count from 1, and a negative index removes
• grep() finds, sub()/gsub() replace, toupper()/trimws() normalise
• NA propagates on purpose — count it, then remove it or na.rm = TRUE
• Vectorise, don’t loop — colMeans(), apply(margin), lapply/sapply, tapply
• Package logic into functions and validate inputs with stop()

🚀 What’s Next?

R Graphics & Genomic Statistics — turning numbers into figures

• Reshaping: melt() / dcast() to move between wide and long tables
• Plotting: base R plot() for a first look, ggplot2 for publication figures
• Statistics: normality checks, choosing the right test per data type
• Real data: volcano plots, heatmaps, clustering on gene matrices