Use Case: Genomic Data Smoothing

Overview

LOESS is well-suited for genomic data such as DNA methylation profiles, ChIP-seq signals, and expression arrays.

Methylation Profile Smoothing

The Challenge

DNA methylation data (from bisulfite sequencing or arrays) shows position-dependent patterns that can be obscured by measurement noise.

Solution

A small fraction = 0.1 lets LOESS follow fine-scale spatial structure without smearing the transitions between methylated and unmethylated regions. intervals = intervals_opts(confidence = 0.95) produces uncertainty bands that naturally widen at positions with sparser CpG coverage, making low-confidence segments immediately apparent in the plot.

library(rfastloess)

set.seed(42)
n <- 1000
positions <- sort(runif(n, 0, 1e6))

true_meth <- 0.5 + 0.3 * sin(positions / 1e5)
observed  <- true_meth + rnorm(n, sd = 0.15)
observed  <- pmax(0, pmin(1, observed))

model <- Loess(
    fraction = 0.1,
    iterations = 3,
    intervals = intervals_opts(confidence = 0.95)
)
result <- fit(model, positions, observed)

plot(positions, observed, pch = ".", col = "gray",
    xlab = "Genomic Position (bp)", ylab = "Methylation Level",
    main = "Methylation Profile Smoothing")
lines(result$x, result$y, col = "blue", lwd = 2)
lines(result$x, result$confidence_lower, col = "blue", lty = 2)
lines(result$x, result$confidence_upper, col = "blue", lty = 2)
legend("topright", c("Observed", "Smoothed", "95% CI"),
        pch = c(1, NA, NA), lty = c(NA, 1, 2),
        col = c("gray", "blue", "blue"))

ChIP-seq Signal Smoothing

Application

ChIP-seq experiments produce sparse, noisy coverage data. LOESS can help identify binding regions.

fraction = 0.05 provides high spatial resolution — important for resolving narrow binding peaks that would otherwise be smeared into the background. The larger iterations = 5 is deliberate: Poisson-distributed read counts produce tall, isolated spikes, and extra robustness iterations progressively down-weight them so the estimated background level is not inflated by a handful of extreme counts.

library(rfastloess)
set.seed(123)

positions <- seq(0, 9990, by = 10)
n <- length(positions)

background <- 10
peak1 <- 50 * exp(-((positions - 2000)^2) / (2 * 200^2))
peak2 <- 80 * exp(-((positions - 5000)^2) / (2 * 300^2))
peak3 <- 40 * exp(-((positions - 8000)^2) / (2 * 150^2))

true_signal <- background + peak1 + peak2 + peak3
observed <- rpois(n, true_signal)

model <- Loess(
    fraction = 0.05,   # Very local smoothing,
    iterations = 5,    # Strong robustness,
    outputs = "residuals"
)
result <- fit(model, positions, observed)

# Identify peaks (smoothed signal significantly above background)
threshold <- quantile(result$y, 0.75)
peaks <- positions[result$y > threshold]
cat(sprintf("Found %d peak positions (first 5): %s\n",
            length(peaks), toString(head(peaks, 5))))
#> Found 250 peak positions (first 5): 1630, 1640, 1650, 1660, 1670

Large Genome Coverage (Streaming)

For whole-genome data that doesn’t fit in memory:

library(rfastloess)
set.seed(42)

positions <- seq(0, 9990, by = 10)
coverage  <- rpois(length(positions), 50)

# Process chromosome-by-chromosome or in chunks
model <- StreamingLoess(
    fraction   = 0.05,
    chunk_size = 100000,   # 100kb chunks,
    overlap    = 10000,    # 10kb overlap,
    merge_strategy = "weighted_average"
)
process_chunk(model, positions, coverage)
#> <LoessResult>
#>   Points:            0 
#>   Fraction Used:     0.05 
#>   Iterations Used:   3
result <- finalize(model)
cat("Smoothed y[0]:", result$y[1], "\n")
#> Smoothed y[0]: 56.83401

Best Practices for Genomic Data

Consideration Recommendation
Fraction 0.05–0.15 (preserve local features)
Iterations 3–5 (handle sequencing outliers)
Large data Use streaming mode
Sparse regions Use boundary_policy = "extend"
Multiple chromosomes Process separately or ensure sorted

See Also

  • vignette("concepts", package = "rfastloess") — How LOESS works
  • ?Loess — All options
  • vignette("robustness", package = "rfastloess") — Outlier downweighting in depth
  • vignette("merge", package = "rfastloess") — Streaming chunk reconciliation
  • vignette("boundary", package = "rfastloess") — Edge handling for sparse regions
  • vignette("use-case-real-time", package = "rfastloess") — For sequencing runs
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] rfastloess_3.0.0 rmarkdown_2.32  
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39    R6_2.6.1         fastmap_1.2.0    xfun_0.61       
#>  [5] maketools_1.3.2  cachem_1.1.0     knitr_1.52       htmltools_0.5.9 
#>  [9] buildtools_1.0.0 lifecycle_1.0.5  cli_3.6.6        sass_0.4.10     
#> [13] jquerylib_0.1.4  compiler_4.6.1   sys_3.4.3        tools_4.6.1     
#> [17] evaluate_1.0.5   bslib_0.12.0     yaml_2.3.12      otel_0.2.0      
#> [21] jsonlite_2.0.0   rlang_1.3.0