Tutorial - differential expression analysis - Biomedical Data Science Laboratory

Introdcution

RNA sequencing (RNA-seq) measures the abundance of RNA molecules (primarily mRNA) in a biological sample. By mapping sequencing reads to a reference genome, we can estimate gene expression levels based on the number of reads aligned to each gene.

Subsequently, differential expression analysis can be performed to identify genes whose expression levels differ significantly between conditions (e.g., healthy vs. disease), providing biological insights into underlying mechanisms.

Installation of STAR

STAR (Spliced Transcripts Alignment to a Reference) is a fast RNA-seq aligner that maps sequencing reads (FASTQ files) to a reference genome.

You can install STAR using conda.

conda create -n star_env -c bioconda star
conda activate star_enva

Generate index file

Before alignment, STAR requires a genome index built from:

Reference genome (FASTA) Gene annotation (GTF)

RNA-seq reads often span exon–exon junctions. STAR addresses this by incorporating splice junction information into the index using the –sjdbGTFfile and –sjdbOverhang options.

sjdbOverhang ≈ read length − 1 (e.g., for 100 bp reads → 99)

STAR --runThreadN 40 \
        --runMode genomeGenerate \
        --genomeDir /home/shared/genome_index_hg38_l100 \
        --genomeFastaFiles /home/shared/ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna \
        --sjdbGTFfile /home/shared/ncbi/genomic.gtf \
        --sjdbOverhang 99

Align pairs to the reference genome

For each sample, paired-end reads are aligned to the reference genome.

STAR outputs:

Sorted BAM file (alignment)
Gene-level counts (optional)

Although STAR can generate gene counts (–quantMode GeneCounts), many pipelines prefer using a dedicated counting tool such as featureCounts for better control and consistency.

The following is the codes that handle two samples.

STAR \
  --runThreadN 8 \
  --genomeDir /home/shared/genome_index_hg38_l100 \
  --readFilesIn inputs/hc_1_1.fastq inputs/hc_1_2.fastq \
  --outFileNamePrefix results/hc_ \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outSAMattributes NH HI AS nM

output
- hc_ReadsPerGene.out.tab
- hc_Aligned.sortedByCoord.out.bam

STAR \
  --runThreadN 8 \
  --genomeDir /home/shared/genome_index_hg38_l100 \
  --readFilesIn inputs/db_1_1.fastq inputs/db_1_2.fastq \
  --outFileNamePrefix results/db_ \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outSAMattributes NH HI AS nM

Gene-level counting with featureCounts

featureCounts assigns aligned reads to genes using the annotation file.

featureCounts \
  -T 8 \
  -p \
  -s 2 \
  -a /home/shared/ncbi/genomic.gtf \
  -o counts.txt \
  results/*.bam

Differential expression analysis with DESeq2

Gene count data can be analyzed using DESeq2 to identify differentially expressed genes.

The following is R script for DESeq2 to detect the differentially aboundant genes, plot expression levels and their statistical significances.

# Load count table
counts <- read.table("counts.txt", header=TRUE, comment.char="#")

# Extract count matrix and set row names
count_matrix <- counts[,7:ncol(counts)]
rownames(count_matrix) <- counts$Geneid

# Set conditions
coldata <- data.frame(
  condition = c("Diabetes", "Diabetes", "Diabetes", "Healthy", "Healthy", "Healthy" )
)
rownames(coldata) <- colnames(count_matrix)

# Differential expression analysis with DESeq2
library(DESeq2)

dds <- DESeqDataSetFromMatrix(
  countData = count_matrix,
  colData = coldata,
  design = ~ condition
)
dds <- DESeq(dds)
res <- results(dds)

# Ploat MA plot
plotMA(res, ylim=c(-5,5))

# Ploat volcano plot
plot(res$log2FoldChange, -log10(res$pvalue))

Tutorial – differential expression analysis