library(tidyverse)Exploratory Analysis
Exploratory analysis will cover basic visualisations and characterisations of the data before we move on to identifying differentially expressed genes.
What is exploratory analysis and why is it necessary?
Exploratory analysis is looking at the data. You are attempting to get a sense of what the data looks like, to notice outliers or issues.
If your data set was a physical problem (like a squeaky wheel or a well wrapped present) you’d physically explore it (pick it up, is it heavy, do I shake it?). Exploratory analysis is the data science equivalent.
Exploratory analysis
From here and for the rest of the workshop, we will be working within R.
In the previous sections we generated counts for 3 samples. We are going to need a few more samples for the rest of this workshop when we explore gene expression, so we are loading in a counts file with 15 samples. We are also going to load in a metadata file that contains information about each sample, which we will load into R as a dataframe called coldata.
We have 3 phenotype groups, each with 5 biological replicates:
- F = female gonad samples
- MT = mid-transitioning gonad samples
- TPM = terminal phase male gonad samples
Run the code below to load the file:
setwd("~/RNA_seq")
coldata <- read.table("DataFiles/coldata.txt", sep="\t", header=TRUE)
counts <- read.table("DataFiles/counts.txt", sep="\t", header=TRUE, row.names=1)
Want to load this dataset on your own computer?
If you are completing this workshop on your own local machine, run the code below to load the file from Github:
coldata <- read.table("https://raw.githubusercontent.com/GenomicsAotearoa/RNA-seq-data-analysis-workflow-NZ/refs/heads/main/data-files/coldata.txt", sep="\t", header=TRUE)
counts <- read.table("https://raw.githubusercontent.com/GenomicsAotearoa/RNA-seq-data-analysis-workflow-NZ/refs/heads/main/data-files/counts.txt", sep="\t", header=TRUE, row.names=1)Now let’s take a look at the counts data.
counts |> head() SI16_21G SI16_23G SI18_30G SI18_48G SI18_52G SI16_22G SI18_15G
LOC117820358 0 0 0 0 0 0 0
LOC117813751 0 0 0 0 0 5 0
ccdc51 55 62 63 46 93 57 83
tma7 1256 2246 1492 1870 1664 2506 1514
LOC117815572 2 1 2 1 1 11 1
LOC117818779 0 1 0 0 0 2 2
SI18_16G SI18_18G SI18_22G SI16_1G SI16_2G SI16_3G SI18_47G
LOC117820358 0 0 0 0 0 0 0
LOC117813751 0 0 0 723 346 482 140
ccdc51 52 61 85 241 181 175 190
tma7 1322 1053 1129 3262 3098 3050 2336
LOC117815572 3 0 4 3 3 22 30
LOC117818779 2 0 1 0 3 0 16
SI18_49G
LOC117820358 0
LOC117813751 33
ccdc51 103
tma7 2245
LOC117815572 15
LOC117818779 0ncol(counts)[1] 15Let’s look at the metadata we have for these samples:
coldata |> head() sample batchname condition histology timepoint adj.timepoint length.mm
1 SI16_1G SI16 T TPM 0 0 225
2 SI16_21G SI16 C F 50 0 172
3 SI16_22G SI16 T MT 50 20 178
4 SI16_23G SI16 C F 50 0 173
5 SI16_2G SI16 T TPM 0 0 220
6 SI16_3G SI16 T TPM 0 0 215
weight.g
1 175.7
2 72.6
3 80.0
4 68.9
5 165.8
6 130.3The first column has all our sample names, which match exactly the column names in our counts table. We then have multiple columns that contain various metadata about our samples, including whether or not samples were part of control (C) or treatment (T) groups, the batch the sample was processed in (these samples come from experiments across two years - 2016 and 2018), the histologically determined phenotype (F = female gonad, MT = mid-transitioning gonad, TPM = terminal phase male gonad), timepoint columns and the original length and weight of the fish that the gonad was retrieved from.
Remove non-expressed genes (optional)
If you are working in a complex organism, many genes will not be expressed in all tissue types. It can be worth checking to see how many genes have zero counts and opting to remove these from your data object.
How many genes have counts of zero? What’s the distribution of our counts?
colSums(counts == 0)SI16_21G SI16_23G SI18_30G SI18_48G SI18_52G SI16_22G SI18_15G SI18_16G
9539 8885 9832 9090 9442 5806 8658 9600
SI18_18G SI18_22G SI16_1G SI16_2G SI16_3G SI18_47G SI18_49G
8763 9704 5714 5625 4951 3678 4823 summary(counts) SI16_21G SI16_23G SI18_30G SI18_48G
Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.0
Median : 7.0 Median : 10.0 Median : 4.0 Median : 6.0
Mean : 316.7 Mean : 339.7 Mean : 324.6 Mean : 331.9
3rd Qu.: 182.0 3rd Qu.: 194.0 3rd Qu.: 116.0 3rd Qu.: 116.0
Max. :90878.0 Max. :128697.0 Max. :201617.0 Max. :192109.0
SI18_52G SI16_22G SI18_15G SI18_16G
Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 0.0 1st Qu.: 1.0 1st Qu.: 0.0 1st Qu.: 0.0
Median : 5.0 Median : 24.0 Median : 6.0 Median : 4.0
Mean : 359.7 Mean : 331.5 Mean : 312.4 Mean : 299.2
3rd Qu.: 133.0 3rd Qu.: 192.0 3rd Qu.: 78.0 3rd Qu.: 47.0
Max. :208754.0 Max. :123405.0 Max. :297486.0 Max. :543596.0
SI18_18G SI18_22G SI16_1G SI16_2G
Min. : 0 Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 0 1st Qu.: 0.0 1st Qu.: 1.0 1st Qu.: 1.0
Median : 7 Median : 5.0 Median : 24.0 Median : 25.0
Mean : 267 Mean : 329.7 Mean : 334.1 Mean : 294.4
3rd Qu.: 55 3rd Qu.: 78.0 3rd Qu.: 222.0 3rd Qu.: 211.0
Max. :454735 Max. :344393.0 Max. :54936.0 Max. :54766.0
SI16_3G SI18_47G SI18_49G
Min. : 0.0 Min. : 0.0 Min. : 0.0
1st Qu.: 2.0 1st Qu.: 6.0 1st Qu.: 2.0
Median : 46.0 Median : 78.0 Median : 36.0
Mean : 340.1 Mean : 305.3 Mean : 193.1
3rd Qu.: 303.0 3rd Qu.: 323.0 3rd Qu.: 180.0
Max. :49808.0 Max. :36463.0 Max. :28195.0 Zero counts sometimes need to be dealt with for differential expression (we’ll get to this in the next section, as different pipelines handle zeros differently), but for visualisation purposes we often do need to filter them out manually. It is common to remove genes with very low counts before any downstream analysis or plotting, but there are different approaches to chosing the threshold. For now, when plotting genes we need to log2 the counts, after adding 1 to all counts so that we don’t log zero (which will result in an error or infinite value).
To remove all rows with 0 counts (genes where no sample has a read mapping to it), you could run something like this:
# Demo code - don't need to run this
counts_filtered <- counts[rowSums(counts) > 0, ]Visualising our data set
In this section, focus on what we are trying to do with the code rather than the code itself. If you can understand what we are trying to achieve by the visualisations, then you find your own way to do this in R. We are not going to focus on how ggplot2 works or the nuances of transforming data specifically for ggplot2.
Count distribution
Let’s first look at the distribution of counts across our samples. Boxplots can give as a quick visualisation of the gene expression distribution across samples. Depending on how different our samples are, we may or may not see major differences between the samples.
# Convert counts data to long format for ggplot
counts_long <- counts |>
as.data.frame() |>
tibble::rownames_to_column("Gene") |>
tidyr::pivot_longer(cols = -Gene, names_to = "sample", values_to = "Counts")
# Create boxplot with ggplot2
ggplot(counts_long, aes(x = sample, y = Counts)) +
geom_boxplot() +
labs(x = "Samples", y = "Counts") +
theme_minimal() 
RNAseq raw count data can range from the 10s to 100,000s of reads mapped to a gene, which is highlighted by this figure. This indicates we should transform our data, and a good option for this is logging our data (for the purpose of visualisations).
Next we’ll log transform the data, adding 1 to avoid log(0) errors. We use log2 for counts for RNAseq - more on why later!
# Create new df that includes counts in long form and all metadata in the coldata
counts_long_metadata <- counts_long |>
left_join(coldata, by = "sample") |>
mutate(sample = fct_reorder(sample, histology))
# Create log2 counts boxplot with ggplot2
ggplot(counts_long_metadata, aes(
x = sample,
y = log2(Counts + 1),
fill = histology)) +
geom_boxplot() +
labs(x = "Samples",
y = "Log Counts +1",
title = "Count distribution") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
The log count boxplots reveal there may be some differences between our three sample groups. The average counts for the TPM (terminal phase males) group appear (we have not done any statistical analysis yet!) to be higher than the other two sample groups, F (female) and MT (mid-transitioning).
This is not a count of how many genes are expressed per sample, this is the distribution of gene expression in each sample.
A common question to ask at this point is, are the differences we can see between our two sample groups a true biological effect, or is this an experimental artifact? Often, we will not know the answer to this question. For now, it is enough to be aware of these differences and to keep them in mind moving forward.
Library size
Library size refers to the total number of reads for the entire sample. If two samples have vastly different library sizes, genes would look erroneously differentially expressed between samples e.g., if one library was twice as large as another, all genes would look to be upregulated in the larger library sample. Differences in library sizes are something that needs to be corrected before differential gene expression can be performed on our sample groups. Often this is included as part of the standard DGE pipeline, but it is worth checking the library sizes before moving forward.
In this code below, we are summing up the total counts in each column, which is how many reads mapped to or are assigned to each sample. This count is the total library size for each sample.
colSums(counts) |>
as.data.frame() |>
setNames("reads") |>
tibble::rownames_to_column("sample") |>
left_join(coldata, by = "sample") |>
mutate(sample = fct_reorder(sample, histology)) |>
ggplot(aes(x = sample, y = reads, fill = histology)) +
geom_col(colour = "black") +
labs(
y = "Reads mapped per sample",
x = "Sample",
title = "Library size") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Library sizes differ a bit across the samples and it appears (still no statistics performed yet!) there is a trend in differences between groups. This indicates it will be important to correct for library size, which can be achieved by calculating a sample-specific value called a size factor. If we take the raw read counts and divide by the size factor for each sample, the resulting values would become comparable for all samples. The way in which size factors are implemented will differ depending on the analysis method used. For some methods the size factor must be manually used to correct for library sizes. For other methods (e.g., DESeq2, which we will use shortly) automatically estimates size factors (and also corrects for library composition) and implements them into the differential gene expression analysis workflow.
Note: Library size and sequencing depth are not quite the same thing. All your samples could start out with the same sequencing depth (e.g., all have 30 M reads), but if one sample has more reads that are filtered out during trimming or have less reads assigned to reads during quantification, it will have a smaller library size than the other samples.
Principal Component Analysis
A principal component analysis (PCA) is linear dimensionality technique that can be used as a way to visualise sample-to-sample distances. In this ordination method, the data points (i.e., here, the samples) are projected onto the 2D plane such that they spread out in the two directions which explain most of the differences in the data. The x-axis is the direction (or principal component) which separates the data points the most (PC1). The amount of the total variance which is contained in the direction is printed in the axis label. The y-axis shows the next most variance (principal component 2). More principal components (PC3, PC4 and so forth) can also be extracted and plotted to investigate trends in the data.
In simple turns, we are making a plot that shows us how similar or different each sample is, taking in to account all of the gene expression. A PCA plot can help us see the biggest sources of biological variation and help determine which has more of an affect on the gene expression – whether that be the treatment type, the tissue type or something else (e.g., an unintended batch effect or outliers).
First, we need to remove all genes with low counts across all samples. We can apply a slightly stricter filtering threshold than the one we used earlier (counts_filtered <- counts[rowSums(counts) > 0, ]) where all genes with zero counts across all samples are removed. Below, this code will keep all genes where at least two columns (samples) have a count of 10 or more. This is arbitrary, you could chose to make it stricter or more lenient depending on your dataset and what you are trying to achieve.
counts_filtered <- counts[rowSums(counts >= 10) >= 2, ]Next the counts object needs to be flipped for prcomp() analysis.
Our counts dataframe currently looks like this:
Sample 1 Sample 2 Sample 3
geneA 1 3 4
geneB 46 986 60
geneC 0 0 2
But we need to transpose it to look like this:
geneA geneB geneC
Sample 1 1 3 4
Sample 2 46 986 60
Sample 3 0 0 2
Use the function t() to transpose the table. Raw RNAseq counts should not be used for PCA due to issues with library size difference (which we saw above) and heteroscedasticity (variance is dependant on the mean). There are multiple transformations we could do (log2, rlog, vst – variance stabilizing transformation) to normalise the counts. We’ll log2 the data here (later we’ll use the vst() function from the DEseq package) and we can do this at the same time as transposing. We still need to add 1 – some zeros will remain in the dataset as only genes with zeros for all samples are removed.
pca <- prcomp(t(log2(counts_filtered +1)))
str(pca)List of 5
$ sdev : num [1:15] 200.1 79.7 62.7 42.4 34.4 ...
$ rotation: num [1:19878, 1:15] 0.01863 0.00325 0.00229 0.00541 0.01484 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:19878] "LOC117813751" "ccdc51" "tma7" "LOC117815572" ...
.. ..$ : chr [1:15] "PC1" "PC2" "PC3" "PC4" ...
$ center : Named num [1:19878] 2.78 6.49 10.87 2.17 5.48 ...
..- attr(*, "names")= chr [1:19878] "LOC117813751" "ccdc51" "tma7" "LOC117815572" ...
$ scale : logi FALSE
$ x : num [1:15, 1:15] -115.4 -88.9 -175.7 -154.6 -159.4 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:15] "SI16_21G" "SI16_23G" "SI18_30G" "SI18_48G" ...
.. ..$ : chr [1:15] "PC1" "PC2" "PC3" "PC4" ...
- attr(*, "class")= chr "prcomp"The PCA object is a list with several attributes we will pull out for plotting.
What we don’t currently have loaded in is our metadata, which we’ll want for plotting. Let’s pull out the PC data, convert rownames to a column so that the sample names can be matched to the coldata dataframe using left join.
pca_data <- as.data.frame(pca$x) |>
tibble::rownames_to_column("sample") |>
left_join(coldata, by = "sample")
head(pca_data) sample PC1 PC2 PC3 PC4 PC5 PC6
1 SI16_21G -115.40201 -133.60413 2.733765 24.14914 49.71906 -13.165591
2 SI16_23G -88.86506 -127.88765 49.296473 67.37976 37.72320 -4.914023
3 SI18_30G -175.67976 -44.87182 -41.570272 -47.88684 -21.39136 12.865041
4 SI18_48G -154.56811 -39.71784 -19.771000 -33.43030 -19.17439 39.083682
5 SI18_52G -159.38555 -73.66377 -36.980099 -52.33758 -19.49532 14.777381
6 SI16_22G 59.85799 -15.38120 81.324431 72.71312 -91.71984 12.409531
PC7 PC8 PC9 PC10 PC11 PC12 PC13
1 0.9259945 18.686625 4.434831 -2.206133 -40.32875 41.052867 -1.382718
2 -2.3950814 -18.807824 -0.351286 -1.981447 37.30643 -35.657086 -3.692596
3 -25.2849768 -42.844168 -24.546734 16.239548 20.94450 30.634983 -23.511774
4 57.5562841 35.857002 16.696809 3.486320 19.49014 -2.120426 -16.625133
5 -18.2533754 6.133761 -19.033976 -14.682713 -22.02857 -31.787392 44.914579
6 -2.4049836 -4.288527 3.801546 8.398512 -17.16222 9.393724 5.035271
PC14 PC15 batchname condition histology timepoint adj.timepoint
1 -5.5988087 2.183614e-13 SI16 C F 50 0
2 2.1084595 2.029630e-13 SI16 C F 50 0
3 -4.8237442 2.397880e-13 SI18 T F 61 61
4 -2.4412137 1.880517e-13 SI18 C F 117 117
5 9.9707917 1.746582e-13 SI18 C F 117 117
6 -0.0457486 2.300565e-13 SI16 T MT 50 20
length.mm weight.g
1 172 72.6
2 173 68.9
3 170 67.0
4 170 71.4
5 161 56.3
6 178 80.0Now we have a data frame with a column each for sample name, each of the PCs, and each of our metadata variables.
We can also pull out the sdev values from our pca object into a vector to use this for plotting in a moment too:
pct_var <- round(100 * pca$sdev^2 / sum(pca$sdev^2),1)
pct_var [1] 68.7 10.9 6.7 3.1 2.0 1.3 1.1 1.0 1.0 0.9 0.9 0.8 0.8 0.7 0.0Now we are ready to plot our data, which we can do using a ggplot scatterplot. There are functions that automatically build PCA plots (e.g., ggfortify::autoplot(pca) or DESeq2::plotPCA(vsd, intgroup="condition") ), rather than building them from scratch with ggplot2, but they have slightly less flexibilty than using ggplot2.
ggplot(pca_data, aes(x = PC1, y = PC2, colour = histology, shape = batchname)) +
geom_point(size = 4) +
labs(x = paste0("PC1: ", pct_var[1], "% variance"),
y = paste0("PC2: ", pct_var[2], "% variance")) +
theme_classic()
DISCUSSION 🧠🏋️♀️ (3 mins)
What’s the next thing we might want to add to our plot from our metadata to investigate trends in our data using PCA?
SOLUTION
The next thing we should investigate is batch effect. We have two different experimental batches that our data was generated from – one in 2016 (SI16) and one in 2018 (SI18). We can change the shape of the points to show batch and we can add labels for the points too:
library(ggrepel)
ggplot(pca_data, aes(x = PC1, y = PC2, colour = histology, shape = batchname)) +
geom_point(size = 4) +
geom_text_repel(aes(label = sample), size = 3, max.overlaps = Inf) +
labs(x = paste0("PC1: ", pct_var[1], "% variance"),
y = paste0("PC2: ", pct_var[2], "% variance")) +
theme_classic()
What patterns do you see?
- Is sample SI16_22G an outlier and should we remove it?
- Do we have enough samples from each batch in each group to determine whether we are seeing batch effects?
Answer: SI16_22G is potentially worth removing. With only 5 samples per histology group, the batch split is at best 2:3, making it difficult to confidently distinguish batch effects from biological signal. The MT group in particular has only 1 SI16 sample and 4 SI18 samples. That said, samples largely cluster by histology group, suggesting histology is the stronger driver of variance – batch effects may be present but appear secondary.
Next up …
After exploring our dataset, removing any outliers and deciding which factors we may need to correct for that are causing batch effects, we are now ready to move onto differentially gene expression analysis.