8. Bulk Transcriptomics (RNA-seq)

Salmon is an ultra-fast RNA-seq quantification engine using quasi-mapping and selective alignment. It incorporates rich statistical models that correct for sequence-specific, fragment-length, and GC-content amplification biases.

# 1. Build Salmon index with genome decoys
salmon index -t gencode.v46.transcripts.fa -d GRCh38.fa -i salmon_index --gencode

# 2. Quantify paired-end RNA-seq with automated library type and GC-bias correction
salmon quant -i salmon_index \
  -l A \
  -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
  --validateMappings \
  --gcBias --seqBias \
  -p 12 \
  -o quant_out/sample_quant

# Output transcript abundance in: quant_out/sample_quant/quant.sf

Kallisto quantifies transcript abundances from RNA-seq reads using pseudoalignment on de Bruijn graphs of the transcriptome, computing millions of reads per minute while assessing inferential variance via bootstrapping.

# 1. Build index from transcript FASTA
kallisto index -i transcripts.idx gencode.v46.transcripts.fa

# 2. Run quantification with 100 bootstrap iterations (for Sleuth differential testing)
kallisto quant -i transcripts.idx \
  -o kallisto_out/sample_quant \
  -b 100 \
  -t 12 \
  sample_R1.fastq.gz sample_R2.fastq.gz

# Output: kallisto_out/sample_quant/abundance.h5 and abundance.tsv

featureCounts is an ultrafast read summarization program written in C that assigns mapped sequencing reads to genomic features (genes, exons, promoters) from coordinate-sorted BAM files.

# 1. Count paired-end reads assigned to genes from coordinate-sorted BAMs
featureCounts \
  -T 12 \
  -p --countReadPairs \
  -B -C \
  -s 2 \
  -a gencode.v46.annotation.gtf \
  -o raw_gene_counts.txt \
  sample1.sorted.bam sample2.sorted.bam sample3.sorted.bam

# 2. View summary counting metrics
cat raw_gene_counts.txt.summary

HTSeq-count is an established Python program for counting aligned reads overlapping annotated genes, supporting union, intersection-strict, and intersection-nonempty overlap resolution modes.

# 1. Run htseq-count on coordinate-sorted BAM
htseq-count \
  --format=bam \
  --order=pos \
  --mode=union \
  --stranded=reverse \
  --minaqual=10 \
  --type=exon \
  --idattr=gene_id \
  sample.sorted.bam \
  gencode.v46.annotation.gtf > sample.htseq_counts.tsv

DESeq2 estimates variance-mean dependence in count data from high-throughput sequencing assays and tests for differential expression based on a model using the negative binomial distribution, with shrinkage of log2 fold changes for gene ranking.

library(DESeq2)
library(EnhancedVolcano)

# 1. Load raw count matrix and metadata
counts <- read.csv("raw_counts.csv", row.names=1)
coldata <- read.csv("metadata.csv", row.names=1)
coldata$condition <- factor(coldata$condition, levels=c("Control", "Treated"))

# 2. Build DESeqDataSet and pre-filter low counts
dds <- DESeqDataSetFromMatrix(countData=counts, colData=coldata, design= ~ condition)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

# 3. Fit negative binomial generalized linear model
dds <- DESeq(dds)

# 4. Extract results with apeglm effect size shrinkage
res <- results(dds, contrast=c("condition", "Treated", "Control"), alpha=0.05)
resLFC <- lfcShrink(dds, coef="condition_Treated_vs_Control", type="apeglm")
write.csv(as.data.frame(resLFC), "deg_results.csv")

# 5. Volcano Plot
EnhancedVolcano(resLFC, lab=rownames(resLFC), x='log2FoldChange', y='pvalue',
                pCutoff=0.05, FCcutoff=1.5, title="Differential Expression")

edgeR performs statistical analysis of replicated RNA-seq count data using generalized linear models and empirical Bayes dispersion estimation. It features Trimmed Mean of M-values (TMM) scaling factors and quasi-likelihood F-tests.

library(edgeR)

# 1. Create DGEList and calculate TMM normalization factors
y <- DGEList(counts = counts_matrix, group = coldata$condition)
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y, method="TMM")

# 2. Estimate dispersions and fit Quasi-Likelihood negative binomial model
design <- model.matrix(~ group, data=y$samples)
y <- estimateDisp(y, design)
fit <- glmQLFit(y, design)

# 3. Test for differential expression
qlf <- glmQLFTest(fit, coef=2)
topTags(qlf, n=20)
write.csv(as.data.frame(topTags(qlf, n=Inf)), "edger_deg_results.csv")

limma-voom transforms RNA-seq count data to log2-counts per million (logCPM) and calculates precision weights based on the empirical mean-variance relationship, allowing standard linear modeling and complex multi-factor designs.

library(limma)
library(edgeR)

# 1. Preprocess with edgeR DGEList
dge <- DGEList(counts = counts_matrix)
dge <- calcNormFactors(dge)

# 2. Apply voom transformation to model mean-variance trend
design <- model.matrix(~ 0 + condition, data=coldata)
v <- voom(dge, design, plot=TRUE)

# 3. Fit linear model and empirical Bayes moderation
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(Treated_vs_Control = conditionTreated - conditionControl, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

# 4. Extract top DE genes
topTable(fit2, sort.by="P", n=20)

Sleuth evaluates differential expression at both the transcript and gene levels for Kallisto pseudoalignment outputs, incorporating technical inferential variance derived from bootstrap resamples into generalized linear models.

library(sleuth)

# 1. Prepare Sleuth object with Kallisto directories and metadata
s2c <- read.table("metadata.tsv", header=TRUE, stringsAsFactors=FALSE)
so <- sleuth_prep(s2c, extra_bootstrap_summary=TRUE, read_bootstrap_tpm=TRUE)

# 2. Fit full and reduced models for Likelihood Ratio Test
so <- sleuth_fit(so, ~ condition, 'full')
so <- sleuth_fit(so, ~ 1, 'reduced')
so <- sleuth_lrt(so, 'reduced', 'full')

# 3. Extract differentially expressed transcripts
sleuth_table <- sleuth_results(so, 'reduced:full', 'lrt', show_all=FALSE)
head(sleuth_table)
write.csv(sleuth_table, "sleuth_deg_transcripts.csv")