13: Advanced Topics in Python for Biology

Chapter 1: Getting Started with R & RStudio for Life Scientists

In computational biology and high-throughput genomics, R is the undisputed gold standard for statistical modeling, exploratory data analysis, and publication graphics. Backed by the Bioconductor project—which hosts over 2,200 specialized packages for RNA-Seq, ChIP-Seq, variant annotation, and single-cell cytometry—mastering R is an essential capability for modern biologists.

1. The R & RStudio Analytical Workspace

RStudio provides an integrated environment for reproducible genomics across four coordinated quadrants:

• Source Editor (Top-Left): Author and save reproducible scripts (.R), R Markdown documents (.Rmd), and Quarto notebooks (.qmd).
• R Console & Terminal (Bottom-Left): Interactive REPL execution with real-time feedback.
• Environment & History (Top-Right): Inspect loaded genomic matrices, S4 objects, factors, and active data frames.
• Files, Plots & Viewer (Bottom-Right): High-resolution vector plot rendering, directory navigation, and documentation manual pages.

2. Package Management: CRAN vs. Bioconductor

General statistical packages are installed from CRAN using install.packages(), whereas genomic packages are version-locked to R releases and managed via BiocManager:

# 1. Install BiocManager from CRAN
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

# 2. Install core genomics packages from Bioconductor
BiocManager::install(c("Biostrings", "GenomicRanges", "DESeq2", "ComplexHeatmap", "clusterProfiler"))

# 3. Verify installation integrity
BiocManager::valid()

3. Calculating GC Content & Nucleotide Composition

Calculating GC percentage across a genomic region in base R without external dependencies:

# Define a sample genomic DNA sequence
sample_dna <- "ATGCGATCGATCGTAGCTAGCTAGCTAGCTAATCGATCGATCGATCGTAAACCCGGGTTT"

# Split string into individual nucleotide characters
bases <- strsplit(sample_dna, "")[[1]]

# Calculate frequency counts and percentages
base_table <- table(factor(bases, levels = c("A", "C", "G", "T")))
cat("Nucleotide Frequencies:\n")
print(base_table)

# GC Percentage formula: ((G + C) / Total Bases) * 100
gc_content <- function(seq_str) {
    b <- strsplit(toupper(seq_str), "")[[1]]
    gc_count <- sum(b %in% c("G", "C"))
    round((gc_count / length(b)) * 100, 2)
}

cat("Calculated GC Content:", gc_content(sample_dna), "%\n")

4. Antiparallel Reverse Complement Generation

DNA double helices run antiparallel (5′ → 3′ vs 3′ → 5′). Generating the reverse complement sequence:

dna_reverse_complement <- function(seq_str) {
    comp_map <- c(A="T", T="A", C="G", G="C", N="N", a="t", t="a", c="g", g="c", n="n")
    b <- strsplit(seq_str, "")[[1]]
    comp <- comp_map[b]
    paste(rev(comp), collapse = "")
}

cat("Forward Strand (5'->3'):", sample_dna, "\n")
cat("Reverse Strand (5'->3'):", dna_reverse_complement(sample_dna), "\n")