Estimates variant-importance scores (PI) based on functional annotations using supervised learning. The function trains an ensemble of prediction models (LASSO or GLM) on case and control annotations and applies them to testing data to estimate the relative importance of each variant for the phenotype of interest.

get_PI(
  training_caseAnnotation,
  training_controlAnnotation,
  training_control_need_N,
  model_need_N,
  modelType,
  testing_annotations
)

Arguments

training_caseAnnotation

Numeric matrix (N_case x M) of functional annotation scores for case variants (known causal variants from training data). Each row represents one variant, each column represents one annotation feature.

training_controlAnnotation

Numeric matrix (N_control x M) of functional annotation scores for control variants (known non-causal variants from training data). Must have the same number of columns as training_caseAnnotation.

training_control_need_N

Integer, number of control variants to sample for each model in the ensemble. Should be less than or equal to nrow(training_controlAnnotation). Typically chosen to balance the training dataset (e.g., equal to or 2-3 times nrow(training_caseAnnotation)).

model_need_N

Integer, number of models to fit in the ensemble. Higher values increase prediction stability at the cost of computation time. Typical values: 5-20.

modelType

Character string specifying the model type: "LASSO" or "GLM".

  • "LASSO": L1-regularized logistic regression with automatic feature selection via cross-validation. Recommended when M is large or features are correlated.

  • "GLM": Standard logistic regression using all features. Recommended when M is small relative to N_case and features are relatively independent.

testing_annotations

Numeric matrix or data frame (N_test x M) of functional annotation scores for the variants being analyzed. Must have the same columns (annotation features) as the training data.

Value

A numeric vector of length N_test containing estimated variant-importance scores (PI values) for each variant in testing_annotations. Values are in the range (0, 1), where higher values indicate higher estimated importance based on the annotation profile.

Details

Methodology:

This function implements annotation-based variant-importance estimation for the GLOW framework. The key idea is to leverage functional annotations (e.g., CADD, PolyPhen, conservation scores) to estimate which variants are more likely to be important, allowing GLOW to upweight likely important variants in the association test.

Algorithm:

  1. Ensemble Model Training:

    • For i = 1 to model_need_N:

      1. Sample training_control_need_N controls randomly from training_controlAnnotation

      2. Combine all cases with sampled controls

      3. Fit a logistic regression model (LASSO or GLM) predicting case/control status (1 for cases, 0 for controls) from annotations

    • This creates model_need_N independent models

  2. Prediction:

    • Apply each of the model_need_N models to testing_annotations

    • Average predictions across models to obtain final PI estimates

Why Ensemble Approach?

Using multiple models with different control samples:

  • Reduces dependence on specific control sample

  • Improves prediction stability and robustness

  • Accounts for variability in control annotations

  • Provides implicit uncertainty quantification

LASSO vs GLM:

Use LASSO when:

  • Number of annotations M is large (e.g., M > 20)

  • Annotations are correlated

  • You want automatic feature selection

  • Sample size is limited relative to M

Use GLM when:

  • Number of annotations M is small (e.g., M < 10)

  • Sample size is large relative to M

  • All annotations are believed to be informative

  • You want interpretable coefficients for all features

Training Data Requirements:

Good training data should have:

  • Cases: Variants known or strongly suspected to be associated with the phenotype (e.g., from previous GWAS hits, Mendelian variants, functional studies)

  • Controls: Variants known or strongly suspected to be neutral (e.g., common synonymous variants, intergenic variants far from genes)

  • Same trait: Training and testing data should be for the same or similar phenotype

  • Sufficient sample size: At least 50-100 cases recommended; controls should be 2-10x more abundant

Annotation Data Format:

Each annotation should be:

  • Numeric (continuous or ordinal)

  • Higher values should generally indicate higher pathogenicity (if opposite, multiply by -1 before input)

  • Reasonably scaled (extreme ranges may require normalization)

  • Available for all variants (no missing values)

Common annotations include:

  • CADD: Combined Annotation Dependent Depletion

  • PolyPhen-2: Prediction of functional effects

  • SIFT: Sorts Intolerant From Tolerant

  • GERP: Genomic Evolutionary Rate Profiling

  • PhyloP/PhastCons: Conservation scores

  • Functional class indicators (missense, LoF, etc.)

Computational Complexity

O(model_need_N * cost_per_model + N_test * model_need_N * M) where:

  • For LASSO: cost_per_model = O(N_train * M * n_lambda) with cross-validation

  • For GLM: cost_per_model = O(N_train * M^2) for matrix inversion

  • Prediction: O(N_test * model_need_N * M)

Typical runtime: seconds to minutes depending on M and model_need_N.

Assumptions

  • training_caseAnnotation and training_controlAnnotation have the same number and types of columns

  • testing_annotations has the same columns as training data

  • training_control_need_N <= nrow(training_controlAnnotation)

  • No missing values in annotation data

  • Annotations are informative for case/control status

Output Interpretation

  • PI values close to 1: High predicted importance (based on annotations similar to known associated variants)

  • PI values close to 0: Low predicted importance (based on annotations similar to known neutral variants)

  • PI values near 0.5: Uncertain prediction (annotations don't strongly favor associated or neutral)

The estimated PI values are then used as input to Optimal_Weights_M to calculate optimal weights that incorporate biological prior knowledge.

Cross-Validation

For LASSO models, the optimal lambda parameter is selected using 10-fold cross-validation (default in cv.glmnet) to minimize binomial deviance. Each model in the ensemble performs its own cross-validation, which can lead to different levels of sparsity across ensemble members.

References

Zhang, H., Liu, M., Landers, J. E., and Wu, Z. (2026). Integrated weighted association test with application to genetic association studies. The Annals of Applied Statistics, 20(3), 2074-2097. doi:10.1214/26-AOAS2223

Ionita-Laza, I., Lee, S., Makarov, V., Buxbaum, J. D., and Lin, X. (2013). Sequence kernel association tests for the combined effect of rare and common variants. The American Journal of Human Genetics, 92(6), 841-853.

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 1-22.

See also

  • model_PI: Generate ensemble of PI prediction models

  • LASSOmodel: Fit single LASSO model (internal)

  • GLMmodel: Fit single GLM model (internal)

  • Optimal_Weights_M: Calculate optimal weights using PI estimates

Examples

# Example 1: Basic usage with LASSO
set.seed(123)
# Simulate training data: cases have higher annotation scores
n_case <- 100
n_control <- 500
n_anno <- 5
training_case <- matrix(rnorm(n_case * n_anno, mean = 0.5, sd = 1),
                       ncol = n_anno)
training_control <- matrix(rnorm(n_control * n_anno, mean = 0, sd = 1),
                          ncol = n_anno)
colnames(training_case) <- colnames(training_control) <- paste0("Anno", 1:n_anno)

# Simulate testing data
n_test <- 200
testing_anno <- matrix(rnorm(n_test * n_anno, mean = 0.2, sd = 1),
                      ncol = n_anno)
colnames(testing_anno) <- paste0("Anno", 1:n_anno)

# Estimate PI using LASSO ensemble
PI_estimates <- get_PI(
  training_caseAnnotation = training_case,
  training_controlAnnotation = training_control,
  training_control_need_N = 100,
  model_need_N = 10,
  modelType = "LASSO",
  testing_annotations = testing_anno
)

# Check output
summary(PI_estimates)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.01749 0.25439 0.44533 0.45396 0.62292 0.95885 
hist(PI_estimates, main = "Distribution of PI Estimates",
     xlab = "Estimated Variant Importance")


if (FALSE) { # \dontrun{
# Example 2: Using GLM instead of LASSO
PI_glm <- get_PI(
  training_caseAnnotation = training_case,
  training_controlAnnotation = training_control,
  training_control_need_N = 100,
  model_need_N = 10,
  modelType = "GLM",
  testing_annotations = testing_anno
)

# Compare LASSO vs GLM predictions
plot(PI_estimates, PI_glm,
     xlab = "LASSO PI", ylab = "GLM PI",
     main = "LASSO vs GLM PI Estimates")
abline(0, 1, col = "red")
cor(PI_estimates, PI_glm)

# Example 3: Using real annotation data (conceptual)
# Load variant annotations (CADD, PolyPhen, etc.)
# case_anno <- read.csv("causal_variants_annotations.csv")
# control_anno <- read.csv("neutral_variants_annotations.csv")
# test_anno <- read.csv("study_variants_annotations.csv")
#
# PI_real <- get_PI(
#   training_caseAnnotation = as.matrix(case_anno),
#   training_controlAnnotation = as.matrix(control_anno),
#   training_control_need_N = 200,
#   model_need_N = 20,
#   modelType = "LASSO",
#   testing_annotations = as.matrix(test_anno)
# )
} # }