TRANSCRIPTOME ANALYSIS
Background, Pipeline, HISAT2 & DESeq2
Elaborative Study Notes
[Link]. Bioinformatics | Practical Session: 07-04-2026
Includes Conceptual Questions & Answers
1. Introduction to Transcriptome & Transcriptomics
1.1 What is the Transcriptome?
The transcriptome is defined as the complete set of all RNA molecules produced by a cell under a
specific set of conditions or at a specific developmental stage. Unlike the genome — which is
largely static and the same in every cell of an organism — the transcriptome is dynamic. It
changes in response to environmental signals, developmental cues, disease states, and
experimental treatments.
Key Concept
Think of the genome as the 'instruction manual' (fixed), while the transcriptome represents
the 'pages currently being read' — these pages change depending on what the cell needs to
do at any given moment.
Understanding the transcriptome is essential for three major reasons:
• Interpreting functional elements of the genome — not all DNA sequences are expressed;
the transcriptome reveals which genes are active.
• Revealing molecular constituents of cells and tissues — different cell types express
different sets of genes, defining their identity and function.
• Understanding development and disease — comparing transcriptomes of healthy vs.
diseased cells, or across developmental time points, reveals what goes wrong in disease
and how cells change during development.
Q: Why is the transcriptome considered more informative than the genome for
understanding cellular function?
Ans: The genome is essentially the same in all cells of an organism and remains relatively
fixed. The transcriptome, however, reflects gene expression — i.e., which genes are actually
being used by a cell at a given time. It changes with environmental conditions,
developmental stage, and disease states, making it a real-time snapshot of cellular activity.
Therefore, transcriptomics directly tells us what a cell is doing, not just what it could
potentially do.
1.2 Transcriptomics — Study of RNA
Transcriptomics is the large-scale study of the entire transcriptome. It is also referred to as
expression profiling because it quantifies expression levels of RNA molecules — most commonly
mRNA — giving insight into which genes are being expressed in a cell under specific conditions
or stages.
The three primary aims of transcriptomics are:
• Cataloguing all species of transcripts — including mRNA, non-coding RNA (ncRNA), and
small RNA.
• Determining transcriptional structure — start and end sites, splicing patterns, alternative
isoforms of genes.
• Quantifying changing expression levels — measuring how transcript abundance changes
under various conditions.
Important Distinction
mRNA is just one type of RNA in the transcriptome. ncRNAs (like lncRNA, miRNA, siRNA)
also play critical regulatory roles. Modern transcriptomics aims to capture all of these.
Q: What is the difference between transcriptomics and genomics?
Ans: Genomics is the study of an organism's entire DNA — all genes and non-coding
sequences — and is largely static. Transcriptomics studies RNA molecules, specifically
which genes are expressed and at what levels, under specific conditions. Genomics tells us
the 'potential' of a cell; transcriptomics tells us the 'reality' of what is happening.
2. Basic Steps of Transcriptome Sequencing & Data Analysis
RNA-Seq (RNA Sequencing) involves a series of carefully ordered steps, from biological sample
collection all the way to bioinformatic interpretation. Each step must be performed carefully, as
errors propagate downstream.
Step 1: Sampling
Biological samples are collected from the organism or tissue of interest. In the example from the
slides, aphid species (Lipaphis erysimi and Myzus persicae) were reared inside clip-cages
attached to leaves of Brassica juncea. After 10 days, adult and nymphal insects were harvested in
liquid nitrogen (N₂) to instantly freeze and preserve RNA integrity.
Why Liquid Nitrogen?
RNA is highly unstable and degrades rapidly by enzymes called RNases. Flash-freezing in
liquid nitrogen halts all enzymatic activity instantly, preserving the transcriptome at the
moment of harvest. This is critical — even a few minutes delay can significantly alter the
gene expression profile.
Step 2: RNA Isolation
Total RNA is extracted from the biological sample. This involves cell lysis, removal of proteins and
DNA, and precipitation of RNA. The quality and purity of RNA is assessed using metrics such as
the RNA Integrity Number (RIN) and absorbance ratios (A260/A280, A260/A230).
Step 3: cDNA Library Preparation
Since RNA is single-stranded and unstable, it must be converted into complementary DNA
(cDNA) — a more stable double-stranded DNA form — before sequencing. This conversion is
done by the enzyme reverse transcriptase. The steps involved are:
• Start with mRNA or total RNA.
• Remove contaminating genomic DNA (using DNase treatment).
• Fragment RNA into smaller pieces (for efficient library construction).
• Reverse transcribe RNA into cDNA using random hexamers or oligo-dT primers.
• Ligate sequencing adaptors to cDNA ends — these short DNA sequences are essential
for the sequencing platform to recognize and process the fragments.
• Sequence the cDNA ends (paired-end sequencing generates reads from both ends of
each fragment).
Why Paired-End Sequencing?
Paired-end (PE) sequencing reads both ends of each cDNA fragment. This provides more
accurate alignment (especially across splice junctions), better detection of structural
variants, and improved assembly of transcripts. PE data requires paired R1 and R2 FASTQ
files.
Step 4: Sequencing
The library is sequenced on a next-generation sequencing (NGS) platform — typically Illumina —
generating millions of short reads (usually 100–150 bp). The output is stored in FASTQ format.
Step 5: Bioinformatic Analyses
The raw sequencing data undergoes a series of computational analyses: quality control, trimming,
alignment, assembly, quantification, normalization, and differential expression analysis. This is the
major focus of this practical.
Q: Why is cDNA made from RNA before sequencing, rather than sequencing RNA
directly?
Ans: Current NGS platforms (e.g., Illumina) are designed to sequence DNA, not RNA. RNA
is also chemically less stable than DNA and is single-stranded, making direct sequencing
difficult. By converting RNA to cDNA via reverse transcription, we get a stable double-
stranded DNA template that can be amplified, fragmented, and sequenced by standard NGS
platforms. (Note: Some newer platforms like Oxford Nanopore can directly sequence RNA,
but Illumina-based RNA-Seq still relies on cDNA.)
3. Key Terminologies
3.1 Single-End vs. Paired-End Sequencing
Single-End (SE): Only one read is generated per DNA fragment (from one end). This is simpler
and cheaper but provides less information.
Paired-End (PE): Two reads are generated per fragment — one from each end (R1 and R2). PE
reads are more expensive but provide significantly better alignment, transcript assembly, and
detection of alternative splicing.
3.2 N50 Statistic
N50 is a measure of assembly quality used particularly in de novo transcriptome assembly. It is
defined as the length of the shortest contig (or transcript) such that 50% of the total assembled
transcriptome length is contained in contigs of that length or longer.
Simple Analogy for N50
Imagine sorting all assembled transcripts from longest to shortest, and walking from the
longest end. N50 is the length of the transcript at which you have 'walked' halfway through
the total assembly. A higher N50 = longer, more contiguous assembly = better quality.
Q: If an assembly has an N50 of 1500 bp, what does it mean?
Ans: It means that 50% of the total assembled sequence is contained in transcripts that are
1500 bp or longer. In other words, if you took all transcripts ≥1500 bp and summed their
lengths, it would account for at least half the total assembly. This indicates reasonably good
assembly contiguity.
3.3 Phred Quality Score
The Phred score (Q) quantifies the confidence of a base call in sequencing. It is calculated as:
Q = -10 × log₁₀(P) where P = probability of an incorrect base
call
Phred Score Accuracy Error Rate Interpretation
Q10 90% 1 in 10 Poor — generally filtered out
Q20 99% 1 in 100 Acceptable minimum
Q30 99.9% 1 in 1,000 Good quality (standard threshold)
Q40 99.99% 1 in 10,000 Excellent quality
Exam Tip
In the pipeline shown in this practical, the fastp tool uses --qualified_quality_phred 30, which
means only bases with Q ≥ 30 are retained. This corresponds to 99.9% base call accuracy
— a very strict and high-quality threshold.
Q: A base has a Phred score of 20. What is the probability that this base call is
incorrect?
Ans: Q = -10 × log₁₀(P) → 20 = -10 × log₁₀(P) → log₁₀(P) = -2 → P = 10 ⁻² = 0.01. So
there is a 1% probability (1 in 100 chance) that this base call is wrong. This corresponds to
99% accuracy.
4. Quality Control (QC): Read Filtering with fastp
Before any analysis, raw FASTQ reads must undergo quality control. This involves:
• Assessing read quality using quality scores.
• Trimming low-quality bases from the 3' ends of reads (where quality typically degrades).
• Removing adapter sequences — short synthetic DNA sequences ligated during library
preparation. If not removed, they interfere with alignment.
• Filtering out reads that are too short after trimming.
4.1 FASTQ File Format
Each read in a FASTQ file is represented by 4 lines:
• Line 1: Header (starts with @, contains read ID and sequencing info)
• Line 2: Nucleotide sequence
• Line 3: + (separator)
• Line 4: Quality score string (ASCII-encoded Phred scores, one character per base)
4.2 fastp Command and Options
fastp is a fast, all-in-one FASTQ pre-processing tool. Below are the key options used in this
pipeline:
Option Purpose / Explanation
-i / -I Input R1 and R2 FASTQ files for paired-end sequencing
-o / -O Output filtered R1 and R2 FASTQ files after processing
--detect_adapter_for_pe Automatically detects and removes adapter sequences specific to
paired-end data. Critical — without this, adapter contamination
persists.
--correction Corrects mismatched bases in the overlapping region of paired
reads using consensus correction.
--qualified_quality_phred 30 Only retains bases with Phred score ≥ 30 (99.9% accuracy). Bases
below Q30 are trimmed.
-w 16 Use 16 CPU threads for parallel processing — speeds up
computation significantly on multi-core systems.
--html / --json Generate HTML and JSON quality control reports for visual
inspection of filtering results.
Q: Why is adapter removal necessary in RNA-Seq? What happens if adapters are not
removed?
Ans: Adapters are short synthetic DNA sequences ligated to cDNA fragments during library
preparation so the sequencer can recognize and process them. However, if the insert
(actual cDNA) is shorter than the read length, the sequencer reads into the adapter
sequence. These adapter sequences are not part of the organism's transcriptome, so during
alignment they would either fail to map (wasting data) or map incorrectly, introducing
alignment artifacts and false expression signals. Removal ensures only true biological
sequence is used.
5. Transcriptome Assembly
After quality-checked reads are obtained, the next step is to reconstruct the original transcripts.
This is called transcriptome assembly and can be done in two ways:
5.1 Reference-Based (Genome-Guided) Assembly
Reads are aligned (mapped) to an existing, known reference genome or transcriptome. The
alignment information is then used to reconstruct transcripts. This approach is faster, more
accurate, and preferred when a high-quality reference genome is available.
Tools used: HISAT2, STAR, TopHat2
When to Use Reference-Based Assembly?
Use this approach when your organism has a well-annotated reference genome. The
assembled transcriptome will be more complete and accurate because you use the genome
as a scaffold.
5.2 De Novo Assembly
Reads are assembled entirely from scratch — without any reference genome — based solely on
overlapping sequences between reads. This approach is necessary when no reference genome
exists for the organism of study.
Tools used: Trinity (most commonly used), SOAPdenovo-Trans
5.3 De Bruijn Graph — The Mathematics Behind De Novo Assembly
Trinity and most de novo assemblers use De Bruijn graph theory to reconstruct transcripts from
millions of short reads. Here is how it works:
• Each read is broken down into overlapping subsequences of length k, called k-mers. For
example, the sequence 'ATGGCGT' with k=4 gives k-mers: ATGG, TGGC, GGCG, GCGT.
• In a De Bruijn graph, nodes represent (k-1)-mers and edges represent k-mers. Each k-mer
corresponds to an edge connecting two (k-1)-mer nodes.
• The assembly problem becomes: find an Eulerian path (a path that traverses every edge
exactly once) through this graph. This is mathematically easier to solve than a Hamiltonian
path (which visits every node once).
• By traversing the Eulerian path, the original sequence is reconstructed.
Why Eulerian Path (not Hamiltonian)?
Hamiltonian paths (visiting each node once) are NP-complete — computationally extremely
hard for large datasets. Eulerian paths (using each edge once) can be solved in linear time
using efficient algorithms, making them practical for genome/transcriptome assembly with
billions of reads.
Q: What is a k-mer, and why is the choice of k value important in de novo assembly?
Ans: A k-mer is a subsequence of length k extracted from a read. For example, with k=4
and the read 'ATGGCGT', the k-mers are ATGG, TGGC, GGCG, GCGT. The choice of k is
critical: a small k leads to more connections in the graph (less specific, higher chance of
misassembly due to repeat sequences), while a large k leads to fewer connections (more
specific, but requires longer reads and deeper coverage). Trinity typically uses multiple k
values (k=25 to 32) and merges results to balance specificity and sensitivity.
6. HISAT2 — Reference Genome Alignment
HISAT2 (Hierarchical Indexing for Spliced Alignment of Transcripts 2) is a fast and sensitive
splice-aware aligner. It maps RNA-Seq reads to a reference genome while correctly handling
exon-intron boundaries (splicing).
Why 'Splice-Aware'?
mRNA in eukaryotes is produced by splicing out introns. When we generate cDNA from
mRNA and sequence it, the reads may span exon-exon junctions. A regular DNA aligner
would fail to map such reads because the intron sequence is absent from the cDNA. Splice-
aware aligners like HISAT2 and STAR recognize and correctly map reads across these
junctions.
6.1 Step 1 — Indexing the Reference Genome
Before alignment, the reference genome must be indexed using hisat2-build. This converts the
genome FASTA file into an efficient data structure (FM-index based on the Burrows-Wheeler
Transform) that allows ultra-fast lookup during alignment.
Command: hisat2-build Saff_A2_final_assembly.fa saff
• Saff_A2_final_assembly.fa — the reference genome FASTA file.
• saff — the base name prefix for the output index files (generates saff.1.ht2, saff.2.ht2 ...
saff.8.ht2).
6.2 Step 2 — Aligning Reads to the Reference
After indexing, each sample's paired-end reads are aligned to the indexed genome:
Option Purpose / Explanation
--dta Downstream Transcriptome Assembly mode — optimizes output for
transcript assemblers like StringTie by reporting soft-clipped reads
and secondary alignments appropriately.
-p 40 Use 40 CPU threads for parallel alignment — greatly speeds up
processing.
-x Prefix path to the genome index files created by hisat2-build.
-1 / -2 Input R1 and R2 FASTQ files for paired-end reads.
--summary-file Save an alignment summary report (alignment rate, error rate) for
each sample.
6.3 SAM/BAM File Conversion with Samtools
HISAT2 outputs alignments in SAM (Sequence Alignment Map) format — a human-readable text
format. For efficiency, this is converted to BAM (Binary Alignment Map) format, which is
compressed and faster to process:
• samtools view -b: Converts SAM to BAM format.
• samtools sort: Sorts BAM file by genomic coordinates — required for all downstream tools
(StringTie, DESeq2, etc.).
Q: What does 'splice-aware alignment' mean and why is it critical for RNA-Seq data?
Ans: Splice-aware alignment means the aligner can map reads that span exon-exon
junctions — sequences that exist in mRNA but not in the genomic DNA (because introns
have been spliced out). In RNA-Seq, reads are derived from processed mRNA. If a read
spans two exons with an intron between them (in the genome), a standard DNA aligner
would fail to map it. HISAT2 uses a hierarchical graph index that includes known splice
sites, allowing it to correctly identify and align junction-spanning reads. Without this, a large
proportion of reads would be unmapped, resulting in massive loss of data.
7. StringTie — Transcript Assembly & Quantification
StringTie takes the sorted BAM files from HISAT2 and assembles transcripts, then quantifies their
expression levels. It can work in two modes: de novo (finding novel transcripts) or reference-
guided (using a known annotation GTF file).
7.1 Key StringTie Options
Option Purpose / Explanation
-G [Link] Mandatory for reference-guided mode. Provides known gene
annotations to guide transcript assembly and improve accuracy.
-o [Link] Output GTF file containing assembled/predicted transcripts.
-A [Link] Output gene-level abundance table (TPM, FPKM, coverage
values).
-B Creates output in a format suitable for Ballgown R package
differential analysis.
-e Strict mode — only quantifies known reference transcripts;
suppresses novel transcript discovery. Faster and recommended
when reference annotation is complete.
--merge Merges GTF files from all samples into one unified, non-redundant
reference annotation.
-G [Link] (merge) During merging, guides the process using the original annotation to
maintain consistency.
-o [Link] Output the merged, consolidated annotation GTF.
The Two-Pass StringTie Strategy
In a rigorous pipeline: (1) Run StringTie on each sample individually to detect transcripts. (2)
Merge all sample GTF files using --merge to create a unified reference. (3) Re-run StringTie
on each sample using the merged GTF with -e flag to get consistent counts across all
samples. This two-pass approach ensures that novel transcripts found in any sample are
captured for all samples.
Q: What is the purpose of the -e flag in StringTie, and when should you use it?
Ans: The -e (expression estimation only) flag restricts StringTie to quantify only transcripts
present in the provided annotation GTF (via -G) and suppresses de novo transcript
discovery. It should be used in the second pass of the pipeline (after merging) to get a
consistent, comparable count matrix across all samples. It is faster and ensures that every
sample is quantified against the same set of transcripts — essential for downstream
statistical analysis with DESeq2.
8. Read Quantification and Normalization
After assembly, we must count how many reads map to each gene or transcript. This raw count is
then normalized to allow fair comparison between genes and between samples.
8.1 Why Normalize?
Raw read counts are biased by:
• Sequencing depth — a sample sequenced more deeply will have higher counts for all
genes, even if expression levels are the same.
• Gene length — longer genes generate more reads than shorter genes, even at equal
expression levels.
• Library composition — if some highly expressed genes dominate, others will appear
under-represented.
8.2 Normalization Methods
Method Description and Use Case
RPKM Reads Per Kilobase per Million mapped reads. Designed for single-
end reads. Normalizes for gene length first, then for sequencing
depth. NOT comparable between samples.
FPKM Fragments Per Kilobase per Million. Like RPKM but for paired-end
reads (one fragment = two reads). Also NOT directly comparable
between samples.
TPM Transcripts Per Million. Normalizes for gene length first, THEN for
sequencing depth. The sum of all TPM values in a sample is always
1 million. TPM IS comparable across samples — the gold standard
for cross-sample comparisons.
Critical Distinction: FPKM vs TPM
Both FPKM and TPM normalize for gene length and sequencing depth, but in a different
order. In FPKM, you normalize for depth first, so the total expression level can vary between
samples (making comparison unreliable). In TPM, you normalize for gene length first, then
scale so the total equals 1 million — meaning the relative proportions are directly
comparable between any two samples.
Q: Why is TPM preferred over FPKM for comparing expression across different
samples?
Ans: In FPKM, the denominator (total mapped reads) varies between samples, making the
values mathematically incomparable across samples — a gene with FPKM=10 in sample A
might not have the same relative expression as FPKM=10 in sample B. TPM corrects this by
first dividing by gene length (to get reads per kb), then dividing each value by the sum of all
length-normalized values in that sample and multiplying by 10⁶. This ensures all TPM
values within a sample sum to exactly 1 million, making the proportions directly comparable
between samples regardless of sequencing depth or library composition.
9. DESeq2 — Differential Expression Analysis
DESeq2 is an R/Bioconductor package used to identify genes that are significantly differentially
expressed between two or more experimental conditions. It uses raw read counts (not normalized
values like TPM) as input and handles normalization internally using a method based on negative
binomial distribution.
9.1 Preparing the Count Matrix: [Link]
The Python script [Link] converts StringTie output into count matrices suitable for DESeq2:
Option Purpose
-i merged_file.txt Input file: lists sample names and paths to their StringTie GTF files.
-l 151 Read length (151 bp in this experiment) — used for accurate count
estimation.
-g gene_count_matrix.csv Output: gene-level count matrix (rows = genes, columns =
samples).
-t transcript_count_matrix.csv Output: transcript-level count matrix (for isoform-level analysis).
-v Verbose mode — prints processing log for debugging.
9.2 DESeq2 Analysis — Step 1: Creating the DESeqDataSet
In R, the DESeq2 analysis begins by creating a DESeqDataSet object:
• DESeqDataSetFromMatrix(): Creates the dataset from the count matrix. Requires: (1)
countData — the gene count matrix; (2) colData — a dataframe with sample metadata
(e.g., condition labels); (3) design — a formula specifying the experimental variable to test
(e.g., ~ condition).
• DESeq(): Runs the full DESeq2 analysis pipeline — estimates size factors (normalizes for
sequencing depth), estimates dispersion (biological variability), fits a negative binomial
model, and performs Wald tests for differential expression.
• saveRDS(): Saves the DESeqDataSet object to disk (.rds file) for future use without re-
running the analysis.
9.3 DESeq2 Analysis — Step 2: Pairwise Comparisons
DESeq2 performs pairwise comparisons between all condition pairs:
• results(dds, contrast = c('condition', conditionB, conditionA)): Extracts differential
expression results comparing conditionB vs. conditionA.
• res[order(res$padj), ]: Sorts results by adjusted p-value (most significant first).
• [Link](): Saves results to CSV files (e.g., DESeq2_conditionB_vs_conditionA.csv).
9.4 Understanding DESeq2 Output Columns
Column Meaning
ENTREZID Official NCBI gene identifier number.
SYMBOL Gene symbol (short name, e.g., Egf).
GENENAME Full descriptive name of the gene.
logFC (log2FC) Log₂ fold change. Positive = upregulated in condition B; Negative =
downregulated.
AveExpr Average expression level across all samples (log₂ scale).
t / stat Test statistic from Wald test.
[Link] Raw (unadjusted) p-value for differential expression.
[Link] (padj) Adjusted p-value (corrected for multiple testing using Benjamini-
Hochberg). Use this for significance calls. Typically padj < 0.05.
Why Adjusted P-value?
In a typical RNA-Seq experiment, we test 10,000–30,000 genes simultaneously. If we use a
5% significance threshold with raw p-values, we'd expect ~500–1500 false positives by
chance alone (multiple testing problem). The Benjamini-Hochberg correction controls the
False Discovery Rate (FDR), ensuring that among all genes called significant, only ~5% are
expected to be false positives.
Q: A gene has logFC = -3 and padj = 0.001. What does this tell you?
Ans: logFC = -3 means the gene is 2³ = 8-fold downregulated in the test condition
compared to the reference condition (i.e., expression is 8 times lower). padj = 0.001 (well
below the standard 0.05 threshold) means this result is highly statistically significant after
correcting for multiple testing. Therefore, this gene is significantly and substantially
downregulated.
10. Visualizing DESeq2 Results
10.1 PCA Plot (Principal Component Analysis)
A PCA plot reduces the high-dimensional expression data (thousands of genes) into 2 principal
components (PC1 and PC2) and visualizes how samples cluster together. Each point represents
one sample:
• Samples from the same condition should cluster tightly together (low biological variability).
• Samples from different conditions should be well-separated.
• PC1 typically explains the most variance (68% in the example), PC2 the second most
(13%).
Interpreting the PCA
In the example, D2-20 and D2-30 DAP (days after pollination) cluster together on the right
side of PC1, while D2-5 and D2-10 cluster on the left. This shows that gene expression
changes significantly across developmental stages. Replicate clustering (DAP1, DAP2,
DAP3) within each group confirms good experimental reproducibility.
10.2 Heatmap
A heatmap displays expression levels of multiple genes across multiple samples simultaneously
as a color-coded matrix:
• Rows = genes; Columns = samples.
• Red = high expression; Blue = low expression (or vice versa, depending on color scheme).
• Genes and samples are often clustered hierarchically to reveal patterns.
• Useful for visualizing the top differentially expressed genes and identifying expression
patterns across conditions.
10.3 Volcano Plot
A volcano plot displays all tested genes simultaneously as a scatter plot:
• X-axis: log₂ fold change (logFC) — effect size.
• Y-axis: -log₁₀(p-value) — statistical significance (higher = more significant).
• Top-right (dark red): significantly upregulated genes (high logFC, low p-value).
• Top-left (blue): significantly downregulated genes (low logFC, low p-value).
• Centre/bottom (grey): non-significant genes.
Reading the Volcano Plot
The 'volcano' shape arises because highly significant genes tend to also have large fold
changes. Genes of interest can be labeled — in the example, Egf (epidermal growth factor)
is a highly significant downregulated gene, and genes like Ctgf and Lif are significantly
upregulated.
Q: In a volcano plot, a gene is located at position (x = 0.3, y = 1.1). Is this gene likely
to be biologically meaningful?
Ans: No. x = 0.3 means only 2^0.3 ≈ 1.23-fold change — almost no change in expression. y
= 1.1 means -log₁₀(p) = 1.1, so p ≈ 0.08 — not statistically significant (above the typical
threshold of p < 0.05, i.e., y < 1.3). This gene would be in the 'not significant' grey zone of
the volcano plot and is unlikely to be biologically meaningful.
11. Complete Pipeline Summary
The complete bioinformatics pipeline for reference-based transcriptome analysis using HISAT2
and DESeq2:
Step Tool / Action Input Output
1. Quality Control fastp Raw FASTQ (R1, R2) Filtered FASTQ (R1, R2) +
QC Report
2. Genome hisat2-build Reference genome FASTA HISAT2 index files (.ht2)
Indexing
3. Alignment HISAT2 + Filtered FASTQ + Index Sorted BAM files
samtools
4. Transcript StringTie (Pass 1) BAM + Annotation GTF Per-sample GTF files
Assembly
5. Merge StringTie --merge All sample GTF files Merged master GTF
Annotations
6. StringTie (Pass 2, BAM + Merged GTF Re-quantified GTF files
Requantification -e)
7. Count Matrix [Link] Sample GTF files gene/transcript count CSV
8. DE Analysis DESeq2 (R) Count matrix + metadata DE gene tables (CSV)
9. Visualization R (ggplot2, DE results + count data PCA plot, Heatmap, Volcano
DESeq2) plot
12. Quick Revision — Important Points for Exams
• The transcriptome is dynamic and condition-dependent; the genome is static.
• Paired-end sequencing reads from both ends of a fragment — better for alignment and
assembly.
• FASTQ files store sequence + per-base Phred quality scores. Q30 = 99.9% accuracy.
• fastp removes adapters, trims low-quality bases, and generates QC reports. --
qualified_quality_phred 30 is used.
• Two assembly strategies: Reference-based (HISAT2) when reference genome exists; De
novo (Trinity) when no reference is available.
• De Bruijn graphs are used by Trinity for de novo assembly. Uses (k-1)-mers as nodes and
k-mers as edges. Solves Eulerian path problem (efficient).
• HISAT2 is splice-aware — correctly maps reads spanning exon-intron boundaries. Uses --
dta flag for downstream StringTie compatibility.
• StringTie assembles transcripts and quantifies expression. --merge combines all sample
annotations. -e flag = strict reference-only quantification.
• FPKM/RPKM = per-gene normalization, NOT cross-sample comparable. TPM = cross-
sample comparable (sums to 1 million).
• DESeq2 uses raw counts (not TPM). Internally normalizes using size factors. Models data
with negative binomial distribution.
• Key DESeq2 output columns: logFC (fold change), AveExpr, [Link], [Link] (use padj
< 0.05 for significance).
• Positive logFC = upregulated; Negative logFC = downregulated.
• Multiple testing correction (Benjamini-Hochberg) is applied to control False Discovery Rate
(FDR).
• PCA plot: shows sample clustering and reproducibility. Volcano plot: shows all genes with
fold-change vs significance. Heatmap: shows expression patterns across samples.
— END OF NOTES —
[Link]. Bioinformatics | Transcriptome Analysis | 07-04-2026