Before you start

0.1 What this lab covers

By the end of this session you should be able to:

Name the main public repositories for RNA-seq data and explain what kind of data each one holds.
Pull a raw count matrix and its sample metadata into R directly from GEO, and assemble them into a SummarizedExperiment.
Convert between gene identifier systems with biomaRt (and know when not to use biomaRt).
Describe a clinical metadata table properly: counts, plots and the right statistical test per variable type.
Compute CPM and TPM by hand, and say precisely which comparisons each one licenses.
Produce the three plots you will make in every single RNA-seq project: per-gene box/dot plots, a gene-set heatmap, and a PCA.

This is an advanced lab: we assume you are comfortable with base R, data.frame/matrix subsetting, and the basics of ggplot2. Where a step is normally done with a one-line wrapper, we do it by hand once, then show the wrapper.

0.2 Packages

# Run once. This takes a while — do it before the session if you can.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")

BiocManager::install(c(
  # data retrieval
  "GEOquery", "SummarizedExperiment","recount3" ,"airway","TCGAbiolinks", 
  # annotation
  "biomaRt", "org.Hs.eg.db", "AnnotationDbi", "ensembldb", "EnsDb.Hsapiens.v86",
  "GenomicRanges",
  # counts handling
  "edgeR", "DESeq2",
  # plotting
  "pheatmap", "ComplexHeatmap"
))

install.packages(c("ggplot2", "dplyr", "tidyr", "tibble", "ggrepel", "matrixStats",
                   "RColorBrewer", "gridExtra"))

library(rmarkdown)
library(GEOquery)
library(SummarizedExperiment)
library(biomaRt)
library(edgeR)
library(ggplot2)
library(dplyr)
library(tidyr)
library(tibble)
library(matrixStats)
library(pheatmap)
library(RColorBrewer)
library(gridExtra)

theme_set(theme_bw(base_size = 12))

# Everything we download is cached here so we never download twice.
setwd("/Users/birzh586/Desktop/RNA-seq course")
data_dir <- "data"
res_dir  <- "results"
dir.create(data_dir, showWarnings = FALSE)
dir.create(res_dir,  showWarnings = FALSE)

1 Where does public RNA-seq data live?

Almost every published RNA-seq experiment is deposited somewhere public, and journals increasingly enforce this. The practical problem is not finding data — it is knowing which level of processing you are downloading, and whether the metadata is good enough to reuse.

1.1 The three levels of an RNA-seq deposit

Level	What it is	Typical file	Where
Raw reads	Unprocessed sequencer output	`.fastq.gz`	SRA, ENA, DDBJ
Alignments	Reads placed on a genome	`.bam`, `.cram`	rarely deposited; EGA/GDC for controlled data
Processed	Gene- or transcript-level counts, FPKM, TPM	`.txt/.tsv/.csv` supplementary files, `.h5`	GEO, ArrayExpress, recount3, ARCHS4

Rule of thumb. If you only need a count matrix to re-analyse, take the processed level — but check what pipeline produced it. If you need to compare against your own data, or the deposited counts are FPKM-only, go back to the raw reads and quantify everything yourself with one pipeline.

1.2 The repositories you will actually use

GEO (Gene Expression Omnibus, NCBI) — the most common landing place for functional genomics. Accessions: GSE (a series/experiment), GSM (one sample), GPL (a platform). GEO stores the metadata and any processed files the authors uploaded as “supplementary files”; the raw reads are pushed to SRA. In R: GEOquery.

SRA (Sequence Read Archive, NCBI) — raw reads. Accessions: SRP/PRJNA (project), SRX (experiment), SRR (run), SRS (sample). Downloading requires sra-toolkit (prefetch, fasterq-dump) and the files are in a compressed SRA format.

ENA (European Nucleotide Archive, EMBL-EBI) — mirrors SRA content but serves plain fastq.gz over FTP/HTTPS with no special toolkit. Accessions: PRJEB/PRJNA, ERR/SRR. For downloading fastq files, ENA is usually the path of least resistance. It also has a clean REST API that returns a tab-separated table of every run in a project, including FTP paths, file sizes and MD5 sums.

ArrayExpress / BioStudies (EMBL-EBI) — accessions E-MTAB-xxxx. Similar role to GEO; well-structured sample-and-data-relationship files (SDRF/IDF), which are often better machine-readable metadata than GEO’s free text.

Controlled-access archives — EGA (Europe) and dbGaP (US) hold human data that cannot be made public (identifiable germline variation). You need a data access agreement and a committee approval. Plan for months, not days.

Pre-quantified resources — these save you a pipeline run:

recount3 — uniformly processed counts for >700,000 human and mouse samples from SRA/GTEx/TCGA, delivered straight into R as a RangedSummarizedExperiment.
ARCHS4 — uniformly processed counts for most public human/mouse RNA-seq, distributed as HDF5.
Expression Atlas (EBI) — curated, re-analysed experiments with differential results already computed.
GDC / TCGA / TARGET — cancer genomics; use TCGAbiolinks or GenomicDataCommons in R.
GTEx — normal tissue reference, bulk RNA-seq across ~50 tissues.

Bioconductor experiment data packages — small, tidy, versioned datasets designed for teaching and method benchmarking: airway, parathyroidSE, pasilla, tweeDEseqCountData. Perfect when you want a dataset that will still install identically in three years.

1.3 Practical: querying ENA from R without leaving R

The ENA “filereport” endpoint returns a TSV you can read straight with read.delim(). This is the cleanest way to build a download manifest for a project.

ena_runs <- function(accession,
                     fields = c("run_accession", "sample_title", "library_layout",
                                "read_count", "fastq_ftp", "fastq_bytes", "fastq_md5")) {
  url <- sprintf(
    "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=%s&result=read_run&fields=%s&format=tsv",
    accession, paste(fields, collapse = ",")
  )
  read.delim(url, stringsAsFactors = FALSE)
}

runs <- ena_runs("PRJNA229998")   # any project accession
head(runs)
sum(as.numeric(unlist(strsplit(runs$fastq_bytes, ";")))) / 1e9  # total GB before you commit

Task. Before downloading anything, always compute the total size. A “small” 40-sample human RNA-seq project is comfortably 300–600 GB of fastq. This single line has saved many full disks.

1.4 Reusing data responsibly

Cite the original accession and the paper.
Check the licence/terms — most GEO data is unrestricted, but consent for human data can restrict re-use.
Record the exact accession, file name, download date and pipeline version in your script. “Downloaded from GEO” is not a method section.
Never assume the metadata is right. We will demonstrate below that the recorded sex of a sample can be checked — and sometimes contradicted — by the data itself.

2 Getting a count matrix and metadata into R

2.1 The dataset

We will use GSE81089, a bulk RNA-seq study of non-small cell lung cancer (NSCLC) from Uppsala: tumour tissue from ~199 patients plus normal lung tissue, sequenced on Illumina HiSeq2500, aligned to GRCh37 and counted with featureCounts against Ensembl 73 annotation.

Why this dataset for teaching:

The authors deposited a raw count matrix, not just FPKM — so we can do a proper count-based analysis.
Row names are Ensembl gene IDs, which forces us to deal with ID mapping.
The sample metadata is real clinical metadata — age, sex, stage, histology, smoking, survival — with all the messiness that implies.
Two histological subtypes (adenocarcinoma vs squamous cell carcinoma) give us a biologically meaningful, marker-verifiable contrast.

2.2 Downloading the metadata (series matrix)

getGEO() retrieves the series matrix: an ExpressionSet whose phenoData holds one row per GSM sample. For RNA-seq series the expression slot is usually empty — the numbers live in the supplementary files.

gse_id <- "GSE81089"

gse <- getGEO(gse_id, GSEMatrix = TRUE, destdir = data_dir, getGPL = FALSE)
length(gse)          # a series can span several platforms

## [1] 1

eset <- gse[[1]]
eset

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 0 features, 218 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: GSM2142443 GSM2142444 ... GSM2142660 (218 total)
##   varLabels: title geo_accession ... vital date:ch1 (57 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
##   pubMedIds: 29282718
## 33576873
## 41540507 
## Annotation: GPL16791

pheno_raw <- pData(eset)
dim(pheno_raw)

## [1] 218  57

colnames(pheno_raw)[1:15]

##  [1] "title"                 "geo_accession"         "status"                "submission_date"      
##  [5] "last_update_date"      "type"                  "channel_count"         "source_name_ch1"      
##  [9] "organism_ch1"          "characteristics_ch1"   "characteristics_ch1.1" "characteristics_ch1.2"
## [13] "characteristics_ch1.3" "characteristics_ch1.4" "characteristics_ch1.5"

GEO stores user-supplied annotation in free-text columns called characteristics_ch1, characteristics_ch1.1, … each formatted as "key: value". Parsing them is the same chore in every project, so write it once:

#' Turn GEO `characteristics_ch1*` columns into a tidy data.frame
#'
#' @param pheno the pData() of a GEO ExpressionSet
#' @return data.frame with one column per distinct characteristic key
parse_geo_characteristics <- function(pheno) {
  ch_cols <- grep("^characteristics_ch1", colnames(pheno), value = TRUE)
  long <- do.call(rbind, lapply(ch_cols, function(cc) {
    v <- as.character(pheno[[cc]])
    key <- sub(":.*$", "", v)
    val <- trimws(sub("^[^:]*:", "", v))
    data.frame(sample = rownames(pheno), key = trimws(key), value = val,
               stringsAsFactors = FALSE)
  }))
  long <- long[!is.na(long$key) & long$key != "" & long$key != long$value, ]
  # make syntactically valid, stable column names
  long$key <- make.names(tolower(gsub("[^A-Za-z0-9]+", "_", long$key)))
  wide <- tidyr::pivot_wider(long, id_cols = sample, names_from = key,
                             values_from = value,
                             values_fn = function(x) paste(unique(x), collapse = "|"))
  as.data.frame(wide)
}

chars <- parse_geo_characteristics(pheno_raw)
str(chars)

## 'data.frame':    218 obs. of  11 variables:
##  $ sample              : chr  "GSM2142443" "GSM2142444" "GSM2142445" "GSM2142446" ...
##  $ tumor_t_or_normal_n_: chr  "L400T" "L401T" "L404T" "L406T" ...
##  $ stage_tnm           : chr  "3" "5" "3" "1" ...
##  $ histology           : chr  "2" "2" "2" "1" ...
##  $ surgery_date        : chr  "2006-10-05" "2006-02-28" "2006-05-16" "2006-12-11" ...
##  $ age                 : chr  "78" "74" "74" "64" ...
##  $ gender              : chr  "male" "male" "female" "male" ...
##  $ vital_date          : chr  "2010-02-19" "2011-05-01" "2012-09-19" "2013-04-28" ...
##  $ dead                : chr  "1" "1" "1" "0" ...
##  $ smoking             : chr  "3" "2" "2" "1" ...
##  $ ps_who              : chr  "1" "1" "0" "0" ...

pheno <- data.frame(
  geo_accession = pheno_raw$geo_accession,
  title         = as.character(pheno_raw$title),
  source        = as.character(pheno_raw$source_name_ch1),
  stringsAsFactors = FALSE
)
pheno <- cbind(pheno, chars[match(pheno$geo_accession, chars$sample), -1, drop = FALSE])
rownames(pheno) <- pheno$geo_accession

head(pheno)

# Which characteristics did the authors actually provide, and how are they coded?
lapply(pheno[, setdiff(colnames(pheno), c("geo_accession", "title")), drop = FALSE],
       function(x) head(sort(table(x), decreasing = TRUE), 6))

## $source
## x
##         Human NSCLC tissue Human non-malignant tissue 
##                        199                         19 
## 
## $tumor_t_or_normal_n_
## x
## L400T L401T L404T L406T L413T L414T 
##     1     1     1     1     1     1 
## 
## $stage_tnm
## x
##  1  2  5  3  4  7 
## 70 45 33 25 23  3 
## 
## $histology
## x
##   2   1   3 
## 108  67  24 
## 
## $surgery_date
## x
## 2008-12-18 2010-04-21 2010-05-28 2010-06-03 2010-08-11 2010-10-11 
##          2          2          2          2          2          2 
## 
## $age
## x
## 65 62 67 74 69 71 
## 13 12 12 12 11 11 
## 
## $gender
## x
## female   male 
##    103     96 
## 
## $vital_date
## x
## 2013-04-28 2013-04-10 2008-03-26 2007-01-04 2007-03-15 2007-05-27 
##         76         29          2          1          1          1 
## 
## $dead
## x
##   0   1 n/a 
## 105  93   1 
## 
## $smoking
## x
##  1  2  3 
## 96 84 19 
## 
## $ps_who
## x
##   0   1   2 
## 120  77   2

Discuss. Look at the codings above. Several clinical fields are stored as bare integers (1, 2, 3) with the key hidden in the paper, not in GEO. Numeric-looking columns are also read as character. This is normal. Never trust a factor level you have not verified; further down we verify sex against Y-chromosome expression and histology against marker genes.

2.3 Downloading the count matrix (supplementary files)

supp_dir <- file.path(data_dir, gse_id)

if (!dir.exists(supp_dir)) {
  getGEOSuppFiles(gse_id, baseDir = data_dir, makeDirectory = TRUE)
}
supp_files <- list.files(supp_dir, full.names = TRUE)
basename(supp_files)

## [1] "GSE81089_FPKM_cufflinks.tsv.gz"           "GSE81089_readcounts_featurecounts.tsv.gz"

# Pick the raw-count file rather than the FPKM one.
count_file <- grep("count", supp_files, ignore.case = TRUE, value = TRUE)[1]
count_file

## [1] "data/GSE81089/GSE81089_readcounts_featurecounts.tsv.gz"

counts_raw <- read.delim(gzfile(count_file), header = TRUE, row.names = 1,
                         check.names = FALSE, stringsAsFactors = FALSE)
dim(counts_raw)

## [1] 63152   218

counts_raw[1:5, 1:5]

Note. check.names = FALSE matters: R would otherwise prepend X to column names that start with a digit, silently breaking the match to the metadata.

2.4 The unglamorous part: matching counts to metadata

This is where most public-data re-analyses go wrong. Column names of the count matrix and row names of the metadata almost never agree, because the submitter used internal sample IDs in the file and GSM accessions in the series matrix.

head(colnames(counts_raw))

## [1] "L400T" "L401T" "L404T" "L406T" "L413T" "L414T"

head(pheno$title)

## [1] "L400T" "L401T" "L404T" "L406T" "L413T" "L414T"

# Strategy: find the metadata column whose values best match the count columns.
match_score <- sapply(pheno, function(x) sum(colnames(counts_raw) %in% as.character(x)))
sort(match_score, decreasing = TRUE)[1:5]

## tumor_t_or_normal_n_                title        geo_accession               source 
##                  216                  197                    0                    0 
##            stage_tnm 
##                    0

key_col <- names(which.max(match_score))
key_col

## [1] "tumor_t_or_normal_n_"

if (max(match_score) < 0.5 * ncol(counts_raw)) {
  # fall back to fuzzy matching: strip prefixes/suffixes and compare
  norm <- function(x) toupper(gsub("[^A-Za-z0-9]", "", x))
  idx  <- match(norm(colnames(counts_raw)), norm(pheno$title))
} else {
  idx <- match(colnames(counts_raw), as.character(pheno[[key_col]]))
}

sum(is.na(idx))                    # count columns with no metadata

## [1] 2

keep <- !is.na(idx)

counts <- as.matrix(counts_raw[, keep, drop = FALSE])
mode(counts) <- "integer"
coldata <- pheno[idx[keep], , drop = FALSE]
rownames(coldata) <- colnames(counts)

stopifnot(identical(colnames(counts), rownames(coldata)))
dim(counts); dim(coldata)

## [1] 63152   216

## [1] 216  13

Task. Report how many samples were dropped and why. In a real project this sentence belongs in your methods: “Of the N samples in GSE81089, n had counts but no usable metadata and were excluded.”

2.5 Wrapping it up in a `SummarizedExperiment`

Keeping counts, sample annotation and gene annotation in one object means a subsetting operation can never desynchronise them. This is the single best defensive habit in Bioconductor.

se <- SummarizedExperiment(
  assays  = list(counts = counts),
  colData = DataFrame(coldata)
)
se

## class: SummarizedExperiment 
## dim: 63152 216 
## metadata(0):
## assays(1): counts
## rownames(63152): ENSG00000000003 ENSG00000000005 ... ENSG00000272545 TC%
## rowData names(0):
## colnames(216): L400T L401T ... L887T L890T
## colData names(13): geo_accession title ... smoking ps_who

# Subsetting keeps everything aligned automatically:
se[, 1:3]

## class: SummarizedExperiment 
## dim: 63152 3 
## metadata(0):
## assays(1): counts
## rownames(63152): ENSG00000000003 ENSG00000000005 ... ENSG00000272545 TC%
## rowData names(0):
## colnames(3): L400T L401T L404T
## colData names(13): geo_accession title ... smoking ps_who

assay(se, "counts")[1:3, 1:3]

##                 L400T L401T L404T
## ENSG00000000003  6364  5953  3179
## ENSG00000000005    17     1     4
## ENSG00000000419  2255  3068  2342

colData(se)[1:3, 1:4]

## DataFrame with 3 rows and 4 columns
##       geo_accession       title             source tumor_t_or_normal_n_
##         <character> <character>        <character>          <character>
## L400T    GSM2142443       L400T Human NSCLC tissue                L400T
## L401T    GSM2142444       L401T Human NSCLC tissue                L401T
## L404T    GSM2142445       L404T Human NSCLC tissue                L404T

2.6 Other one-liner routes to public counts

2.6.1 recount3

# library(BiocFileCache)
# bfc <- BiocFileCache()
# removeCache(bfc)
# library(recount3)
# projects <- available_projects()                       # ~ 8,600 human projects
# subset(projects, project == "SRP009615")
# rse <- create_rse(subset(projects, project == "SRP009615" & project_type == "data_sources"))
# assay(rse, "raw_counts") <- transform_counts(rse)      # base-pair coverage -> read counts

2.6.2 A Bioconductor data package

library(airway)
data(airway)
airway                     # 8 samples, human airway smooth muscle, +/- dexamethasone
head(assay(airway))
colData(airway)

2.6.3 TCGA

library(TCGAbiolinks)
query <- GDCquery(project = "TCGA-LUAD",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "STAR - Counts",
                  sample.type = c("Primary Tumor", "Solid Tissue Normal")
)

getResults(query) |> dim()
nrow(getResults(query))
# GDCdownload(query)
# luad <- GDCprepare(query)   # a RangedSummarizedExperiment

3 From Ensembl IDs to gene symbols

3.1 Why identifiers are hard

Ensembl gene IDs (ENSG00000141510) are stable and unambiguous, but unreadable. They carry a version suffix (ENSG00000141510.12) which changes when the model changes — and which breaks every match() if one of your two tables has it and the other does not.
Entrez Gene IDs (7157) are stable integers used by NCBI tools and most pathway databases (KEGG, and internally by many enrichment packages).
HGNC symbols (TP53) are readable and unstable: symbols get renamed, and the mapping to Ensembl IDs is many-to-many. Symbols are also famously corrupted by Excel (SEPT2 → 2-Sep, MARCH1 → 1-Mar). HGNC eventually renamed those genes (SEPTIN2, MARCHF1) because of it.
Annotation is versioned. Our counts were generated against Ensembl 73 / GRCh37. If we annotate with the current Ensembl (GRCh38) some IDs will be retired and return NA.

The rule: keep the stable ID as your row identifier all the way through the analysis, and attach the symbol as an extra column used only for display. Never re-key your matrix on symbols.

3.2 biomaRt: the general-purpose route

biomaRt queries the live Ensembl BioMart web service, so it works for any species Ensembl hosts and can return almost any per-gene attribute.

library(biomaRt)
ensembl <- useEnsembl(biomart = "genes")
datasets <- listDatasets(ensembl)
head(datasets)
# datasets[grep("sapiens", datasets$dataset), ]

ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")

Attributes are the columns you can ask for; filters are the columns you can query on.

attributes <- listAttributes(ensembl)
attributes[1:15, ]
nrow(attributes)                             # thousands of them
attributes[grep("symbol|external_gene", attributes$name), ][1:10, ]

filters <- listFilters(ensembl)
head(filters)

3.2.1 Matching the annotation version

listEnsemblArchives()          # every past release with its URL

# Our counts came from Ensembl 73 / GRCh37, so the honest thing to do is:
ensembl_grch37 <- useEnsembl(biomart = "genes",
                             dataset = "hsapiens_gene_ensembl",
                             GRCh = 37)

For teaching we will annotate against the current release and simply quantify how many IDs are lost — that number is itself informative.

3.2.2 Fetching gene information

# gene_ids <- sub("\\..*$", "", rownames(se))     # strip version suffixes if present
# head(gene_ids)
# 
# ann_file <- file.path(data_dir, "biomart_gene_info.rds")
# if (file.exists(ann_file)) {
#   biomart_gene_info <- readRDS(ann_file)
# } else {
#   biomart_gene_info <- getBM(
#     attributes = c("ensembl_gene_id", "external_gene_name", "description",
#                    "gene_biotype", "chromosome_name", "start_position",
#                    "end_position", "strand"),
#     filters    = "ensembl_gene_id",
#     values     = gene_ids,
#     mart       = ensembl_grch37
#   )
#   saveRDS(biomart_gene_info, ann_file)
# }
# 
# head(biomart_gene_info)
# nrow(biomart_gene_info)
# mean(gene_ids %in% biomart_gene_info$ensembl_gene_id)   # retrieval rate

# Download the correct GTF annotation from the Ensembl release-73 archive and put it in  data
# https://ftp.ensembl.org/pub/release-73/gtf/homo_sapiens/Homo_sapiens.GRCh37.73.gtf.gz


# Read the GTF
# BiocManager::install("rtracklayer")
library(rtracklayer)
gtf_file <- file.path(
  data_dir,
  "Homo_sapiens.GRCh37.73.gtf.gz"
)
gtf <- import(gtf_file)
colnames(mcols(gtf))

##  [1] "source"          "type"            "score"           "phase"           "gene_id"        
##  [6] "transcript_id"   "exon_number"     "gene_name"       "gene_biotype"    "transcript_name"
## [11] "exon_id"         "protein_id"

table(gtf$type)

## 
##        exon         CDS start_codon  stop_codon 
##     1302592      789588       92706       83203

# A GTF file contains genomic annotation such as: - chromosome, genomic start, genomic end, strand, gene ID, gene name...
# Extract gene-level annotation
# Keep exon records
gtf_exon <- gtf[gtf$type == "exon"]
head(gtf_exon)

## GRanges object with 6 ranges and 12 metadata columns:
##       seqnames      ranges strand |                 source     type     score     phase
##          <Rle>   <IRanges>  <Rle> |               <factor> <factor> <numeric> <integer>
##   [1]        1 11869-12227      + | processed_transcript       exon        NA      <NA>
##   [2]        1 12613-12721      + | processed_transcript       exon        NA      <NA>
##   [3]        1 13221-14409      + | processed_transcript       exon        NA      <NA>
##   [4]        1 11872-12227      + | unprocessed_pseudogene     exon        NA      <NA>
##   [5]        1 12613-12721      + | unprocessed_pseudogene     exon        NA      <NA>
##   [6]        1 13225-14412      + | unprocessed_pseudogene     exon        NA      <NA>
##               gene_id   transcript_id exon_number   gene_name gene_biotype transcript_name
##           <character>     <character> <character> <character>  <character>     <character>
##   [1] ENSG00000223972 ENST00000456328           1     DDX11L1   pseudogene     DDX11L1-002
##   [2] ENSG00000223972 ENST00000456328           2     DDX11L1   pseudogene     DDX11L1-002
##   [3] ENSG00000223972 ENST00000456328           3     DDX11L1   pseudogene     DDX11L1-002
##   [4] ENSG00000223972 ENST00000515242           1     DDX11L1   pseudogene     DDX11L1-201
##   [5] ENSG00000223972 ENST00000515242           2     DDX11L1   pseudogene     DDX11L1-201
##   [6] ENSG00000223972 ENST00000515242           3     DDX11L1   pseudogene     DDX11L1-201
##               exon_id  protein_id
##           <character> <character>
##   [1] ENSE00002234944        <NA>
##   [2] ENSE00003582793        <NA>
##   [3] ENSE00002312635        <NA>
##   [4] ENSE00002234632        <NA>
##   [5] ENSE00003608237        <NA>
##   [6] ENSE00002306041        <NA>
##   -------
##   seqinfo: 255 sequences from an unspecified genome; no seqlengths

# Build exon-level annotation table
gene_annotation <- data.frame(
  ensembl_gene_id = as.character(gtf_exon$gene_id),
  external_gene_name = as.character(gtf_exon$gene_name),
  gene_biotype = as.character(gtf_exon$gene_biotype),
  chromosome_name = as.character(seqnames(gtf_exon)),
  start_position = start(gtf_exon),
  end_position = end(gtf_exon),
  strand = as.character(strand(gtf_exon)),
  stringsAsFactors = FALSE
)

# Convert to one row per gene
gene_annotation <- gene_annotation %>%
  group_by(ensembl_gene_id) %>%
  summarise(
    external_gene_name = first(external_gene_name),
    gene_biotype = first(gene_biotype),
    chromosome_name = first(chromosome_name),
    start_position = min(start_position),
    end_position = max(end_position),
    strand = first(strand),
    .groups = "drop"
  )

head(gene_annotation)

dim(gene_annotation)

## [1] 63152     7

gene_ids <- rownames(counts)
gene_annotation <- gene_annotation[gene_annotation$ensembl_gene_id %in% gene_ids,]
mean(gene_ids %in% gene_annotation$ensembl_gene_id)   # retrieval rate

## [1] 0.9999842

Task. What fraction of IDs failed to retrieve? Those are genes retired between Ensembl 73 and the current release (mostly clone-based and predicted models). If that fraction is large for your data, re-quantify rather than patch.

3.2.3 Searching the annotation

# # Genes whose description mentions an oncogene:
# head(biomart_gene_info[grep("oncogene", biomart_gene_info$description), c(2, 3)])
# 
# # Genes whose description mentions tumor protein:
# head(biomart_gene_info[grep("tumor protein", biomart_gene_info$description), c(2, 3)])

# Genes whose symbol contains TP53:
gene_annotation[grep("^TP53", gene_annotation$external_gene_name), c(1, 2, 4)]

# What biotypes are in this matrix?
sort(table(gene_annotation$gene_biotype), decreasing = TRUE)[1:10]

## 
##       protein_coding           pseudogene              lincRNA            antisense 
##                22695                15533                 6969                 5241 
##                miRNA             misc_RNA                snRNA               snoRNA 
##                 3352                 2171                 2063                 1537 
## processed_transcript       sense_intronic 
##                 1124                  742

3.3 Attaching annotation to the object, safely

ann <- gene_annotation[match(gene_ids, gene_annotation$ensembl_gene_id), ]
rownames(ann) <- rownames(se)
ann$ensembl_gene_id[is.na(ann$ensembl_gene_id)] <- gene_ids[is.na(ann$ensembl_gene_id)]

# Display label: symbol where we have one, otherwise fall back to the Ensembl ID.
ann$symbol <- ifelse(is.na(ann$external_gene_name) | ann$external_gene_name == "",
                     ann$ensembl_gene_id, ann$external_gene_name)

rowData(se) <- DataFrame(ann)
head(rowData(se))

## DataFrame with 6 rows and 8 columns
##                 ensembl_gene_id external_gene_name   gene_biotype chromosome_name start_position
##                     <character>        <character>    <character>     <character>      <integer>
## ENSG00000000003 ENSG00000000003             TSPAN6 protein_coding               X       99883667
## ENSG00000000005 ENSG00000000005               TNMD protein_coding               X       99839799
## ENSG00000000419 ENSG00000000419               DPM1 protein_coding              20       49551404
## ENSG00000000457 ENSG00000000457              SCYL3 protein_coding               1      169821804
## ENSG00000000460 ENSG00000000460           C1orf112 protein_coding               1      169631245
## ENSG00000000938 ENSG00000000938                FGR protein_coding               1       27938575
##                 end_position      strand      symbol
##                    <integer> <character> <character>
## ENSG00000000003     99894988           -      TSPAN6
## ENSG00000000005     99854882           +        TNMD
## ENSG00000000419     49575092           -        DPM1
## ENSG00000000457    169863408           -       SCYL3
## ENSG00000000460    169823221           +    C1orf112
## ENSG00000000938     27961788           -         FGR

# How much duplication is there in symbols?
sum(duplicated(rowData(se)$symbol))

## [1] 6910

head(sort(table(rowData(se)$symbol), decreasing = TRUE), 5)

## 
##    Y_RNA   snoU13       U3       U6 SNORD112 
##      881      476       87       61       50

length(rowData(se)$symbol)

## [1] 63152

length(unique(rowData(se)$symbol))

## [1] 56242

nrow(se)

## [1] 63152

length(unique(rowData(se)$symbol))

## [1] 56242

length(names(which(table(rowData(se)$symbol) > 1)))

## [1] 2942

#' Row indices of a set of gene symbols, robust to duplicates and misses
#' @return named integer vector (names = symbol)
rows_for_symbols <- function(se, symbols) {
  hit <- match(symbols, rowData(se)$symbol)
  names(hit) <- symbols
  missing <- symbols[is.na(hit)]
  if (length(missing)) message("Not found: ", paste(missing, collapse = ", "))
  hit[!is.na(hit)]
}

rows_for_symbols(se, c("TP53", "EGFR", "NAPSA", "KRT5", "NOT_A_GENE"))

##  TP53  EGFR NAPSA  KRT5 
##  8214  8960  6503 16453

3.4 Alternatives to biomaRt

biomaRt needs the network and the answer can change between sessions. For reproducible pipelines prefer a versioned local annotation package.

3.4.1 org.Hs.eg.db

library(org.Hs.eg.db)
columns(org.Hs.eg.db)

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"    
##  [7] "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"     "GENETYPE"     "GO"          
## [13] "GOALL"        "IPI"          "MAP"          "OMIM"         "ONTOLOGY"     "ONTOLOGYALL" 
## [19] "PATH"         "PFAM"         "PMID"         "PROSITE"      "REFSEQ"       "SYMBOL"      
## [25] "UCSCKG"       "UNIPROT"

sym <- AnnotationDbi::mapIds(org.Hs.eg.db,
                             keys = head(gene_ids, 10),
                             column = "SYMBOL", keytype = "ENSEMBL",
                             multiVals = "first")
sym

## ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000938 
##        "TSPAN6"          "TNMD"          "DPM1"         "SCYL3"         "FIRRM"           "FGR" 
## ENSG00000000971 ENSG00000001036 ENSG00000001084 ENSG00000001167 
##           "CFH"         "FUCA2"          "GCLC"          "NFYA"

# Ensembl -> Entrez, which most pathway tools want:
AnnotationDbi::mapIds(org.Hs.eg.db, keys = head(gene_ids, 10),
                      column = "ENTREZID", keytype = "ENSEMBL", multiVals = "first")

## ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000938 
##          "7105"         "64102"          "8813"         "57147"         "55732"          "2268" 
## ENSG00000000971 ENSG00000001036 ENSG00000001084 ENSG00000001167 
##          "3075"          "2519"          "2729"          "4800"

multiVals deserves attention: "first" silently discards alternatives, "list" keeps them, "filter" drops ambiguous keys entirely. Choose deliberately.

3.4.2 EnsDb (Ensembl, versioned, offline)

library(EnsDb.Hsapiens.v86)
edb <- EnsDb.Hsapiens.v86

head(ensembldb::genes(edb, filter = GeneIdFilter(head(gene_ids, 5)),
                      return.type = "data.frame")[, c("gene_id", "gene_name", "gene_biotype")])

3.4.3 Cache your biomaRt result

# Minimum viable reproducibility: freeze the query result next to the analysis.
saveRDS(gene_annotation, file.path(data_dir, "gene_annotation.rds"))
writeLines(capture.output(sessionInfo()), file.path(res_dir, "sessionInfo_1.txt"))

3.5 Exercise 1

3.5.1 Question

How many genes in se are protein_coding? How many are miRNA?
Build a data frame of all mitochondrial genes (chromosome MT) with their symbols.
TP53 has one row here. Find a gene symbol that maps to more than one Ensembl ID in this object, and explain what could cause that.

3.5.2 Solution

# 1
sort(table(rowData(se)$gene_biotype), decreasing = TRUE)[c("protein_coding", "miRNA")]

## 
## protein_coding          miRNA 
##          22695           3352

# 2
mt <- as.data.frame(rowData(se))[which(rowData(se)$chromosome_name == "MT"),
                                 c("ensembl_gene_id", "symbol", "gene_biotype")]
nrow(mt)

## [1] 37

head(mt, 15)

# 3
dup_sym <- names(which(table(rowData(se)$symbol) > 1))
head(dup_sym)

## [1] "5S_rRNA"     "7SK"         "AADACL2"     "AB019438.63" "AB019438.66" "AB019439.68"

example <- dup_sym[1]
as.data.frame(rowData(se))[rowData(se)$symbol %in% example,
                           c("ensembl_gene_id", "symbol", "chromosome_name", "gene_biotype")]

Causes: Duplicated gene symbols arise because gene symbols are not always unique across the genome annotation. They may occur for genes located on patch or haplotype scaffolds, pseudoautosomal genes present on both X and Y chromosomes, read-through gene models, or related paralogous genes. As a result, multiple Ensembl gene IDs may correspond to the same gene symbol.

4 Understanding the sample metadata

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

4.1 Cleaning the clinical variables

cd <- as.data.frame(colData(se))

# Coerce the obviously-numeric fields
num_candidates <- grep("age|survival|days|time", colnames(cd), value = TRUE)
num_candidates

## [1] "stage_tnm" "age"

for (v in num_candidates) cd[[v]] <- suppressWarnings(as.numeric(cd[[v]]))

# Everything else that has few levels becomes a factor
fct_candidates <- setdiff(colnames(cd), c("geo_accession", "title", num_candidates))
for (v in fct_candidates) {
  if (length(unique(cd[[v]])) <= 12) cd[[v]] <- factor(cd[[v]])
}
str(cd)

## 'data.frame':    216 obs. of  13 variables:
##  $ geo_accession       : chr  "GSM2142443" "GSM2142444" "GSM2142445" "GSM2142446" ...
##  $ title               : chr  "L400T" "L401T" "L404T" "L406T" ...
##  $ source              : Factor w/ 2 levels "Human non-malignant tissue",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ tumor_t_or_normal_n_: chr  "L400T" "L401T" "L404T" "L406T" ...
##  $ stage_tnm           : num  3 5 3 1 5 2 2 4 4 4 ...
##  $ histology           : Factor w/ 3 levels "1","2","3": 2 2 2 1 2 3 3 1 2 2 ...
##  $ surgery_date        : chr  "2006-10-05" "2006-02-28" "2006-05-16" "2006-12-11" ...
##  $ age                 : num  78 74 74 64 72 71 65 57 80 80 ...
##  $ gender              : Factor w/ 2 levels "female","male": 2 2 1 2 1 2 2 1 2 2 ...
##  $ vital_date          : chr  "2010-02-19" "2011-05-01" "2012-09-19" "2013-04-28" ...
##  $ dead                : Factor w/ 3 levels "0","1","n/a": 2 2 2 1 2 2 2 2 1 2 ...
##  $ smoking             : Factor w/ 3 levels "1","2","3": 3 2 2 1 2 1 1 1 2 2 ...
##  $ ps_who              : Factor w/ 3 levels "0","1","2": 2 2 1 1 1 1 1 2 2 1 ...

summary(cd)

##  geo_accession         title                                  source    tumor_t_or_normal_n_
##  Length:216         Length:216         Human non-malignant tissue: 19   Length:216          
##  Class :character   Class :character   Human NSCLC tissue        :197   Class :character    
##  Mode  :character   Mode  :character                                    Mode  :character    
##                                                                                             
##                                                                                             
##                                                                                             
##                                                                                             
##    stage_tnm     histology  surgery_date            age           gender     vital_date       
##  Min.   :1.000   1   : 67   Length:216         Min.   :45.00   female:102   Length:216        
##  1st Qu.:1.000   2   :106   Class :character   1st Qu.:63.00   male  : 95   Class :character  
##  Median :2.000   3   : 24   Mode  :character   Median :69.00   NA's  : 19   Mode  :character  
##  Mean   :2.589   NA's: 19                      Mean   :67.91                                  
##  3rd Qu.:4.000                                 3rd Qu.:74.00                                  
##  Max.   :7.000                                 Max.   :84.00                                  
##  NA's   :19                                    NA's   :19                                     
##    dead     smoking    ps_who   
##  0   :103   1   :95   0   :118  
##  1   : 93   2   :83   1   : 77  
##  n/a :  1   3   :19   2   :  2  
##  NA's: 19   NA's:19   NA's: 19  
##                                 
##                                 
##

# Identify the variables we care about. Adjust the names to the printed output above.
age_var    <- grep("^age",       colnames(cd), value = TRUE)[1]
sex_var    <- grep("gender|sex", colnames(cd), value = TRUE)[1]
hist_var   <- grep("histo",      colnames(cd), value = TRUE)[1]
stage_var  <- grep("stage",      colnames(cd), value = TRUE)[1]
smoke_var  <- grep("smok",       colnames(cd), value = TRUE)[1]
c(age = age_var, sex = sex_var, histology = hist_var, stage = stage_var, smoking = smoke_var)

##         age         sex   histology       stage     smoking 
##       "age"    "gender" "histology" "stage_tnm"   "smoking"

table(cd[[sex_var]],  useNA = "ifany")

## 
## female   male   <NA> 
##    102     95     19

table(cd[[hist_var]], useNA = "ifany")

## 
##    1    2    3 <NA> 
##   67  106   24   19

4.1.1 Verifying the sex label against the data

We read sex off the transcriptome: XIST is expressed almost exclusively in samples with two X chromosomes, and RPS4Y1/DDX3Y/UTY only from a Y chromosome.

cpm_quick <- edgeR::cpm(assay(se, "counts"), log = TRUE, prior.count = 1)

sex_genes <- rows_for_symbols(se, c("XIST", "RPS4Y1", "DDX3Y", "UTY", "KDM5D"))
sex_expr  <- as.data.frame(t(cpm_quick[sex_genes, , drop = FALSE]))
colnames(sex_expr) <- names(sex_genes)

sex_expr$y_score  <- rowMeans(sex_expr[, intersect(c("RPS4Y1", "DDX3Y", "UTY", "KDM5D"),
                                                   colnames(sex_expr)), drop = FALSE])
sex_expr$label    <- cd[[sex_var]]
sex_expr$inferred <- ifelse(sex_expr$y_score > median(range(sex_expr$y_score)), "male", "female")

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

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

##         inferred
## recorded female male
##   female    102    0
##   male        0   95

# Read the direction off the cross-tabulation above, then recode explicitly.
lvl_male <- names(which.max(tapply(sex_expr$y_score, sex_expr$label, mean)))
cd$sex   <- factor(ifelse(as.character(cd[[sex_var]]) == lvl_male, "male", "female"),
                   levels = c("female", "male"))
table(cd$sex)

## 
## female   male 
##    102     95

# Samples where label and biology disagree — flag, do not silently drop.
discordant <- rownames(cd)[as.character(cd$sex) != sex_expr$inferred]
length(discordant); head(discordant)

## [1] 19

## [1] NA NA NA NA NA NA

Discuss. Two clusters should be cleanly separated. Any sample sitting between them is a candidate for sample swap or contamination. What would you do with a discordant sample in your own study?

4.1.2 Verifying histology against marker genes

# Histology diagnosis spring 2013 HB:
# 1=squamous cell cancer (Lung squamous cell carcinoma (LUSC)); 2=AC unspecified (Lung Adenocarcinoma (LUAD)); 3=Large cell/ NOS

mk <- rows_for_symbols(se, c("NAPSA", "NKX2-1", "SFTPC",      # Lung Adenocarcinoma (LUAD)
                             "KRT5", "TP63", "SOX2"))         # Lung squamous cell carcinoma (LUSC)
mk_df <- as.data.frame(t(cpm_quick[mk, , drop = FALSE]))
colnames(mk_df) <- names(mk)
mk_df$histology <- cd[[hist_var]]
table(mk_df$histology)

## 
##   1   2   3 
##  67 106  24

#  1   2   3 
#  67 106  24 

mk_long <- pivot_longer(mk_df, -histology, names_to = "gene", values_to = "log2cpm")
mk_long$gene <- factor(mk_long$gene, levels = c("NAPSA", "NKX2-1", "SFTPC", "KRT5", "TP63", "SOX2"))

ggplot(mk_long, aes(x = factor(histology), y = log2cpm, fill = factor(histology))) +
  geom_boxplot(outlier.size = 0.4) +
  facet_wrap(~ gene, scales = "free_y", nrow = 2) +
  labs(x = "recorded histology code", y = "log2 CPM", fill = "code",
       title = "Which histology code is adenocarcinoma, and which is squamous?") +
  theme(legend.position = "none")

ad_score <- rowMeans(mk_df[, intersect(c("NAPSA", "NKX2-1", "SFTPC"), colnames(mk_df)), drop = FALSE])
sq_score <- rowMeans(mk_df[, intersect(c("KRT5", "TP63", "SOX2"),    colnames(mk_df)), drop = FALSE])
by_code  <- tapply(ad_score - sq_score, cd[[hist_var]], mean)
sort(by_code, decreasing = TRUE)

##         2         3         1 
##  5.389745  1.800665 -3.603480

# Highest adeno-minus-squamous score => adenocarcinoma, lowest => squamous.
code_ad <- names(sort(by_code, decreasing = TRUE))[1]
code_sq <- names(sort(by_code))[1]

cd$histology <- factor(dplyr::case_when(
  as.character(cd[[hist_var]]) == code_ad ~ "LUAD",
  as.character(cd[[hist_var]]) == code_sq ~ "LUSC",
  TRUE ~ "other"), levels = c("LUAD", "LUSC", "other"))
table(cd$histology, useNA = "ifany")

## 
##  LUAD  LUSC other 
##   106    67    43

colData(se) <- DataFrame(cd)

4.2 Describing the cohort

4.2.1 Categorical variables: counts and a (reluctant) pie chart

sex_tab <- as.data.frame(table(sex = colData(se)$sex))
sex_tab$prop <- sex_tab$Freq / sum(sex_tab$Freq)

ggplot(sex_tab, aes(x = "", y = prop, fill = sex)) +
  geom_col(width = 1, colour = "white") +
  coord_polar(theta = "y") +
  geom_text(aes(label = sprintf("%s\n%d (%.0f%%)", sex, Freq, 100 * prop)),
            position = position_stack(vjust = 0.5), size = 4) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Sex distribution of the cohort") +
  theme_void() + theme(legend.position = "none")

Note. Pie charts are hard to read because humans compare angles badly. Reviewers like them; a stacked or grouped bar chart is almost always clearer. Here is the same information, better:

ggplot(as.data.frame(colData(se)), aes(x = histology, fill = sex)) +
  geom_bar(position = position_dodge(preserve = "single")) +
  geom_text(stat = "count", aes(label = after_stat(count)),
            position = position_dodge(width = 0.9), vjust = -0.3, size = 3.2) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "number of samples", title = "Samples per histology and sex")

4.2.2 Continuous variables: box plot with the individual points

cdf <- as.data.frame(colData(se))

ggplot(cdf[!is.na(cdf[[age_var]]), ], aes(x = histology, y = .data[[age_var]], fill = histology)) +
  geom_boxplot(width = 0.5, outlier.shape = NA, alpha = 0.7) +
  geom_jitter(width = 0.15, size = 1.4, alpha = 0.6) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "age at diagnosis (years)",
       title = "Age distribution by histological subtype") +
  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 the sample size, the spread and any bimodality immediately visible.

4.2.3 The right test for each variable type

grp <- droplevels(cdf$histology[cdf$histology %in% c("LUAD", "LUSC")])
sub <- cdf[cdf$histology %in% c("LUAD", "LUSC"), ]

# Categorical vs categorical
tab <- table(sub$sex, droplevels(sub$histology))
tab

##         
##          LUAD LUSC
##   female   68   24
##   male     38   43

chisq.test(tab)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 12.119, df = 1, p-value = 0.000499

fisher.test(tab)          # use when any expected count < 5

## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab
## p-value = 0.0003236
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.617089 6.391729
## sample estimates:
## odds ratio 
##   3.183382

# Continuous vs categorical
by(sub[[age_var]], droplevels(sub$histology), function(x) summary(x))

## droplevels(sub$histology): LUAD
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   45.00   62.00   68.50   67.72   74.00   83.00 
## --------------------------------------------------------------------------- 
## droplevels(sub$histology): LUSC
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   46.00   62.50   67.00   67.64   73.50   84.00

shapiro.test(sub[[age_var]][sub$histology == "LUAD"])   # normality check first

## 
##  Shapiro-Wilk normality test
## 
## data:  sub[[age_var]][sub$histology == "LUAD"]
## W = 0.97605, p-value = 0.052

t.test(sub[[age_var]] ~ droplevels(sub$histology))

## 
##  Welch Two Sample t-test
## 
## data:  sub[[age_var]] by droplevels(sub$histology)
## t = 0.05992, df = 137.93, p-value = 0.9523
## alternative hypothesis: true difference in means between group LUAD and group LUSC is not equal to 0
## 95 percent confidence interval:
##  -2.406035  2.556415
## sample estimates:
## mean in group LUAD mean in group LUSC 
##           67.71698           67.64179

wilcox.test(sub[[age_var]] ~ droplevels(sub$histology))

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  sub[[age_var]] by droplevels(sub$histology)
## W = 3568, p-value = 0.959
## alternative hypothesis: true location shift is not equal to 0

#' Minimal "Table 1": one row per variable, correct test chosen by type
make_table_one <- function(df, group, vars) {
  g <- droplevels(factor(df[[group]]))
  do.call(rbind, lapply(vars, function(v) {
    x <- df[[v]]
    if (is.numeric(x)) {
      p <- tryCatch(wilcox.test(x ~ g)$p.value, error = function(e) NA_real_)
      summ <- paste(tapply(x, g, function(z)
        sprintf("%.1f (%.1f-%.1f)", median(z, na.rm = TRUE),
                quantile(z, .25, na.rm = TRUE), quantile(z, .75, na.rm = TRUE))),
        collapse = " | ")
      data.frame(variable = v, type = "median (IQR)", summary = summ, p = p)
    } else {
      tb <- table(x, g)
      p <- tryCatch(fisher.test(tb, simulate.p.value = TRUE)$p.value,
                    error = function(e) NA_real_)
      summ <- paste(apply(tb, 1, paste, collapse = "/"), collapse = "; ")
      data.frame(variable = v, type = "n per level", summary = summ, p = p)
    }
  }))
}

vars <- c(age_var, "sex", stage_var, smoke_var)
vars <- vars[!is.na(vars) & vars %in% colnames(sub)]
make_table_one(sub, "histology", vars)

Discuss. If sex is unbalanced between LUAD and LUSC (it usually is, because smoking history differs), then sex and histology are confounded. On Day 2 we handle that by putting sex in the design matrix — but only a table like this tells you that you need to.

5 Normalisation: CPM and TPM

Raw counts are not comparable, for two reasons that are easy to confuse:

Library size (sequencing depth). A sample sequenced to 40 M reads gives roughly twice the count of the same sample at 20 M. This affects every comparison across samples.
Gene length. A 10 kb transcript generates more fragments than a 1 kb transcript at the same molar concentration. This affects comparisons between genes, not between samples.

5.1 CPM — counts per million

\[\text{CPM}_{ij} = \frac{c_{ij}}{\sum_k c_{kj}} \times 10^6\]

CPM corrects for library size only. Gene length cancels out when you compare the same gene across samples, which is exactly what differential expression does — hence CPM is the natural scale for cross-sample views.

counts_mat <- assay(se, "counts")
lib_size   <- colSums(counts_mat)

summary(lib_size / 1e6)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   23.99   49.08   54.27   55.78   60.63   95.87

hist(lib_size / 1e6, breaks = 40, col = "grey80",
     xlab = "library size (million reads)", main = "Sequencing depth per sample")

cpm_manual <- t(t(counts_mat) / lib_size) * 1e6
cpm_edger  <- edgeR::cpm(counts_mat)
all.equal(cpm_manual[1:100, 1:5], cpm_edger[1:100, 1:5])

## [1] TRUE

colSums(cpm_manual)[1:5]         # 1e6 by construction

## L400T L401T L404T L406T L413T 
## 1e+06 1e+06 1e+06 1e+06 1e+06

5.1.1 Log transformation and the prior count

log2(CPM) makes RNA-seq expression data more manageable, less skewed, and easier to interpret. We use log2(CPM) because it: Compresses very large expression values, reducing skew. Makes the expression distribution easier to visualize and compare. Makes fold changes easier to interpret: a 2-fold change corresponds to a difference of 1 on the log2 scale.

## 
log_cpm <- edgeR::cpm(counts_mat, log = TRUE, prior.count = 2)

par(mfrow = c(1, 2))
hist(cpm_manual[, 1], breaks = 100, main = "CPM (sample 1)", xlab = "CPM", col = "grey80")
hist(log_cpm[, 1],    breaks = 100, main = "log2 CPM (sample 1)", xlab = "log2 CPM", col = "grey80")

par(mfrow = c(1, 1))

Adding a prior.count before the log does two things: it avoids log(0) = -Inf, and it shrinks the wild variance of low-count genes. Without it, a gene going 0 → 1 count looks like an infinite fold change.

5.1.2 Filtering low-expressed genes

We filter out lowly expressed genes because they provide little reliable information for differential expression analysis and can add noise.

dge <- DGEList(counts = counts_mat, samples = as.data.frame(colData(se)))
summary(rowSums(cpm(dge) > 1))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00   56.76  129.25  216.00

# filterByExpr() is specifically designed to remove genes with too little expression to provide reliable statistical evidence in differential expression analysis.
keep <- filterByExpr(dge, group = colData(se)$histology) 
table(keep)

## keep
## FALSE  TRUE 
## 38777 24375

dge_filt <- dge[keep, , keep.lib.sizes = FALSE]
se_filt  <- se[keep, ]

nrow(se); nrow(se_filt)

## [1] 63152

## [1] 24375

plot_density <- function(m, main) {
  plot(density(m[, 1]), ylim = c(0, 0.35), xlab = "log2 CPM", main = main, lwd = 1.5,
       col = "grey30")
  for (i in 2:min(30, ncol(m))) lines(density(m[, i]), col = alpha_col(i))
}
alpha_col <- function(i) adjustcolor(RColorBrewer::brewer.pal(8, "Dark2")[(i %% 8) + 1], 0.5)

par(mfrow = c(1, 2))

plot_density(
  edgeR::cpm(dge, log = TRUE),
  "Before filtering"
)

plot_density(
  edgeR::cpm(dge_filt, log = TRUE),
  "After filtering"
)

par(mfrow = c(1, 1))

The “before filtering” panel has a large peak at low log-CPM: genes that are essentially never detected. They add nothing but multiple-testing burden and noise.

Discuss. filterByExpr keeps genes with enough counts in at least as many samples as the smallest group. Why is a fixed cutoff like “CPM > 1 in all samples” a bad idea when group sizes are unequal?

5.2 TPM — transcripts per million

TPM corrects for both length and depth, and does the divisions in the other order:

\[r_{ij} = \frac{c_{ij}}{\ell_i}, \qquad \text{TPM}_{ij} = \frac{r_{ij}}{\sum_k r_{kj}} \times 10^6\]

Dividing by length first makes the per-sample values proportional to molar concentration; the columns then sum to exactly 10^6, so the value is a genuine proportion of the transcript pool. (FPKM normalises by depth first and by length second, which is why FPKM columns do not sum to a constant and FPKM is not comparable across samples.)

5.2.1 Getting a gene length

The correct length for gene-level counts is the union exon length: merge all exons of all transcripts of a gene and sum the widths.

library(GenomicRanges)

len_file <- file.path(data_dir, "gene_lengths.rds")
if (file.exists(len_file)) {
  gene_len <- readRDS(len_file)
} else {
  ex <- ensembldb::exonsBy(EnsDb.Hsapiens.v86, by = "gene")
  gene_len <- sum(width(GenomicRanges::reduce(ex)))       # named integer vector
  saveRDS(gene_len, len_file)
}
head(gene_len)

## ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000938 
##            4535            1610            1207            6883            5967            3474

summary(as.numeric(gene_len))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       8     415     987    2320    3114  205012

The same thing with biomaRt (works for any Ensembl species)

# ex_bm <- getBM(attributes = c("ensembl_gene_id", "chromosome_name",
#                               "exon_chrom_start", "exon_chrom_end"),
#                filters = "ensembl_gene_id",
#                values  = gene_ids,      # ~1.2 M rows for a whole transcriptome: be patient
#                mart    = ensembl)
# 
# gr <- GRanges(seqnames = ex_bm$chromosome_name,
#               ranges   = IRanges(ex_bm$exon_chrom_start, ex_bm$exon_chrom_end),
#               gene     = ex_bm$ensembl_gene_id)
# gene_len_bm <- sum(width(GenomicRanges::reduce(split(gr, gr$gene))))

#' Counts -> TPM
#' @param counts integer matrix, genes x samples
#' @param len    numeric vector of gene lengths in bp, same order as rows of counts
counts_to_tpm <- function(counts, len) {
  stopifnot(nrow(counts) == length(len), all(len > 0))
  rate <- counts / (len / 1000)                    # reads per kilobase
  t(t(rate) / colSums(rate)) * 1e6
}

ids_filt <- sub("\\..*$", "", rownames(se_filt))
have_len <- ids_filt %in% names(gene_len)
mean(have_len)

## [1] 0.9693949

se_tpm <- se_filt[have_len, ]
len_vec <- as.numeric(gene_len[sub("\\..*$", "", rownames(se_tpm))])

tpm <- counts_to_tpm(assay(se_tpm, "counts"), len_vec)
assay(se_tpm, "tpm") <- tpm

colSums(tpm)[1:5]        # exactly 1e6

## L400T L401T L404T L406T L413T 
## 1e+06 1e+06 1e+06 1e+06 1e+06

5.2.2 CPM and TPM disagree — and that is the point

cpm_sub <- edgeR::cpm(assay(se_tpm, "counts"))
s <- 1

comp <- data.frame(
  cpm    = cpm_sub[, s],
  tpm    = tpm[, s],
  length = len_vec
)
comp <- comp[comp$cpm > 1, ]

ggplot(comp, aes(x = log2(cpm), y = log2(tpm), colour = log10(length))) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  scale_colour_viridis_c() +
  labs(x = "log2 CPM", y = "log2 TPM", colour = "log10\ngene length",
       title = sprintf("CPM vs TPM, sample %s", colnames(tpm)[s]),
       subtitle = "Long genes move down, short genes move up")

5.3 When to use what

Task	Use	Why
Differential expression	raw counts + DESeq2/edgeR size factors	The models need the mean–variance relationship of counts; pre-normalised values break it
Clustering, PCA, heatmaps, quality check	log2 CPM, or `vst`/`rlog`	Variance stabilised, library size corrected
Comparing one gene across samples	CPM or TPM	Both fine; length is constant within the gene
Comparing genes within one sample (“is A higher than B here?”)	TPM	Only TPM removes the length effect
Deconvolution, signature scoring, cross-study comparison	TPM	Method expects concentration-like, column-normalised values
Anything	avoid Fragments Per Kilobase of transcript per Million mapped reads(FPKM)	Less suitable for cross-sample comparison; mainly encountered in older or legacy pipelines

Caveat on TPM across samples. TPM columns sum to a constant, so if a handful of genes dominate one sample (haemoglobin in blood, immunoglobulin in a plasma-cell-rich tumour), every other gene’s TPM is depressed. This composition effect is exactly what TMM (edgeR) and median-of-ratios (DESeq2) are designed to correct, and neither CPM nor TPM does it.

5.4 Exercise 2

5.4.1 Question

Which sample has the smallest library size, and how many times smaller is it than the largest?
Pick a long gene and a short gene that have similar CPM in sample 1. Compare their TPM. Explain the direction of the change.
If your collaborator sends you a matrix of FPKM values, can you convert them to TPM? Write the function.

5.4.2 Solution

# 1
ls_mb <- sort(colSums(assay(se, "counts"))) / 1e6
c(smallest = ls_mb[1], largest = ls_mb[length(ls_mb)], ratio = ls_mb[length(ls_mb)] / ls_mb[1])

## smallest.L563T  largest.L563N    ratio.L563N 
##      23.986003      95.873366       3.997055

# 2
df <- data.frame(cpm = cpm_sub[, 1], tpm = tpm[, 1], len = len_vec, symbol = rowData(se_tpm)$symbol)
df <- df[df$cpm > 5 & df$cpm < 6, ]
df <- df[order(df$len), ]
rbind(head(df, 2), tail(df, 2))

# 3  FPKM is already length-normalised; you only need to rescale the columns.
fpkm_to_tpm <- function(fpkm) t(t(fpkm) / colSums(fpkm)) * 1e6

TPM corrects for gene length, so when genes have similar CPM, short genes get higher TPM values while long genes get lower TPM values. This is because the same number of reads from a short gene indicates a higher RNA molecule concentration than the same number of reads from a long gene.

6 Visualising genes of interest

Two plots cover most of what you will ever need to show: a per-gene box/dot plot for a handful of genes, and a heatmap for a gene list.

6.1 Reshaping to long format once

#' Extract genes by symbol as a tidy data frame
#' @param se     SummarizedExperiment
#' @param genes  character vector of symbols
#' @param assay  which assay to pull ("counts", "tpm", ...)
#' @param log    log2-transform (with pseudocount) after extraction
gene_long <- function(se, genes, assay = "logcpm", log = FALSE) {
  idx <- rows_for_symbols(se, genes)
  m   <- assays(se)[[assay]][idx, , drop = FALSE]
  rownames(m) <- names(idx)
  if (log) m <- log2(m + 1)
  df <- as.data.frame(t(m))
  df$sample <- rownames(df)
  out <- pivot_longer(df, -sample, names_to = "gene", values_to = "expression")
  cbind(out, as.data.frame(colData(se))[out$sample, , drop = FALSE])
}

assay(se_tpm, "logcpm") <- edgeR::cpm(assay(se_tpm, "counts"), log = TRUE, prior.count = 2)

6.2 Box plot with individual points

nsclc_markers <- c("NAPSA", "NKX2-1", "SFTPC",     # adeno
                   "KRT5", "TP63", "SOX2",         # squamous
                   "MKI67", "EGFR", "CD274")  # proliferation, growth signaling, immune checkpoint

dat <- gene_long(se_tpm, nsclc_markers, assay = "logcpm")
dat <- dat[dat$histology %in% c("LUAD", "LUSC"), ]
dat$histology <- droplevels(dat$histology)

pvals <- dat %>%
  group_by(gene) %>%
  summarise(p = wilcox.test(expression ~ histology)$p.value, .groups = "drop") %>%
  mutate(padj = p.adjust(p, method = "BH"),
         label = ifelse(padj < 0.001, "FDR < 0.001", sprintf("FDR = %.3g", padj)))
pvals

ggplot(dat, aes(x = histology, y = expression, fill = histology)) +
  geom_boxplot(width = 0.55, outlier.shape = NA, alpha = 0.75) +
  geom_jitter(width = 0.16, size = 0.7, alpha = 0.45) +
  geom_text(
    data = pvals,
    aes(x = 1.5, y = Inf, label = label),
    inherit.aes = FALSE,
    vjust = 1.4,
    size = 3
  ) +
  facet_wrap(~ factor(gene,levels = c("NAPSA", "NKX2-1", "SFTPC","KRT5", "TP63", "SOX2","MKI67", "EGFR", "CD274")), 
                      scales = "free_y", ncol = 3) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "log2 CPM", title = "Classic NSCLC lineage markers") +
  theme(legend.position = "none")

Note. These p-values are for description, not discovery: we picked the genes because we expected them to differ. Genome-wide, properly modelled testing is Day 2’s job. Reporting a Wilcoxon p-value on hand-picked genes as if it were a finding is a common and avoidable error.

6.3 Heatmap of a gene list

gene_set <- c(nsclc_markers,
              "KRT6A", "KRT14", "DSG3", "PKP1",           # squamous programme
              "SFTPB", "SFTPA1", "MUC1", "CEACAM6",       # alveolar / adeno programme
              "TP53", "CDKN2A", "RB1", "KEAP1", "STK11",  # commonly mutated
              "PTPRC", "CD3E", "CD8A", "FOXP3", "COL1A1") # microenvironment

idx <- rows_for_symbols(se_tpm, gene_set)
mat <- assay(se_tpm, "logcpm")[idx, ]
rownames(mat) <- names(idx)

# Z-score per gene: without this, highly expressed genes dominate the colour scale
mat_z <- t(scale(t(mat)))

keep_samp <- colData(se_tpm)$histology %in% c("LUAD", "LUSC")
mat_z <- mat_z[, keep_samp]

ann_col <- data.frame(
  histology = droplevels(colData(se_tpm)$histology[keep_samp]),
  sex       = colData(se_tpm)$sex[keep_samp],
  row.names = colnames(mat_z)
)
ann_colors <- list(
  histology = c(LUAD = "#66C2A5", LUSC = "#FC8D62"),
  sex       = c(female = "#8DA0CB", male = "#E78AC3")
)

pheatmap(mat_z,
         annotation_col = ann_col,
         annotation_colors = ann_colors,
         show_colnames = FALSE,
         clustering_distance_rows = "correlation",
         clustering_distance_cols = "euclidean",
         clustering_method = "ward.D2",
         breaks = seq(-3, 3, length.out = 101),
         color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100),
         fontsize_row = 8,
         main = "Selected NSCLC genes (row z-scored log2 CPM)")

Choices that change what a heatmap says — make them deliberately:

Scaling. Row z-scores show relative patterns; unscaled values show absolute level. Say which you used in the legend.
Colour saturation. Capping the scale at ±3 SD stops two outlier samples from washing out the whole map.
Distance and linkage. Correlation distance groups genes with a similar shape; Euclidean groups by magnitude. ward.D2 gives compact clusters; complete is more chain-prone.
Diverging palette, centred at zero — for z-scores, always. A sequential palette on centred data hides the sign.

The same figure with ComplexHeatmap, for publication-grade control

library(ComplexHeatmap)
library(circlize)

ha <- HeatmapAnnotation(
  histology = ann_col$histology,
  sex       = ann_col$sex,
  col = ann_colors
)

Heatmap(mat_z,
        name = "z-score",
        top_annotation = ha,
        col = colorRamp2(c(-3, 0, 3), c("#2166AC", "white", "#B2182B")),
        show_column_names = FALSE,
        clustering_method_rows = "ward.D2",
        row_names_gp = gpar(fontsize = 8))

7 Principal component analysis

PCA answers one question before you do any statistics: what is the dominant structure in my data, and does it correspond to my biological variable or to something else?

7.1 Doing it properly

log_cpm_f <- assay(se_tpm, "logcpm")

# 1. Use the most variable genes: uninformative genes add noise, not signal.
vars_gene <- rowVars(log_cpm_f)
top_n     <- 1000
top_genes <- order(vars_gene, decreasing = TRUE)[seq_len(top_n)]

# 2. prcomp expects samples in rows. Center always; scale only if you want every
#    gene to contribute equally regardless of its variance.
pca <- prcomp(t(log_cpm_f[top_genes, ]), center = TRUE, scale. = FALSE)

pct <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1)
head(pct, 10)

##  [1] 22.3  7.1  6.1  4.6  3.8  3.3  2.5  2.1  1.8  1.7

scree <- data.frame(PC = seq_along(pct)[1:15], var = pct[1:15])
ggplot(scree, aes(PC, var)) +
  geom_col(fill = "steelblue") +
  geom_line(colour = "grey30") + geom_point() +
  scale_x_continuous(breaks = 1:15) +
  labs(y = "% variance explained", title = "Scree plot")

scores <- as.data.frame(pca$x[, 1:5])
scores <- cbind(scores, as.data.frame(colData(se_tpm)))

p1 <- ggplot(scores, aes(PC1, PC2, colour = histology)) +
  geom_point(size = 2, alpha = 0.85) +
  stat_ellipse(level = 0.68) +
  scale_colour_brewer(palette = "Set2") +
  labs(x = sprintf("PC1 (%.1f%%)", pct[1]), y = sprintf("PC2 (%.1f%%)", pct[2]))

p2 <- ggplot(scores, aes(PC1, PC2, colour = sex)) +
  geom_point(size = 2, alpha = 0.85) +
  scale_colour_brewer(palette = "Dark2") +
  labs(x = sprintf("PC1 (%.1f%%)", pct[1]), y = sprintf("PC2 (%.1f%%)", pct[2]))

gridExtra::grid.arrange(p1, p2, nrow = 1)

7.2 Which metadata variable drives which component?

Eyeballing colours does not scale past three variables. Test each PC against each variable:

test_vars <- intersect(c("histology", "sex", age_var, stage_var, smoke_var),colnames(colData(se_tpm)))
pcs <- 1:6

assoc <- sapply(test_vars, function(v) {
  x <- colData(se_tpm)[[v]]
  sapply(pcs, function(i) {
    ok <- !is.na(x)
    if (length(unique(x[ok])) < 2) return(NA_real_)
    fit <- lm(pca$x[ok, i] ~ x[ok])
    a   <- anova(fit)
    a$`Pr(>F)`[1]
  })
})
rownames(assoc) <- paste0("PC", pcs)
signif(assoc, 3)

##     histology      sex    age stage_tnm smoking
## PC1  8.32e-72 5.63e-05 0.8760     0.591  0.2250
## PC2  5.08e-03 5.42e-03 0.3940     0.121  0.8590
## PC3  1.87e-12 1.17e-14 0.7220     0.718  0.3920
## PC4  1.50e-03 1.28e-06 0.8890     0.770  0.4790
## PC5  4.99e-01 1.23e-21 0.5130     0.166  0.3050
## PC6  3.29e-01 7.32e-07 0.0647     0.882  0.0708

pheatmap(-log10(assoc), cluster_rows = FALSE, cluster_cols = FALSE,
         display_numbers = matrix(sprintf("%.1f", -log10(assoc)), nrow = nrow(assoc)),
         color = colorRampPalette(c("white", "firebrick"))(50),
         main = "-log10 p, association of each PC with each covariate")

Discuss. Suppose PC1 is strongly associated with histology — good, your biology dominates. But if PC1 were associated with library size or with a processing date, you would have a technical driver larger than your effect of interest, and you would need to model it. Always add library size to this table:

libsize <- colSums(assay(se_tpm, "counts"))
sapply(pcs, function(i) cor(pca$x[, i], log10(libsize), method = "spearman"))

## [1] -0.10957263  0.07505924 -0.04916467 -0.12211862  0.12458591  0.02952048

7.3 Sample–sample distances as a complement

set.seed(1)
sel <- sample(ncol(log_cpm_f), min(60, ncol(log_cpm_f)))   # subsample for legibility
d <- dist(t(log_cpm_f[top_genes, sel]))
mat_d <- as.matrix(d)

pheatmap(mat_d,
         clustering_distance_rows = d,
         clustering_distance_cols = d,
         annotation_col = data.frame(histology = colData(se_tpm)$histology[sel],
                                     row.names = colnames(log_cpm_f)[sel]),
         show_rownames = FALSE, show_colnames = FALSE,
         color = colorRampPalette(brewer.pal(8, "Set2"))(100),
         main = "Euclidean distance between samples")

An outlying sample shows up here as a row/column that is far from everything, including its own group. Combined with the PCA and the sex check above, this is your sample-level QC.

7.4 Exercise 3

7.4.1 Question

Re-run the PCA using all genes instead of the top 1000. Does the variance explained by PC1 go up or down, and why?
Plot the same PCA, coloring samples by source, smoking, or another variable in the scores dataframe.

7.4.2 Solution

# 1
pca_all <- prcomp(t(log_cpm_f), center = TRUE, scale. = FALSE)
pct_all <- round(100 * pca_all$sdev^2 / sum(pca_all$sdev^2), 1)
c(top1000_PC1 = pct[1], allgenes_PC1 = pct_all[1])

##  top1000_PC1 allgenes_PC1 
##         22.3         13.4

#Variable-gene selection makes the biological structure easier to see.

# 2
ggplot(scores, aes(PC1, PC2, colour =smoking)) +
  geom_point(size = 2, alpha = 0.85) +
  scale_colour_brewer(palette = "Dark2") +
  labs(x = sprintf("PC1 (%.1f%%)", pct[1]), y = sprintf("PC2 (%.1f%%)", pct[2]))

The percentage usually drops: adding thousands of low-variance genes increases total variance without adding structure, so PC1’s share shrinks. The biological interpretation of PC1 typically stays the same.

8 Save the objects for Day 2

saveRDS(se,      file.path(data_dir, "se_raw.rds"))
saveRDS(se_filt, file.path(data_dir, "se_filtered.rds"))
saveRDS(se_tpm,  file.path(data_dir, "se_tpm.rds"))

writeLines(capture.output(sessionInfo()), file.path(res_dir, "sessionInfo_day1.txt"))

Take-home messages

Public data is easy to download and hard to trust. Budget most of your time for reconciling counts with metadata, not for the statistics.
Keep stable identifiers as keys; symbols are for humans, and only at the end.
Every categorical level you did not personally define should be verified against the data (sex via XIST/Y genes, subtype via markers, tissue via lineage genes).
CPM corrects depth; TPM corrects depth and length; neither corrects composition, and neither belongs inside a differential expression model.
PCA is a diagnostic, not a result. If the largest component is not your biology, find out what it is before going further.

Försöksdjursvetenskap (teoretisk) och avancerad databearbetning

Finding public data, annotating genes, and exploring counts and metadata

Birong Zhang

07 September 2026