The number of different genes detected when sequencing a transcriptome library increases slower and slower as reads are added: by definition the first read is always new; the second has high chances to be different from the first, etc., but after millions of them each extra read will have much higer chances to be identical to another read already sequenced. The consequence of this is that analysing the same library sequenced deeply, for instance on HiSeq, or shallowly, for instance on MiSeq, will not yield the same result in terms of gene discovery, etc. See Dave Tang's blog for a longer discussion on sequencing depth.
This tutorial is about comparing libraries on properties that are not invariant
with sequencing depth. The solution presented here is to normalise the
libraries to the same depth, that is, to use the same number of reads for
each library. There are mainly two ways: either input a fixed number of reads
in the alignment pipeline, or sub-sample a fixed number of alignments after
using all the reads. This tutorial shows how do the second solution for CAGE
data using R.
Busy people familiar with CAGE and R can skip the tutorial and read the manual
of the rrarefy command of the vegan package.
See the main README for general recommendations on how or what to prepare before running this tutorial.
The data downloaded here is the count of CAGE tags in the FANTOM5 CAGE peaks for all the Phase 1 libraries of the FANTOM 5 project. See the README file for more information on that file
wget --quiet --timestamping https://fantom.gsc.riken.jp/5/datafiles/latest/extra/CAGE_peaks/hg19.cage_peak_phase1and2combined_counts.osc.txt.gz
echo "0b288555ef51d1f1f9f04a2536d51a1d hg19.cage_peak_phase1and2combined_counts.osc.txt.gz" | md5sum -c## hg19.cage_peak_phase1and2combined_counts.osc.txt.gz: OK
The table is in Order-Switchable Columns format. The read.table
command in R will automatically discard its comment lines, set up column names
with the header=TRUE option and row names with the row.names=1 option.
The table is large: so loading the table will take time…
osc <- read.table('hg19.cage_peak_phase1and2combined_counts.osc.txt.gz', row.names=1, header=TRUE)
dim(osc)## [1] 201803 1829
The name of the libraries are long because they contain a plain English
description of their contents. We will shorten them to their identifier. For
example, counts.Adipocyte%20-%20breast%2c%20donor1.CNhs11051.11376-118A8
becomes CNhs11051. The association can be re-made using FANTOM5 SDRF
files.
colnames(osc) <- regmatches(colnames(osc), regexpr('CNhs.....', colnames(osc)))The first line, 01STAT:MAPPED, is special and contains the total number of
tags for each library. The sum of all tags in all peaks is lower than this,
because some tags are not in peaks. We will change the value of
01STAT:MAPPED so that the sum of the columns in our table is the total number
of reads for each library.
summary(t(osc["01STAT:MAPPED",]))## 01STAT:MAPPED
## Min. : 6602
## 1st Qu.: 1787336
## Median : 3229346
## Mean : 3808352
## 3rd Qu.: 5254150
## Max. :16059002
osc['01STAT:MAPPED',] <- osc['01STAT:MAPPED',] - colSums(osc[-grep('01STAT:MAPPED', rownames(osc)),])
summary(colSums(osc), digits=10)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6602 1787336 3229346 3808352 5254150 16059002
Let's see now how many CAGE peaks are detected per library. The command osc > 0 produces a data frame containing TRUE where a tag count was higher than
zero, and FALSE otherwise. In R, since TRUE equals 1 and FALSE equals
0, the colSums command applied on this data fame of TRUE/FALSE values will
then count the detected peaks. (Note: the proper name of TRUE/FALSE values is
boolean.)
summary(colSums(osc > 0), digits=10)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 360 44966 55314 54276 64576 121973
There are large variations in the number of peaks detected. Let's take the example of subcutaneous adipocytes (libraries CNhs12494, CNhs11371 and CNhs12017).
numberOfPeaks <- colSums(osc > 0)
numberOfTags <- colSums(osc)
adipocytes <- c('CNhs12494', 'CNhs11371', 'CNhs12017')
numberOfPeaks[adipocytes]## CNhs12494 CNhs11371 CNhs12017
## 27728 41048 54730
numberOfTags[adipocytes]## CNhs12494 CNhs11371 CNhs12017
## 535250 1384056 4224299
Do we see more peaks just because there were more tags ?
Here, we will remove tags from the data until each library has the same number of tags, that is, we will normalise the sequencing depth on the most shallow library.
We will use the rrarefy command of the vegan package. Unlike most R
commands that work on data frames, this command uses rows by default, so we
will transpose the table, sub-sample it, and transpose it again.
The sub-sampling removes tags randomly, so each time it is run it will never
produce exactly the same result. However, highly expressed peaks will stay
highly expressed, etc. For this tutorial, the computation is made identical
across runs by setting the random number seed with the set.seed command (and
resetting it with rm(.Random.seed)). Note: do not use set.seed() in your
project if you do not understand well the consequences.
library(vegan)## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
set.seed(1)
minNumberOfTags <- 500000
# Let's discard the libraries that do not have enough tags.
summary(numberOfTags > minNumberOfTags)## Mode FALSE TRUE
## logical 105 1724
osc.sub <- t(rrarefy(t(osc[,numberOfTags > minNumberOfTags]), minNumberOfTags))
rm(.Random.seed)
summary(colSums(osc.sub), digits=10)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5e+05 5e+05 5e+05 5e+05 5e+05 5e+05
That is all ! Now all the libraries contain the same number of tags.
normalisedNumberOfPeaks <- colSums(osc.sub > 0)
normalisedNumberOfPeaks[adipocytes]## CNhs12494 CNhs11371 CNhs12017
## 27011 28705 27507
colSums(osc.sub)[adipocytes]## CNhs12494 CNhs11371 CNhs12017
## 5e+05 5e+05 5e+05
barplot( log10( cbind( normalisedNumberOfPeaks[adipocytes]
, numberOfPeaks[adipocytes]
, numberOfTags[adipocytes])
)
, beside=TRUE
, legend=TRUE
, args.legend=list( x="topleft"
, title="library ID")
, main="Effect of sub-sampling normalisation on number of detected peaks"
, ylab="log10"
)
The number of detected peaks is now very similar between the three biological replicates !
One can also use sub-sampling to compare the expression profile of CAGE peaks, normalising for the fact that low-expressed peaks will be found in less libraries. The solution would be asking a question like “what would be the profile of a promoter if only 100 reads had aligned to it ?”. This question is related to the supplementary note number 4 of the FANTOM 5 paper (Forrest et al., 2014), on ubiquitous and tissue-restricted expression, and will be the topic of a future tutorial.
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 12 (bookworm)
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.11.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.11.0
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8
## [6] LC_MESSAGES=en_GB.UTF-8 LC_PAPER=en_GB.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] vegan_2.6-4 lattice_0.20-45 permute_0.9-7
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-162 cli_3.6.1 knitr_1.43 rlang_1.1.1 xfun_0.39 highr_0.10 jsonlite_1.8.5 htmltools_0.5.5 sass_0.4.6
## [10] rmarkdown_2.22 grid_4.3.0 evaluate_0.21 jquerylib_0.1.4 MASS_7.3-58.2 fastmap_1.1.1 yaml_2.3.7 cluster_2.1.4 compiler_4.3.0
## [19] mgcv_1.8-41 rstudioapi_0.14 digest_0.6.31 R6_2.5.1 parallel_4.3.0 splines_4.3.0 bslib_0.5.0 Matrix_1.5-3 tools_4.3.0
## [28] cachem_1.0.8