SummarizedExperiment objectWe’ve finally reached the point in which we can load data into
SummarizedExperiment objects.
We are going to simplify the analysis by only using the 22 cerebellum samples. Expression quantification was done using STAR to align to the mouse genome and then counting reads that map to genes. In addition to the counts per gene per sample, we also need information on which sample belongs to which Sex/Time point/Replicate. And for the genes, it is helpful to have extra information called annotation.
Let’s read in the data files that we downloaded in the last episode and start to explore them:
One of the things Bioconductor requires: our counts matrix needs to
have row names. The tidyverse does not like row names, so
we’ll use read.csv() (instead of
readr::read_csv()).
Because we need the rownames, we need to use a different process for
loading in our data from CSV format. We’ll use read.csv()
with the row.names argument to specify which columns are
the rownames.
counts <- read.csv(system.file("extdata/GSE96870_counts_cerebellum.csv", package="bioc2026Intro"),
row.names = 1)
dim(counts)## [1] 41786 22
head(counts)## GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
## Xkr4 1891 2410 2159 1980 1977 1945
## LOC105243853 0 0 1 4 0 0
## LOC105242387 204 121 110 120 172 173
## LOC105242467 12 5 5 5 2 6
## Rp1 2 2 0 3 2 1
## Sox17 251 239 218 220 261 232
## GSM2545342 GSM2545343 GSM2545344 GSM2545345 GSM2545346 GSM2545347
## Xkr4 1757 2235 1779 1528 1644 1585
## LOC105243853 1 3 3 0 1 3
## LOC105242387 177 130 131 160 180 176
## LOC105242467 3 2 2 2 1 2
## Rp1 3 1 1 2 2 2
## Sox17 179 296 233 271 205 230
## GSM2545348 GSM2545349 GSM2545350 GSM2545351 GSM2545352 GSM2545353
## Xkr4 2275 1881 2584 1837 1890 1910
## LOC105243853 1 0 0 1 1 0
## LOC105242387 161 154 124 221 272 214
## LOC105242467 2 4 7 1 3 1
## Rp1 3 6 5 3 5 1
## Sox17 302 286 325 201 267 322
## GSM2545354 GSM2545362 GSM2545363 GSM2545380
## Xkr4 1771 2315 1645 1723
## LOC105243853 0 1 0 1
## LOC105242387 124 189 223 251
## LOC105242467 4 2 1 4
## Rp1 3 3 1 0
## Sox17 273 197 310 246We could also do the following:
counts2 <- read.csv(system.file("extdata/GSE96870_counts_cerebellum.csv", package="bioc2026Intro"))
rownames(counts2) <- counts2[,1]
counts2 <- counts2[,-1]
head(counts2)## GSM2545336 GSM2545337 GSM2545338 GSM2545339 GSM2545340 GSM2545341
## Xkr4 1891 2410 2159 1980 1977 1945
## LOC105243853 0 0 1 4 0 0
## LOC105242387 204 121 110 120 172 173
## LOC105242467 12 5 5 5 2 6
## Rp1 2 2 0 3 2 1
## Sox17 251 239 218 220 261 232
## GSM2545342 GSM2545343 GSM2545344 GSM2545345 GSM2545346 GSM2545347
## Xkr4 1757 2235 1779 1528 1644 1585
## LOC105243853 1 3 3 0 1 3
## LOC105242387 177 130 131 160 180 176
## LOC105242467 3 2 2 2 1 2
## Rp1 3 1 1 2 2 2
## Sox17 179 296 233 271 205 230
## GSM2545348 GSM2545349 GSM2545350 GSM2545351 GSM2545352 GSM2545353
## Xkr4 2275 1881 2584 1837 1890 1910
## LOC105243853 1 0 0 1 1 0
## LOC105242387 161 154 124 221 272 214
## LOC105242467 2 4 7 1 3 1
## Rp1 3 6 5 3 5 1
## Sox17 302 286 325 201 267 322
## GSM2545354 GSM2545362 GSM2545363 GSM2545380
## Xkr4 1771 2315 1645 1723
## LOC105243853 0 1 0 1
## LOC105242387 124 189 223 251
## LOC105242467 4 2 1 4
## Rp1 3 3 1 0
## Sox17 273 197 310 246Genes are in rows and samples are in columns, so we have counts for 41,786 genes and 22 samples.
Next read in the sample annotations. Because samples are in columns
in the count matrix, we will name the object coldata:
coldata <- read.csv(system.file("extdata/GSE96870_coldata_cerebellum.csv", package="bioc2026Intro"),
row.names = 1)
dim(coldata)## [1] 22 10
head(coldata)## title geo_accession organism age sex
## GSM2545336 CNS_RNA-seq_10C GSM2545336 Mus musculus 8 weeks Female
## GSM2545337 CNS_RNA-seq_11C GSM2545337 Mus musculus 8 weeks Female
## GSM2545338 CNS_RNA-seq_12C GSM2545338 Mus musculus 8 weeks Female
## GSM2545339 CNS_RNA-seq_13C GSM2545339 Mus musculus 8 weeks Female
## GSM2545340 CNS_RNA-seq_14C GSM2545340 Mus musculus 8 weeks Male
## GSM2545341 CNS_RNA-seq_17C GSM2545341 Mus musculus 8 weeks Male
## infection strain time tissue mouse
## GSM2545336 InfluenzaA C57BL/6 Day8 Cerebellum 14
## GSM2545337 NonInfected C57BL/6 Day0 Cerebellum 9
## GSM2545338 NonInfected C57BL/6 Day0 Cerebellum 10
## GSM2545339 InfluenzaA C57BL/6 Day4 Cerebellum 15
## GSM2545340 InfluenzaA C57BL/6 Day4 Cerebellum 18
## GSM2545341 InfluenzaA C57BL/6 Day8 Cerebellum 6Now samples are in rows with the GEO sample IDs as the rownames, and
we have 10 columns of information. The columns that are the most useful
for this workshop are geo_accession (GEO sample IDs again),
sex and time.
The counts only have gene symbols, which while short and somewhat recognizable to the human brain, are not always good absolute identifiers for exactly what gene was measured. For this we need additional gene annotations that were provided by the authors.
The count and coldata files were in comma
separated value (.csv) format, but we cannot use that for our gene
annotation file because the descriptions can contain commas that would
prevent a .csv file from being read in correctly.
Instead the gene annotation file is in tab separated value (.tsv)
format. Likewise, the descriptions can contain the single quote
' (e.g., 5’), which by default R assumes indicates a
character entry. So we have to use a more generic function
read.delim() with extra arguments to specify that we have
tab-separated data (sep = "\t") with no quotes used
(quote = "").
We also put in other arguments to specify that the first row contains
our column names (header = TRUE), the gene symbols that
should be our row.names are in the 5th column
(row.names = 5), and that NCBI’s species-specific gene ID
(i.e., ENTREZID) should be read in as character data even though they
look like numbers (colClasses argument). You can look up
this details on available arguments by simply entering the function name
starting with question mark. (e.g., ?read.delim)
rowranges <- read.delim(system.file("extdata/GSE96870_rowranges.tsv", package="bioc2026Intro"),
sep = "\t",
colClasses = c(ENTREZID = "character"),
header = TRUE,
quote = "",
row.names = 5)
dim(rowranges)## [1] 41786 7Let’s look at the first few rows of rowranges
head(rowranges)## seqnames start end strand ENTREZID
## Xkr4 1 3670552 3671742 - 497097
## LOC105243853 1 3357323 3366505 + 105243853
## LOC105242387 1 3658847 3670456 - 105242387
## LOC105242467 1 4233436 4233728 - 105242467
## Rp1 1 4409170 4409241 - 19888
## Sox17 1 4496291 4497354 - 20671
## product
## Xkr4 X Kell blood group precursor related family member 4, transcript variant X1
## LOC105243853 uncharacterized LOC105243853, transcript variant X2
## LOC105242387 uncharacterized LOC105242387
## LOC105242467 lipoxygenase homology domain-containing protein 1, transcript variant X2
## Rp1 retinitis pigmentosa 1 (human), transcript variant 2
## Sox17 SRY (sex determining region Y)-box 17, transcript variant 5
## gbkey
## Xkr4 mRNA
## LOC105243853 ncRNA
## LOC105242387 ncRNA
## LOC105242467 mRNA
## Rp1 mRNA
## Sox17 mRNAFor each of the 41,786 genes, we have the seqnames
(e.g., chromosome number), start and end
positions, strand, ENTREZID, gene product
description (product) and the feature type
(gbkey). These gene-level metadata are useful for the
downstream analysis. For example, from the gbkey column, we
can check what types of genes and how many of them are in our
dataset:
table(rowranges$gbkey)##
## C_region D_segment exon J_segment misc_RNA
## 20 23 4008 94 1988
## mRNA ncRNA precursor_RNA rRNA tRNA
## 21198 12285 1187 35 413
## V_segment
## 535We will create a SummarizedExperiment from these
objects:
count object will be saved in assays
slotcoldata object with sample information will be
stored in colData slot (sample
metadata)rowranges object describing the genes will be
stored in rowRanges slot (features
metadata)Before we put them together, you ABSOLUTELY MUST MAKE SURE THE
SAMPLES AND GENES ARE IN THE SAME ORDER! Even though we saw that
count and coldata had the same number of
samples and count and rowranges had the same
number of genes, we never explicitly checked to see if they were in the
same order. One quick way to check:
## [1] TRUE## [1] TRUE
# If the first is not TRUE, you can match up the samples/columns in
# counts with the samples/rows in coldata like this (which is fine
# to run even if the first was TRUE):
tempindex <- match(colnames(counts), rownames(coldata))
coldata <- coldata[tempindex, ]
# Check again:
all.equal(colnames(counts), rownames(coldata)) ## [1] TRUEOnce we have verified that samples and genes are in the same order,
we can then create our SummarizedExperiment object. Then we
can also create our DESeqDataSet object as well:
# One final check:
stopifnot(rownames(rowranges) == rownames(counts), # features
rownames(coldata) == colnames(counts)) # samples
se <- SummarizedExperiment(
assays = list(counts = as.matrix(counts)),
rowRanges = as(rowranges, "GRanges"),
colData = coldata
)
se## class: RangedSummarizedExperiment
## dim: 41786 22
## metadata(0):
## assays(1): counts
## rownames(41786): Xkr4 LOC105243853 ... TrnT TrnP
## rowData names(3): ENTREZID product gbkey
## colnames(22): GSM2545336 GSM2545337 ... GSM2545363 GSM2545380
## colData names(10): title geo_accession ... tissue mouseNow we can move forward as usual.