Package {RAS}

RAS: Regional Association Score for Genome-Wide Association Studies

Description

The RAS package implements the Regional Association Score method for genome-wide association studies. It converts per-SNP effect sizes into a genomic -\log_{10}(p) profile and locates association regions on that profile with one of two detectors. The method supports both continuous and binary traits.

Option A, one call (recommended)

ras runs the complete pipeline and returns a "ras" object. Its first argument is either an in-memory genotype matrix or the path to a chunked on-disk .rasbin file (geno_to_rasbin, bed_to_rasbin); genotypes are streamed from the file in compiled code, so peak memory is set by chunk_snps rather than by the chromosome. ras_original is the pure-R in-memory implementation of 1.0.x, kept for reference. Use plot.ras to visualise results.

result <- ras(geno, phenotype, covariates,
              covariate_cols = c("age", "sex", paste0("pc", 1:10)),
              is_continuous  = TRUE,
              chrom = 1, save_dir = "results/")

bed_to_rasbin("cohort_chr1.bed", "cohort_chr1.rasbin")
result <- ras("cohort_chr1.rasbin", phenotype, covariates,
              covariate_cols = c("age", "sex", paste0("pc", 1:10)),
              is_continuous  = TRUE,
              chrom = 1, save_dir = "results/")

print(result)             # detected positions / regions
plot(result)              # full-chromosome scan profile
plot(result, zoom = TRUE) # zoomed view around each detection

Option B, step by step (advanced)

ras_scan \to ras_detect \to ras_validate \to plot.ras, or ras_scan \to ras_box_detect \to plot.ras.

## Step 1: averaged -log10(p) profile
scan     <- ras_scan(geno, phenotype, covariates,
                     covariate_cols = c("age", "sex", paste0("pc", 1:10)),
                     is_continuous = TRUE, chrom = 1, save_dir = "results/")

## Step 2 + 3, changepoint detector
detected <- ras_detect(scan$x, scan$y,
                       window_size = 3000,
                       slope.p.values.threshold.left  = 1e-10,
                       slope.p.values.threshold.right = 1e-20)
final    <- ras_validate(detected, x = scan$x, y = scan$y,
                         this.skip = 10, p.value.threshold = 1e-10)

## Step 2 + 3, box-scan detector (threshold from ras_box_calibrate())
final    <- ras_box_detect(scan$x, scan$y, calibration = cal)

## Step 4: plot
result   <- structure(
  list(scan = scan, detection = final, chrom = 1, save_dir = "results/"),
  class = "ras")
plot(result)
plot(result, zoom = TRUE)

Stage 1, scan

ras_scan performs num_rep independent 50/50 train/hold-out splits. In each repetition it fits per-SNP regressions on the training half (compute_gwas_weights) and runs a forward scan on the hold-out half (screen_forward_max_region) that slides an expanding window and records the minimum p-value at each grid position, streaming chunk_snps SNP columns at a time from the .rasbin file. The -\log_{10}(p) vectors are averaged across repetitions. For continuous traits, covariates are residualised from the training phenotype each repetition (no hold-out leakage); for binary traits, a Rao score test is used in the scan step. The pure-R in-memory implementation (ras_scan_original, compute_gwas_weights_original, compute_pgs_matrix, screen_forward_max_region_original) computes the same profile and additionally offers the per-window logistic regression (scan_test = "glm") for binary traits.

When an independent GWAS of the same trait is available, its effect sizes can replace the training split: ras_harmonize_sumstats aligns the external summary statistics to the target variant map and ras_scan_external scans all samples once with those fixed weights.

Stage 2 and 3, changepoint detector (detector = "changepoint")

ras_detect (pure-R: ras_detect_original) slides a window of width cp_window_size across the profile. At each position get_break_points fits a segmented linear model and applies the Davies test. A candidate is retained only when three conditions all hold: Davies p < cp_p_threshold, left-side slope significantly positive (slope_test, p < cp_slope_left), and right-side slope significantly negative (p < cp_slope_right). Each accepted candidate is snapped to the nearest local peak via get_local_maximum. ras_validate then re-tests each candidate with local Davies tests on both sides within a window of half-width second_window_size; a candidate is accepted if either side passes p < second_p_threshold, and candidates with -\log_{10}(p) \le min_signal at the estimated position are discarded regardless.

Stage 2 and 3, box-scan detector (detector = "box")

ras_box_detect scores every window of a set of widths by how far its mean stands above both adjacent flanks and above the profile background, in the raw -\log_{10}(p) units of the profile (ras_box_stat), and reports the intervals [\tau_L, \tau_R] whose score clears a threshold. A window inside a plateau or on the shoulder of a peak scores about zero, so broad regions are delimited at their edges and neighbouring peaks stay separate. The threshold is calibrated on null profiles with ras_box_calibrate.

Stage 4, plots

plot.ras with zoom = FALSE produces two files: the RAS scan profile with detected changepoints or intervals highlighted, and a dual-axis overlay with the Davies test significance (or the box-scan score) at every candidate. plot.ras with zoom = TRUE produces a multi-panel figure with one zoomed panel per detection.

Author(s)

Jiahe Jin jiahe.jin@yale.edu, Yiran Jiang yiran.jiang@uky.edu, Heping Zhang heping.zhang@yale.edu (maintainer)

See Also

• ras, ras_scan, ras_detect, ras_validate, ras_scan_external

• ras_original, ras_scan_original, ras_detect_original (pure-R, in-memory)

• ras_box_detect, ras_box_calibrate, ras_box_stat

• geno_to_rasbin, bed_to_rasbin, rasbin_header, rasbin_read_chunk

• ras_harmonize_sumstats, ras_weights_from_sumstats, ras_sumstats_report

• compute_gwas_weights_original, compute_pgs_matrix, screen_forward_max_region_original

• get_break_points, get_local_maximum, slope_test

• plot.ras, print.ras

• ras_memory, release_memory

Convert a PLINK 1 .bed/.bim/.fam Fileset to .rasbin

Description

Converts a PLINK 1 binary fileset directly into the chunked .rasbin format, without ever materializing a dense genotype matrix in R. Both the read from .bed and the write to .rasbin are done in fixed-size SNP chunks in C (see src/plink_io.c), so peak memory is O(\text{n\_samples} \times \text{chunk\_snps}) regardless of how many samples or SNPs the fileset contains.

Usage

bed_to_rasbin(bed_path, out_path, chunk_snps = 2000)

Arguments

Details

Only SNP-major .bed files are supported (the only mode any current PLINK version writes). Genotype dosage counts copies of the A1 allele (the fourth column of .bim), matching PLINK's own --recode A / plink2's --export A convention: 2 = homozygous A1, 1 = heterozygous, 0 = homozygous A2, NA = missing.

Unlike geno_to_rasbin (dtype 0, double, 8 bytes/genotype), this writes dtype 1 (int8, 1 byte/genotype) – 8x smaller on disk, though still ~4x larger than .bed's 2-bit packing (see src/rasbin.h for why byte-per-genotype was chosen over bit-packing). Both dtypes are read transparently by rasbin_read_chunk and every ⁠_fast⁠ pipeline function.

Value

Invisibly, the out_path that was written.

Per-SNP GWAS Effect Size Weights

Description

Fits the per-SNP regressions of RAS Stage 1 on the training half of the sample: for each variant, phenotype1 ~ dosage for a continuous trait or phenotype1 ~ dosage + covariates for a binary trait. SNP columns are streamed from the .rasbin file in chunks and only the training rows are kept, so the full genotype matrix is never held in memory. The pure-R, in-memory implementation is compute_gwas_weights_original.

Usage

compute_gwas_weights(
  rasbin_path,
  phenotype1,
  this.sample,
  this.df,
  is_continuous,
  covariate_formula = NULL,
  chunk_snps = 5000
)

Arguments

Details

Binary traits are fitted by the Frisch-Waugh-Lovell identity: the covariate projection is computed once and each SNP costs one residual regression. SNPs with a missing dosage in the training rows are refitted on their complete rows by downdating the shared covariate cross-products in blocks, so post-QC sequencing data with sparse missingness stays fast.

Value

Numeric matrix with N rows (variants) and columns Estimate, Std. Error, t value and Pr(>|t|). The Estimate column is the weight vector used by screen_forward_max_region.

See Also

ras_scan, which calls this function once per repetition; compute_gwas_weights_original.

Per-SNP GWAS Effect Size Weights (Pure-R, In-Memory)

Description

This is the pure-R, in-memory implementation that RAS 1.0.x shipped under the plain name. Since 1.1.0 the plain name (compute_gwas_weights) is the compiled, disk-backed implementation; this function is kept for reference and gives the same results. Runs a per-SNP regression on the training split to obtain effect-size estimates used as polygenic score (PGS) weights in the forward scan.

Usage

compute_gwas_weights_original(
  geno,
  phenotype1,
  this.sample,
  this.df,
  is_continuous,
  covariate_formula = NULL
)

Arguments

Details

For continuous traits (is_continuous = TRUE) the function fits lm(phenotype1 ~ this.x) for each SNP, excluding samples with missing dosage via na.omit. The covariates are assumed to have already been residualised from phenotype1 before this function is called (as done in ras_scan_original).

For binary traits (is_continuous = FALSE) the function builds a data frame combining the training-split rows of this.df with the dosage column (this.x) and phenotype (phenotype1), removes incomplete cases with na.omit, and fits lm(phenotype1 ~ covariate_formula). Note that lm rather than glm(family = binomial) is used at the GWAS weight stage to obtain a linear probability approximation to the effect size; this is consistent with standard GWAS practice for weight derivation.

Both branches are vectorised: instead of fitting one lm per SNP, all NA-free SNPs are solved in closed form with a single BLAS cross-product (continuous) or a batched Frisch-Waugh-Lovell residualisation (binary). SNP columns that contain missing dosages are handled individually so the result matches lm + na.omit exactly. A single “Starting GWAS ...” message is printed.

Value

A numeric matrix with N rows and 4 columns:

Estimate

Per-SNP effect-size estimate (slope of the dosage term). The first column is used as PGS weights by compute_pgs_matrix.

Std. Error

Standard error of the estimate.

t value

t-statistic.

Pr(>|t|)

Two-sided p-value.

See Also

compute_pgs_matrix for the next step, which multiplies these weights by the hold-out genotype matrix. ras_scan_original for the recommended high-level entry point that calls this function automatically.

Examples

set.seed(1)
n_samp <- 60; n_snp <- 10
geno <- matrix(sample(0:2, n_samp * n_snp, replace = TRUE,
                      prob = c(0.6, 0.3, 0.1)),
               nrow = n_samp, ncol = n_snp)
pheno_train <- rnorm(n_samp)

coef_mat <- compute_gwas_weights_original(
  geno              = geno,
  phenotype1        = pheno_train,
  this.sample       = seq_len(n_samp),
  this.df           = data.frame(row.names = seq_len(n_samp)),
  is_continuous     = TRUE
)
dim(coef_mat)        # n_snp x 4
head(coef_mat, 3)

Build the Per-Individual PGS Contribution Matrix

Description

Multiplies each hold-out individual's genotype dosages by the per-SNP PGS weights to produce a matrix of weighted dosage contributions used in the forward scan.

Usage

compute_pgs_matrix(geno, this.leftout, pgs.weights)

Arguments

Details

Missing genotype values (NA) in geno are imputed to zero before computing the product. This is a simple mean-imputation equivalent under a centred dosage scale and is sufficient for the forward scan, where the primary goal is to aggregate regional signals rather than obtain individual-level accuracy.

The caller's geno object is not modified: the NA-replacement assignment triggers R's copy-on-modify semantics, so a local copy of geno is made inside the function. Callers working with very large genotype matrices should be aware that this temporarily doubles the memory footprint of geno during the call; use ras_memory to check feasibility before running.

Value

A numeric matrix of dimensions |\texttt{this.leftout}| \times N, where entry [i, j] is the weighted dosage contribution of SNP j for hold-out individual i (i.e., geno[this.leftout[i], j] * pgs.weights[j]).

See Also

compute_gwas_weights_original for the step that produces pgs.weights. screen_forward_max_region_original for the step that consumes this matrix. ras_memory for pre-flight memory estimation.

Examples

set.seed(2)
geno    <- matrix(sample(0:2, 50 * 20, replace = TRUE), nrow = 50, ncol = 20)
weights <- rnorm(20)
leftout <- 1:15

pgs_mat <- compute_pgs_matrix(geno, leftout, weights)
dim(pgs_mat)   # 15 x 20

## Entry [i, j] equals geno[leftout[i], j] * weights[j]
stopifnot(pgs_mat[1, 1] == geno[leftout[1], 1] * weights[1])

Write a Genotype Matrix to the Chunked .rasbin Format

Description

Converts an in-memory genotype dosage matrix (or a matrix stored in an .rds/.csv/.tsv file) into .rasbin: a flat, column-major binary file that supports reading an arbitrary SNP-column range with a single seek, without loading the whole matrix into memory.

Usage

geno_to_rasbin(geno, out_path, chunk_write_size = 2000)

Arguments

Details

File layout (see src/rasbin.h for the authoritative spec):

• Header (32 bytes): 8-byte magic "RASBIN01", 8-byte n_samples (int64), 8-byte n_snps (int64), 4-byte dtype (int32; 0 = double), 4 reserved bytes.

• Body: n_samples * n_snps doubles, column-major, so column j (0-based) occupies bytes [32 + j*n_samples*8, 32 + (j+1)*n_samples*8).

NA dosages round-trip exactly: they are ordinary IEEE-754 doubles at the bit level, so a raw byte copy preserves R's NA_real_ sentinel without any special-casing.

Value

Invisibly, the out_path that was written.

Detect a Single Changepoint via Segmented Regression

Description

Fits a segmented linear model to detect one breakpoint in the relationship between x and y up to index t, and tests its significance using the Davies test.

Usage

get_break_points(x, y, t)

Arguments

Details

The function first fits an ordinary linear model y ~ x over the first t observations, then calls segmented to estimate a single breakpoint. If the segmented fit succeeds, the Davies test (davies.test) is applied to the base linear model to assess whether the slope change is significant.

Both the segmented call and the Davies test are wrapped in try(..., silent = TRUE): if either fails (e.g., due to collinear data or a degenerate window), the function returns p.values = 1 and NULL for the breakpoint and slopes rather than propagating an error.

The breakpoint position is returned as an index into the full x vector (not just the first t elements), found by locating the element of x[1:t] closest to the estimated breakpoint coordinate.

Value

A named list with four elements:

break.points

Integer. Index of the estimated breakpoint in x[1:t], or NULL if no breakpoint was found.

p.values

Numeric. Davies test p-value for the breakpoint, or 1 if the fit failed or no breakpoint was found.

slope.left

Numeric. Estimated slope to the left of the breakpoint, or NULL if not found.

slope.right

Numeric. Estimated slope to the right of the breakpoint, or NULL if not found.

See Also

ras_detect_original which calls this function repeatedly in a sliding-window loop. slope_test for the one-tailed slope verification step. segmented, davies.test for the underlying segmented-regression routines.

Examples

set.seed(1)
x <- 1:60
y <- c(seq(0, 6, length.out = 30),
       seq(6, 2, length.out = 30)) + rnorm(60, sd = 0.4)
result <- get_break_points(x, y, t = 60)
cat("Breakpoint index:", result$break.points, "\n")
cat("Davies p-value:  ", result$p.values,    "\n")
cat("Left slope:      ", result$slope.left,   "\n")
cat("Right slope:     ", result$slope.right,  "\n")

Find the Local Maximum Within a Window

Description

Returns the index of the maximum value of y within a symmetric window of half-width window.size centred at x0.

Usage

get_local_maximum(y, x0, window.size = 50)

Arguments

Details

This function is used by ras_detect (and ras_detect_original) after the sliding-window loop to snap each accepted changepoint index to the nearest local peak in the -\log_{10}(p) profile. Snapping to the peak ensures that reported positions correspond to the most significant SNP in the association region rather than to the mathematical breakpoint of the piecewise linear fit, which may be slightly offset.

Value

Integer. The index (into y) of the local maximum within the window centred at x0. If multiple positions tie for the maximum, which.max returns the first.

See Also

ras_detect which calls this function to refine detected changepoint positions.

Examples

y <- c(1, 3, 7, 5, 2, 8, 4, 1)

## Peak within a window of half-width 2 centred at index 3
get_local_maximum(y, x0 = 3, window.size = 2)  # returns 3 (value = 7)

## Peak within a window of half-width 3 centred at index 3
## covers indices 1:6; max is at index 6 (value = 8)
get_local_maximum(y, x0 = 3, window.size = 3)  # returns 6

Plot a RAS Result Object

Description

S3 plot method for objects of class "ras" returned by ras or ras_original. When zoom = FALSE (default) produces the full-chromosome scan profile with detected changepoints marked. When zoom = TRUE produces the zoomed multi-panel figure around each detected changepoint. For a result obtained with detector = "box" the shaded areas are the detected intervals [\tau_L, \tau_R], the vertical lines mark each interval's anchor, and the right-hand axis of the overlay plot shows the box-scan score instead of the Davies test significance.

Usage

## S3 method for class 'ras'
plot(
  x,
  zoom = FALSE,
  device = "pdf",
  p.threshold = 8,
  y_cap = NULL,
  min_display_p = 1,
  xlim = NULL,
  zoom_half_width = 3000,
  ncol = 3,
  min_signal = 2.5,
  ...
)

Arguments

Value

Invisibly returns x. Called for its side effect of writing plot files or rendering to the active graphics device.

See Also

ras for the function that produces the "ras" object. print.ras for the console summary.

Examples

## Build a minimal "ras" object by hand (see ras() for the full pipeline)
set.seed(7)
xg <- seq(1, 2000, by = 10)
yg <- c(seq(0, 9, length.out = length(xg) %/% 2),
        seq(9, 1, length.out = length(xg) - length(xg) %/% 2)) +
      rnorm(length(xg), sd = 0.4)
detection <- list(
  tau_hats         = xg[100],
  all.changepoints = xg[c(95, 100, 105)],
  all.p.values     = c(5, 12, 4),
  left.slopes      = 0.3,
  right.slopes     = -0.3
)
result <- structure(
  list(scan = list(x = xg, y = yg), detection = detection,
       chrom = 1, save_dir = tempdir()),
  class = "ras")

plot(result, device = "screen")
plot(result, zoom = TRUE, device = "screen")

Plot RAS Scan Profile with Detected Changepoints

Description

Produces two diagnostic plots for one chromosome: the RAS scan profile with detected changepoints marked, and a dual-axis overlay of the profile with Davies test significance at each candidate position.

Usage

plot_ras_scan(
  x,
  y,
  detection.result,
  this_chrom,
  save.directory,
  p.threshold = 8,
  device = "pdf",
  y_cap = NULL,
  min_display_p = 1,
  xlim = NULL,
  min_signal = 2.5
)

Arguments

Details

Plot 1: scan profile. The scan line is drawn segment-by-segment with colour determined by the local -\log_{10}(p) value:

• Red (#d73027): \ge p.threshold

• Orange (#fc8d59): \ge p.threshold / 2

• Yellow (#fee090): \ge min_signal

• Blue (#91bfdb): below min_signal

Detected changepoints (tau_hats) are marked with vertical dashed red lines, shaded bands, and labelled with their position and Davies -\log_{10}(p). Left and right slope lines (dark green and purple) are overlaid at each changepoint; their half-width is fixed at 2.5\ the total x-range. When y_cap is set, positions above the cap are shown as arrows with their true value annotated.

Plot 2: dual-axis overlay. The scan profile (left y-axis, blue line) is overlaid with triangle markers showing the Davies -\log_{10}(p) for each candidate (right y-axis, coloured by whether it exceeds p.threshold). Vertical dashed segments connect markers to the baseline.

Output files are named: chr-<chrom>-cp-plot.<ext> and chr-<chrom>-cp-p-values-plot.<ext>.

Value

Invisibly returns NULL. Called for its side effect of writing plot files or rendering to the active graphics device.

See Also

plot_ras_zoom_regions for zoomed panels around each changepoint. ras which calls this function automatically. ras_validate whose output is passed as detection.result.

Zoom-In Plots Around Each Detected Changepoint

Description

Generates a multi-panel figure with one panel per detected changepoint, each showing a zoomed view of the RAS scan profile centred on that position.

Usage

plot_ras_zoom_regions(
  x,
  y,
  detection.result,
  this_chrom,
  save.directory,
  p.threshold = 8,
  device = "pdf",
  zoom_half_width = 3000,
  ncol = 3,
  min_signal = 2.5
)

Arguments

Details

Each panel shows a [\tau - \texttt{zoom\_half\_width},\; \tau + \texttt{zoom\_half\_width}] window around changepoint \tau. The scan line is drawn in blue with colour-coded overlays for regions above the significance thresholds (min_signal, p.threshold / 2, p.threshold). Left and right slope lines are plotted in dark green and purple respectively, with their horizontal extent scaled to span approximately 45\ zoom_half_width).

If detection.result$tau_hats is empty the function returns invisibly with a message. Blank panels are added to complete the grid if the number of changepoints is not a multiple of ncol.

The output file is named chr-<chrom>-zoom.<ext>.

Value

Invisibly returns NULL. Called for its side effect of writing a multi-panel plot file or rendering to the active graphics device.

See Also

plot_ras_scan for the full-chromosome overview plots. ras which calls this function automatically.

Print a RAS Result Object

Description

Prints a one-line summary of a "ras" object showing the chromosome and detected changepoint positions.

Usage

## S3 method for class 'ras'
print(x, ...)

Arguments

Value

Invisibly returns x.

See Also

ras, plot.ras.

RAS: Regional Association Score Analysis

Description

One-call entry point of the package. Runs the RAS pipeline on one chromosome: the Stage-1 scan (ras_scan), region detection on the resulting profile with the box-scan detector (ras_box_detect, default) or the changepoint detector (ras_detect followed by ras_validate), and the diagnostic plots (plot.ras). Genotypes are streamed from a .rasbin file; an in-memory genotype matrix is accepted as well.

Usage

ras(
  geno,
  phenotype,
  covariates,
  covariate_cols,
  is_continuous,
  num_rep = 5,
  skip1 = 10,
  skip2 = 20,
  chrom = 1,
  save_dir = file.path(tempdir(), "RAS"),
  min_window_size = 5,
  max_window_size = 100,
  scan_test = c("score", "glm"),
  chunk_snps = 5000,
  cores = 1,
  keep_reps = FALSE,
  cp_p_threshold = 0.01,
  cp_window_size = 3000,
  cp_min_length = 10,
  cp_slope_check_window = 30,
  cp_slope_left = 1e-10,
  cp_slope_right = 1e-20,
  second_window_size = 50,
  second_p_threshold = 1e-10,
  min_signal = 2.5,
  detector = c("box", "changepoint"),
  box_threshold = NULL,
  box_calibration = NULL,
  box_level0 = NA_real_,
  box_null = 100,
  box_alpha = 0.05,
  box_null_n = NULL,
  box_method = c("order", "gumbel"),
  run_plots = TRUE,
  plot_device = "pdf",
  plot_p_threshold = 8,
  plot_y_cap = NULL
)

Arguments

Details

The box-scan threshold is calibrated on the data: unless box_calibration or box_threshold is supplied, the scan is repeated on box_null permuted phenotypes and the threshold is the order statistic of their maximum scores at level box_alpha (ras_box_calibrate), which controls the family-wise error rate without distributional assumptions. This multiplies the run time by about box_null; for large cohorts calibrate once with ras_box_calibrate and pass the result.

The pure-R, in-memory implementation that RAS 1.0.x shipped as ras() is available as ras_original; it defaults to the changepoint detector, and with detector = "changepoint" the two functions give the same scan profile to machine precision and the same detected positions.

Value

Invisibly, an object of class "ras": a list with scan (the ras_scan result), detection (for detector = "box" the ras_box_detect result, whose regions table has one row per interval; for detector = "changepoint" the ras_validate result, whose tau_hats are the detected positions), box_calibration (the calibration used, or NULL), chrom and save_dir. Use print() and plot() on it.

See Also

ras_scan, ras_detect, ras_validate, ras_box_detect, plot.ras; geno_to_rasbin and bed_to_rasbin to prepare the genotype file; ras_original for the pure-R implementation.

Examples

set.seed(3)
n_samp <- 120; n_snp <- 400
geno <- matrix(sample(0:2, n_samp * n_snp, replace = TRUE,
                      prob = c(0.6, 0.3, 0.1)), n_samp, n_snp)
causal <- 181:220                       # 40 causal SNPs, together explaining
pheno  <- 2 * as.numeric(scale(rowSums(geno[, causal]))) + rnorm(n_samp)  # ~80% of the variance
cov_df <- data.frame(age = rnorm(n_samp), sex = rbinom(n_samp, 1, 0.5))

## default: box-scan detector, threshold calibrated on permuted phenotypes
## (20 permutations with the Gumbel fit keep the example short; the
## default is the exact order statistic on 100 permutations)
res <- ras(geno, pheno, cov_df, covariate_cols = c("age", "sex"),
           is_continuous = TRUE, num_rep = 2, skip1 = 2, skip2 = 5,
           min_window_size = 2, max_window_size = 20,
           box_null = 20, box_method = "gumbel",
           save_dir = tempdir(), run_plots = FALSE)
res$detection$regions
res$box_calibration

## large data: convert once and pass the file; here with the changepoint
## detector of RAS 1.0.x (settings scaled down to the 200-point profile)
rb <- tempfile(fileext = ".rasbin")
geno_to_rasbin(geno, rb)
res_cp <- ras(rb, pheno, cov_df, covariate_cols = c("age", "sex"),
              is_continuous = TRUE, num_rep = 2, skip1 = 2, skip2 = 5,
              min_window_size = 2, max_window_size = 20,
              detector = "changepoint",
              cp_window_size = 100, cp_slope_check_window = 10,
              cp_slope_left = 1e-2, cp_slope_right = 1e-2,
              second_window_size = 20, second_p_threshold = 1e-2,
              save_dir = tempdir(), run_plots = FALSE)
res_cp$detection$tau_hats
unlink(c(rb, paste0(rb, ".meta.rds")))

Aliases Kept from the Development Versions

Description

In the development versions of this package the disk-backed C implementations carried a _fast suffix. Since RAS 1.1.0 they are the default implementations and carry the plain names, while the original pure-R, in-memory implementations carry an _original suffix. The _fast names are kept as thin aliases so that scripts written against the development versions keep working; new code should use the plain names.

Usage

ras_fast(...)

ras_scan_fast(...)

ras_detect_fast(...)

compute_gwas_weights_fast(...)

screen_forward_max_region_fast(...)

ras_scan_external_fast(...)

Arguments

Value

Whatever the aliased function returns.

See Also

ras, ras_scan, ras_detect, compute_gwas_weights, screen_forward_max_region, ras_scan_external.

Examples

identical(formals(ras_fast), formals(ras))   # FALSE only because of `...`

Calibrate the Box-Scan Threshold on Null Profiles

Description

Turns a set of null scan profiles into the threshold that ras_box_detect needs. Null profiles are obtained by running the same scan (ras_scan or ras_scan_external, with the same settings) on the same genotypes with the phenotype permuted, or with a phenotype simulated under the null.

Usage

ras_box_calibrate(
  null_profiles,
  alpha = 0.05,
  widths = .ras_box_default_widths(),
  gamma = 0.5,
  length_ratio = 1,
  calib = NULL,
  edge = c("open", "strict"),
  method = c("order", "gumbel")
)

Arguments

Details

level0 is the median of the null profiles' medians. Both thresholds are always computed and returned (threshold_order, threshold_gumbel); threshold is the one selected by method.

Value

An object of class "ras_box_calibration": a list with level0, threshold, threshold_order, threshold_gumbel, threshold_scaled, method, length_ratio, maxT (the maximum score of every null profile), calib, alpha, gumbel (fitted location and scale), null_exceed_scaled, fwer_holdout (NA when no profile was held out), widths, gamma, edge.

See Also

ras_box_detect, which consumes the result.

Examples

## 60 short synthetic null profiles (a real calibration uses several
## hundred profiles from the actual scan on permuted phenotypes)
set.seed(2)
nulls <- lapply(1:60, function(k) 0.7 + abs(rnorm(200, sd = 0.3)))
cal <- ras_box_calibrate(nulls, alpha = 0.05, calib = 1:40)
cal$threshold
cal$fwer_holdout           # error rate on the 20 held-out profiles

y <- nulls[[41]]; y[91:110] <- y[91:110] + 6
ras_box_detect(seq_along(y), y, calibration = cal)$regions

Detect Elevated Regions with the Box Scan

Description

The box-scan region detector, the alternative to the changepoint detector (ras_detect + ras_validate) selected with detector = "box" in ras. It reports intervals [\tau_L, \tau_R] of the scan profile whose mean stands above both adjacent flanks and above the profile background, which lets it delimit broad plateau-shaped association regions as well as sharp peaks.

Usage

ras_box_detect(
  x,
  y,
  threshold = NULL,
  calibration = NULL,
  scaled = FALSE,
  level0 = NULL,
  widths = NULL,
  gamma = NULL,
  edge = NULL,
  floor = 0.2,
  max_regions = Inf,
  keep_stat = FALSE
)

Arguments

Details

Windows are scored with ras_box_stat and taken greedily by score; a window is dropped when it touches an already accepted one expanded by 10 + max(5, width/2) grid points, which removes side lobes. Each accepted window becomes one region; its anchor is the highest point of the profile inside it.

The threshold is a property of the null distribution of the maximum score along a profile and therefore depends on the profile's length and on the scan settings. Obtain it with ras_box_calibrate from null profiles produced by the same scan on permuted phenotypes; the order statistic it returns controls the family-wise error rate at the requested level without any distributional assumption.

Value

A list that plot.ras and print.ras accept as a detection element:

detector

"box".

regions

Data frame with one row per region, ordered by decreasing score: tau_L, tau_R (grid indices), pos_L, pos_R (the same boundaries in x units), width, T_box, contrast, height, anchor_idx, anchor_pos, anchor_y. NULL when nothing reaches the threshold.

tau_hats

Numeric vector. Anchor positions in x units (the same role as in ras_validate's result).

all.changepoints, all.p.values

Anchor positions and their T_box scores, so that the candidate overlay of plot.ras works.

left.slopes, right.slopes

NULL; the box scan estimates no slopes.

threshold, level0, level, level_factor

The threshold used, the null level, the median of y and the resulting level factor.

stat

The ras_box_stat table when keep_stat = TRUE, otherwise NULL.

See Also

ras_box_calibrate for the threshold, ras_box_stat for the statistic, ras_detect for the changepoint detector, ras which calls this function when detector = "box".

Examples

set.seed(1)
x <- seq(1, 6000, by = 10)                 # 600 grid points
y <- 0.7 + abs(rnorm(600, sd = 0.3))       # null-like background
y[301:342] <- y[301:342] + 8               # a 42-point plateau
det <- ras_box_detect(x, y, threshold = 2)
det$regions[, c("tau_L", "tau_R", "pos_L", "pos_R", "T_box")]

## a window inside the plateau scores about 0: only its edges are reported
S <- ras_box_stat(y)
max(S$T[S$w == 20 & S$i >= 305 & S$i + 19 <= 338])

Box-Scan Statistic for Every Window Start and Width

Description

Computes the box-scan score T(i, w) of every window of every width in widths along a RAS scan profile. This is the building block of ras_box_detect; most users will call that function instead.

Usage

ras_box_stat(
  y,
  widths = .ras_box_default_widths(),
  level0 = NA_real_,
  gamma = 0.5,
  edge = c("open", "strict")
)

Arguments

Details

For a window of w grid points starting at i, with flanks of f = \max(5, w/2) points on each side,

contrast = mean(window) - max(mean(left flank), mean(right flank))

height = mean(window) - median(y)

T = min(contrast, height) / max(1, median(y)/level0)^{gamma}.

All quantities are in the raw units of the profile; no noise scale is estimated. A window inside a plateau, on the shoulder of a peak, or spanning a cluster of neighbouring peaks scores about 0 because at least one flank is as high as the window.

Value

A data frame with one row per (start, width) pair and columns i (window start, grid index), w (width), contrast, height and T.

See Also

ras_box_detect, ras_box_calibrate.

Examples

y <- rep(1, 400); y[201:230] <- 8
S <- ras_box_stat(y)
S[which.max(S$T), ]          # the 30-point window starting at 201

First-Pass Changepoint Detection via Sliding Window

Description

Slides a window of window_size grid points across the scan profile. In each window a single-breakpoint segmented regression of y on x is fitted (Muggeo's algorithm with bootstrap restarts) and the Davies test for the existence of a breakpoint is applied. A breakpoint is kept as a candidate when the Davies p-value is below p.values.threshold, the fitted slope rises on its left and falls on its right, and one-tailed slope tests in a local window of half-width slope_check_window_size confirm both directions. Accepted candidates are snapped to the nearest local maximum of y. The result is passed to ras_validate for the second pass.

Usage

ras_detect(
  x,
  y,
  p.values.threshold = 0.01,
  min.length = 10,
  skip = 1,
  window_size = 3000,
  slope_check_window_size = 30,
  slope.p.values.threshold.left = 1e-10,
  slope.p.values.threshold.right = 1e-20
)

Arguments

Details

The segmented fit and the Davies test are a port to compiled code of segmented and davies.test. The bootstrap restarts draw from R's random-number stream, so set.seed() makes a run reproducible, while results are statistically rather than bit-for-bit equal to the pure-R implementation ras_detect_original, which calls the segmented package.

Value

A list with tau_hats (accepted candidate positions, as indices into x), p.values, slope.left, slope.right and slope.angle for those candidates, all.changepoints and all.p.values for every window examined, and previous_tau_hats (a copy of tau_hats used by ras_validate).

See Also

ras_validate for the second pass; ras_box_detect for the alternative detector; ras, which runs both passes; ras_detect_original.

Examples

set.seed(42)
x <- 1:100
y <- c(seq(0, 8, length.out = 50), seq(8, 1, length.out = 50)) +
     rnorm(100, sd = 0.5)
cp <- ras_detect(x, y, window_size = 50, slope_check_window_size = 10,
                 slope.p.values.threshold.left  = 1e-3,
                 slope.p.values.threshold.right = 1e-3)
cp$tau_hats

First-Pass Changepoint Detection (Pure-R)

Description

This is the pure-R, in-memory implementation that RAS 1.0.x shipped under the plain name. Since 1.1.0 the plain name (ras_detect) is the compiled, disk-backed implementation; this function is kept for reference and gives the same results. Scans a -\log_{10}(p)-value sequence using a sliding window to detect positions where the slope changes significantly from positive to negative.

Usage

ras_detect_original(
  x,
  y,
  p.values.threshold = 0.01,
  min.length = 10,
  skip = 1,
  window_size = 3000,
  slope_check_window_size = 30,
  slope.p.values.threshold = 1e-08,
  slope.p.values.threshold.left = 1e-10,
  slope.p.values.threshold.right = 1e-20
)

Arguments

Details

At each window start position the function calls get_break_points on the window to estimate and test a single breakpoint. A candidate is retained when three conditions are all met:

1. The Davies test p-value is below p.values.threshold.

2. The left-side slope (estimated by segmented regression) is positive.

3. The right-side slope is negative.

Candidates passing these three filters are then subjected to one-tailed slope tests (slope_test) using centred sub-sequences of half-width slope_check_window_size. Only candidates that also pass both slope p-value thresholds are recorded.

After the sliding-window loop, each accepted changepoint is refined to the nearest local peak in y via get_local_maximum.

The function records all candidates examined (all.changepoints, all.p.values) in addition to the accepted ones, so that downstream functions can display the full candidate landscape.

Value

A named list with eight elements:

tau_hats

Integer vector. Accepted changepoint indices (refined to local peaks).

p.values

Numeric vector. Davies test p-values at accepted changepoints.

slope.left

Numeric vector. Left-side slopes at accepted changepoints.

slope.right

Numeric vector. Right-side slopes at accepted changepoints.

all.changepoints

Integer vector. All candidate positions examined, including those that failed the slope tests.

all.p.values

Numeric vector. Davies p-values for all examined candidates; set to 1 for candidates that failed the slope direction or slope significance filters.

slope.angle

Numeric vector. Interior angle (degrees) at each accepted changepoint, computed from the left and right slope estimates.

previous_tau_hats

Integer vector. Copy of tau_hats before any downstream modification; used by ras_validate.

See Also

ras_validate for the second-pass validation step. get_break_points for the per-window segmented regression. slope_test for the one-tailed slope verification. get_local_maximum for the peak-refinement step. ras_original for the recommended end-to-end entry point.

Examples

set.seed(42)
x <- 1:100
y <- c(seq(0, 8, length.out = 50),
       seq(8, 1, length.out = 50)) + rnorm(100, sd = 0.5)
result <- ras_detect_original(
  x, y,
  window_size             = 50,
  skip                    = 2,
  slope_check_window_size = 10,
  slope.p.values.threshold.left  = 1e-3,
  slope.p.values.threshold.right = 1e-3
)
cat("Detected changepoints:", result$tau_hats, "\n")

Harmonise External GWAS Summary Statistics into RAS Weights

Description

The entry point for summary-informed RAS. Takes a raw external summary-statistics file and a local genotype map, chooses the naming/build route from the data rather than from an assumption, and returns a verified weight vector together with a reproducible QC table.

Usage

ras_harmonize_sumstats(
  sumstats,
  map,
  dict = NULL,
  chain = NULL,
  chain_back = NULL,
  cup = NULL,
  converter = .liftover_converter_ucsc,
  check_reverse = TRUE,
  assume_same_build = FALSE,
  cols = list(),
  ...
)

Arguments

Details

The driver never guesses about the genome build. It inspects both sides' ID conventions and then either matches directly, translates the summary statistics through a dictionary, lifts their coordinates over, or refuses – listing exactly what would unblock it. Every path ends at ras_weights_from_sumstats with build_check = TRUE unless the caller has explicitly asserted the builds agree or a liftOver round trip has already established them.

Value

A list with

See Also

ras_scan_external to scan with these weights; ras_sumstats_report to inspect an alignment without committing to it.

LiftOver External Summary Statistics onto the Map's Build (Route B2)

Description

Converts the summary statistics' coordinates to another genome build and verifies every conversion by round-tripping it (from -> to -> from); anything that cannot return to its own chromosome and position is dropped. This is the LAST-RESORT route. liftOver has real error rates (about 1.7% of positions discordant GRCh37 to GRCh38, 3.1% the other way), so prefer a dictionary lookup (ras_read_variant_dictionary) when one is available: a lookup can only fail to find a variant, it cannot find the wrong one.

Usage

ras_liftover_sumstats(
  ss,
  chain,
  chain_back = NULL,
  cols = list(),
  converter = .liftover_converter_ucsc,
  check_reverse = TRUE,
  cup = NULL,
  ...
)

Arguments

Details

Requires the UCSC liftOver binary on PATH unless a replacement converter is supplied.

Value

A list with sumstats (only successfully converted rows, with the pre-conversion position kept in BP_source) and report (a per-status count table).

See Also

ras_harmonize_sumstats

Rename a Local Genotype Map's Variants to rsIDs (Harmonisation Route B0)

Description

Gives the LOCAL map rsIDs using a dictionary on the SAME genome build, so that afterwards both sides are keyed on rsID, the build gate becomes evaluable, and nothing about positions has to be trusted. Preferred when the map is ours – inside All of Us the dictionary is the Variant Annotation Table.

Usage

ras_map_ids_to_rsid(map, dict, update_name_file = NULL)

Arguments

Details

Matching is chr:pos plus the unordered allele pair, so multi-allelic sites cannot be confused. Variants with no dictionary entry keep their original ID (they simply will not match the summary statistics afterwards).

Value

A list with map (the renamed map, carrying the original IDs in SNP_original), report (counts), and update_name (the rename pairs as a data frame).

See Also

ras_harmonize_sumstats

Estimate Memory Requirements and Check System Readiness

Description

Detects the current machine's RAM and CPU configuration, computes per-stage peak memory estimates for the RAS pipeline, and issues a go / no-go verdict.

Usage

ras_memory(
  n_total,
  n_train,
  n_holdout,
  n_snps,
  bytes_per_element = 8,
  abort = FALSE
)

Arguments

Details

Peak memory is estimated for three pipeline stages:

Stage 1: compute_gwas_weights_original

Holds the full geno matrix plus a small N \times 4 coefficient matrix. Peak \approx geno_mb + coefmat_mb.

Stage 2: compute_pgs_matrix

Worst-case holds two copies of geno (R copy-on-modify triggered by geno[is.na(geno)] <- 0) plus the output pgs.mat. Peak \approx 2 * geno_mb + pgsmat_mb.

Stage 3: screen_forward_max_region_original

Holds geno and pgs.mat simultaneously. Peak \approx geno_mb + pgsmat_mb.

The overall estimated peak is the maximum across the three stages.

Available RAM is queried via wmic on Windows and /proc/meminfo on Linux; the function degrades gracefully (prints a caution message and returns can_proceed = TRUE) if the query fails.

The function prints a formatted report to the console and invisibly returns the numeric estimates for programmatic use.

Value

Invisibly returns a named list with two elements:

memory

Named numeric vector with elements geno_mb, pgsmat_mb, coefmat_mb, stage1_peak_mb, stage2_peak_mb, stage3_peak_mb, overall_peak_mb, available_mb, total_mb.

system

Named list with elements cores_physical, cores_logical, r_version, platform, and can_proceed (logical: TRUE if RAM is sufficient or cannot be measured).

See Also

ras_scan_original, ras_original for the functions whose memory use is being estimated. release_memory to reclaim heap memory after each repetition.

Examples

## Check feasibility for 5,000 samples and 500,000 SNPs
ras_memory(
  n_total   = 5000,
  n_train   = 2500,
  n_holdout = 2500,
  n_snps    = 500000
)

## Abort if memory is insufficient (wrapped in try() so the example runs)
try(ras_memory(
  n_total   = 100000,
  n_train   = 50000,
  n_holdout = 50000,
  n_snps    = 1000000,
  abort     = TRUE
))

RAS Pipeline (Pure-R, In-Memory)

Description

This is the pure-R, in-memory implementation that RAS 1.0.x shipped under the plain name. Since 1.1.0 the plain name (ras) is the compiled, disk-backed implementation; this function is kept for reference and gives the same results. Detect significant association regions in GWAS data. Returns an object of class "ras".

Usage

ras_original(
  geno,
  phenotype,
  covariates,
  covariate_cols,
  is_continuous,
  num_rep = 5,
  skip1 = 10,
  skip2 = 20,
  chrom = 1,
  save_dir = file.path(tempdir(), "RAS"),
  min_window_size = 5,
  max_window_size = 100,
  scan_test = c("glm", "score"),
  cp_p_threshold = 0.01,
  cp_window_size = 3000,
  cp_min_length = 10,
  cp_slope_check_window = 30,
  cp_slope_left = 1e-10,
  cp_slope_right = 1e-20,
  second_window_size = 50,
  second_p_threshold = 1e-10,
  min_signal = 2.5,
  detector = c("changepoint", "box"),
  box_threshold = NULL,
  box_calibration = NULL,
  box_level0 = NA_real_,
  run_plots = TRUE,
  plot_device = "pdf",
  plot_p_threshold = 8,
  plot_y_cap = NULL
)

Arguments

Details

For a stage-by-stage description of the pipeline see RAS.

Two region detectors are available for Stage 2. The default detector = "changepoint" reproduces the published method and the results of earlier versions of this package exactly. detector = "box" runs ras_box_detect instead; it needs a threshold, which is a property of the null distribution of the scan and is obtained with ras_box_calibrate from null profiles (the same scan on permuted phenotypes).

Value

Invisibly returns a named list with two elements:

scan

A list with elements:

x

Integer vector. SNP position index grid seq(1, N, by = skip1).

y

Numeric vector (same length as x). Averaged -\log_{10}(p)-value profile across all repetitions.

detection

For detector = "changepoint", the list returned by ras_validate with elements:

tau_hats

Integer vector. Validated changepoint positions (re-mapped to genomic coordinates via this.start and this.skip).

all.changepoints

Integer vector. All candidate positions examined during the first pass, re-mapped to genomic coordinates.

all.p.values

Numeric vector. -\log_{10}(p) Davies test values for every candidate in all.changepoints.

left.slopes

Numeric vector. Estimated left-side slopes at each validated changepoint.

right.slopes

Numeric vector. Estimated right-side slopes at each validated changepoint.

detector

"changepoint".

For detector = "box", the list returned by ras_box_detect: regions (one row per interval with pos_L, pos_R, T_box, ...), tau_hats (the anchor of each interval) and the threshold used.

Plot files (if run_plots = TRUE) are written to save_dir with names chr-<chrom>-cp-plot.<ext>, chr-<chrom>-cp-p-values-plot.<ext>, and chr-<chrom>-zoom.<ext>.

Simple (recommended)

## one call does everything
result <- ras_original(geno, phenotype, covariates,
              covariate_cols = c("age", "sex", paste0("pc", 1:10)),
              is_continuous  = TRUE,
              chrom = 1, save_dir = "results/")

print(result)             # detected changepoint positions
plot(result)              # full-chromosome scan profile
plot(result, zoom = TRUE) # zoomed view around each changepoint

For step-by-step control see ras_scan_original, ras_detect_original, ras_validate.

See Also

ras_scan_original for the scan-only step (useful when tuning changepoint parameters separately). ras_detect_original, ras_validate for the changepoint detection steps. plot.ras for the plotting step. ras_memory to check memory requirements before loading large genotype data.

Examples

set.seed(42)
n_samp <- 80; n_snp <- 60
geno <- matrix(
  sample(0:2, n_samp * n_snp, replace = TRUE, prob = c(0.6, 0.3, 0.1)),
  nrow = n_samp, ncol = n_snp)
pheno <- rnorm(n_samp)
cov_df <- data.frame(age = rnorm(n_samp), sex = rbinom(n_samp, 1, 0.5))
cov_df$age_squared <- cov_df$age^2
cov_df$age_sex     <- cov_df$age * cov_df$sex
for (i in 1:10) cov_df[[paste0("pc", i)]] <- rnorm(n_samp)

result <- ras_original(
  geno              = geno,
  phenotype         = pheno,
  covariates        = cov_df,
  covariate_cols    = c("age", "sex", "age_squared", "age_sex",
                        paste0("pc", 1:10)),
  is_continuous     = TRUE,
  num_rep           = 2,
  skip1             = 1,        # every SNP on the grid: a 60-point profile
  skip2             = 2,
  min_window_size   = 2,
  max_window_size   = 10,
  # detector settings scaled down to the 60-point toy profile (the defaults
  # are tuned for genome-scale grids with tens of thousands of points)
  cp_window_size        = 30,
  cp_min_length         = 5,
  cp_slope_check_window = 5,
  cp_slope_left         = 1e-3,
  cp_slope_right        = 1e-3,
  second_window_size    = 10,
  second_p_threshold    = 1e-3,
  chrom             = 1,
  save_dir          = tempdir(),
  run_plots         = TRUE,
  plot_device       = "pdf"
)

cat("Detected changepoints:", result$detection$tau_hats, "\n")
plot(result$scan$x, result$scan$y, type = "l",
     xlab = "SNP index", ylab = expression(-log[10](p)),
     main = "RAS scan profile")

## The box-scan region detector on the same scan profile. In the one-call
## form this is ras_original(..., detector = "box", box_threshold = 2); the threshold
## would normally come from ras_box_calibrate() on null profiles.
result_box <- result
result_box$detection <- ras_box_detect(result$scan$x, result$scan$y,
                                       threshold = 2)
print(result_box)

Read a Variant Dictionary for Summary-Statistics Harmonisation

Description

Reads a table mapping rsID to (chromosome, position, REF, ALT) on ONE genome build, for use by ras_harmonize_sumstats's dictionary route and by ras_map_ids_to_rsid. Sources include a dbSNP slice, the All of Us Variant Annotation Table, or any delimited file with those columns.

Usage

ras_read_variant_dictionary(src, cols = list(), build = NA_character_)

Arguments

Details

A plain VCF is accepted: lines beginning ## are skipped and #CHROM is used as the header. An ALT field carrying several comma-separated alleles is expanded to one row per alternate allele, so multi-allelic sites resolve by their alleles rather than by row order.

Value

A data frame with columns RSID, CHR, BP, REF, ALT, KEY (the unordered allele-pair key), with the build label attached as an attribute.

See Also

ras_harmonize_sumstats, ras_map_ids_to_rsid

RAS Stage 1: Averaged Regional Association Profile

Description

Computes the -\log_{10}(p) profile that the RAS detectors work on. The sample is split at random into a training half and a hold-out half num_rep times; in each repetition per-SNP regression weights are fitted on the training half (compute_gwas_weights) and a forward regional scan is run on the hold-out half (screen_forward_max_region); the profiles are averaged over repetitions. Genotypes are streamed from a .rasbin file in chunks of chunk_snps SNP columns, so peak memory does not grow with the chromosome. For a genotype matrix small enough to sit in memory the pure-R implementation ras_scan_original gives the same profile.

Usage

ras_scan(
  geno,
  phenotype,
  covariates,
  covariate_cols,
  is_continuous,
  num_rep = 5,
  skip1 = 10,
  skip2 = 20,
  chrom = 1,
  save_dir = file.path(tempdir(), "RAS"),
  min_window_size = 5,
  max_window_size = 100,
  scan_test = c("score", "glm"),
  chunk_snps = 5000,
  cores = 1,
  keep_reps = FALSE,
  rows = NULL
)

Arguments

Value

Invisibly, a list with x (the SNP index of each grid point, seq(1, N, by = skip1)), y (the averaged -\log_{10}(p) profile) and reps (a length(x) by num_rep matrix of per-repetition profiles when keep_reps = TRUE, otherwise NULL). The averaged profile is also written to save_dir as mean_p_values_chr<chrom>_reps1-<num_rep>.rds.

See Also

ras for the full pipeline; ras_detect and ras_box_detect for the detectors that consume the profile; ras_scan_external for a scan with external weights; ras_scan_original for the pure-R implementation.

RAS Scan with External Weights on All Samples

Description

Runs the Stage-1 forward scan with weights taken from an independent external GWAS instead of a within-sample training split. Because the weights are independent of the target cohort, no 50/50 split is drawn and no repetitions are averaged: every sample is scanned once. Genotypes are streamed from a .rasbin file in chunks. The in-memory implementation is ras_scan_external_original.

Usage

ras_scan_external(
  geno,
  phenotype,
  covariates,
  covariate_cols,
  weights,
  is_continuous = TRUE,
  skip1 = 10,
  skip2 = 20,
  min_window_size = 5,
  max_window_size = 100,
  chunk_snps = 5000,
  chrom = 1,
  save_dir = NULL
)

Arguments

Details

The thresholds of the changepoint detector were set on averaged profiles, whose noise floor is lower than that of a single-pass profile; calibrate before treating a detection from this scan as a finding (for the box-scan detector, ras_box_calibrate on permuted phenotypes).

Value

A list with x (the SNP index of each grid point) and y (the -\log_{10}(p) profile).

See Also

ras_harmonize_sumstats to build weights; ras_scan for the split-based scan; ras_box_detect and ras_detect to detect regions on the profile.

Summary-Informed RAS Scan on All Samples (In-Memory)

Description

The in-memory counterpart of ras_scan_external, for cohorts small enough that a dense n \times m genotype matrix and its PGS matrix fit in RAM. At biobank scale use ras_scan_external instead: at 453,698 samples by 332,690 variants the dense double alone is over a petabyte.

Usage

ras_scan_external_original(
  geno,
  phenotype,
  covariates,
  covariate_cols,
  weights,
  is_continuous = TRUE,
  skip1 = 10,
  skip2 = 20,
  min_window_size = 5,
  max_window_size = 100,
  scan_test = c("glm", "score"),
  chrom = 1,
  save_dir = NULL
)

Arguments

Value

A list with x (grid point variant indices) and y (the profile).

See Also

ras_scan_external, ras_harmonize_sumstats

RAS Stage 1 Scan (Pure-R, In-Memory)

Description

This is the pure-R, in-memory implementation that RAS 1.0.x shipped under the plain name. Since 1.1.0 the plain name (ras_scan) is the compiled, disk-backed implementation; this function is kept for reference and gives the same results. Computes the averaged -\log_{10}(p)-value profile over num_rep independent 50/50 train/hold-out splits.

Usage

ras_scan_original(
  geno,
  phenotype,
  covariates,
  covariate_cols,
  is_continuous,
  num_rep = 5,
  skip1 = 10,
  skip2 = 20,
  chrom = 1,
  save_dir = file.path(tempdir(), "RAS"),
  min_window_size = 5,
  max_window_size = 100,
  scan_test = c("glm", "score")
)

Arguments

Details

For each of the num_rep repetitions the function performs three steps:

1. GWAS weights. A fresh 50/50 random split of all n samples is drawn. compute_gwas_weights_original fits a per-SNP regression on the training half and returns an N \times 4 coefficient matrix; the first column (effect-size estimates) is used as PGS weights.

2. PGS contribution matrix. compute_pgs_matrix multiplies each hold-out individual's dosage vector by the per-SNP weights, producing an n_{\text{holdout}} \times N matrix of weighted contributions.

3. Forward scan. screen_forward_max_region_original slides an expanding window across the genome. At each position it accumulates weighted dosages, regresses them against the hold-out phenotype, and records the minimum p-value over window sizes in [min_window_size, max_window_size]. The result is a vector of -\log_{10}(p) values on the grid seq(1, N, by = skip1).

The -\log_{10}(p) vectors from all repetitions are summed and divided by num_rep to form the final profile.

For continuous traits (is_continuous = TRUE), covariates are residualised from the training phenotype using only training individuals before the GWAS step, so no hold-out information leaks into the effect-size estimates. The residualisation is repeated from the original phenotype vector each repetition. For binary traits (is_continuous = FALSE) the raw phenotype is passed directly to the GWAS step, and the scan step fits a logistic model via glm(family = binomial()).

After each repetition the large intermediate matrices (coef.mat, pgs.mat) are removed from the R session and release_memory is called to return free heap pages to the OS. On Linux/glibc this calls malloc_trim(0) and can substantially reduce RSS between repetitions.

Value

Invisibly returns a named list with two elements:

x

Integer vector of length \lceil N / \texttt{skip1} \rceil. The SNP position index grid seq(1, N, by = skip1).

y

Numeric vector, same length as x. Averaged -\log_{10}(p)-value profile across all num_rep repetitions.

The following .rds files are written to save_dir:

chr-<chrom>_coef_mat-<rep>.rds

The N \times 4 GWAS coefficient matrix from repetition rep.

mean_p_values_chr<chrom>_reps1-<num_rep>.rds

The averaged -\log_{10}(p) vector y.

See Also

ras_original for a single-call wrapper that also runs changepoint detection and plotting. compute_gwas_weights_original, compute_pgs_matrix, screen_forward_max_region_original for the constituent steps. ras_detect_original, ras_validate for downstream changepoint analysis. ras_memory to check memory requirements before loading large genotype data. release_memory for OS-level heap reclamation.

Examples

set.seed(42)
n_samp <- 80; n_snp <- 60
geno <- matrix(
  sample(0:2, n_samp * n_snp, replace = TRUE, prob = c(0.6, 0.3, 0.1)),
  nrow = n_samp, ncol = n_snp)
pheno <- rnorm(n_samp)
cov_df <- data.frame(age = rnorm(n_samp), sex = rbinom(n_samp, 1, 0.5))
cov_df$age_squared <- cov_df$age^2
cov_df$age_sex     <- cov_df$age * cov_df$sex
for (i in 1:10) cov_df[[paste0("pc", i)]] <- rnorm(n_samp)

scan <- ras_scan_original(
  geno            = geno,
  phenotype       = pheno,
  covariates      = cov_df,
  covariate_cols  = c("age", "sex", "age_squared", "age_sex",
                      paste0("pc", 1:10)),
  is_continuous   = TRUE,
  num_rep         = 2,
  skip1           = 5,
  skip2           = 5,
  min_window_size = 2,
  max_window_size = 10,
  chrom           = 1,
  save_dir        = tempdir()
)
plot(scan$x, scan$y, type = "l",
     xlab = "SNP index", ylab = expression(-log[10](p)),
     main = "RAS scan profile")

Inspect an External-Summary Alignment Without Committing to It

Description

A non-fatal report on how an external summary-statistics file lines up with a local genotype map: which ID convention each side uses, the suggested route, how many variants overlap by ID, what the build gate would say, and how the overlap is distributed along the chromosome. Run this before ras_harmonize_sumstats to see what the gates will do without triggering them.

Usage

ras_sumstats_report(ss, map, cols = list(), block_size = 500L)

Arguments

Details

The coverage line matters more here than in a PRS setting: RAS's statistic is regional, so weights that are missing in one stretch of the chromosome remove the scan's ability to detect there, rather than merely weakening it overall.

Value

Invisibly, a list with plan, n_id, build and blocks. Called for the report it prints.

See Also

ras_harmonize_sumstats

Second-Pass Changepoint Validation

Description

Validates candidate changepoints from ras_detect_original by re-running local Davies tests in windows around each candidate.

Usage

ras_validate(
  this.result,
  x,
  y,
  this.start = 1,
  this.skip = 30,
  second_window_size = 50,
  p.value.threshold = 1e-10,
  min_signal = 2.5
)

Arguments

Details

For each candidate \hat{\tau} in this.result$tau_hats the function fits two local linear models and applies the Davies test to each: one on the window [\hat{\tau},\; \min(\hat{\tau} + \texttt{second\_window\_size},\, n)] and one on [\max(\hat{\tau} - \texttt{second\_window\_size},\, 1),\; \hat{\tau}]. If either Davies p-value is below p.value.threshold, the candidate is accepted; otherwise its all.p.values entry is set to zero to suppress it in downstream plots.

Windows with fewer than four observations are not tested (Davies test requires at least four points) and their p-value is set to 1.0.

After filtering, any remaining candidate with y[tau_hat] <= min_signal is removed, as such positions are below the minimum signal threshold.

Accepted changepoint indices and all candidate indices are re-mapped from the scan grid to genomic coordinates via \texttt{this.start} + (\text{index} - 1) \times \texttt{this.skip}.

Value

A named list with five elements:

tau_hats

Integer vector. Validated changepoint positions, re-mapped to genomic coordinates.

all.changepoints

Integer vector. All candidate positions examined, re-mapped to genomic coordinates.

all.p.values

Numeric vector. -\log_{10}(p) Davies values for each candidate; set to 0 for rejected candidates.

left.slopes

Numeric vector. Left-side slopes at accepted changepoints (carried over from first-pass output).

right.slopes

Numeric vector. Right-side slopes at accepted changepoints.

See Also

ras_detect_original for the first-pass detection step whose output this function takes as input. plot.ras for visualising the validated changepoints. ras_original for the recommended end-to-end entry point.

Examples

set.seed(42)
x <- 1:100
y <- c(seq(0, 8, length.out = 50),
       seq(8, 1, length.out = 50)) + rnorm(100, sd = 0.5)

cp_result <- ras_detect_original(
  x, y,
  window_size             = 50,
  skip                    = 2,
  slope_check_window_size = 10,
  slope.p.values.threshold.left  = 1e-3,
  slope.p.values.threshold.right = 1e-3
)

final <- ras_validate(
  cp_result, x = x, y = y,
  this.skip          = 1,
  second_window_size = 20,
  p.value.threshold  = 1e-3
)
cat("Validated changepoints:", final$tau_hats, "\n")

Align External GWAS Summary Statistics to a Local Genotype Map

Description

Turns an external GWAS's effect estimates into a weight vector aligned to a local genotype matrix, one weight per genotype column, in order. This is the matcher underneath ras_harmonize_sumstats; call that instead unless you have already settled the naming and build questions yourself.

Usage

ras_weights_from_sumstats(
  sumstats,
  map,
  cols = list(),
  drop_ambiguous = TRUE,
  on_or = TRUE,
  match_by = c("id", "pos"),
  match_min_prop = 0.2,
  build_check = TRUE,
  target_af = NULL,
  af_tol = 0.2,
  block_size = 500L,
  pos_fallback = NULL
)

Arguments

Details

The alignment rule is deliberately conservative:

• w_j = +\beta_j when the summary statistics' effect allele is the allele the local dosage counts (map$A1);

• w_j = -\beta_j when it is the other allele (map$A2);

• w_j = 0 when the variant is absent, the alleles cannot be reconciled, or (with drop_ambiguous = TRUE) the variant is strand-ambiguous.

Strand-ambiguous (A/T, C/G) variants are removed rather than rescued. Rescuing them requires a frequency rule, which is a modelling decision this layer does not make.

Value

A list with

See Also

ras_harmonize_sumstats for the driver that chooses the naming/build route first; ras_scan_external to scan with the resulting weights.

Read the .rasbin File Header

Description

Read the .rasbin File Header

Usage

rasbin_header(path)

Arguments

Value

Named numeric vector with elements n_samples, n_snps, dtype.

Read a Column Range from a .rasbin File

Description

Read a Column Range from a .rasbin File

Usage

rasbin_read_chunk(path, col_start, col_end)

Arguments

Value

An n_samples x (col_end - col_start + 1) numeric matrix.

Release Memory Back to the OS

Description

Runs full garbage collection and, on Linux with the GNU C library (glibc), calls malloc_trim(0) to return free heap pages to the operating system. Useful after large temporary matrices are removed in memory-heavy RAS pipeline stages.

Usage

release_memory(verbose = TRUE)

Arguments

Details

malloc_trim() is a glibc extension. On every other platform, including Linux systems with another C library such as musl (Alpine Linux), Windows and macOS, the call is not compiled in and NA is returned silently; no error is raised.

Value

Invisibly returns the malloc_trim(0) result: 1 if heap pages were returned to the OS, 0 if nothing was returned, NA_integer_ where malloc_trim() is unavailable (all platforms other than glibc Linux).

Forward Scan of the RAS Profile

Description

The Stage-1 forward scan on the hold-out half of the sample. At every grid position (every skip1 SNPs) windows of min_window_size, min_window_size + skip2, ... up to max_window_size SNPs are formed, each window's weighted dosage score is tested against the phenotype, and the smallest p-value is recorded as the profile value. SNP columns are streamed from the .rasbin file and the weighted scores are accumulated on the fly, so neither the genotype matrix nor a score matrix is held in memory. The pure-R, in-memory implementation is screen_forward_max_region_original together with compute_pgs_matrix.

Usage

screen_forward_max_region(
  rasbin_path,
  weights,
  this.leftout,
  this.df,
  is_continuous,
  covariate_formula,
  skip1 = 10,
  skip2 = 20,
  min_window_size = 5,
  max_window_size = 100,
  chunk_snps = 5000
)

Arguments

Value

Numeric vector with one -\log_{10}(p) value per grid position, ceiling(N / skip1) values in total.

See Also

ras_scan, which averages this scan over repetitions; ras_scan_external; screen_forward_max_region_original.

Forward Scan of the RAS Profile (Pure-R, In-Memory)

Description

This is the pure-R, in-memory implementation that RAS 1.0.x shipped under the plain name. Since 1.1.0 the plain name (screen_forward_max_region) is the compiled, disk-backed implementation; this function is kept for reference and gives the same results. Slides an expanding window across the genome, accumulating weighted genotype contributions and testing association with the hold-out phenotype at each position.

Usage

screen_forward_max_region_original(
  geno,
  pgs.mat,
  this.df,
  num_signals,
  start.point = 1,
  save.directory = tempdir(),
  this.chrome = 1,
  min_window_size = 5,
  max_window_size = 100,
  isSimulation = TRUE,
  this.repetition = 1,
  screening_round = 1,
  isPlot = FALSE,
  skip1 = 100,
  skip2 = 5,
  is_continuous,
  covariate_formula = NULL,
  scan_test = c("glm", "score"),
  signal.starts = NULL,
  signal.window.size = NULL
)

Arguments

Details

At each grid position j \in \{1, 1 + \texttt{skip1}, \ldots, N\} the function builds a polygenic score (PGS) by accumulating columns of pgs.mat from a growing symmetric window [j - w + 1,\; j + w - 1] for w \in \{0, \texttt{min\_window\_size}, \ldots, \texttt{max\_window\_size}\}. The accumulation is incremental: only the newly added columns are summed at each step, avoiding redundant computation.

At window size w = 0, only column j of pgs.mat is used (the single-SNP PGS). For each window size the PGS is tested against phenotype2 via lm (continuous) or glm(family = binomial()) (binary). The minimum p-value over all window sizes is stored for position j.

The final return value is -\log_{10} of these per-position minimum p-values.

Value

Numeric vector of length \lceil N / \texttt{skip1} \rceil. Each element is the -\log_{10}(p)-value at the corresponding grid position, where the p-value is the minimum over all tested window sizes.

See Also

compute_pgs_matrix for the step that produces pgs.mat. ras_scan_original for the recommended high-level entry point that calls this function automatically for each repetition.

Examples

set.seed(3)
n_samp <- 80; n_snp <- 50
geno <- matrix(sample(0:2, n_samp * n_snp, replace = TRUE,
                      prob = c(0.6, 0.3, 0.1)),
               nrow = n_samp, ncol = n_snp)
weights  <- rnorm(n_snp)
leftout  <- 31:80
pgs_mat  <- compute_pgs_matrix(geno, leftout, weights)

scan_df           <- data.frame(matrix(rnorm(50 * 12), 50, 12))
colnames(scan_df) <- c("age", "sex", "age_squared", "age_sex",
                       paste0("pc", 1:8))
scan_df$phenotype2 <- rnorm(50)

p_vec <- screen_forward_max_region_original(
  geno              = geno,
  pgs.mat           = pgs_mat,
  this.df           = scan_df,
  num_signals       = -1,
  is_continuous     = TRUE,
  covariate_formula = paste(c("age", "sex"), collapse = " + "),
  skip1             = 5,
  skip2             = 5,
  min_window_size   = 2,
  max_window_size   = 8,
  isPlot            = FALSE
)
plot(seq(1, n_snp, by = 5), p_vec, type = "l",
     xlab = "SNP index", ylab = expression(-log[10](p)))

One-Tailed Slope Test Through the Origin

Description

Fits a no-intercept linear model y ~ x - 1 and performs a one-tailed t-test on the slope coefficient.

Usage

slope_test(x, y, lower.tail)

Arguments

Details

The model lm(y ~ x - 1) forces the regression line through the origin, which is appropriate after centring both x and y at the candidate changepoint. Under this parameterisation a positive left slope and a negative right slope correspond to the peak shape expected at a true association region.

The function is called by ras_detect_original on both sides of each Davies-significant candidate, using asymmetric thresholds (cp_slope_left = 1e-10, cp_slope_right = 1e-20 by default) to require a steeper descending edge than ascending edge.

Value

Numeric scalar. The one-tailed p-value for the slope coefficient.

See Also

get_break_points for the Davies test step that precedes this slope check. ras_detect_original for the full detection workflow.

Examples

# Left-side slope: x goes from negative to 0, y should be rising (positive slope)
x_left <- -10:0
y_left <- x_left * 0.8 + rnorm(11, sd = 0.3)
slope_test(x_left, y_left, lower.tail = FALSE)  # expect small p (slope > 0)

# Right-side slope: x goes from 0 to positive, y should be falling (negative slope)
x_right <- 0:10
y_right <- x_right * (-0.8) + rnorm(11, sd = 0.3)
slope_test(x_right, y_right, lower.tail = TRUE)  # expect small p (slope < 0)