Försöksdjursvetenskap (teoretisk) och avancerad databearbetning
Lab 6 — Public RNA-seq data in R: download, gene annotation and cohort description
Birong Zhang
16 October 2026
Before we start
0.1 Introduction
Thousands of RNA-seq experiments on human patients are already public. Re-using them is one of the most direct ways to answer a biological question without a new animal experiment. Today you will take one such study from the public database GEO and turn it into clean objects you can analyse.
The study is GSE161731: blood samples from patients with COVID_19, influenza, a common-cold coronavirus, or a bacterial infection, plus healthy volunteers. We will
- Get the data — download the table of gene counts and the table that describes the patients, and make sure the two tables match;
- Annotate genes — replace cryptic Ensembl IDs such
as
ENSG00000165949with gene names you recognise, such as IFI27; - Describe the cohort — clean the patient information, check it against the RNA itself, and summarise who is in each group.
At the end we save everything, so Lab 7 can start straight away with PCA and differential expression.
How each part works. Every part ends with an Exercise (questions 1–3, then a Solution tab — try before you peek) followed by an Assignment. The assignment is almost the same as the exercise: you copy the solution code and change a small thing (a gene name, a colour, a group). If you can do the exercise, you can do the assignment.
You do not need to be a programmer. Every line of code has a comment. If a line does not make sense, run it alone and look at what comes out — that is the fastest way to learn.
| Part | Topic | Ends with |
|---|---|---|
| 1 | Get the data: download GSE161731 into R | Exercise 1 → Assignment 12 |
| 2 | Annotate genes: from Ensembl IDs to gene symbols | Exercise 2 → Assignment 13 |
| 3 | Describe the cohort: understanding the sample metadata | Exercise 3 → Assignment 14 |
| — | Save the objects for Lab 7 |
0.2 Packages
# Run this ONCE, before the session. It takes a while.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c(
"GEOquery", # download from GEO
"org.Hs.eg.db", # human gene names (works offline)
"AnnotationDbi", # the tool to look things up in org.Hs.eg.db
"edgeR" # today only used for a quick sex check
))
install.packages(c("ggplot2", "dplyr", "RColorBrewer"))library(GEOquery) # getGEOSuppFiles()
library(org.Hs.eg.db) # human gene annotation database
library(AnnotationDbi) # mapIds() to look up gene names
library(edgeR) # cpm()
library(ggplot2) # all plots
library(dplyr) # small table summaries
theme_set(theme_bw(base_size = 12)) # one plot style for the whole document1 Get the data: download GSE161731 into R
1.1 Where does public RNA-seq data live?
The same experiment is stored at three levels of processing. Knowing which one you download is the whole game.
| Level | What it is | Where | Size |
|---|---|---|---|
| Raw reads | FASTQ files: the sequences from the machine | SRA / ENA | 1–10 GB per sample |
| Aligned reads | BAM files: reads placed on the genome | rarely public | large |
| Processed data | count table: how many reads each gene got in each sample | GEO supplementary files | a few MB |
Today we take the processed level: a gene × sample table of raw counts. Each number says “this many RNA fragments from this gene were found in this blood sample”. That is exactly what DESeq2 needs in Lab 7.
The databases you will meet:
- GEO (NCBI, USA):
GSE…= a study,GSM…= one sample. - ArrayExpress / BioStudies (EMBL-EBI, Europe): the European equivalent.
- SRA / ENA: the raw reads.
- recount3, ARCHS4: many studies re-processed in the same way.
- TCGA / GTEx: large human cohorts (tumours, and normal tissues).
Re-use data responsibly: cite the accession number and the paper, respect patient consent (some human data is controlled access), and read the methods.
1.2 Why GSE161731?
McClain et al. (2021), Nature Communications: whole-blood RNA-seq from patients who came to the emergency department at Duke University with an acute respiratory infection, plus healthy controls.
- It is human blood. Blood can be sampled again and again in patients — a realistic alternative to an animal infection model.
- It has five groups: COVID_19, Influenza, other (seasonal) coronavirus, bacterial infection, and healthy. The real clinical question is which infection, not just “sick or healthy”.
- The expected biology is well known. Viruses switch on interferon genes (IFI27, IFI44L, ISG15, RSAD2); bacteria switch on neutrophil genes (S100A8, S100A9, MMP8, CD177). If our analysis is right, these genes must show up. That is a positive control.
- The patient table is realistically messy: age
written as
>89, some people sampled twice, groups of different sizes.
1.3 Downloading the files
GEO keeps the authors’ processed files as “supplementary files”. For GSE161731 there are two: the count table and a count key (the patient table).
# Change this to YOUR folder. Everything is downloaded once and kept here.
setwd("/Users/birzh586/Desktop/RNA-seq course")
data_dir <- "data" # downloads and saved objects go here
dir.create(data_dir, showWarnings = FALSE) # showWarnings = FALSE: no message if it exists
# Download every supplementary file of the study into data/GSE161731/
supp <- getGEOSuppFiles("GSE161731", baseDir = data_dir, makeDirectory = TRUE)
list.files(file.path(data_dir, "GSE161731")) # which files did we get?## [1] "GSE161731_counts_key.csv.gz" "GSE161731_counts.csv.gz" "GSE161731_key.csv.gz"
## [4] "GSE161731_xpr_nlcpm.csv.gz" "GSE161731_xpr_tpm_geo.txt.gz"# The two file paths. file.path() glues folder and file names together.
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 can read .gz without unzipping).
# row.names = 1 -> the first column (gene ID / sample ID) becomes the row names
# check.names = FALSE -> keep the column names exactly as they are
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) # rows = genes, columns = samples## [1] 60675 201dim(key_raw) # rows = samples, columns = patient information## [1] 198 8counts_raw[1:5, 1:4] # top-left corner of the count tablehead(key_raw) # first rows of the patient tableTwo quick checks to do with every table you download:
# 1. Are these whole numbers? Raw counts must be integers (you cannot find half a read).
all(counts_raw == round(counts_raw)) # TRUE = raw counts## [1] TRUE# 2. How many reads does each sample have in total? This is the "library size".
head(colSums(counts_raw))## 94189 DU09-03S0000604 DU09-03S0000611 105920 DU09-03S0000774 DU09-03S0000775
## 52431553 50351794 52172232 53276313 43639732 56439096# 3. What do the gene names look like? These are 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 table and the rows of the patient table must be the same samples, in the same order. If they are shifted by one, every result afterwards is wrong — and R will not warn you. So we check.
# The sample names are written slightly differently in the two files.
# We make them identical: add an "S" in front, and turn "-" and "_" into ".".
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))
# intersect(): the sample names found in BOTH tables.
common <- intersect(colnames(counts_raw), rownames(key_raw))
length(common)## [1] 198# Keep only those samples, sorted by cohort, and put the counts in the SAME order.
meta <- key_raw[common, ]
meta <- meta[order(meta$cohort), ]
counts <- counts_raw[, rownames(meta)]
# The most important check of the day. It must print TRUE.
all(colnames(counts) == rownames(meta))## [1] TRUEstr(meta) # which patient variables do we have, and of which type?## '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 group##
## Bacterial CoV other COVID-19 healthy Influenza
## 24 61 77 19 171.5 Exercise 1
1.5.1 Question
- How many genes and how many samples are in
counts? - Which sample has the largest library size?
- How many samples were in the count file but not in the patient table?
1.5.2 Solution
# 1 dim() gives: number of rows (genes), number of columns (samples)
dim(counts)## [1] 60675 198# 2
lib <- colSums(counts) # total reads per sample
names(which.max(lib)) # which.max() finds the position; names() gives the sample## [1] "S434325"# 3 setdiff(a, b): things in a that are not in b
length(setdiff(colnames(counts_raw), rownames(key_raw))) # in counts, not in the key## [1] 31.6 Assignment 12 — Get the data
Hint: start from the code of Exercise 1 and change one small thing.
- Report the number of samples and variables in
meta. - Which sample has the smallest library size, and how
many reads does it have? (Hint:
which.min(),min()) - How many samples were in the patient table but not in the count file?
2 Annotate genes: from Ensembl IDs to gene symbols
2.1 Why do we need to translate gene names?
ENSG00000165949 is precise but unreadable.
IFI27 is readable, but gene symbols sometimes change, and a few
symbols belong to more than one Ensembl ID. Think of Ensembl IDs as
personal ID numbers and symbols as first
names: the number is unique, the name is what people
recognise.
We translate once, early, so that later we can simply write
counts["IFI27", ].
2.2 Looking up gene symbols
org.Hs.eg.db is a human gene “dictionary” stored on your
computer. mapIds() looks up each Ensembl ID and gives back
one answer per gene.
gene_ids <- rownames(counts) # the Ensembl IDs we have
head(gene_ids)## [1] "ENSG00000223972" "ENSG00000227232" "ENSG00000278267" "ENSG00000243485" "ENSG00000274890"
## [6] "ENSG00000237613"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"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 there are several answers, take the first
)
# The same lookup, now asking for the gene type (protein-coding, ncRNA, ...)
gene_biotype <- mapIds(org.Hs.eg.db,
keys = gene_ids,
keytype = "ENSEMBL",
column = "GENETYPE",
multiVals = "first")
# How well did it work?
sum(!is.na(gene_symbol)) # IDs that got a symbol## [1] 35529sum(is.na(gene_symbol)) # IDs with no symbol (NA)## [1] 25146sum(duplicated(gene_symbol[!is.na(gene_symbol)])) # symbols used more than once## [1] 82Notes Some symbols appear on more than one row. For each symbol we keep the row with the most reads on average — the version of the gene that is really expressed in blood — and drop the others. Genes with no symbol are dropped too.
2.3 One row per gene symbol
Some symbols appear on more than one row. For each symbol we keep the row with the most reads on average — the version of the gene that is really expressed in blood — and drop the others. Genes with no symbol are dropped too.
# 1. Sort all rows from highest to lowest average count.
o <- order(rowMeans(counts), decreasing = TRUE)
counts <- counts[o, ]
gene_ids <- gene_ids[o]
gene_symbol <- gene_symbol[o]
gene_biotype <- gene_biotype[o]
# 2. Keep a row if it HAS a symbol AND that symbol was not already seen higher up.
keep_row <- !is.na(gene_symbol) & !duplicated(gene_symbol)
sum(keep_row) # number of genes we keep## [1] 35447counts <- 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 gene symbols
# 3. A small annotation table: one row per gene, same order as counts.
ann <- data.frame(symbol = rownames(counts),
ensembl = gene_ids,
biotype = gene_biotype,
row.names = rownames(counts))
head(ann)# Which kinds of genes are in the table?
sort(table(ann$biotype), decreasing = TRUE)##
## protein-coding pseudo ncRNA snoRNA other snRNA
## 19208 8893 6155 726 368 69
## rRNA scRNA
## 24 4# Are the genes we care about present?
c("IFI27", "ISG15", "S100A8", "XIST", "RPS4Y1") %in% ann$symbol## [1] TRUE TRUE TRUE TRUE TRUEBiology note. Most blood studies focus on protein-coding genes: they are the best annotated and the easiest to interpret. Non-coding RNAs are real, but harder to link to a function.
2.4 Exercise 2
2.4.1 Question
- How many genes are protein-coding, and what proportion of all genes is that?
- Find all gene 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 == compares every value; sum() counts the TRUEs; mean() gives the proportion
sum(ann$biotype == "protein-coding", na.rm = TRUE)## [1] 19208mean(ann$biotype == "protein-coding", na.rm = TRUE)## [1] 0.5418794# 2 grep(): find names matching a pattern; "^" means "starts with"
grep("^IFI", ann$symbol, 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 rowSums(): total count of each gene over all samples
g <- c("IFI6", "IFI27", "IFI44L", "ISG15")
g %in% ann$symbol## [1] TRUE TRUE TRUE TRUE2.5 Assignment 13 — Annotate genes
Hint: start from the code of Exercise 2 and change the names.
- How many genes are of the type
"ncRNA"(non-coding RNA), and what proportion of all genes is that? - Find all gene names starting with
"S100"(a family of neutrophil proteins). - Which of these neutrophil genes has the highest total count across all samples: S100A8, S100A9, S100A12, MMP8? Is that what you expect in blood?
3 Describe the cohort: understanding the sample metadata
Before any analysis, get to know your patients. This is where you catch mislabelled samples, repeated patients, and groups that differ in age or sex. Finding these now is much cheaper than after the analysis.
3.1 Cleaning the patient information
# Look at the raw values of each variable before changing anything.
table(meta$cohort, useNA = "ifany")##
## Bacterial CoV other COVID-19 healthy Influenza
## 24 61 77 19 17table(meta$gender, useNA = "ifany")##
## Female Male <NA>
## 84 93 21table(meta$hospitalized, useNA = "ifany")##
## No Yes <NA>
## 65 12 121head(sort(unique(meta$age))) # note the text value ">89"## [1] ">89" "0" "14" "15" "16" "18"# 1. Delete cohort=NA and group names contain spaces and dashes ("COVID_19"). Replace them with "_".
unique(meta$cohort)## [1] "Bacterial" "CoV other" "COVID-19" "healthy" "Influenza"meta$cohort <- gsub("[^A-Za-z0-9]+", "_", meta$cohort)
meta <- meta[meta$cohort %in% c("healthy", "Bacterial", "CoV_other", "COVID_19", "Influenza"), ]
table(meta$cohort)##
## Bacterial CoV_other COVID_19 healthy Influenza
## 24 61 77 19 17# 2. Age is stored as text because of ">89" (very old patients are not given an
# exact age, to protect their privacy). Replace it with 90 and make it a number.
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# 3. Delete gender=NA
meta <- meta[meta$gender %in% c("Female", "Male"), ]
table(meta$gender)##
## Female Male
## 84 93# 4. Make grouping variables "factors" (R's word for a category).
# We choose the order of the groups ourselves: healthy first.
meta$cohort <- factor(meta$cohort,
levels = c("healthy", "Bacterial", "CoV_other", "COVID_19", "Influenza"))
meta$gender <- factor(meta$gender)
meta$race <- factor(meta$race)
str(meta)## 'data.frame': 177 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 "healthy","Bacterial",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ 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 ...counts <- counts[, rownames(meta)] # remove those samples from counts
all(colnames(counts) == rownames(meta)) # check again — always## [1] TRUE3.2 Describing the cohort
Which plot and which test? It depends on the type of variable.
| Variable | Example | Plot | Test (2 groups) |
|---|---|---|---|
| Category | sex, hospitalised | bar chart | chi-square test, or Fisher’s test for small groups |
| Number | age | box plot with the points | Wilcoxon test (or t-test if the data look bell-shaped) |
3.2.1 Categories: bar chart (and a pie chart)
ggplot(meta, aes(x = cohort, fill = gender)) +
geom_bar(position = position_dodge(preserve = "single"), na.rm = TRUE) + # bars side by side
geom_text(stat = "count", aes(label = after_stat(count)), # write n on top
position = position_dodge(width = 0.9), vjust = -0.5, size = 3.2) +
scale_fill_brewer(palette = "Set2") + # the colours
labs(x = NULL, y = "number of samples", title = "Samples per cohort and sex")
# The same information as a pie chart.
tab <- as.data.frame(table(cohort = meta$cohort)) # counts per cohort
tab$prop <- tab$Freq / sum(tab$Freq) # as a proportion
ggplot(tab, aes(x = "", y = prop, fill = cohort)) +
geom_col(width = 1, colour = "white") +
coord_polar(theta = "y") + # this bends the bar into a circle
geom_text(aes(label = scales::percent(prop, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3.5) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Cohort composition") +
theme_void()
Note. People compare angles badly, so a bar chart is almost always clearer than a pie chart.
3.2.2 Numbers: 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) + # the box
geom_jitter(width = 0.15, size = 1.4, alpha = 0.6) + # one dot per patient
scale_fill_brewer(palette = "Set2") +
labs(x = NULL, y = "age (years)", title = "Age by cohort") +
theme(legend.position = "none")
Always show the points. A box made of 4 patients looks the same as a box made of 400. The dots show how many patients there really are.
3.2.3 Are COVID_19 patients and healthy people comparable?
In Lab 7 we compare COVID_19 with healthy. If the two groups differ in age or sex, some gene differences may be caused by age or sex rather than by the virus.
# Keep only the two groups we will compare in Lab 7.
sub <- meta[meta$cohort %in% c("COVID_19", "healthy"), ]
sub$cohort <- droplevels(sub$cohort)
table(sub$cohort)##
## healthy COVID_19
## 16 45# Sex (category) vs group -> chi-square, or Fisher when groups are small
tab2 <- table(sub$gender, sub$cohort)
tab2##
## healthy COVID_19
## Female 6 22
## Male 10 23chisq.test(tab2)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tab2
## X-squared = 0.24319, df = 1, p-value = 0.6219fisher.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.1596541 2.3137555
## sample estimates:
## odds ratio
## 0.6320615# Age (number) vs group -> Wilcoxon test
tapply(sub$age, sub$cohort, median) # median age in each group## healthy COVID_19
## 18 33wilcox.test(age ~ cohort, data = sub)##
## Wilcoxon rank sum test with continuity correction
##
## data: age by cohort
## W = 0.5, p-value = 3.471e-09
## alternative hypothesis: true location shift is not equal to 0# A small "Table 1", like in a clinical paper.
Table <- meta %>%
group_by(cohort) %>%
summarise(
n = n(),
median_age = median(age, na.rm = TRUE),
pct_female = round(100 * mean(gender == "Female"), 1)
)
TableDiscuss. Healthy volunteers are often younger than hospital patients. If age differs between the groups, age and disease are mixed up (confounded). Only a table like this tells you so.
3.3 Exercise 3
3.3.1 Question
- How many samples are left in each cohort, and how many samples did the sex check remove?
- Draw the age box plot again, but with the colour palette
"Pastel1". - Compare the age of Influenza patients with healthy controls. Which test did you use, and what is the p-value?
3.3.2 Solution
# 1
table(meta$cohort)##
## healthy Bacterial CoV_other COVID_19 Influenza
## 16 24 35 45 17length(exclude_ids)## [1] 18# 2 Only the palette name changed.
ggplot(meta, aes(x = cohort, y = age, fill = cohort)) +
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 = "Pastel1") +
labs(x = NULL, y = "age (years)", title = "Age by cohort") +
theme(legend.position = "none")
# 3 Age is a number and there are two groups -> Wilcoxon test.
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 03.4 Assignment 14 — Describe the cohort
Hint: start from the code of Exercise 3 and change one small thing.
- How many female and male samples
are left? (Hint:
table(meta$gender)) - Draw the bar chart “Samples per cohort and sex” again, but with the
colour palette
"Dark2". - Compare the age of Bacterial patients with
healthy controls with the right test and report the
p-value. In one sentence: are COVID_19 patients and healthy controls
comparable in age and sex (use the results of the
statschunk)?
4 Extra — Checking the recorded sex against the RNA (optional)
The RNA itself tells us the sex of each sample:
- XIST is switched on in cells with two X chromosomes (female);
- RPS4Y1, DDX3Y, UTY, KDM5D sit on the Y chromosome (male only).
If the label in the table disagrees with the RNA, the sample may have been swapped or mislabelled.
# A quick normalisation so samples with more reads do not look "higher".
# (log2 CPM is explained properly in Lab 7.)
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
# One row per sample: XIST level, average Y-gene level, and the recorded sex.
sex_df <- data.frame(
sample = colnames(logcpm_quick),
XIST = logcpm_quick["XIST", ],
y_score = colMeans(logcpm_quick[y_genes, ]),
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 sex seen in the RNA")
Discuss. You should see two clear clouds: females (high XIST, no Y) and males (Y genes, low XIST). A dot of the “wrong” colour in a cloud is a suspicious sample.
# A simple rule: high Y-gene expression = male. The cut-off is halfway
# between the lowest and highest Y score.
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# Samples where the recorded sex disagrees with the RNA
discordant <- sex_df$sample[!is.na(sex_df$label) &
as.character(sex_df$label) != sex_df$inferred]
# Samples with no recorded sex
label_missing <- sex_df$sample[is.na(sex_df$label)]
exclude_ids <- c(discordant, label_missing)
exclude_ids## [1] "S434766" "S434777" "SDU18.02S0011675" "SDU09.02S0000101" "S434393"
## [6] "SDU14.02S0000003" "S434399" "SDU14.02S0000001" "S434486" "S434570"
## [11] "S434532" "S434596" "SDU14.02S0000007" "S434537" "S434482"
## [16] "S434517" "S434518" "SDU14.02S0000005"# Remove them from meta AND counts, and check the order again.
meta_clean <- meta[!rownames(meta) %in% exclude_ids, ]
meta_clean$gender <- droplevels(meta_clean$gender) # forget empty categories
counts_clean <- counts[, rownames(meta_clean)]
all(colnames(counts_clean) == rownames(meta_clean))## [1] TRUE5 Save the objects for Lab 7
Lab 7 starts from exactly these three objects.
saveRDS(counts, file.path(data_dir, "gse161731_counts.rds")) # raw counts, gene symbols
saveRDS(meta, file.path(data_dir, "gse161731_meta.rds")) # clean patient table
saveRDS(ann, file.path(data_dir, "gse161731_gene_annotation.rds")) # symbol / ensembl / biotype
sessionInfo() # the software versions: 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] dplyr_1.2.1 ggplot2_4.0.3 edgeR_4.6.3 limma_3.64.3
## [5] org.Hs.eg.db_3.21.0 AnnotationDbi_1.72.0 IRanges_2.44.0 S4Vectors_0.48.1
## [9] GEOquery_2.76.0 Biobase_2.70.0 BiocGenerics_0.56.0 generics_0.1.4
##
## loaded via a namespace (and not attached):
## [1] KEGGREST_1.50.0 SummarizedExperiment_1.40.0 gtable_0.3.6
## [4] xfun_0.60 bslib_0.12.0 lattice_0.23-1
## [7] tzdb_0.5.0 vctrs_0.7.3 tools_4.5.0
## [10] curl_8.0.0 tibble_3.3.1 RSQLite_3.53.3
## [13] blob_1.3.0 pkgconfig_2.0.3 Matrix_1.7-6
## [16] data.table_1.18.6.1 RColorBrewer_1.1-3 rentrez_1.2.4
## [19] S7_0.2.2 lifecycle_1.0.5 GenomeInfoDbData_1.2.14
## [22] farver_2.1.2 compiler_4.5.0 Biostrings_2.78.0
## [25] statmod_1.5.2 codetools_0.2-20 Seqinfo_1.0.0
## [28] GenomeInfoDb_1.44.3 htmltools_0.5.9 sass_0.4.10
## [31] yaml_2.3.12 crayon_1.5.3 pillar_1.11.1
## [34] jquerylib_0.1.4 tidyr_1.3.2 DelayedArray_0.36.1
## [37] cachem_1.1.0 abind_1.4-8 locfit_1.5-9.12
## [40] tidyselect_1.2.1 digest_0.6.39 purrr_1.2.2
## [43] labeling_0.4.3 fastmap_1.2.0 grid_4.5.0
## [46] cli_3.6.6 SparseArray_1.10.10 magrittr_2.0.5
## [49] S4Arrays_1.10.1 XML_3.99-0.24 withr_3.0.3
## [52] readr_2.2.0 scales_1.4.0 UCSC.utils_1.4.0
## [55] bit64_4.8.6 rmarkdown_2.32 XVector_0.50.0
## [58] httr_1.4.9 matrixStats_1.5.0 bit_4.6.0
## [61] otel_0.2.0 png_0.1-9 hms_1.1.4
## [64] memoise_2.0.1 evaluate_1.0.5 knitr_1.51
## [67] GenomicRanges_1.62.1 rlang_1.3.0 glue_1.8.1
## [70] DBI_1.3.0 xml2_1.6.0 rstudioapi_0.19.0
## [73] jsonlite_2.0.0 R6_2.6.1 MatrixGenerics_1.22.0Take-home messages
- For differential expression, download raw counts.
- Check that counts and metadata are in the same order — and check again after every time you remove samples.
- Analyse with stable Ensembl IDs, show results with gene symbols, and say how you handled missing and repeated symbols.
- Never trust a label blindly: the RNA can confirm (or contradict) the recorded sex.
- Describe your groups before comparing them; differences in age or sex can pretend to be disease effects.
Further reading
- McClain MT et al. (2021) Dysregulated transcriptional responses to SARS-CoV-2 in the periphery. Nat Commun 12:1079. (The GSE161731 paper.)
- Conesa A et al. (2016) A survey of best practices for RNA-seq data analysis. Genome Biology 17:13.