The bulk RNA-Seq assay measuring immune blood cells is a component of the CMI-Flu prediction challenge. An overview of all the data, including links to more detailed descriptions of the other data, may be found at CMI-x .

Experimental protocol

RNA was extracted from PBMCs using the miRNeasy Mini Kit (Qiagen). Libraries were prepped using Watchmaker mRNA Library Prep Kit (Watchmaker Genomics). Libraries were sequenced on a NovaSeq6000 (Illumina) system, targeting 20 million reads per sample. Processing was performed using the default nf-core/rna-seq pipeline , aligning to the ensembl reference genome v113. We provide a mapping file of all genes (as ensembl IDs and gene symbols) in the “reference files” folder.

Data Standardization

RNA-Seq data generated for CMI-Flu was combined with other publicly-available data.

  • Where possible, public data was rerun through the nf-core/rna-seq pipeline .
  • Data is provided (where possible) as raw counts (integers) and transcripts per million (TPM) which normalizes each gene’s read count as a fraction of transcripts in the sample, then scaled to 1e06.

The Data

RNA-Seq data from each study has its own “rnaseq” table. This allows contestants to download and explore one table at a time before deciding which data they wish to incorporate into their models. While these data are not directly linked to a specific prediction task, you can leverage them to build more robust models and improve predictions for other data types within the 2025 LJI prediction challenge cohort. There are 5 publicly curated datasets with pre- and post-vaccination values, and one in-house generated set (2025LJI prediction/challenge cohort) with pre-vaccination data only. Note that in 2025LJI, 3 participants are missing Day-14 measures, and 2 donors are missing Day0 measures.

Data tables

list.files('../datasets/260512/train/', pattern='rnaseq')
## [1] "2025LJI_rnaseq-BATCH-CORRECTED.tsv"  "2025LJI_rnaseq.tsv"                 
## [3] "publicData_rnaseq_2019_UGA.tsv"      "publicData_rnaseq_2020_UGA.tsv"     
## [5] "publicData_rnaseq_2024_UGA.tsv"      "publicData_rnaseq_SDY224_Bcells.tsv"
## [7] "publicData_rnaseq_SDY2867.tsv"       "publicData_rnaseq_SDY2941.tsv"
RNA-Seq Assay
Data dictionary
Column Description
participant_id Donor ID, one for each subject in each study, in the format studyID.subjectID (links to participants.tsv)
timepoint Days relative to influenza vaccination.
ensembl_gene_id Ensembl gene identifier (ENSG...), unversioned.
counts raw counts (integers).
tpm transcripts per million (library size normalization).
material Sample type: PBMCs, B cells, or other.
study_accession Study identifier. Links to the investigations.tsv file.
subject subject ID.

Data exploration

Data preview:

rnaseq = read_tsv('../datasets/260512/train/2025LJI_rnaseq.tsv')
## Rows: 5919900 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): participant_id, ensembl_gene_id, material, comments, study_accessio...
## dbl (3): timepoint, counts, tpm
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(rnaseq)
## # A tibble: 6 × 9
##   participant_id  timepoint ensembl_gene_id counts    tpm material comments     
##   <chr>               <dbl> <chr>            <dbl>  <dbl> <chr>    <chr>        
## 1 2025LJI.SUB1829       -14 ENSG00000000003      4 0.0932 PBMCs    Batch : Plat…
## 2 2025LJI.SUB1829         0 ENSG00000000003     15 0.699  PBMCs    Batch : Plat…
## 3 2025LJI.SUB2491       -14 ENSG00000000003     14 0.694  PBMCs    Batch : Plat…
## 4 2025LJI.SUB2491         0 ENSG00000000003     14 0.211  PBMCs    Batch : Plat…
## 5 2025LJI.SUB3691       -14 ENSG00000000003      9 0.198  PBMCs    Batch : Plat…
## 6 2025LJI.SUB3691         0 ENSG00000000003     13 0.103  PBMCs    Batch : Plat…
## # ℹ 2 more variables: study_accession <chr>, subject <chr>

Batch correction

In the 2025LJI Prediction/challenge cohort data generated in house, we detected batch effects across our different plates when examining both pre- and post-vaccination data. To account for this, we applied limma’s removeBatchEffect() function across paired samples in the pre-vaccination data (D-14 and D0), and then subtracted the calculated correction to the remaining data (unpaired samples). This batch-corrected data is supplied in 2025LJI_rnaseq_BATCH-CORRECTED.tsv. The code is below.

library(tidyverse)
library(DESeq2)
library(limma)


rnaseq = read_tsv('../datasets/260512/train/2025LJI_rnaseq.tsv')

# extract counts into wide format
counts = rnaseq %>%
  mutate(sample = paste(participant_id, timepoint, sep = '_')) %>%
  select(ensembl_gene_id, sample, counts) %>%
  pivot_wider(names_from = sample, values_from = counts)
dim(counts)
## [1] 78932    76
counts[1:5, 1:5]
## # A tibble: 5 × 5
##   ensembl_gene_id `2025LJI.SUB1829_-14` `2025LJI.SUB1829_0`
##   <chr>                           <dbl>               <dbl>
## 1 ENSG00000000003                     4                  15
## 2 ENSG00000000005                     0                   0
## 3 ENSG00000000419                    51                 296
## 4 ENSG00000000457                    75                 261
## 5 ENSG00000000460                    26                  81
## # ℹ 2 more variables: `2025LJI.SUB2491_-14` <dbl>, `2025LJI.SUB2491_0` <dbl>
# create metadata table
meta = rnaseq %>%
  mutate(sample = paste(participant_id, timepoint, sep = '_')) %>%
  distinct(sample, subject, timepoint, comments) %>%
  mutate(plate = if_else(timepoint %in% c(0, 1), 'plate2', 'plate1'))

## to start, only select subjects with both timepoints (paired data)
meta = meta %>%
  group_by(subject) %>%
  filter(n() == 2) %>%
  ungroup()
dim(meta)
## [1] 70  5
head(meta)
## # A tibble: 6 × 5
##   sample              subject timepoint comments        plate 
##   <chr>               <chr>       <dbl> <chr>           <chr> 
## 1 2025LJI.SUB1829_-14 SUB1829       -14 Batch : Plate 1 plate1
## 2 2025LJI.SUB1829_0   SUB1829         0 Batch : Plate 2 plate2
## 3 2025LJI.SUB2491_-14 SUB2491       -14 Batch : Plate 1 plate1
## 4 2025LJI.SUB2491_0   SUB2491         0 Batch : Plate 2 plate2
## 5 2025LJI.SUB3691_-14 SUB3691       -14 Batch : Plate 1 plate1
## 6 2025LJI.SUB3691_0   SUB3691         0 Batch : Plate 2 plate2
counts = counts %>%
  select(ensembl_gene_id, all_of(meta$sample)) %>%
    column_to_rownames('ensembl_gene_id')

## variance-stabilizing transform
dds_vst <- DESeqDataSetFromMatrix(
    countData = counts,
    colData   = meta,
    design    = ~ 1
)
dds_vst     <- estimateSizeFactors(dds_vst)
vst_mat     <- log2(counts(dds_vst, normalized=TRUE) + 1)

# batch coefficients on pre-vaccination data

pre_samples <- meta$sample[
    meta$timepoint %in% c("-14","0")]
meta_pre    <- meta %>% filter(sample %in% pre_samples)
table(meta_pre$timepoint)
## 
## -14   0 
##  35  35
vst_pre     <- vst_mat[, pre_samples]

# Fit plate effect within pre-vaccination samples, blocking on donor
plate_pre        <- as.numeric(factor(meta_pre$plate)) - 1  # 0=Plate1, 1=Plate2
table(plate_pre)
## plate_pre
##  0  1 
## 35 35
design_batch     <- model.matrix(~ subject + plate_pre, data=meta_pre)
fit_batch        <- lmFit(vst_pre, design_batch)
fit_batch        <- eBayes(fit_batch)
batch_coef       <- fit_batch$coefficients[, "plate_pre"]

cat("Batch coefficients estimated for", length(batch_coef), "genes\n")
## Batch coefficients estimated for 78932 genes
cat("Batch coefficient range:", round(range(batch_coef), 3), "\n")
## Batch coefficient range: -2.767 3.889
# Subtract plate effect from every sample
plate_all     <- as.numeric(factor(meta$plate)) - 1
vst_corrected <- vst_mat - outer(batch_coef, plate_all)

# Remove zero-variance genes
gene_vars          <- apply(vst_corrected, 1, var)
vst_corrected_filt <- vst_corrected[gene_vars > 1e-10, ]
cat("Genes after zero-variance filter:", nrow(vst_corrected_filt), "\n")
## Genes after zero-variance filter: 43495
# select the top 2000 genes by variance
top_genes          <- names(sort(gene_vars[gene_vars > 1e-10], decreasing=TRUE))[1:2000]
vst_corrected_filt <- vst_corrected_filt[top_genes, ]
cat("Genes after top-variance filter:", nrow(vst_corrected_filt), "\n")
## Genes after top-variance filter: 2000
# plots

pca_corr <- prcomp(t(vst_corrected_filt), scale.=TRUE)
var_exp  <- summary(pca_corr)$importance[2,] * 100

pca_df <- as.data.frame(pca_corr$x[,1:3])
pca_df$sample <- rownames(pca_df)
pca_df <- left_join(pca_df,
                    meta %>% dplyr::select(sample, subject, timepoint, plate),
                    by="sample") %>%
    mutate(timepoint = as.factor(timepoint))

# R² of PC1 for plate — want this to be low after correction
fit_pc1 <- lm(pca_corr$x[,1] ~ meta$plate)
fit_pc2 <- lm(pca_corr$x[,2] ~ meta$plate)
cat("PC1 R² for Plate:", round(summary(fit_pc1)$r.squared, 3), "\n")
## PC1 R² for Plate: 0
cat("PC2 R² for Plate:", round(summary(fit_pc2)$r.squared, 3), "\n")
## PC2 R² for Plate: 0
timepoint_cols <- c("-14"="#756BB1","0"="#31A354")

p_time <- ggplot(pca_df, aes(PC1, PC2, color=timepoint, shape=plate)) +
    geom_point(size=2.5, alpha=0.85) +
    scale_color_manual(values=timepoint_cols) +
    labs(title="Pre-vacc limma-corrected PCA",
         x=paste0("PC1 (",round(var_exp[1],1),"%)"),
         y=paste0("PC2 (",round(var_exp[2],1),"%)")) +
    theme_bw()

p_time

# long format for competition (not run)

vst_corrected_long = vst_corrected %>%
  as.data.frame() %>%
  rownames_to_column('ensembl_gene_id') %>%
  pivot_longer(cols = -ensembl_gene_id,
               names_to = 'ID',
               values_to = 'batch_corrected_expression') %>%
  separate(ID, into = c('participant_id', 'timepoint'), sep = '_') %>%
  mutate(study_accession = '2025LJI',
         subject         = str_remove(participant_id, '2025LJI\\.'),
         comments        = 'limma batch correction',
         material        = 'PBMCs') %>%
  select(participant_id, timepoint, ensembl_gene_id,
         batch_corrected_expression, material, comments,
         study_accession, subject)
head(vst_corrected_long)

table(unique(vst_corrected_long$participant_id) %in% unique(rnaseq$participant_id))
table(unique(vst_corrected_long$ensembl_gene_id) %in% unique(rnaseq$ensembl_gene_id))