How today is organised

Four hours, five short blocks. Every block ends with a small exercise; the solution is in the tab next to the question, but please try it before you open it. The five exercises are exactly the practice version of the five assignment tasks listed in section.

	Block	Time	Ends with
1	Download public RNA-seq data and get GSE161731 into R	45 min	Exercise 1
2	From Ensembl IDs to gene symbols	35 min	Exercise 2
	break	15 min
3	Understanding the sample metadata	60 min	Exercise 3
4	Normalisation: CPM and TPM	40 min	Exercise 4
	break	10 min
5	Visualising genes of interest + saving objects for Lab 7	35 min	Exercise 5

Nothing today is difficult R. Almost every line is base R (read.csv, [ ], table, colSums) plus a little ggplot2. Every line is commented. If a line does not make sense, run it alone in the console and look at the result — that is the fastest way to learn.

0.1 Introduction

This lab introduces a practical workflow for working with public RNA-seq data in R. Using the GSE161731 whole-blood dataset from the Gene Expression Omnibus (GEO), we retrieve the raw count matrix and its sample metadata, align the two, and convert Ensembl gene identifiers into gene symbols. The lab then covers cleaning and verification of clinical metadata, descriptive statistics and the appropriate test for each variable type, the calculation and interpretation of CPM and TPM, and the gene-level expression plots and gene-set heatmaps used to inspect biologically selected marker genes. The dataset compares patients with COVID-19, influenza, seasonal coronavirus and bacterial infection against healthy controls, so the expected interferon and neutrophil signatures serve throughout as a positive control for the analysis.

0.2 Learning objectives

By the end of Lab 6 you should be able to:

LO12 — Get data. Download a raw count matrix and its sample table from GEO, and align the two so that column i of the counts is row i of the metadata.
LO13 — Annotate genes. Translate Ensembl gene IDs into gene symbols, and deal honestly with genes that have no symbol or a duplicated symbol.
LO14 — Describe a cohort. Clean clinical variables, check the recorded sex against the transcriptome, and summarise the cohort with the right plot and the right statistical test for each variable type.
LO15 — Normalise. Compute CPM and TPM by hand, filter low-expressed genes, and say which comparison each measure allows.
LO16 — Visualise genes. Make a per-gene box/dot plot and a gene-set heatmap for biologically chosen marker genes, and explain what they show.

0.3 Assignment

Assignment 12 — Get the data (LO12)

Download the GSE161731 count matrix and count key, and report the dimensions of both.
Keep only the samples present in both files, and show with one line of code that the columns of the counts and the rows of the metadata are in the same order.
In two sentences: what does one row and one column of this matrix represent, and what processing level are these numbers (raw counts, CPM, or TPM)?

Assignment 13 — Annotate genes (LO13)

Map the Ensembl IDs to gene symbols and report how many IDs got a symbol and how many did not.
Report how many symbols occur more than once, and state the rule you used to keep one row per symbol.
Show the number of genes per biotype and say why a blood study is often restricted to protein-coding genes.

Assignment 14 — Describe the cohort (LO14)

Produce a clean metadata table: numeric age, one sample per subject, tidy cohort labels. Report how many samples you removed and why.
Check the recorded sex against XIST and the Y-chromosome genes, and report any discordant sample.
Compare age and sex between the COVID-19 and healthy groups with the correct test for each variable type, and state in one sentence whether the two groups are comparable.

Assignment 15 — Normalise (LO15)

Filter low-expressed genes, stating your rule and how many genes survive.
Compute CPM by hand and confirm that every column sums to 1e6.
Compute TPM and show one gene where CPM and TPM rank the samples differently; explain in one sentence why that can happen.

Assignment 16 — Visualise genes (LO16)

Make a box plot with individual points of IFI27 across the five cohorts.
Make a heatmap of the interferon-stimulated gene set, annotated by cohort.
In two or three sentences: which cohort has the strongest interferon signature, and does that match what you expected from the biology?

0.4 Packages

# Run this ONCE, ideally before the session. It takes a while.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")

BiocManager::install(c(
  "GEOquery",              # download from GEO
  "biomaRt",               # gene annotation from Ensembl (online)
  "org.Hs.eg.db",          # gene annotation (offline, human)
  "AnnotationDbi",         # the query interface for org.Hs.eg.db
  "EnsDb.Hsapiens.v86",    # Ensembl annotation incl. exon coordinates (offline)
  "edgeR"                  # we only use its cpm() and filterByExpr() helpers today
))

install.packages(c("ggplot2", "dplyr", "tidyr", "pheatmap", "RColorBrewer"))

library(GEOquery)        # getGEOSuppFiles(), getGEO()
library(org.Hs.eg.db)    # human gene annotation database
library(AnnotationDbi)   # mapIds() to query that database
library(edgeR)           # cpm(), filterByExpr()
library(ggplot2)         # all plots
library(dplyr)           # small data.frame manipulations
library(tidyr)           # pivot_longer(): wide -> long for ggplot
library(pheatmap)        # heatmaps
library(RColorBrewer)    # colour palettes

theme_set(theme_bw(base_size = 12))   # one plot theme for the whole document

# Change this to YOUR folder. Everything is downloaded once and cached here.
# setwd("/Users/birzh586/Desktop/RNA-seq course")

data_dir <- "data"       # raw downloads go here
res_dir  <- "results"    # tables and figures we produce go here
dir.create(data_dir, showWarnings = FALSE)   # showWarnings = FALSE: silent if it exists
dir.create(res_dir,  showWarnings = FALSE)

1 Block 1 — Download public RNA-seq data

1.1 Where does public RNA-seq data live?

1.1.1 The three levels of an RNA-seq deposit

The same experiment is stored three times, at three levels of processing. Knowing which one you are downloading is the whole game.

Level	What it is	Where	Size
Raw reads	FASTQ files, one per run	SRA / ENA	1–10 GB per sample
Aligned data	BAM files	rarely public (controlled access)	large
Processed data	count matrix, FPKM/TPM tables	GEO / ArrayExpress supplementary files	a few MB

Today we take the processed level: a gene × sample matrix of raw counts. That is the correct input for DESeq2 in Lab 7. FPKM or TPM tables are not.

1.1.2 The repositories you will actually use

GEO (NCBI): accession GSE… for a study, GSM… for a sample. Metadata plus author-supplied processed files.
ArrayExpress / BioStudies (EMBL-EBI): the European equivalent.
SRA / ENA: the raw reads themselves. ENA has friendly FTP links.
recount3, ARCHS4: thousands of studies already re-quantified with one uniform pipeline — useful when you want counts that are comparable across studies.
TCGA / GTEx: large, deeply curated human cohorts (tumour, and normal tissue).

# Optional: ENA can be queried with a plain URL, and read.delim() reads it directly.
# Replace the accession with the BioProject shown on the GEO page of your study.
ena_url <- paste0("https://www.ebi.ac.uk/ena/portal/api/filereport",
                  "?accession=PRJNA324260&result=read_run",
                  "&fields=run_accession,sample_title,fastq_ftp&format=tsv")
runs <- read.delim(ena_url)   # a data.frame with one row per sequencing run
head(runs)

1.1.3 Reusing data responsibly

Cite the original accession and the paper. Check the licence and the consent (controlled-access human data cannot simply be downloaded). And read the paper’s methods: library prep, read length, and quantification tool all change what you can conclude.

1.2 why do we chose GSE161731 dataset

GSE161731 — McClain et al. (2021), Nature Communications, “Dysregulated transcriptional responses to SARS-CoV-2 in the periphery”. Whole blood (PAXgene) RNA-seq from patients presenting with acute respiratory illness at Duke University, plus healthy controls.

Why this dataset for a course on animal-free / human data analysis?

It is human, and it is blood. Blood is the tissue you can actually sample repeatedly in patients — the realistic alternative to an animal infection model.
It has five clinically defined groups: COVID-19, Influenza, other seasonal coronavirus, bacterial infection, and healthy. So the interesting question is not the trivial “sick vs healthy” but “which infection”, which is exactly the diagnostic problem in an emergency department.
The expected biology is strong and well documented. Viral infection triggers type I interferon signalling (IFI27, IFI44L, ISG15, RSAD2, OAS1, MX1, SIGLEC1); bacterial infection triggers a neutrophil/inflammatory programme (S100A8, S100A9, S100A12, MMP8, CD177, DEFA4, OLFM4). If our pipeline is correct, these genes must come out. That gives you a positive control: a way to know the analysis is working, before you trust anything new.
It has real, messy clinical metadata: age recorded as >89, some subjects sampled more than once, unbalanced groups. That is what real data looks like.

Whole blood also has a specific complication we will meet again in Lab 7: the biggest source of variation is often which cells are in the tube (neutrophils vs lymphocytes), not gene regulation inside one cell type.

1.3 Downloading the supplementary files

GEO stores the authors’ processed files as “supplementary files”. For GSE161731 there are two: the count matrix, and a “count key” that is the sample table.

# getGEOSuppFiles() downloads every supplementary file of the series into
# <baseDir>/GSE161731/ and returns a data.frame whose row names are the file paths.
supp <- getGEOSuppFiles("GSE161731", baseDir = data_dir, makeDirectory = TRUE)
rownames(supp)

## [1] "data/GSE161731/GSE161731_counts.csv.gz"      "data/GSE161731/GSE161731_counts_key.csv.gz" 
## [3] "data/GSE161731/GSE161731_key.csv.gz"         "data/GSE161731/GSE161731_xpr_nlcpm.csv.gz"  
## [5] "data/GSE161731/GSE161731_xpr_tpm_geo.txt.gz"

# supp <- tryCatch(
#   getGEOSuppFiles(
#     "GSE161731",
#     baseDir = data_dir,
#     makeDirectory = TRUE
#   ),
#   error = function(e) {
#     message("Some supplementary files could not be downloaded.")
#     NULL
#   }
# )

files <- list.files(file.path(data_dir, "GSE161731"), full.names = TRUE)
files # look at what we actually got

## [1] "data/GSE161731/GSE161731_counts_key.csv.gz"  "data/GSE161731/GSE161731_counts.csv.gz"     
## [3] "data/GSE161731/GSE161731_key.csv.gz"         "data/GSE161731/GSE161731_xpr_nlcpm.csv.gz"  
## [5] "data/GSE161731/GSE161731_xpr_tpm_geo.txt.gz"

# Build the two file paths. file.path() glues folder names with the correct slash.
counts_file <- file.path(data_dir, "GSE161731", "GSE161731_counts.csv.gz")
key_file    <- file.path(data_dir, "GSE161731", "GSE161731_counts_key.csv.gz")

# read.csv() reads a comma-separated file; it understands .gz without unzipping.
#   row.names = 1     -> use the first column (gene ID / sample ID) as row names
#   check.names=FALSE -> do NOT rewrite column names into "valid" R names
counts_raw <- read.csv(counts_file, row.names = 1, check.names=FALSE )
key_raw    <- read.csv(key_file,    row.names = 1, check.names=FALSE )

dim(counts_raw)   # 60675 genes x 201 samples

## [1] 60675   201

dim(key_raw)      # 198 samples x 8 clinical variables

## [1] 198   8

counts_raw[1:5, 1:4]   # top-left corner of the matrix

head(key_raw)          # first rows of the sample table

Two sanity checks that you should do with every matrix you download:

# 1. Are these integers? Raw counts must be whole numbers.
head(colSums(counts_raw))            # library size (total reads assigned) per sample

##           94189 DU09-03S0000604 DU09-03S0000611          105920 DU09-03S0000774 DU09-03S0000775 
##        52431553        50351794        52172232        53276313        43639732        56439096

all(counts_raw == round(counts_raw)) # TRUE means integer counts, not FPKM/TPM

## [1] TRUE

# 2. What do the row names look like? These should be Ensembl gene IDs.
head(rownames(counts_raw))

## [1] "ENSG00000223972" "ENSG00000227232" "ENSG00000278267" "ENSG00000243485" "ENSG00000274890"
## [6] "ENSG00000237613"

1.4 Matching counts to metadata

The columns of the count matrix and the rows of the sample table must describe the same samples, in the same order. Never assume it — check it.

# intersect() returns the IDs present in BOTH objects.
colnames(counts_raw) <- paste0("S",colnames(counts_raw))
rownames(key_raw) <- paste0("S",rownames(key_raw))

colnames(counts_raw) <- gsub("[-_]", ".", colnames(counts_raw))
rownames(key_raw) <- gsub("[-_]", ".", rownames(key_raw))

common <- intersect(colnames(counts_raw), rownames(key_raw))
length(common)   # how many samples survive

## [1] 198

# Subset both objects with the same vector, so both end up in the same order.
meta   <- key_raw[common, ] %>% arrange(cohort)       # select these rows, keep all columns
counts <- counts_raw[, rownames(meta)]     # keep all rows (genes), select these columns
# The single most important check of the day:
all(colnames(counts) == rownames(meta))   # must print TRUE

## [1] TRUE

str(meta)                    # variable names and their types

## 'data.frame':    198 obs. of  8 variables:
##  $ subject_id      : chr  "A1BD46" "ADC430" "AB87B2" "896282" ...
##  $ age             : chr  "57" "69" "72" "68" ...
##  $ gender          : chr  "Female" "Male" "Female" "Female" ...
##  $ race            : chr  "Black/African American" "Unknown/Not reported" "Black/African American" "Black/African American" ...
##  $ cohort          : chr  "Bacterial" "Bacterial" "Bacterial" "Bacterial" ...
##  $ time_since_onset: chr  NA NA NA NA ...
##  $ hospitalized    : chr  NA NA NA NA ...
##  $ batch           : int  1 1 1 2 1 1 1 1 1 1 ...

table(meta$cohort)           # how many samples per clinical group

## 
## Bacterial CoV other  COVID-19   healthy Influenza 
##        24        61        77        19        17

1.5 Exercise 1

1.5.1 Question

How many genes and how many samples are in counts?
Which sample has the largest library size, and how many times larger is it than the smallest?
key_raw and counts_raw might not contain exactly the same samples. How many samples were in one file but not the other?

1.5.2 Solution

# 1
dim(counts)

## [1] 60675   198

# 2
lib <- colSums(counts)                     # total counts per sample
names(which.max(lib))                      # which.max() gives the position; names() the ID

## [1] "S434325"

round(max(lib) / min(lib), 1)              # ratio largest / smallest

## [1] 79.6

# 3
length(setdiff(colnames(counts_raw), rownames(key_raw)))  # in counts, not in key

## [1] 3

length(setdiff(rownames(key_raw), colnames(counts_raw)))  # in key, not in counts

## [1] 0

2 Block 2 — From Ensembl IDs to gene symbols

2.1 Why identifiers are hard

ENSG00000126709 is stable and unambiguous; IFI6 is readable but not unique over time (symbols get renamed, and several Ensembl IDs can share one symbol). The rule: analyse with Ensembl IDs, display with symbols. Today we will translate once, early, so that every later line of code can be written as logcpm["IFI27", ].

2.2 Mapping with org.Hs.eg.db

org.Hs.eg.db is an offline annotation database for human. mapIds() looks up each key and returns one value per key.

gene_ids <- rownames(counts)
head(gene_ids)

## [1] "ENSG00000223972" "ENSG00000227232" "ENSG00000278267" "ENSG00000243485" "ENSG00000274890"
## [6] "ENSG00000237613"

gene_symbol <- mapIds(
  org.Hs.eg.db,
  keys      = gene_ids,     # what we have
  keytype   = "ENSEMBL",    # what type it is
  column    = "SYMBOL",     # what we want back
  multiVals = "first"       # if an ID maps to several symbols, take the first
)

gene_biotype <- mapIds(org.Hs.eg.db, 
                       keys = gene_ids, 
                       keytype = "ENSEMBL",
                       column = "GENETYPE", 
                       multiVals = "first")

# How well did it work?
length(gene_symbol)                 # one entry per gene, in the same order as the rows

## [1] 60675

sum(is.na(gene_symbol))             # IDs with no symbol

## [1] 25146

mean(!is.na(gene_symbol))           # proportion successfully mapped

## [1] 0.5855624

sum(duplicated(gene_symbol[!is.na(gene_symbol)]))   # symbols used more than once

## [1] 82

# The general-purpose alternative is biomaRt: online, versioned, more attributes.
# Use it when you need a specific Ensembl release, or extra fields like description.
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
ann_bm <- getBM(
  attributes = c("ensembl_gene_id", "external_gene_name", "gene_biotype", "description"),
  filters    = "ensembl_gene_id",
  values     = gene_ids,
  mart       = ensembl
) #,"entrezgene_id"
# Always save the result so your analysis does not change when Ensembl updates:
saveRDS(ann_bm, file.path(data_dir, "biomart_gene_info.rds"))
head(ann_bm)

2.3 One row per gene symbol

Duplicated symbols happen because of patch scaffolds, X/Y pseudoautosomal genes, read-through models and paralogues. We keep, for each symbol, the row with the highest average count — the most informative copy — and drop the rest.

# Put the rows in order of decreasing average expression.
# Why this matters: Sometimes the same gene (same symbol) maps to multiple Ensembl IDs (e.g. different annotation versions or isoforms). We don't want to just pick # one at random — we want to keep the version that's actually expressed the most, since that's usually the more reliable/meaningful measurement.
o <- order(rowMeans(counts), decreasing = TRUE)
counts      <- counts[o, ]
gene_ids    <- gene_ids[o]
gene_symbol <- gene_symbol[o]
gene_biotype<- gene_biotype[o]

# Keep a row if: it has a symbol, AND this symbol has not been seen higher up.
keep_row <- !is.na(gene_symbol) & !duplicated(gene_symbol)
sum(keep_row)   # genes we keep

## [1] 35447

counts       <- counts[keep_row, ]
gene_ids     <- gene_ids[keep_row]
gene_biotype <- gene_biotype[keep_row]
rownames(counts) <- gene_symbol[keep_row]   # from now on, rows are named by symbol

# A small annotation table, one row per gene, in the same order as counts.
ann <- data.frame(symbol  = rownames(counts),
                  ensembl = gene_ids,
                  biotype = gene_biotype,
                  row.names = rownames(counts))
head(ann)

# What kinds of genes are in the matrix?
sort(table(ann$biotype), decreasing = TRUE)

## 
## protein-coding         pseudo          ncRNA         snoRNA          other          snRNA 
##          19208           8893           6155            726            368             69 
##           rRNA          scRNA 
##             24              4

# Check that a few genes we care about are present:
c("IFI27", "ISG15", "S100A8", "XIST", "RPS4Y1") %in% rownames(counts)

## [1] TRUE TRUE TRUE TRUE TRUE

2.4 Exercise 2

2.4.1 Question

How many genes are protein-coding, and what proportion of the matrix is that?
Find all row names starting with "IFI" (interferon-induced genes).
Which of these interferon genes has the highest total count across all samples: IFI6, IFI27, IFI44L, ISG15?

2.4.2 Solution

# 1
sum(ann$biotype == "protein-coding", na.rm = TRUE)

## [1] 19208

mean(ann$biotype == "protein-coding", na.rm = TRUE)

## [1] 0.5418794

# 2  grep() with value = TRUE returns the matching names themselves.
grep("^IFI", rownames(counts), value = TRUE)

##  [1] "IFITM2"  "IFITM3"  "IFIT3"   "IFI6"    "IFIT2"   "IFI16"   "IFITM1"  "IFIT1"   "IFI44L" 
## [10] "IFI27"   "IFI35"   "IFI44"   "IFIH1"   "IFIT5"   "IFIT1B"  "IFI27L2" "IFI30"   "IFI27L1"
## [19] "IFITM10" "IFITM9P" "IFITM4P" "IFITM5"

# 3
g <- c("IFI6", "IFI27", "IFI44L", "ISG15")
sort(rowSums(counts[g, ]), decreasing = TRUE)

##    IFI6   ISG15  IFI44L   IFI27 
## 5113654 4214768 1711951 1575791

3 Block 3 — Understanding the sample metadata

Before any modelling, describe your samples. This step catches confounding, batch structure and mislabelled samples — all far cheaper to find now than after the differential expression.

3.1 Cleaning the clinical variables

# Look at the raw values of each variable before touching anything.
table(meta$cohort,       useNA = "ifany")

## 
## Bacterial CoV other  COVID-19   healthy Influenza 
##        24        61        77        19        17

table(meta$gender,       useNA = "ifany")

## 
## Female   Male   <NA> 
##     84     93     21

table(meta$hospitalized, useNA = "ifany")

## 
##   No  Yes <NA> 
##   65   12  121

table(meta$time_since_onset, useNA = "ifany")

## 
##  early   late middle   <NA> 
##     19     22     36    121

head(sort(unique(meta$age)))   # note the character value ">89"

## [1] ">89" "0"   "14"  "15"  "16"  "18"

# 1. Age is stored as text because of ">89". Replace that string, then convert.
meta$age <- as.numeric(gsub(">89", "90", meta$age))
summary(meta$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   18.00   27.50   33.93   50.00   90.00

# 2. Cohort labels contain spaces and dashes, which break formulas and file names.
#    gsub() replaces every character that is not a letter or digit with "_".
meta$cohort <- gsub("[^A-Za-z0-9]+", "_", meta$cohort)
table(meta$cohort)

## 
## Bacterial CoV_other  COVID_19   healthy Influenza 
##        24        61        77        19        17

# 3. Turn the categorical variables into factors (R's type for a grouping variable).
meta$cohort <- factor(meta$cohort)
meta$gender <- factor(meta$gender)
meta$race   <- factor(meta$race)
str(meta)

## 'data.frame':    198 obs. of  8 variables:
##  $ subject_id      : chr  "A1BD46" "ADC430" "AB87B2" "896282" ...
##  $ age             : num  57 69 72 68 75 38 76 55 38 87 ...
##  $ gender          : Factor w/ 2 levels "Female","Male": 1 2 1 1 1 1 2 2 2 1 ...
##  $ race            : Factor w/ 7 levels "American Indian/Alaska Native",..: 3 6 3 3 3 3 3 3 3 3 ...
##  $ cohort          : Factor w/ 5 levels "Bacterial","CoV_other",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ time_since_onset: chr  NA NA NA NA ...
##  $ hospitalized    : chr  NA NA NA NA ...
##  $ batch           : int  1 1 1 2 1 1 1 1 1 1 ...

# 4. One row per person. Some subjects were sampled more than once; two samples
#    from the same person are not independent observations, and most statistical
#    tests assume independence.
sum(duplicated(meta$subject_id))          # how many repeated samples are there?

## [1] 43

keep_subj <- !duplicated(meta$subject_id) # TRUE for the first sample of each subject
meta   <- meta[keep_subj, ]
counts <- counts[, rownames(meta)]             # subset the counts the SAME way

all(colnames(counts) == rownames(meta))   # check again — always

## [1] TRUE

dim(counts); table(meta$cohort)

## [1] 35447   155

## 
## Bacterial CoV_other  COVID_19   healthy Influenza 
##        24        51        46        17        17

3.2 Verifying the sex label against the data

We can read sex off the transcriptome: XIST is expressed almost only in samples with two X chromosomes, and RPS4Y1 / DDX3Y / UTY / KDM5D only from a Y chromosome. If the recorded label disagrees, something is wrong with the sample.

# A quick log2-CPM matrix, only to look at a few genes (proper version comes later).
logcpm_quick <- cpm(counts, log = TRUE, prior.count = 1)

y_genes <- c("RPS4Y1", "DDX3Y", "UTY", "KDM5D")
y_genes <- y_genes[y_genes %in% rownames(logcpm_quick)]   # keep the ones present

# Build a small data.frame with one row per sample and three columns.
sex_df <- data.frame(
  sample  = colnames(logcpm_quick),
  XIST    = logcpm_quick["XIST", ],                  # one gene, all samples
  y_score = colMeans(logcpm_quick[y_genes, ]),      # average of the Y gene      
  label   = meta$gender
)

ggplot(sex_df, aes(x = XIST, y = y_score, colour = label)) +
  geom_point(size = 2, alpha = 0.8) +
  labs(x = "XIST (log2 CPM)", y = "mean Y-gene expression (log2 CPM)",
       colour = "recorded sex", title = "Recorded sex vs transcriptome-inferred sex")

Discuss. Two clusters should be cleanly separated. A sample sitting between them is a candidate for a sample swap or a contaminated library. What would you do with a discordant sample in your own study — drop it, keep it, or ask the data manager?

# A simple rule: high Y-gene expression = male. Pick the cut half-way between the
# two obvious clusters, which you can read off the plot above.
cutoff <- mean(range(sex_df$y_score))
sex_df$inferred <- ifelse(sex_df$y_score > cutoff, "Male", "Female")

table(recorded = sex_df$label, inferred = sex_df$inferred)

##         inferred
## recorded Female Male
##   Female     69    3
##   Male        1   68

# Category 1: truly discordant samples (label exists but disagrees with inferred sex)
discordant <- sex_df$sample[!is.na(sex_df$label) & 
                             as.character(sex_df$label) != sex_df$inferred]

# Category 2: samples with missing sex label
label_missing <- sex_df$sample[is.na(sex_df$label)]

# Combine both into one exclusion list
exclude_ids <- c(discordant, label_missing)
exclude_ids

##  [1] "SDU18.02S0011675" "S434766"          "S434777"          "SDU09.02S0000101" "S434393"         
##  [6] "SDU14.02S0000003" "S434399"          "SDU14.02S0000001" "S434486"          "S434570"         
## [11] "S434532"          "S434596"          "SDU14.02S0000007" "S434537"          "S434482"         
## [16] "S434517"          "S434518"          "SDU14.02S0000005"

# Remove these samples from sex_df
sex_df_clean <- sex_df[!sex_df$sample %in% exclude_ids, ]

# If you need to clean meta and counts too, keep them in sync:
meta   <- meta[!rownames(meta) %in% exclude_ids, ]   # adjust column name as needed
meta <- meta %>% arrange(cohort)
counts <- counts[, rownames(meta)]
all(colnames(counts) ==  rownames(meta))

## [1] TRUE

3.3 Describing the cohort

3.3.1 Categorical variables: counts, a bar chart, and a (reluctant) pie chart

ggplot(meta, aes(x = cohort, fill = gender)) +
  geom_bar(position = position_dodge(preserve = "single")) +      # bars side by side
  geom_text(stat = "count",
            aes(label = after_stat(count)),       # print the n on top
            position = position_dodge(width = 0.9), vjust = -0.5, hjust = 0.5, size = 3.2) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "number of samples", title = "Samples per cohort and sex")

# The same information as a pie chart: a stacked bar bent into a circle.
tab <- as.data.frame(table(cohort = meta$cohort))   # a data.frame with Freq column
tab$prop <- tab$Freq / sum(tab$Freq)

ggplot(tab, aes(x = "", y = prop, fill = cohort)) +
  geom_col(width = 1, colour = "white") +
  coord_polar(theta = "y") +                        # this is what makes it a pie
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Cohort composition") +
    geom_text(aes(label = scales::percent(prop, accuracy = 1)),
            position = position_stack(vjust = 0.5),
            size = 3.5, colour = "black")+
  theme_void()

Note. Pie charts are hard to read because people compare angles badly. Reviewers like them; a bar chart is almost always clearer.

3.3.2 Continuous variables: box plot with the individual points

ggplot(meta, aes(x = cohort, y = age, fill = cohort)) +
  geom_boxplot(width = 0.5, outlier.shape = NA, alpha = 0.7) +  # hide outliers…
  geom_jitter(width = 0.15, size = 1.4, alpha = 0.6) +          # …because we draw all points
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "age (years)", title = "Age distribution by cohort") +
  theme(legend.position = "none")

Always show the points. A box plot of n = 4 looks exactly like a box plot of n = 400. Overlaying the observations makes sample size, spread and bimodality visible.

3.3.3 The right test for each variable type

# Restrict to the two groups we will compare in Lab 7.
sub <- meta[meta$cohort %in% c("COVID_19", "healthy"), ]
sub$cohort <- droplevels(sub$cohort)     # forget the unused factor levels
table(sub$cohort)

## 
## COVID_19  healthy 
##       45       16

# Categorical vs categorical  -> chi-square, or Fisher when any expected count < 5
tab2 <- table(sub$gender, sub$cohort)
tab2

##         
##          COVID_19 healthy
##   Female       22       6
##   Male         23      10

chisq.test(tab2)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab2
## X-squared = 0.24319, df = 1, p-value = 0.6219

fisher.test(tab2)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab2
## p-value = 0.562
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.4321978 6.2635420
## sample estimates:
## odds ratio 
##   1.582124

# Continuous vs categorical -> t-test if roughly normal, Wilcoxon if not
tapply(sub$age, sub$cohort, summary)          # summary of age within each group

## $COVID_19
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   20.00   29.00   33.00   44.38   60.00   90.00 
## 
## $healthy
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   18.00   18.00   18.38   19.00   20.00

shapiro.test(sub$age[sub$cohort == "COVID_19"])  # normality check first

## 
##  Shapiro-Wilk normality test
## 
## data:  sub$age[sub$cohort == "COVID_19"]
## W = 0.87202, p-value = 0.0001422

# t.test(age ~ cohort, data = sub)
wilcox.test(age ~ cohort, data = sub)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  age by cohort
## W = 719.5, p-value = 3.471e-09
## alternative hypothesis: true location shift is not equal to 0

# A minimal "Table 1", built with dplyr: n, median age, and % female per cohort.
meta %>%
  group_by(cohort) %>%
  summarise(
    n           = n(),
    median_age  = median(age, na.rm = TRUE),
    age_IQR     = paste0(quantile(age, .25, na.rm = TRUE), "-",
                         quantile(age, .75, na.rm = TRUE)),
    pct_female  = round(100 * mean(gender == "Female"), 1)
  )

Discuss. If age differs between COVID-19 patients and healthy controls (it usually does — controls are volunteers, patients are older), then age and cohort are confounded. In Lab 7 we handle that by putting age in the design formula. But only a table like this tells you that you need to.

3.4 Exercise 3

3.4.1 Question

How many samples are there per cohort after removing repeated subjects?
Is sex balanced across all five cohorts? Test it with one line.
Compare the age of Influenza patients with healthy controls. Which test did you choose, and what is the p-value?

3.4.2 Solution

# 1
table(meta$cohort)

## 
## Bacterial CoV_other  COVID_19   healthy Influenza 
##        24        35        45        16        17

# 2
chisq.test(table(meta$gender, meta$cohort))

## 
##  Pearson's Chi-squared test
## 
## data:  table(meta$gender, meta$cohort)
## X-squared = 1.9977, df = 4, p-value = 0.7362

# 3  Small groups and skewed age -> Wilcoxon is the safe choice.
sub_flu <- meta[meta$cohort %in% c("Influenza", "healthy"), ]
wilcox.test(age ~ cohort, data = sub_flu)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  age by cohort
## W = 151.5, p-value = 0.5652
## alternative hypothesis: true location shift is not equal to 0

4 Block 4 — Normalisation: CPM and TPM

Raw counts are not comparable across samples, for two separate reasons:

Library size (sequencing depth). A sample with 40 M reads gives roughly twice the count of a sample with 20 M reads, for the same biology. → fixed by CPM.
Gene length. A 10 kb gene collects more reads than a 1 kb gene at the same molar concentration. → fixed by TPM.

4.1 Filtering low-expressed genes

Genes with almost no reads carry no statistical information, but they cost you multiple-testing power. Remove them before normalising.

# The manual rule: keep a gene if it has more than 1 CPM in at least as many
# samples as the smallest group. min(table(...)) is that smallest group size.
cpm_all  <- cpm(counts)                     # CPM, not logged
min_grp  <- min(table(meta$cohort))
keep_gene <- rowSums(cpm_all > 1) >= min_grp
table(keep_gene)                            # FALSE = dropped, TRUE = kept

## keep_gene
## FALSE  TRUE 
## 21100 14347

# edgeR's automatic version of the same idea, group-aware:
keep_auto <- filterByExpr(counts, group = meta$cohort)
table(keep_auto)

## keep_auto
## FALSE  TRUE 
## 17898 17549

counts_f <- counts[keep_auto, ]             # the filtered count matrix
ann_f    <- ann[rownames(counts_f), ]       # keep the annotation table aligned
dim(counts_f)

## [1] 17549   137

4.2 CPM — counts per million

\[\text{CPM}_{gi} = \frac{\text{count}_{gi}}{\text{library size}_i} \times 10^6\]

lib_size <- colSums(counts_f)          # one total per sample

# t() transposes. Dividing a matrix by a vector works row-wise in R, so we
# transpose, divide by the per-sample totals, and transpose back.
cpm_manual <- t(t(counts_f) / lib_size) * 1e6
head(colSums(cpm_manual))

##         S94189         S97389 SDU09.03S19498         S94478         S97392         S95967 
##          1e+06          1e+06          1e+06          1e+06          1e+06          1e+06

4.2.1 Log transformation and the prior count

# Counts are very skewed: a few genes are huge, most are small. log2 fixes that.
# We add a small "prior count" because log2(0) is -Inf.
logcpm <- cpm(counts_f, log = TRUE, prior.count = 1) %>% data.frame() # = log2(CPM + 1), roughly
dim(logcpm)

## [1] 17549   137

logcpm[1:5,1:5]

par(mfrow = c(1, 2))                                   # two plots side by side
hist(cpm_manual[, 1], breaks = 50, main = "CPM", xlab = "CPM")
hist(logcpm[, 1],     breaks = 50, main = "log2 CPM", xlab = "log2 CPM")

par(mfrow = c(1, 1))                                   # reset the plotting layout

4.3 TPM — transcripts per million

CPM divides by library size only. TPM divides by gene length first, then by the sample total. The order matters.

\[\text{RPK}_{gi} = \frac{\text{count}_{gi}}{\text{length}_g/1000}, \qquad \text{TPM}_{gi} = \frac{\text{RPK}_{gi}}{\sum_g \text{RPK}_{gi}} \times 10^6\]

4.3.1 Getting a gene length

The correct length is the union of the exons of the gene, not the distance from start to end of the gene body.

library(EnsDb.Hsapiens.v86)

# exonsBy() returns the exons of every gene; reduce() merges overlapping exons so
# that we do not count the same base twice; width() gives their lengths in bases.
exons_by_gene <- exonsBy(EnsDb.Hsapiens.v86, by = "gene")
gene_len_all  <- sum(width(reduce(exons_by_gene)))   # one number per Ensembl gene ID

# Look up the length of each gene in our matrix, using the Ensembl ID we kept in ann_f.
gene_len <- gene_len_all[ann_f$ensembl]
names(gene_len) <- rownames(counts_f)

sum(is.na(gene_len))                 # genes with no length available

## [1] 0

head(gene_len)

##   ACTB MT-CO1  HLA-B  HLA-C  CSF3R MT-ND4 
##   3256   1542   3012   2507   8429   1378

# Keep only genes with a known length, in both objects.
has_len   <- !is.na(gene_len)
counts_t  <- counts_f[has_len, ]
len_kb    <- gene_len[has_len] / 1000            # length in kilobases

rpk <- counts_t / len_kb                          # step 1: divide each ROW by its length
tpm <- t(t(rpk) / colSums(rpk)) * 1e6             # step 2: rescale each COLUMN to 1e6

head(round(colSums(tpm)))                         # every column sums to 1e6

##         S94189         S97389 SDU09.03S19498         S94478         S97392         S95967 
##          1e+06          1e+06          1e+06          1e+06          1e+06          1e+06

4.3.2 CPM and TPM disagree — and that is the point

g <- "IFI27"
plot(log2(cpm_manual[g, ] + 1), log2(tpm[g, ] + 1),
     xlab = "log2 CPM + 1", ylab = "log2 TPM + 1", pch = 16, cex = 0.7,
     main = paste(g, ": CPM vs TPM across samples"))

For one gene across samples, CPM and TPM are almost perfectly correlated (length is constant, so it cancels). For different genes within one sample, only TPM is meaningful, because a long gene inflates its own CPM.

4.4 When to use what

Question	Use
Differential expression (DESeq2, edgeR)	raw counts — never CPM or TPM
Compare one gene across samples	CPM (or log2 CPM)
Compare different genes inside one sample	TPM
PCA, clustering, heatmaps	log2 CPM (or VST from DESeq2)
Deconvolution, marker-score signatures	TPM

4.5 Exercise 4

4.5.1 Question

How many genes did filtering remove, in absolute numbers and as a percentage?
Confirm that column 5 of tpm sums to 1e6, and that column 5 of cpm_manual does too.
Take S100A8 (short gene, extremely high in neutrophils) and MX1. Which of the two has the larger TPM/CPM ratio, and why?

4.5.2 Solution

# 1
nrow(counts) - nrow(counts_f)

## [1] 17898

(nrow(counts) - nrow(counts_f))/nrow(counts) *100

## [1] 50.49228

# 2
sum(tpm[, 5]); sum(cpm_manual[, 5])

## [1] 1e+06

## [1] 1e+06

# 3  TPM/CPM is driven by gene length: short genes gain, long genes lose.
ratio <- rowMeans(tpm[c("S100A8", "MX1"), ]) / rowMeans(cpm_manual[c("S100A8", "MX1"), ])
ratio

##    S100A8       MX1 
## 2.9967129 0.3686365

gene_len[c("S100A8", "MX1")]     # the short gene has the larger ratio

## S100A8    MX1 
##   1006   8616

5 Block 5 — Visualising genes of interest

Now the biology. Three gene sets, chosen before looking at the data:

Interferon-stimulated genes (ISGs) — the type I interferon programme that any viral infection switches on: IFI27, IFI44L, ISG15, RSAD2, OAS1, MX1, IFIT1, SIGLEC1.
Neutrophil / inflammatory genes — classically higher in bacterial infection: S100A8, S100A9, MMP8, CD177.

isg    <- c("IFI27", "IFI44L", "ISG15", "RSAD2")
neutro <- c("S100A8", "S100A9", "MMP8", "CD177")

# Only keep the genes that survived filtering (some may have been dropped).
isg    <- isg[isg %in% rownames(logcpm)]
neutro <- neutro[neutro %in% rownames(logcpm)]

marker <- c(isg, neutro)

5.1 Reshaping to long format once

ggplot2 wants one row per observation (“long” format), while our matrix has one row per gene and one column per sample (“wide” format).

# Take the ISG rows, transpose so samples are rows, and make it a data.frame.
expr_df <- as.data.frame(t(logcpm[marker, rownames(meta)]))
expr_df$sample <- rownames(expr_df)      # keep the sample ID as a column
expr_df$cohort <- meta$cohort            # attach the group (same order — we checked)

# pivot_longer(): collapse the gene columns into two columns, gene and log2cpm.
expr_long <- pivot_longer(expr_df,
                          cols      = all_of(marker),
                          names_to  = "gene",
                          values_to = "log2cpm")
expr_long$gene <- factor(expr_long$gene, levels = c(marker))
expr_long$cohort <- factor(expr_long$cohort, levels = c("healthy","Bacterial", "CoV_other", "COVID_19",  "Influenza"))
head(expr_long)

5.2 Box plot with individual points

ggplot(expr_long, aes(x = cohort, y = log2cpm, fill = cohort)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.15, size = 0.6, alpha = 0.5) +
  facet_wrap(~ gene, scales = "free_y", nrow = 2) +   # one panel per gene
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "log2 CPM", title = "Interferon-stimulated genes across cohorts") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

# # A single gene, the way you would put it in a figure panel.
one <- data.frame(cohort = meta$cohort, expr = as.numeric(logcpm["IFI27", ]))
one$cohort <- factor(one$cohort, levels = c("healthy","Bacterial", "CoV_other", "COVID_19",  "Influenza"))


ggplot(one, aes(x = cohort, y = expr, fill = cohort)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.15, size = 1.2, alpha = 0.6) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "IFI27 (log2 CPM)",
       title = "IFI27: one of the strongest single-gene markers of viral infection") +
  theme(legend.position = "none")

5.3 Heatmap of a gene list

genes_hm <- c(isg, neutro)          # all three sets in one figure
mat      <- logcpm[genes_hm, rownames(meta)]

# Sample annotation bar: a data.frame whose row names match the matrix columns.
ann_col <- data.frame(cohort = meta$cohort, row.names = colnames(mat))

# Gene annotation bar: which set each gene belongs to.
ann_row <- data.frame(
  set = c(rep("interferon", length(isg)),
          rep("neutrophil", length(neutro))),
  row.names = genes_hm
)

pheatmap(mat,
         scale             = "row",     # z-score each gene, or big genes dominate the colours
         annotation_col    = ann_col,
         annotation_row    = ann_row,
         show_colnames     = FALSE,     # 130+ sample names would be unreadable
         cluster_rows      = FALSE,    
         cluster_cols      = F,      # keep our biological ordering
         color             = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         main              = "Marker gene sets across the GSE161731 cohorts")

Read the heatmap biologically. Do the viral cohorts (COVID-19, Influenza, seasonal coronavirus) cluster together on the interferon block? Does the bacterial cohort separate on the neutrophil block? This is the positive control that tells you the data and your pipeline are behaving.

5.4 Exercise 5

5.4.1 Question

Make a box plot of S100A8 across cohorts. Which group is highest?
Build a simple “interferon score”: the mean log2 CPM of the eight ISGs per sample. Plot it per cohort.
Redo the heatmap using only the T-cell genes, and describe in one sentence what you see for COVID-19.

5.4.2 Solution

# 1
d1 <- data.frame(cohort = meta$cohort, expr = as.numeric(logcpm["S100A8",]))
ggplot(d1, aes(cohort, expr, fill = cohort)) +
  geom_boxplot(alpha = .7, outlier.shape = NA) + geom_jitter(width = .15, size = .8) +
  labs(y = "S100A8 (log2 CPM)", x = NULL) + theme(legend.position = "none")

# 2  colMeans() of the ISG rows = one score per sample.
d2 <- data.frame(cohort = meta$cohort, score = colMeans(logcpm[isg,]))
ggplot(d2, aes(cohort, score, fill = cohort)) +
  geom_boxplot(alpha = .7, outlier.shape = NA) + geom_jitter(width = .15, size = .8) +
  labs(y = "interferon score (mean log2 CPM)", x = NULL) + theme(legend.position = "none")

# 3
pheatmap(logcpm[isg, ], scale = "row", show_colnames = FALSE,
         annotation_col = data.frame(cohort = meta$cohort, row.names = colnames(logcpm)),
          cluster_cols      = F,      # let the samples group themselves
         main = "Interferon-stimulated genes (ISGs)")

6 Save the objects for Lab 7

Lab 7 starts from exactly these four objects, so save them now.

saveRDS(counts_f, file.path(data_dir, "gse161731_counts_filtered.rds"))  # raw counts, filtered
saveRDS(counts,   file.path(data_dir, "gse161731_counts_all.rds"))       # raw counts, unfiltered
saveRDS(meta,     file.path(data_dir, "gse161731_meta.rds"))             # clean sample table
saveRDS(logcpm,   file.path(data_dir, "gse161731_logcpm.rds"))           # log2 CPM for plots
saveRDS(ann_f,    file.path(data_dir, "gse161731_gene_annotation.rds"))  # symbol / ensembl / biotype


sessionInfo()   # paste this into any report: it is your reproducibility record

## R version 4.5.0 (2025-04-11)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Stockholm
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.32.0          AnnotationFilter_1.32.0  
##  [4] GenomicFeatures_1.60.0    GenomeInfoDb_1.44.3       GenomicRanges_1.62.1     
##  [7] Seqinfo_1.0.0             RColorBrewer_1.1-3        pheatmap_1.0.13          
## [10] tidyr_1.3.2               dplyr_1.2.1               ggplot2_4.0.3            
## [13] edgeR_4.6.3               limma_3.64.3              org.Hs.eg.db_3.21.0      
## [16] AnnotationDbi_1.72.0      IRanges_2.44.0            S4Vectors_0.48.1         
## [19] GEOquery_2.76.0           Biobase_2.70.0            BiocGenerics_0.56.0      
## [22] generics_0.1.4           
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.3.0                   bitops_1.1-0                rlang_1.3.0                
##  [4] magrittr_2.0.5              otel_0.2.0                  matrixStats_1.5.0          
##  [7] compiler_4.5.0              RSQLite_3.53.3              png_0.1-9                  
## [10] vctrs_0.7.3                 ProtGenerics_1.40.0         pkgconfig_2.0.3            
## [13] crayon_1.5.3                fastmap_1.2.0               XVector_0.50.0             
## [16] labeling_0.4.3              Rsamtools_2.24.1            rmarkdown_2.32             
## [19] tzdb_0.5.0                  UCSC.utils_1.4.0            purrr_1.2.2                
## [22] bit_4.6.0                   xfun_0.60                   cachem_1.1.0               
## [25] jsonlite_2.0.0              blob_1.3.0                  DelayedArray_0.36.1        
## [28] BiocParallel_1.44.0         parallel_4.5.0              R6_2.6.1                   
## [31] bslib_0.12.0                rtracklayer_1.68.0          jquerylib_0.1.4            
## [34] SummarizedExperiment_1.40.0 knitr_1.51                  readr_2.2.0                
## [37] rentrez_1.2.4               Matrix_1.7-6                tidyselect_1.2.1           
## [40] rstudioapi_0.19.0           abind_1.4-8                 yaml_2.3.12                
## [43] codetools_0.2-20            curl_8.0.0                  lattice_0.23-1             
## [46] tibble_3.3.1                withr_3.0.3                 KEGGREST_1.50.0            
## [49] S7_0.2.2                    evaluate_1.0.5              xml2_1.6.0                 
## [52] Biostrings_2.78.0           pillar_1.11.1               MatrixGenerics_1.22.0      
## [55] RCurl_1.98-1.20             hms_1.1.4                   scales_1.4.0               
## [58] glue_1.8.1                  lazyeval_0.2.3              tools_4.5.0                
## [61] BiocIO_1.18.0               data.table_1.18.6.1         locfit_1.5-9.12            
## [64] GenomicAlignments_1.44.0    XML_3.99-0.24               grid_4.5.0                 
## [67] GenomeInfoDbData_1.2.14     restfulr_0.0.17             cli_3.6.6                  
## [70] S4Arrays_1.10.1             gtable_0.3.6                sass_0.4.10                
## [73] digest_0.6.39               SparseArray_1.10.10         rjson_0.2.23               
## [76] farver_2.1.2                memoise_2.0.1               htmltools_0.5.9            
## [79] lifecycle_1.0.5             httr_1.4.9                  statmod_1.5.2              
## [82] bit64_4.8.6

Take-home messages

Download raw counts, not FPKM, if you intend to do differential expression.
Check that counts and metadata are aligned. Then check again after every subset.
Analyse with stable IDs, display with symbols, and be explicit about how you resolved duplicated symbols.
Clean the clinical variables and verify them against the data (sex genes, marker genes) before you trust a single label.
CPM corrects depth; TPM corrects depth and length; neither belongs in DESeq2.
Choose your marker genes from the biology before you look at the result. They are your positive control.

Forsoksdjursvetenskap (teoretisk) och avancerad databearbetning

Lab 6 - Finding public RNA-seq data, annotating genes, and exploring counts and metadata

Birong Zhang

16 October 2026