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 ENSG00000165949 with 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 document

1 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.

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

There are many public datasets available.

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).

We are going to use Gene Expression Omnibus today. 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/RNASeq_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")
counts_file;key_file

## [1] "data/GSE161731/GSE161731_counts.csv.gz"

## [1] "data/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) becomes the row names
# check.names = FALSE -> keep the column names exactly as they are => 
# Why? In R, when check.names = TRUE (the default in functions such as read.csv()), column names that begin with numbers may be automatically modified to make them syntactically valid, for example, 123 may become X123.
counts_raw <- read.csv(counts_file, row.names = 1, check.names = FALSE)
key_raw    <- read.csv(key_file,    row.names = 1, check.names = FALSE)

# If you encounter internet connectivity or download issues, please use `read.csv()` to load the data you downloaded from Studium.
counts_raw <- read.csv("data/gse161731_counts_raw.csv", row.names = 1, check.names = FALSE)
key_raw    <- read.csv("data/gse161731_key_raw.csv",    row.names = 1, check.names = FALSE)
dim(counts_raw)   # rows = genes, columns = samples

## [1] 60675   201

dim(key_raw)      # rows = samples, columns = patient information

## [1] 198   8

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

head(key_raw)          # first rows of the patient table

Three 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))

##           SS94189 SSDU09.03S0000604 SSDU09.03S0000611          SS105920 SSDU09.03S0000774 
##          52431553          50351794          52172232          53276313          43639732 
## SSDU09.03S0000775 
##          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 ".".
# We added the prefix `"S"` to sample and row identifiers to prevent R from automatically modifying names that start with numbers, ensuring consistent identifiers throughout the analysis.
head(colnames(counts_raw))

## [1] "SS94189"           "SSDU09.03S0000604" "SSDU09.03S0000611" "SS105920"         
## [5] "SSDU09.03S0000774" "SSDU09.03S0000775"

head(rownames(key_raw))

## [1] "SS94189"           "SSDU09.03S0000604" "SSDU09.03S0000611" "SS105920"         
## [5] "SSDU09.03S0000774" "SSDU09.03S0000775"

# paste0() combines strings or character values together without adding any separator.
# paste0("S","a")
colnames(counts_raw) <- paste0("S", colnames(counts_raw))
rownames(key_raw)    <- paste0("S", rownames(key_raw))
# gsub() searches for all occurrences of a specified pattern in a string and replaces them with a given value.
# gsub("S","s", "Star")
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, ]
# order(): arranges the data in ascending order based on the specified variable.
meta   <- meta[order(meta$cohort), ] 
counts <- counts_raw[, rownames(meta)]

# checks whether the column names of counts exactly match the row names of meta in the same order, returning TRUE if all match and FALSE otherwise.
all(colnames(counts) == rownames(meta))

## [1] TRUE

# The most important check of the day. Why is this important? It ensures that each sample's gene expression counts correspond to the correct metadata, preventing samples from being incorrectly assigned to cohorts or other experimental groups.

str(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        17

1.5 Exercise 1

1.5.1 Question

How many genes and how many samples are in counts?
Which sample has the largest library size?
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
#colSums() calculates the sum of each column in a matrix or data frame
lib <- colSums(counts)              # total reads per sample
names(which.max(lib))               # which.max() finds the position; names() gives the sample ID/name

## [1] "SSS434325"

# 3  setdiff(a, b): things in a that are not in b
setdiff(colnames(counts_raw), rownames(key_raw))

## [1] "SSSDU09.02S0000150.batch2" "SSSDU09.02S0000154.batch2" "SSSDU09.02S0000158.batch2"

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

## [1] 3

1.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: 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) lists all the types of annotation information available in the org.Hs.eg.db database, such as gene symbols, Entrez Gene IDs, Ensembl IDs, and gene names.
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"

# Retrieves the gene name for each Ensembl gene ID
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
)

# Retrieves the gene type associated with each Ensembl gene ID, such as protein-coding or non-coding RNA
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] 35529

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

## [1] 25146

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

## [1] 82

Notes 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] 35447

counts       <- counts[keep_row, ]
gene_ids     <- gene_ids[keep_row]
gene_biotype <- gene_biotype[keep_row]
rownames(counts) <- gene_symbol[keep_row]   # from now on, rows are 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))
dim(ann)

## [1] 35447     3

head(ann)

# If you find `One row per gene symbol` diffcult to understand, please use `read.csv()` to load the annotation file you downloaded from Studium.
ann <- read.csv("data/gse161731_gene_annotation.csv")
# Which kinds of genes are in the table?
table(ann$biotype)

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

sort(table(ann$biotype), decreasing = TRUE)

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

sort(table(ann$biotype), decreasing = F)

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

# Are the genes we care about present?
c("IFI27", "ISG15", "S100A8", "XIST", "RPS4Y1") %in% ann$symbol

## [1] TRUE TRUE TRUE TRUE TRUE

Biology 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?
Find all gene names starting with "IFI" (interferon-induced genes).
Are the genes we care about present: IFI6, IFI27, IFI44L, ISG15?

2.4.2 Solution

# 1  == compares every value; sum() counts the TRUEs;
sum(ann$biotype == "protein-coding", na.rm = TRUE)

## [1] 19208

# 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  
g <- c("IFI6", "IFI27", "IFI44L", "ISG15")
g %in% ann$symbol # %in% checks whether each element in one vector is present in another vector, returning TRUE for matches and FALSE otherwise.

## [1] TRUE TRUE TRUE TRUE

2.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)?
Find all gene names starting with "S100" (a family of neutrophil proteins).
Are the genes we care about present: S100A8, S100A9, S100A12, MMP8? ___

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.
# useNA = "ifany": displays the number of missing (NA) values if they exist
table(meta$cohort,useNA = "ifany")

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

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

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

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

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

# sort() arranges the elements of a vector in ascending order by default.
# unique() returns the distinct elements of a vector, removing duplicates.
head(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"

# [^A-Za-z0-9]+ is a regular expression (regex) that matches one or more consecutive characters that are not letters or numbers.
# gsub("[^A-Za-z0-9]+", "_", "ENSG-00001.2 (gene)")
meta$cohort <- gsub("[^A-Za-z0-9]+", "_", meta$cohort)
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 <- gsub(">89", "90", meta$age)
meta$age <- as.numeric(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).
# factor() converts a variable into a categorical variable and optionally sets the category order.
# 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] TRUE

3.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) +   
  # position_dodge(preserve = "single"): places bars side by side and preserves the width of each individual bar.
  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"), ]
table(sub$cohort)

## 
##   healthy Bacterial CoV_other  COVID_19 Influenza 
##        19         0         0        77         0

sub$cohort <- droplevels(sub$cohort) # removes unused factor levels after filtering data.
table(sub$cohort)

## 
##  healthy COVID_19 
##       19       77

# Age (number) vs group -> Wilcoxon test
tapply(sub$age, sub$cohort, median)   # median age in each group

##  healthy COVID_19 
##       18       33

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

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

# Sex (category) vs group -> Fisher when groups are small (Chi-squared test; suitable for sufficiently large expected counts)
tab2 <- table(sub$gender, sub$cohort)
tab2

##         
##          healthy COVID_19
##   Female       8       32
##   Male        11       45

fisher.test(tab2) # Fisher’s exact test; suitable for small expected counts.

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

Discuss. 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?
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  how many samples are left in each cohort
table(meta$cohort)

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

# 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 = 180, p-value = 0.541
## alternative hypothesis: true location shift is not equal to 0

3.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$))
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 stats chunk)?

4 Save the objects for Lab 7

Lab 7 starts from exactly these three objects.

write.csv(counts, file.path(data_dir, "gse161731_counts.csv"))           # raw counts, gene symbols
write.csv(meta,   file.path(data_dir, "gse161731_meta.csv"))             # clean patient table
write.csv(ann,    file.path(data_dir, "gse161731_gene_annotation.csv"))  # symbol / ensembl / biotype

5 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     80    4
##   Male        1   92

# 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] "SSS434766"          "SSS434777"          "SSSDU18.02S0011674" "SSSDU18.02S0011675"
## [5] "SSSDU09.02S0000101"

# 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] TRUE

# 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)
  )
Table

Take-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.

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