1 of 35

The math of RNA-seq.

Applied Computational Genomics

https://github.com/quinlan-lab/applied-computational-genomics

Aaron Quinlan

Departments of Human Genetics and Biomedical Informatics

USTAR Center for Genetic Discovery

University of Utah

quinlanlab.org

2 of 35

Acknowledgements: Michele Busby

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

3 of 35

RNA-seq. Simple right?

Sequence RNA transcripts (usually cDNA) to understand how gene regulation varies in different cell types, genetic backgrounds, or conditions.

4 of 35

RNA-seq analysis is not as simple as it might seem

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0881-8

5 of 35

RNA-seq molecular prep

http://www.nature.com/nrg/journal/v12/n10/fig_tab/nrg3068_F1.html

6 of 35

RNA-seq alignment

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0881-8

7 of 35

RNA-seq read counting

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

The number of reads that align to each gene is an estimate of how many RNA transcripts were present in the sample for that gene.

Problems:

Paralogs (reads don't always map uniquely)

Reads that span 2 or more exons

Sampling error

8 of 35

RNA-seq alignment tools and performance

http://www.nature.com/nmeth/journal/v10/n12/pdf/nmeth.2722.pdf

9 of 35

Raw result. Counts per gene or transcript.

FeatureCounts is a useful tool for generating counts: https://www.ncbi.nlm.nih.gov/pubmed/24227677

Gene	Control	Experimental Condition
Gene A	27	27
Gene B	90	270
Gene C	280	640
Gene D	1003	3021
Gene E	3050	7342

Oooh! The experimental condition caused genes B,C,D, and E to be overexpressed!

Wow! Gene E is expressed at a 3X higher rate than Gene D!

10 of 35

Wrong. Why?

11 of 35

Raw result. Counts per gene or transcript.

Gene	Control	Experimental Condition
Gene A	27	27
Gene B	90	270
Gene C	280	640
Gene D	1003	3021
Gene E	3100	3342

Would you reach the same conclusion if you knew that the control experiment used 10 million reads while the experimental condition used 30 million reads?

Oooh! The experimental condition caused genes B,C, and D to be overexpressed!

12 of 35

Raw result. Counts per gene or transcript.

Gene	Control	Experimental Condition
Gene A	27	27
Gene B	90	270
Gene C	280	640
Gene D	1003	3021
Gene E	3100	3342

Possible explanations:

3 times as many transcripts are expressed for Gene E than Gene D.
Gene E is 3X as long as Gene D.
Gene D and E are the same length and produce the same number of transcripts, but Gene D has a close paralog that has acted as a "sponge" for ⅔ of its alignments.

13 of 35

RPKM (or FPKM): reads per kilobase of exon model per million reads

Gene	Gene Size	Control	Experimental Condition
Gene A	1 kb	27	27
Gene B	2 kb	90	270
Gene C	6 kb	280	640
Gene D	30 kb	1003	3021
Gene E	40 kb	3100	3342

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

14 of 35

RPKM: reads per kilobase of exon model per million reads

Gene	Gene Size	Control	Experimental Condition
Gene A	1 kb	27	27
Gene B	2 kb	90	270
Gene C	6 kb	280	640
Gene D	30 kb	1003	3021
Gene E	40 kb	3100	3342
Total (w/ counts for ~19000 other genes)		20,000,000	40,000,000

Step 1: Normalize (i.e., adjust) gene counts by the total amount of sequences in the experiment

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

15 of 35

RPKM: reads per kilobase of exon model per million reads

Gene	Gene Size	Control	Experimental Condition
Gene A	1 kb	27	27
Gene B	2 kb	90	270
Gene C	6 kb	280	640
Gene D	30 kb	1003	3021
Gene E	40 kb	3100	3342
Total		20,000,000	40,000,000
Millions of reads		20	40

Step 1: Normalize (i.e., adjust) gene counts by the total amount of sequences in the experiment

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

16 of 35

RPKM: reads per kilobase of exon model per million reads

Gene	Gene Size	Control (RPM)	Experimental Condition (RPM)
Gene A	1 kb	1.34	0.675
Gene B	2 kb	4.5	6.75
Gene C	6 kb	14.0	16.0
Gene D	30 kb	50.15	75.525
Gene E	40 kb	155.0	83.55
Total		20,000,000	40,000,000
Millions of reads		20	40

Step 1: Normalize gene counts the total amount of sequences in the experiment

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

17 of 35

RPKM: reads per kilobase of exon model per million reads

Gene	Gene Size	Control (RPM)	Experimental Condition (RPM)
Gene A	1 kb	1.34	0.675
Gene B	2 kb	4.5	6.75
Gene C	6 kb	14.0	16.0
Gene D	30 kb	50.15	75.525
Gene E	40 kb	155.0	83.55
Total		20,000,000	40,000,000
Millions of reads		20	40

Step 2: Normalize gene counts RPM by gene length

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

18 of 35

RPKM: reads per kilobase of exon model per million reads

Gene	Gene Size	Control (RPKM)	Experimental Condition (RPKM)
Gene A	1 kb	1.34	0.675
Gene B	2 kb	2.25	3.375
Gene C	6 kb	2.33	2.67
Gene D	30 kb	1.67	2.52
Gene E	40 kb	3.875	2.09
Total		20,000,000	40,000,000
Millions of reads		20	40

Step 2: Normalize gene counts RPM by gene length

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

19 of 35

Which gene expressed the most transcripts in each condition?

20 of 35

RPKM: reads per kilobase of exon model per million reads

Gene	Gene Size	Control (RPKM)	Experimental Condition (RPKM)
Gene A	1 kb	1.34	0.675
Gene B	2 kb	2.25	3.375
Gene C	6 kb	2.33	2.67
Gene D	30 kb	1.67	2.52
Gene E	40 kb	3.875	2.09
Total		20,000,000	40,000,000
Millions of reads		20	40

RPKM is best for within-sample comparisons of gene expression. TPM is best for inter-sample comparisons

. The measure RPKM (reads per kilobase of exon model per million reads) [1] is a within-sample normalization method that will remove the feature-length and library-size effects. This measure and its subsequent derivatives FPKM (fragments per kilobase of exon model per million mapped reads), a within-sample normalized transcript expression measure analogous to RPKs, and TPM (transcripts per million) are the most frequently reported RNA-seq gene expression values. It should be noted that RPKM and FPKM are equivalent for SE reads and that FPKM can be converted into TPM using a simple formula [38]. The dichotomy of within-sample and between-sample comparisons has led to a lot of confusion in the literature. Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads. Furthermore, programs such as Cufflinks that estimate gene length from the data can find significant differences in gene length between samples that cannot be ignored. TPMs, which effectively normalize for the differences in composition of the transcripts in the denominator rather than simply dividing by the number of reads in the library, are considered more comparable between samples of different origins and composition but can still suffer some biases. These must be addressed with normalization techniques such as TMM

21 of 35

Differential expression between samples

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4728800/pdf/13059_2016_Article_881.pdf

Limma

22 of 35

Poisson in differential (bulk) RNA expression

https://www.youtube.com/watch?v=KuxRYdG1EdQ

# Simulate 10,000 genes

N_genes = 10000

# Simulate a range of 10000 lambdas to reflect mean

# expression of each of the 10000 genes.

lambdas = seq(1,1024,len=N_genes) #creates 10000 evenly-spaced numbers

# between 1 and 1024

# Simulate two technical replicates for an RNA-seq

# experiment. The observed expression for each of the

# 10,000 genes will be sampled from a Poisson distribution

# using the same lambda assigned to the gene for each

# replicate.

rep1 = rpois(N_genes, lambdas)

rep2 = rpois(N_genes, lambdas)

# remove genes where the expression was 0 in either rep.

non_zero = which(rep1>0 & rep2>0)

lambdas = lambdas[non_zero]

rep1 = rep1[non_zero]

rep2 = rep2[non_zero]

# make a data frame of lambdas and expression from replicates

rna_sim = data.frame(lambdas, rep1, rep2)

23 of 35

Sources of variance. Poisson (counting) noise

Number of events

Expected frequency

https://www.youtube.com/watch?v=KuxRYdG1EdQ

# plot the expression ratio from the two replicates as a function of # the mean expression for the gene (lambda)

# plot in log2 space for better separation of the data

library(ggplot2)

ggplot(rna_sim, aes(x=lambdas, y=log2(rep1/rep2))) + geom_point(size=0.5)

Far greater random or "technical" variation in genes with low expression

24 of 35

The need for replicates.

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

Library prep sequencing

Technical replicate #1

Technical replicate #2

Biological replicate #1

Biological replicate #2

Library prep sequencing

A true technical replicate repeats sample and library prep.

https://hbctraining.github.io/Intro-to-rnaseq-hpc-salmon/lessons/experimental_planning_considerations.html#:~:text=Technical%20replicates%3A%20use%20the%20same,the%20biological%20variation%20between%20samples.

25 of 35

Sources of variance. Poisson (counting) noise

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

Gene A

Gene B

Gene C

In each case, there is an apparent 2X difference in the mean read counts between control and experimental.

Yet poisson variance is higher relative to the total count when counts are low versus when they are high. For example, the difference in expression of a gene measured with one read versus two reads is inherently less certain than the differences in expression of a gene measured with 100 reads versus 200 reads, even though both differences are nominally a 2X fold change.

26 of 35

Sources of variance. Poisson (counting) noise

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

27 of 35

Sources of variance. Non-poisson technical variance

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

"Non-Poisson Technical Variance is measurement imprecision that stems from the inability of RNA-Seq measurements to measure expression perfectly.

This imprecision is seen when expression from the same sample is measured twice. The expression measurements will not match exactly, and have error greater than what is expected from Poisson noise alone.

Sources of measurement imprecision may include PCR amplification errors during library preparation or machine errors."

Log(read counts) Technical Replicate 1

Log(read counts) Technical Replicate 2

28 of 35

Sources of variance. Biological variance

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

Biological variance is variance that naturally occurs within the samples under investigation. This variance stems from the fact that the expression of any given gene is likely to naturally fluctuate within the cells themselves, and between samples of the same condition. Sources of biological variance include genetic differences among samples and gene expression responses to the environment.

Green dots: biological replicates from S. cerevisiae

Blue dots: technical replicates from S. cerevisiae

Red dots: simulated replicates with only Poisson noise

Uncertainty in counting

Add'l uncertainty from experiment (e.g., PCR)

Add'l biological uncertainty

29 of 35

Relative contribution of sources of variance.

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

Most uncertainty is driven by biological variance -->

30 of 35

How to decrease uncertainty? 1. Increase # of reads

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

Theoretical doubling of reads

31 of 35

How to decrease uncertainty? 2. Replicates

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

32 of 35

Which is better? Depth or replicates?

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

33 of 35

Which is better? Depth or replicates?

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

34 of 35

Multiple testing

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

If you test 20,000 genes for differential gene expression, and you use a significance cut off of p<0.05, then you should expect to call approximately 1000 (i.e., 5% of 20000) genes to exhibit differential expressed solely random chance.
Thus, if your list of differentially expressed genes at p<0.05 is about 1000 genes long, then either there are no genes differentially expressed between the two conditions, or your experiment is underpowered.

35 of 35

False discovery rate

http://michelebusby.tumblr.com/post/26913184737/thinking-about-designing-rna-seq-experiments-to

It is informative to report a False Discovery Rate (FDR) as well as a p-value. This value is the expected proportion of false positives among all of the significant results (100% in the previous example).

Usually we strive for a FDR at 10% or lower. Importantly, the FDR can only be calculated after an experiment is run, because it requires knowing how many genes were called differentially expressed.