1 of 39

Systems Biology

RNA-seq data analysis

III. Quantification &

Differential expression

April 19, 2022

Bioinformatics & Genomics Lab

Department of Life Science

Hanyang University

Sang-Ho Yoon

2 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq analytic pipeline

3 of 39

Bioinformatics & Genomics, the BIG lab.

Gene/isoform quantification

- Two different levels of expression in RNA-seq analysis

- Gene-level quantification; assign all reads to a gene ignoring isoforms

🡪 FeatureCounts, htseq-count, etc.

- Transcript/isoform-level quantification; assign reads to individual RNA transcripts

🡪 RSEM, BitSeq, Cufflinks, etc.

uses probabilistic models

4 of 39

Bioinformatics & Genomics, the BIG lab.

Gene/isoform quantification

Gene-level quantification

- Gene locus is a union structure of its isoforms

- Directly count all reads aligned to exons of a gene 🡪 simple task w/ some issues

Isoform-level quantification

- Isoforms often share multiple exonic structures with each other

- Reads should be assigned to one of isoforms 🡪 probabilistic models

Incomplete annotation

MAPQ or randomly assign

Nascent transcript

Read strandedness

5 of 39

Bioinformatics & Genomics, the BIG lab.

Gene/isoform quantification

Quantification metrics

- The number of assigned reads are proxy of expression level

- Read counts of a gene/transcript are influenced by numerous factors

*Overall sequencing depth (total mapped reads)

*length of genes/transcripts

expression should be normalized with total mapped reads

and length of gene / transcript

6 of 39

Bioinformatics & Genomics, the BIG lab.

Gene/isoform quantification

Quantification metrics

- RPKM (Reads Per Kilobase of exons per Million mapped reads) for single-end

- FPKM (Fragments Per Kilobase of exons per Million mapped reads) for paired-end

1. Count the total reads in a sample and divide by 1,000,000

🡪 “per million” scaling factor

2. Divide the read counts by the “per million” scaling factor

🡪 Reads per Million (RPM); counts normalized by sequencing depth

3. Divide RPM by the length of a gene in kilobases

🡪 RPKM (1 RPKM ~ 1 copy in a cell)

*The sum of RPKM-/FPKM-normalized reads in each sample is different

- Inter-sample quantile normalization should be performed to compare expressions

7 of 39

Bioinformatics & Genomics, the BIG lab.

Gene/isoform quantification

Quantification metrics

- TPM (Transcripts Per kilobase per Million mapped reads)

- Based on the observation that approx. 1 million transcripts are found in a cell

- TPM calculates relative expression values comparable across different samples

1. Divide the read counts by the length of each transcript in kilobases

🡪 RPK, reads per kilobase

2. Count all RPK values in a sample and divide by 1,000,000

🡪 “per million” scaling factor

3. Divide RPK values by the “per million” scaling factor

🡪 TPM

Pop-up quiz:

Fill in these values to get an extra point ☺

8 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

Connect VPN and open your Shell

호스트: 166.104.118.161

포트 번호: 22

User: biguser

Password: biglab2428

9 of 39

Bioinformatics & Genomics, the BIG lab.

Connect biglab1 and navigate into your working space

$ ssh biglab1

$ cd systems_biology_2022/YOUR_DIRECTORY

$ mkdir session_07

$ cd session_07

mkdir A = make directory named A

cd A = enter the directory A

RNA-seq data analysis: Gene-level quantification

10 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

featureCounts: an ultrafast and accurate read summarization program.

featureCounts is a gene-level quantifier.

Usage in Linux command line:

ex) featureCounts -a genes.gtf -p -o output.txt input.bam

(reference or assembled GTF)

(output)

(input)

(-p: Input is paired-end reads. Count fragments rather than reads)

11 of 39

Bioinformatics & Genomics, the BIG lab.

Run featureCounts to quantify gene-level read count

RNA-seq data analysis: Gene-level quantification

ln –s ~/systems_biology_2022/TA/__human_transcriptome__/gencode.v34.gtf

ln –s ~/systems_biology_2022/TA/session_7/Naive_CD8_T_cells_rep1/star/Naive_CD8_T_cells_rep1_1_Aligned.sortedByCoord.out.bam

12 of 39

Bioinformatics & Genomics, the BIG lab.

Run featureCounts to quantify gene-level read count

mkdir 7_featureCounts

cd 7_featureCounts

featureCounts -p -a ../gencode.v34.gtf \

-o Naive_CD8_T_cells_rep1.txt \

../Naive_CD8_T_cells_rep1_1_Aligned.sortedByCoord.out.bam

RNA-seq data analysis: Gene-level quantification

13 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

14 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

Running FeatureCounts result in two output.

output.txt : main output with # of reads
output.txt.summary : summary file.

less Naive_CD8_T_cells_rep1.txt.summary

Reads mapped to genes

15 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

Running FeatureCounts result in two output.

output.txt : main output with # of reads
output.txt.summary : summary file.

less Naive_CD8_T_cells_rep1.txt

Gene ID

Chr.

(2 exons)

Start

(2 exons)

End

(2 exons)

strand

Gene length

Read counts

We want to collect this values for every genes in this file

16 of 39

Bioinformatics & Genomics, the BIG lab.

Run featureCounts for replicates and the other cell type

RNA-seq data analysis: Gene-level quantification

ln –s ~/systems_biology_2022/TA/session_7/Naive_CD8_T_cells_rep2/star/Naive_CD8_T_cells_rep2_1_Aligned.sortedByCoord.out.bam

ln –s ~/systems_biology_2022/TA/session_7/Naive_CD8_T_cells_rep3/star/Naive_CD8_T_cells_rep3_1_Aligned.sortedByCoord.out.bam

ln –s ~/systems_biology_2022/TA/session_7/Switched_memory_B_cells_rep1/star/Switched_memory_B_cells_rep1_1_Aligned.sortedByCoord.out.bam

ln –s ~/systems_biology_2022/TA/session_7/Switched_memory_B_cells_rep2/star/Switched_memory_B_cells_rep2_1_Aligned.sortedByCoord.out.bam

ln –s ~/systems_biology_2022/TA/session_7/Switched_memory_B_cells_rep3/star/Switched_memory_B_cells_rep3_1_Aligned.sortedByCoord.out.bam

17 of 39

Bioinformatics & Genomics, the BIG lab.

Run featureCounts for other replicates in Naive CD8 T cells

#replicate2

featureCounts -p -a ../gencode.v34.gtf \

-o Naive_CD8_T_cells_rep2.txt \

../Naive_CD8_T_cells_rep2_1_Aligned.sortedByCoord.out.bam

#replicate3

featureCounts -p -a ../gencode.v34.gtf \

-o Naive_CD8_T_cells_rep3.txt \

../Naive_CD8_T_cells_rep3_1_Aligned.sortedByCoord.out.bam

RNA-seq data analysis: Gene-level quantification

18 of 39

Bioinformatics & Genomics, the BIG lab.

And three replicates for Memory B cells

#replicate1

featureCounts -p -a ../gencode.v34.gtf \

-o Switched_memory_B_cells_rep1.txt \

../Switched_memory_B_cells_rep1_1_Aligned.sortedByCoord.out.bam

#replicate2

featureCounts -p -a ../gencode.v34.gtf \

-o Switched_memory_B_cells_rep2.txt \

../Switched_memory_B_cells_rep2_1_Aligned.sortedByCoord.out.bam

#replicate3

featureCounts -p -a ../gencode.v34.gtf \

-o Switched_memory_B_cells_rep3.txt \

../Switched_memory_B_cells_rep3_1_Aligned.sortedByCoord.out.bam

RNA-seq data analysis: Gene-level quantification

19 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

cd ..

cp /home/biguser/systems_biology_2022/TA/src/1.readcount.featureCounts.py .

An in-house python script summarizing read counts

20 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

python 1.readcount.featureCounts.py \

7_featureCounts/Naive_CD8_T_cells_rep1.txt,7_featureCounts/Naive_CD8_T_cells_rep2.txt,7_featureCounts/Naive_CD8_T_cells_rep3.txt,Switched_memory_B_cells_rep1.txt,Switched_memory_B_cells_rep2.txt,Switched_memory_B_cells_rep3.txt \

immune_cells.txt

Comma-separated input

21 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

cp /home/biguser/systems_biology_2022/TA/src/2.add_gene_name.py .

Gene ID > Name dictionary

Exclude “turned-off” genes

Match gene ID and gene name and filter not-expressing genes

22 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: Gene-level quantification

python 2.add_gene_name.py immune_cells.txt immune_cells.symbols.txt

input

output

Gene name & ID

23 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq analytic pipeline

24 of 39

Bioinformatics & Genomics, the BIG lab.

Differential Expression Analysis

DE analysis

- Statistical analysis to discover quantitative changes in gene expression levels

- Different experimental groups have distinct gene expression landscape

e.g., diseased vs. healthy blood, vaccinated vs. unvaccinated, etc.

- DE analysis tries to identify genes expressed at different levels between groups

*statistically

- Differentially expressed gene (DEG) analysis

🡪 DESeq2, EdgeR, etc.

- Differentially expressed transcript (DET) analysis

🡪 Bitseq, Cuffdiff, etc.

1) Normalization of sequencing libraries

2) Statistical test

25 of 39

Bioinformatics & Genomics, the BIG lab.

Differential Expression Analysis

1) Normalization of sequencing libraries

- However, these values violate statistical assumptions in DEG analysis tools

- FPKM/TPM describes expression of genes in general

FPKM/TPM handle this issue

🡪 read count values should be normalized with a different approach

Within condition variance

between condition variance

Condition A

Condition B

26 of 39

Bioinformatics & Genomics, the BIG lab.

DEG analysis using DESeq2

1) Normalization of sequencing libraries

An outlier

Averaged expression of Gene C = (33+55+200)/3 = 96

Averaged expression of Gene C = e^4.3 = 73.7

*Log-transformed data is insensitive to outliers

These average values are called

“Geometric averages”

27 of 39

Bioinformatics & Genomics, the BIG lab.

DEG analysis using DESeq2

1) Normalization of sequencing libraries

Un-expressed genes are removed

Exclude infinity values

Step4. Subtract the average log value from the log(counts)

-1.7

-4.3

-1.0

0.1

0.8

-0.8

-0.3

1.0

Step5. Calculate median values for each sample

-0.9

-0.1

0.9

Step6. Convert them into “normal numbers”

e^-0.9= 0.4

e^-0.1= 0.9

e^0.9= 2.5

Final scaling factors

*NOTE: median also mitigates outliers

28 of 39

Bioinformatics & Genomics, the BIG lab.

DEG analysis using DESeq2

1) Normalization of sequencing libraries

Step6. Convert them into “normal numbers”

e^-0.9= 0.4

e^-0.1= 0.9

e^0.9= 2.5

Final scaling factors

/0.4

/0.9

/2.5

/0.4

/0.9

/2.5

/0.4

/0.9

/2.5

29 of 39

Bioinformatics & Genomics, the BIG lab.

DEG analysis using DESeq2

2) Statistical test

- Shrinkage for dispersion estimates is performed to denoise

Shrunken variances

Condition A

Condition B

Within condition variance

between condition variance

- Data fitting using a generalized linear model

Shrink dispersion for each gene

30 of 39

Bioinformatics & Genomics, the BIG lab.

DEG analysis using DESeq2

2) Statistical test

- RNA-seq follows negative binomial distribution (overdispersion problem)

Microarray data follows Poisson distribution;

mean = variance

In RNA-seq, variance is much higher than mean

- DESeq2 uses a negative binomial (NB) model to fit the observed read counts

- Calculate P-values for log₂ fold changes using Wald test

- Multiple testing correction using Benjamini-Hochberg procedure

https://youtu.be/0b24mpzM_5M

31 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

DESeq2 Input

Read count table (parsed from FeatureCounts output)
Label table

Read count table

Label table

Sample Group

Sample1 cond1

Sample2 cond1

Sample3 cond2

Sample4 cond2

Sample1 Sample2 Sample3 Sample4

Gene1 1 5 7 4

Gene2 0 3 0 0

*Should be in same order

32 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

Or you can make using Excel

cp /home/biguser/systems_biology_2022/TA/src/3.mklabel_DESeq2.py .

33 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

python 3.mklabel_DESeq2.py immune_cells.symbols.txt

34 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

DESeq2, R code

Rscript DESeq2.R <count table> <label table> <outprefix> \

labelTable.txt

Group of interest

Background group

Rscript DESeq2.R readCountTable.txt labelTable.txt output cond1 cond2

Sample Group

Sample1 cond1

Sample2 cond1

Sample3 cond2

Sample4 cond2

/home/biguser/systems_biology_2022/TA/src/DESeq2.R

35 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

cp /home/biguser/systems_biology_2022/TA/src/DESeq2.R .

1. Read tables

2. Run DESeq2

3. Draw MA plot

4. Write DEG stats

36 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

Rscript DESeq2.R immune_cells.symbols.txt \

immune_cells.symbols.label.txt \

immune_cells_deseq2.txt \

Naive_CD8_T_cells \

Switched_memory_B_cells

ex)

Rscript DESeq2.R <count table> <label table> <outprefix> \

- Run this command line to perform DESeq2 analysis

DESeq2 output

output.pdf : MA-plot for DEG
output.txt : Significantly differentially expressed gene� (significant: padj <= 0.05)

37 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

immune_cells_deseq2.txt

Mean expression of all samples

Standard error

Wald statistic

Nominal

p-value

BH adjusted p-value

Multiple testing correction

T cell-specific genes

B cell-specific genes

...

38 of 39

Bioinformatics & Genomics, the BIG lab.

RNA-seq data analysis: DEG analysis using DESeq2

immune_cells_deseq2.pdf

MA-plot:�log fold change (M) on Y-axis�mean expression (A) on X-axis

Red dots = adjusted p-value ≤ 0.05
Blue lines = |log₂ fold-change| ≥ 2

39 of 39

Bioinformatics & Genomics, the BIG lab.

Homework

- From the results we have produced today

1) Summarize each step and explain briefly

2) Command lines & your working directory (use “pwd” command)

3) Screenshot of expression table (top n lines)

4) Screenshot of the final DESeq2 result table and MA plot

5) Select one marker gene for each cell type and briefly summarize its functions

e.g., CD8B: T cell surface marker, its function is ~

IGHA2: immunoglobulin heavy chain usu. expressed in plasma cells ~

Send Microsoft word file to sangho1130@hanyang.ac.kr by 4/24 (sun) midnight
File name; 2016251005_W5_sangho.docx

CAUTION: Never run those code on BIGLAB-MASTER.
Use below command line to log-on to biglab1 server

ssh biglab1