RNA-seq and Transcriptomics:

How to experimentally determine gene expression?

​

​

​

​

Why do we sequence RNA?

Functional studies

• Genome is the same for the organism but different cells have different genes expressed;
• Experiment conditions may have a pronounced effect on gene expression (e.g. Drug treated vs. untreated cell line; wild type versus knock out)

Predicting transcript sequence from genome sequence is impossible at the current stage

• Alternative splicing, RNA editing, etc.

Some molecular features exist only for RNA

• Alternative isoforms, fusion transcripts, RNA editing, etc.

RNA-seq

• RNA (quite often poly A+) converted to a library of cDNA fragments with adaptors
• Molecules sequenced from one end (Single End) or both ends (Pair End)
• Reads of 30-400bp depending on sequence technology

Challenges

• Samples
• Purity, quantity, quality
• RNA is fragile compared to DNA (easily degraded)
• Most eukaryptic RNAs are spliced (consist of exons separated by introns in DNA)
• Mapping reads to genome is a problem
• The relative fraction of particular RNAs vary wildly
• 1 – 10^7 molecules per cell (lncRNA vs ribosomal RNA)
• As sequencing is based on sampling, highly expressed genes covered by majority of reads (ribosomal genes)
• Cannot estimate the overall decrease/increase of transcription
• RNAs come in a wide range of sizes
• Small RNAs must be captured separately

rRNA depletion methods

Типичный параметр: RIN (RQI), DV200

Common aims of RNA-Seq analysis (what can you ask of the data?)

• Gene expression and differential expression
• Isoforms expression, alternative splicing
• Novel transcripts discovery and annotation
• Allele specific expression
• Link to known SNPs or mutations
• SNP/mutation discovery
• Fusion detection
• RNA editing

Какие требования к эксперименту?

Сколько нам нужно чтений?

​

Expression Profiling / Differential expression 5-10 Million

Alternative splicing, quantifying cSNPs 50-100 Million

De Novo Transcriptome Assembly 100-1000 Million

​

​

Какие могут быть ошибки и как от них избавиться?

Replicates

Multiple isolations of cells showing the same phenotype, stage or other experimental condition (Environmental Factors, Growth Conditions, Time)

Correlation Coefficient 0.9

Типы реплик

DOI:10.11613/BM.2021.010502

Что дальше делаем с чтениями?

Three RNA-seq mapping strategies

Diagrams from Cloonan & Grimmond, Nature Methods 2010

https://www.nature.com/articles/nmeth.2714

De novo assembly

Align to transcriptome

Align to reference genome

Alignment

1

15

Module 2

• Uses a reference genome/transcriptome to map reads
• Capable of some novel transcript inference
• Relatively fast runtime
• Tools: HISAT2, STAR, GSNAP, .... and many others

Kim et al. 2015. Nat Methods 12:357–360

rnabio.org

Assembly

• Infer transcript structure directly from the data
• Useful when you do not have a reference sequence
• Other uses – highly rearranged genomes (some cancers)
• Computationally expensive
• Tools: Trinity, Velvet, SPAdes

1

16

Module 2

Haas, et al (2013) doi: https://www.nature.com/articles/nprot.2013.084

Splicing

Module 5

17

rnabio.org

1

Splicing

Module 5

18

rnabio.org

1

Pseudoalignment

1

19

Module 2

• Does not determine where in the genome a read lies, only which transcripts it is compatible with
• Does not produce a bam by default (though pseudo-bams can be created), not useful for variant detection.
• Tools:
• Kallisto, Sailfish, Salmon (one of modes)
• The comparison of the sequencing reads to the transcripts is done using a transcriptome de Bruijn graph
• the graph is constructed from the k-mers present in an input transcriptome as opposed to reads which is done normally for genome/transcriptome assembly.

Bray, 2016 doi:10.1038/nbt.3519 https://tinyheero.github.io/2015/09/02/pseudoalignments-kallisto.html

rnabio.org

Which alignment strategy is best?

• De novo assembly
• If a reference genome does not exist for the species being studied
• If complex polymorphisms/mutations/haplotypes might be missed by comparing to the reference genome

​

• Align to transcriptome
• If you have short reads (< 50bp)
• Relies on known transcripts

​

• Align to reference genome
• All other cases
• Does not rely on known transcripts – allows for discovery

​

• Each strategy involves different alignment/assembly tools

Should I use a splice-aware or unspliced mapper?

• The fragments being sequenced in RNA-seq represent mRNA - introns are removed

​

• But we are usually aligning these reads back to the reference genome

​

• Unless your reads are short (<50bp) you should use a splice-aware aligner
• HISAT2, STAR, MapSplice, etc.

HISAT2

• HISAT is a ‘splice-aware’ RNA-seq read aligner
• HISAT = Hierarchical Indexing for Spliced Alignments of Transcripts
• Requires a reference genome
• Very fast

​

• Uses an indexing scheme based on the Burrows-Wheeler transform and the Ferragina-Manzini (FM) index
• Multiple types of indexes for alignment
• a whole-genome FM index to anchor each alignment
• numerous local FM indexes for very rapid extensions of these alignments.
• Whole-genome indices with SNPs and known transcript structures accounted for

​

https://www.nature.com/articles/s41587-019-0201-4 https://pmc.ncbi.nlm.nih.gov/articles/PMC4655817/

HISAT2 algorithm

• Uses a hierarchical indexing algorithm + several adaptive strategies
• based on the position of a read with respect to splice sites
• Find candidate locations across the whole genome first
• mapping part of each read using the global FM index
• Generally identifies one or a small number of candidates.
• Do local alignment
• selects one of ~48,000 local indexes for each candidate
• uses it to align the remainder of the read.
• For paired reads, each mate is separately aligned
• If a read fails to align, then the alignments of its mate are used as anchors to map the unaligned mate

Read types

Output of HISAT2

• A SAM/BAM file
• SAM stands for Sequence Alignment/Map format
• BAM is the binary version of a SAM file

​

• Remember, compressed files require special handling compared to plain text files

​

• How can I convert BAM to SAM?
• Samtools

https://pmc.ncbi.nlm.nih.gov/articles/PMC5600148/

RNA-seq data analysis

Quality control

Good vs bad quality

Duplicated reads

Bowtie/Tophat/Cufflinks/DESeq2 �RNA-seq Pipeline

Gene annotation

(.gtf file)

Reference genome

(.fa file)

Raw sequence data

(.fastq files)

Inputs

Анализ данных RNA-seq

https://nf-co.re/rnaseq/3.21.0

1. Merge re-sequenced FastQ files (cat)
2. Auto-infer strandedness by subsampling and pseudoalignment (fq, Salmon)
3. Read QC (FastQC)
4. UMI extraction (UMI-tools)
5. Adapter and quality trimming (Trim Galore!)
6. Removal of genome contaminants (BBSplit)
7. Removal of ribosomal RNA (SortMeRNA)
8. Choice of multiple alignment and quantification routes (For STAR the sentieon implementation can be chosen):
1. STAR -> Salmon
2. STAR -> RSEM
3. HiSAT2 -> NO QUANTIFICATION
9. Sort and index alignments (SAMtools)
10. UMI-based deduplication (UMI-tools)
11. Duplicate read marking (picard MarkDuplicates)
12. Transcript assembly and quantification (StringTie)
13. Create bigWig coverage files (BEDTools, bedGraphToBigWig)
14. Extensive quality control:
• RSeQC
• Qualimap
• dupRadar
• Preseq
• DESeq2
• Kraken2 -> Bracken on unaligned sequences; optional
15. Pseudoalignment and quantification (Salmon or ‘Kallisto’; optional)
16. Present QC for raw read, alignment, gene biotype, sample similarity, and strand-specificity checks (MultiQC, R)

Bioinformatics challenges

• Splice junctions
• Several overlapping transcripts from one gene
• Non-uniquely mapped reads
• Transcripts of different length

Alignment (Tophat)

Tophat results

How do we quantify expression from RNA-seq?

RPKM: Reads per Kb million (Mortazavi et al. Nature Methods 2008)

• Longer and more highly expressed transcripts are more likely be represented among RNA-seq reads
• RPKM normalizes by transcript length and the total number of reads captured and mapped in the experiment
• Sequencing depth can alter RPKM values

What is FPKM (RPKM)?

Module 3

39

rnabio.org

1

• RPKM:reads.

Reads Per Kilobase of transcript per Million mapped

• FPKM: Fragments Per Kilobase of transcript per Million mapped reads.

​

• No essential difference - Just a terminology change to better describe paired-end reads!

What is FPKM?

Module 3

40

rnabio.org

• Why not just count reads in my RNAseq data?

• The relative expression of a transcript is proportional to the number of cDNA fragments that originate from it. However:

• # fragments is biased towards larger genes

• # fragments is related to total library depth

Fragments

Per Kilobase of transcript

per Million mapped reads.

1

What is FPKM?

Module 3

41

rnabio.org

1

• FPKM attempts to normalize for gene size and library depth
• remember – RPKM is essentially the same!

​

• C = number of mappable fragments for a gene (transcript)
• N = total number of mappable fragments in the library
• L = number of base pairs in the gene (transcript)
• FPKM = (C / (N x L) ) x 1,000 x 1,000,000
• FPKM = (1,000,000,000 x C) / (N x L)
• FPKM = (C / (N / 1,000,000)) / (L/1000)

How do FPKM and TPM differ?

Module 3

42

rnabio.org

1

• TPM: Transcript per Kilobase Million
• The difference is in the order of operations:

• The sum of all TPMs in each sample is the same. Easier to compare across samples!
• http://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/
• https://www.ncbi.nlm.nih.gov/pubmed/22872506

FPKM

1. Determine total fragment count, divide by 1,000,000 (per Million)
2. Divide each gene/transcript fragment count by #1 (Fragments Per Million)
3. Divide each FPM by length of each gene/transcript in kilobases (FPKM)

TPM

1. Divide each gene/transcript fragment count by length of the transcript in kilobases (Fragments Per Kilobase)

​

• Sum all FPK values for the sample and divide by 1,000,000 (per Million)
• Divide #1 by #2 (TPM)

https://www.youtube.com/watch?v=TTUrtCY2k-w

Differential Gene Expression Analysis

CuffDiff:

-t-test (one can set a threshold)

-replicates encouraged but not needed

-can provide differential splicing and promoter usage

​

DESeq:

- counts reads as following a negative binomial distribution

- fit a generalized linear model (GLM): more then two groups can be tested; Wald test

- models variability between replicates

https://nf-co.re/differentialabundance/1.5.0