Lecture #3

Expression quantification

«Single-cell genomics»

Sergey Isaev

PhD student at the Medical University of Vienna

2025

Roadmap of RNA-Seq analysis

2

Reads’ QC

Alignment

Expression quantification

Pseudoalignment

Normalization

Differential expression and downstream analysis

Random primers on QC

3

“Random” sequences prime DNA not really random

High rRNA / mtDNA content

The presence of rRNA in a sample can be seen by the high number of repeated reads

4

Contamination in sample

A sample is most likely contaminated if there are several (two or more) modes in the distribution of GC-content

5

Roadmap of RNA-Seq analysis

6

Reads’ QC

Alignment

Expression quantification

Pseudoalignment

Normalization

Differential expression and downstream analysis

Mapping, not alignment

Alignment is the process of finding the best (with the best weight) arrange of two or more sequences

Sequence alignment is a bit computationally complex procedure (for fun, you can try aligning 1000 reads per human genome and see how long it takes)

Mapping is the procedure of finding of the exact position of subsequence in some sequence. Mapping is much faster than alignment, and therefore mapping algorithms are used for NGS experiments

For simplicity, we will call this alignment, but it is important to understand that the procedures are different

7

What should I use as reference?

Reference genome: long, but we don’t even need reference

Reference transcriptome: fast, but needs good reference

8

Most commonly used software for alignments

To align reads without large gaps, use Bowtie or BWA (this is relevant when working with WES and WGS)

To align reads with large gaps (the result of splicing), STAR (suffix array) and HISAT2 (bwt) are used

9

STAR

1. Works well even with large differences from the reference and is easy to use
2. Works relatively quickly
3. Requires a large amount of RAM (tens of GB — you can’t run it on a laptop)

10

STAR algorithm

First, the largest subword of the read that matches the genome is searched

​

The search is performed by representing the genome as an uncompressed suffix array (so the search is quite fast, but requires quite a lot of memory)

11

STAR algorithm

STAR then performs a similar procedure for the remained part of the reading

12

STAR algorithm

If after the second iteration the read is not completely aligned, then the presence of SNV or insertion / deletion in the unaligned parts of reads will be allowed

13

STAR algorithm

If the unaligned part is at the beginning or end of the read and does not provide a good alignment, then it will be considered a sequencing adapter and simply cut off (soft clipping procedure)

14

BAM / SAM format

15

Two-pass alignment

STAR has a --twopassMode in which alignment is done twice

​

The logic is that during the first alignment, new splice junctions may not be detected if the coverage is low enough, and some degree of reinforcement is necessary to detect new splicing events

16

Veeneman et al., Bioinf, 2016

SNP calling in RNA-Seq

In general, SNP calling from RNA-Seq can be done; various approaches are recommended for this, which are not always the best for WES, but it is better for this to identify single nucleotide polymorphisms, use WGS or WES

​

The difficulty in determining replacements lies in several details:

1. allele-specific expression is present,
2. the coverage of different positions varies greatly

17

RSeQC

There is a separate quality control stage for alignments, this is provided using the RSeQC package

​

It is very important to pay attention to the quality of the alignment, because it speaks about the quality of the sample (not the quality of the sequencing!) — as a result, FastQC can show good metrics, but the sample itself will be of very low quality (with degraded RNA, etc.)

​

The RSeQC report can be included in the general MultiQC QC report

18

Read distribution across genomic elements

The samples you analyze should have a comparable distribution of reads mapping to similar genomic elements (exons/introns/...)

19

Clipping profile

You should also be confused by large sections of adapter sequences (perhaps you have some problems with preparing the library — for example, the insert size is too small)

20

Gene body coverage

This metric has been discussed before

​

The illustration on the right shows the distribution of coverage along the length of the gene. Guess which of the samples are FFPE and which are FF?

21

Roadmap of RNA-Seq analysis

22

Reads’ QC

Alignment

Expression quantification

Pseudoalignment

Normalization

Differential expression and downstream analysis

HTSeq / featureCounts

HTSeq and featureCounts are the simplest programs used to count expression

​

Their work is based on simple logic: if the read falls on a gene, then we give +1 to gene expression

23

Liao et al., Bioinformatics, 2014 and�Andres et al., Bioinformatics, 2015

Different isoforms

24

Example: gene MS4A1 (protein CD20)

The figure above is an example of the MS4A1 gene, which encodes the CD20 protein. How to understand which isoform the reads came from?

​

You can try to find such a ratio of the number of transcript isoforms that would most likely generate the observed distribution of reads over different parts of the transcriptome

25

Example: gene PTPRC

For example, let's take the PTPRC gene, which encodes the CD45 protein. This gene has several splice isoforms, and expression of a particular isoform is an important distinguishing feature of naive and memory T cells

26

RSEM

RSEM (RNA-Seq by Expectation Maximization) optimizes isoform ratio likelihood based on gene coverage — Li and Dewey, BMC Bioinf, 2011

​

Is a standard for determining different isoforms

27

Garber et al., Nat Meth, 2011

Roadmap of RNA-Seq analysis

28

Reads’ QC

Alignment

Expression quantification

Pseudoalignment

Normalization

Differential expression and downstream analysis

kallisto

kallisto constructs a reference colored de Bruijn graph from transcriptome k-mers

​

Experimental reads are split into k-mers with the same length

​

Using the found paths in the graph, the most plausible transcript counts are found using Maximum Likelihood

29

Kallisto computational efficiency

kallisto is many times faster than most other expression counting approaches

30

Bray et al., Nat Biotechnol, 2016

Disadvantages of kallisto

Does not return a .bam file with alignment, and therefore

1. there is no possibility to search for SNVs,
2. there is no way to search for new splice isoforms

​

To work with kallisto we must be very sure of the correctness of the reference

31

Expression quantification in scRNA-Seq

Reads in .fastq format

Alignment

Expression quantification per cell�(demultiplexing is a procedure of an assignment of each read to one or another cell)

Usually it’s performed by the same tool

32

Demultiplexing problem: dropEst

For the correct reconstruction of UMI and barcode sequence sometimes (specifically with older protocols — Drop-Seq or so on) might be used dropEst

33

Cell Ranger — your safest option

• A 10x-supported tool that is suitable only for 10x libraries
• Automatically detect 10x chemistry version ⇒ more error-proof analysis
• STAR-based, so� ⇒ very computationally heavy (1 Tb disk space, 128 Gb RAM, 16 cores — recommendations)� ⇒ returns alignment in .bam-format� ⇒ works pretty slow (a few hours per sample)
• Works good with all of the 10x protocols (CITE-Seq, 10x Multiome, 10x VDJ and so on)

34

Cell Ranger

• Very simple to run:

cellranger count \� --id={id of the experiment} \� --transcriptome={path to reference genome}� --fastqs={path to directory with reads} \� --sample={sample name} \� --localcores={number of cores} \� ...

• Reference genomes are available on an official web-site with the tool

35

Cell Ranger

Important elements of the output:

1. raw_feature_bc_matrix.h5 — barcode (read: “cell”)-per-feature (read: “gene”) matrix,
2. filtered_feature_bc_matrix.h5 — filtered matrix from (1),
3. web_summary.html — html-report with read, alignment, cell and other QC metrics.

36

kallisto | bustools (bus = barcode, UMI, sequence)

• Works with any UMI- and barcode-based methods (not only 10x)
• Based on kallisto, so� ⇒ works fast and with minimal hardware requirements,� ⇒ doesn’t create .bam-file with alignment,� ⇒ works pretty fast.
• Can’t detect 10x chemistry version automatically, so you have to put coordinates of barcode and UMI by yourself
• Also works with some of the modified protocols (e. g. CITE-Seq)

37

kallisto | bustools

• Command to launch:

kb count \� -i {index file} \� -g {file that links transcripts and genes} \� -x {10x chemistry version or BC and UMI coordinates} \� {reads}

• Index might be downloaded from kallisto website ot be built by yourself
• Doesn’t automatically perform cell filtering!

38

kallisto | bustools

In the simplest case an output contains only one folder and one additional file:

• counts_unfiltered — expression count matrix
1. cells_x_genes.barcodes.txt — cells’ barcodes
2. cells_x_genes.genes.txt — gene names
3. cells_x_genes.mtx — sparse matrix with counts per each gene per each cell
• inspect.json — .json-file with short summary of pseudoalignment QC

39

kallisto | bustools and paralogs

Multimapped reads aren’t included in STAR’s expression matrices

kallisto can detect an expression of such genes with higher accuracy

40

CellRanger

Shainer and Stemmer, 2021

kallisto | bustools and paralogs

Multimapped reads aren’t included in STAR’s expression matrices

kallisto can detect an expression of such genes with higher accuracy

41

kallisto | bustools

Shainer and Stemmer, 2021

Pipelines comparison

Cell Ranger is the most computationally heavy and time consuming tool, but it’s the most convenient one, and now it’s the “golden standard” for 10x scRNA-Seq data preprocessing

42

Melsted et al. 2019