R for Genomics

intermediate

R/Bioconductor Ecosystem

Bioconductor Core Packages

🔬Bioconductor

📏GenomicRangesGRanges objects

📈DESeq2Differential expression

🔍VariantAnnotationVCF read/write

🗃ComplexHeatmapPublication heatmaps

🗺clusterProfilerPathway enrichment

🔬SeuratSingle-cell analysis

📊ggplot2Publication plots

Bioconductor Ecosystem

Bioconductor (bioconductor.org) - 2,300+ packages; quarterly releases; strict quality control
GenomicRanges - GRanges objects; findOverlaps(), subsetByOverlaps(), coverage()
BSgenome - full reference genomes as R objects; getSeq() extracts sequences
Biostrings - DNA/RNA/AA sequence manipulation; pattern matching, alignment
VariantAnnotation - read/write VCF; predictCoding(), locateVariants()
rtracklayer - import/export BED, WIG, BigWig, GTF; liftOver() for coordinate conversion
ComplexHeatmap - publication-quality heatmaps; used in TCGA publications

GenomicRanges Operations

code

library(GenomicRanges)
library(rtracklayer)

# Create GRanges
gr <- GRanges(
  seqnames = c("chr1", "chr1", "chr17"),
  ranges = IRanges(start=c(1000,5000,43044294),
                   end=c(2000,6000,43125364)),
  strand = c("+","-","+")
)

# Import BED file
bed <- import("capture.bed", format="BED")

# Find overlapping variants and exons
exons <- import("gencode.v45.gtf", feature.type="exon")
hits <- findOverlaps(gr, exons)

# Coverage from BAM
library(GenomicAlignments)
bam_reads <- readGAlignments("sample.bam")
cov <- coverage(bam_reads)
mean_cov <- mean(cov[["chr17"]][43044294:43125364])

ggplot2 for Genomics Visualisation

code

library(ggplot2)
library(dplyr)

# Variant distribution across chromosomes
df %>%
  count(CHROM) %>%
  ggplot(aes(x=reorder(CHROM, -n), y=n, fill=CHROM)) +
  geom_bar(stat="identity") +
  labs(title="Variant count per chromosome", x="Chromosome", y="Count") +
  theme_bw() + theme(legend.position="none")

# CADD score distribution
df %>%
  ggplot(aes(x=CADD_phred, fill=ClinVar_class)) +
  geom_density(alpha=0.6) +
  geom_vline(xintercept=20, linetype="dashed") +
  labs(title="CADD score by ClinVar classification") +
  theme_classic()

# Heatmap with ComplexHeatmap
library(ComplexHeatmap)
Heatmap(expr_matrix,
  name = "log2(TPM+1)",
  show_row_names = FALSE,
  col = circlize::colorRamp2(c(0,6,12), c("blue","white","red")))

Python for Genomics Cloud & HPC Genomics