3: Control Flow, Logic & Functions in R

Chapter 3: Control Flow, Logic & Functions in R

While standard procedural scripts use explicit loops, idiomatic R prioritizes vectorization and functional programming (the apply family). Mastering vectorized logic speeds up genome-wide computations by orders of magnitude.

1. Vectorized Conditionals with ifelse()

Classifying Phred quality scores across sequencing bases in a single vectorized pass:

phred_scores <- c(38, 22, 14, 40, 28, 35, 18, 33)

quality_tier <- ifelse(phred_scores >= 30, "Q30+ (High Accuracy, <= 0.1% error)",
                ifelse(phred_scores >= 20, "Q20+ (Acceptable, <= 1.0% error)",
                                           "Low Quality (Needs Trimming)"))

for (i in seq_along(phred_scores)) {
    cat(sprintf("Base %d [Q%02d]: %s\n", i, phred_scores[i], quality_tier[i]))
}

2. Functional Programming: The Apply Family in Genomics

Operating across rows (genes) or columns (samples) in count matrices without manual loops:

# Calculate gene variance across 100 samples
gene_vars <- apply(counts_matrix, 1, var)

# Calculate sample median read counts
sample_medians <- apply(counts_matrix, 2, median)

# Normalize sample columns with sweep
size_factors <- c(1.1, 0.95)
norm_counts <- sweep(counts_matrix, 2, size_factors, "/")

3. Custom Bioinformatics Functions: Genetic Code Translator

Building reusable functions with input validation:

translate_dna <- function(dna_sequence, to_stop = TRUE) {
    codon_map <- c(
        "TTT"="F", "TTC"="F", "TTA"="L", "TTG"="L",
        "TCT"="S", "TCC"="S", "TCA"="S", "TCG"="S",
        "TAT"="Y", "TAC"="Y", "TAA"="*", "TAG"="*",
        "TGT"="C", "TGC"="C", "TGA"="*", "TGG"="W",
        "CTT"="L", "CTC"="L", "CTA"="L", "CTG"="L",
        "ATG"="M", "ACT"="T", "AAT"="N", "AAA"="K",
        "GTT"="V", "GCT"="A", "GAT"="D", "GGA"="G"
    )
    
    seq_clean <- toupper(gsub("[^ACGT]", "", dna_sequence))
    n <- nchar(seq_clean)
    num_codons <- floor(n / 3)
    codons <- substring(seq_clean, seq(1, by = 3, length.out = num_codons),
                                   seq(3, by = 3, length.out = num_codons))
    
    aa <- character(num_codons)
    for (i in seq_along(codons)) {
        res <- codon_map[codons[i]]
        if (is.na(res)) res <- "X"
        if (to_stop && res == "*") break
        aa[i] <- res
    }
    paste(aa[aa != ""], collapse = "")
}

cat("Peptide:", translate_dna("ATGGCTTTTAAATGA"), "\n")