Skip to content

Feature Request: Add pgen support to Step1 #199

Open
Open
Feature Request: Add pgen support to Step1#199

Description

@jjfarrell
jjfarrell
opened on Feb 13, 2026

Adding Plink2 PGEN Support to SAIGE Step 1

Step2 already has pgen support. However, Step1 does not( see #181) . It would be great if plink2 input could be also added to step1 to allow a single format to be used for both steps. Plink2 format is much better for WGS over plink2 and very fast. I am not an expert on C and would not be the best person to produce a pull request. So I asked Claude to come up the modifications need. Hope this helps with moving forward on adding this feature.

Overview

This document outlines the modifications needed to add plink2 pgen/pvar/psam support to SAIGE step 1 (fitNULLGLMM). Step 2 already has pgen support, so we can leverage similar patterns.

Files to Modify

1. R Script: extdata/step1_fitNULLGLMM.R

Add new command-line options for pgen files:

#!/usr/bin/env Rscript
options(stringsAsFactors=F)

library(SAIGE)
require(optparse)

print(sessionInfo())

## set list of cmd line arguments
option_list <- list(
  make_option("--plinkFile", type="character",default="",
    help="Path to plink file for creating the genetic relationship matrix (GRM). minMAFforGRM can be used to specify the minimum MAF of markers in the plink file to be used for constructing GRM. Genetic markers are also randomly selected from the plink file to estimate the variance ratios"),
  
  # NEW: Add pgen support options
  make_option("--pgenFile", type="character",default="",
    help="Path to pgen file prefix (without .pgen/.pvar/.psam extensions) for creating the genetic relationship matrix (GRM). Alternative to --plinkFile for plink2 format files."),
  
  make_option("--phenoFile", type="character", default="",
    help="Path to the phenotype file..."),
  
  # ... rest of existing options ...
)

## Parse arguments
parser <- OptionParser(usage="%prog [options]", option_list=option_list)
args <- parse_args(parser, positional_arguments = 0)
opt <- args$options
print(opt)

# NEW: Validation - ensure only one genotype file format is specified
if(opt$plinkFile != "" && opt$pgenFile != "") {
  stop("Please specify either --plinkFile or --pgenFile, not both\n")
}

if(opt$plinkFile == "" && opt$pgenFile == "") {
  stop("Please specify either --plinkFile or --pgenFile for genotype data\n")
}

# Parse covariates
covars <- strsplit(opt$covarColList,",")[[1]]
qcovars <- strsplit(opt$qCovarColList,",")[[1]]

# ... rest of existing code for parsing numeric options ...

# Set seed
set.seed(1)

# Call fitNULLGLMM with appropriate genotype file
fitNULLGLMM(plinkFile=opt$plinkFile,
            pgenFile=opt$pgenFile,  # NEW: Add pgenFile parameter
            phenoFile = opt$phenoFile,
            phenoCol = opt$phenoCol,
            traitType = opt$traitType,
            invNormalize = opt$invNormalize,
            covarColList = covars,
            qCovarCol = qcovars,
            sampleIDColinphenoFile = opt$sampleIDColinphenoFile,
            tol=opt$tol,
            maxiter=opt$maxiter,
            tolPCG=opt$tolPCG,
            maxiterPCG=opt$maxiterPCG,
            nThreads = opt$nThreads,
            SPAcutoff = opt$SPAcutoff,
            numMarkers = opt$numRandomMarkerforVarianceRatio,
            skipModelFitting = opt$skipModelFitting,
            memoryChunk = opt$memoryChunk,
            tauInit = tauInit,
            LOCO = opt$LOCO,
            traceCVcutoff = opt$traceCVcutoff,
            ratioCVcutoff = opt$ratioCVcutoff,
            outputPrefix = opt$outputPrefix,
            outputPrefix_varRatio = opt$outputPrefix_varRatio,
            IsOverwriteVarianceRatioFile = opt$IsOverwriteVarianceRatioFile,
            useSparseGRMtoFitNULL = opt$useSparseGRMtoFitNULL,
            useSparseGRMforVarRatio = opt$useSparseGRMforVarRatio,
            IsSparseKin = opt$IsSparseKin,
            sparseGRMFile = opt$sparseGRMFile,
            sparseGRMSampleIDFile = opt$sparseGRMSampleIDFile,
            relatednessCutoff = opt$relatednessCutoff,
            isCateVarianceRatio = opt$isCateVarianceRatio,
            cateVarRatioMinMACVecExclude = cateVarRatioMinMACVecExclude,
            cateVarRatioMaxMACVecInclude = cateVarRatioMaxMACVecInclude,
            isCovariateTransform = opt$isCovariateTransform,
            isDiagofKinSetAsOne = opt$isDiagofKinSetAsOne,
            skipVarianceRatioEstimation = opt$skipVarianceRatioEstimation,
            minMAFforGRM = opt$minMAFforGRM)

2. R Function: R/SAIGE_fitGLMM_fast.R

Modify the fitNULLGLMM function signature and implementation:

#' @param plinkFile character. Path to PLINK bed/bim/fam file prefix
#' @param pgenFile character. Path to PLINK2 pgen/pvar/psam file prefix (NEW)
#' @param phenoFile character. Path to phenotype file
#' ... other parameters ...
#' @export
fitNULLGLMM = function(
  plinkFile = "",
  pgenFile = "",  # NEW parameter
  phenoFile = "",
  phenoCol = "",
  traitType = "binary",
  invNormalize = FALSE,
  covarColList = NULL,
  qCovarCol = NULL,
  sampleIDColinphenoFile = "",
  tol=0.02,
  maxiter=20,
  tolPCG=1e-5,
  maxiterPCG=500,
  nThreads = 1,
  SPAcutoff = 2,
  numMarkers = 30,
  skipModelFitting = FALSE,
  memoryChunk = 2,
  tauInit = c(0,0),
  LOCO = TRUE,
  traceCVcutoff = 0.0025,
  ratioCVcutoff = 0.001,
  outputPrefix = "",
  outputPrefix_varRatio = NULL,
  IsOverwriteVarianceRatioFile = FALSE,
  IsSparseKin = FALSE,
  sparseGRMFile = NULL,
  sparseGRMSampleIDFile = NULL,
  numRandomMarkerforSparseKin = 1000,
  relatednessCutoff = 0.125,
  isCateVarianceRatio = FALSE,
  cateVarRatioIndexVec = NULL,
  cateVarRatioMinMACVecExclude = c(0.5,1.5,2.5,3.5,4.5,5.5,10.5,20.5),
  cateVarRatioMaxMACVecInclude = c(1.5,2.5,3.5,4.5,5.5,10.5,20.5),
  isCovariateTransform = TRUE,
  isDiagofKinSetAsOne = FALSE,
  skipVarianceRatioEstimation = FALSE,
  minMAFforGRM = 0.01
) {

  # NEW: Determine which genotype file format is being used
  isPgen <- (pgenFile != "")
  isPlinkBed <- (plinkFile != "")
  
  if(!isPgen && !isPlinkBed) {
    stop("Must specify either plinkFile or pgenFile")
  }
  
  if(isPgen && isPlinkBed) {
    stop("Cannot specify both plinkFile and pgenFile")
  }
  
  # Set genoFile to whichever is provided
  genoFile <- ifelse(isPgen, pgenFile, plinkFile)
  genoType <- ifelse(isPgen, "pgen", "bed")  # NEW: track file type
  
  cat("Using genotype file:", genoFile, "format:", genoType, "\n")
  
  # ... Read phenotype file and process covariates (unchanged) ...
  
  # NEW: Initialize genotype reading based on file type
  if(genoType == "pgen") {
    # Use pgen reading functions (see C++ modifications below)
    setgeno_pgen(paste0(genoFile, ".pgen"),
                 paste0(genoFile, ".pvar"),
                 paste0(genoFile, ".psam"),
                 subSampleInGeno,
                 memoryChunk,
                 isDiagofKinSetAsOne)
  } else {
    # Existing bed file reading
    setgeno(paste0(genoFile, ".bed"),
            subSampleInGeno,
            memoryChunk,
            isDiagofKinSetAsOne)
  }
  
  # Get number of markers - works for both formats
  M = GetNumMarkers()
  
  # ... rest of function continues using the initialized genotype data ...
  # The downstream code doesn't need to know the file format since
  # both pgen and bed readers expose the same interface
  
  # Calculate GRM, fit null model, estimate variance ratios, etc.
  # (existing code unchanged)
  
  closegeno()  # Works for both formats
}

3. C++ Source: src/SAIGE_MAIN.cpp (or similar)

Add pgen reading functionality. Based on step 2's implementation, add these functions:

#include <RcppArmadillo.h>
#include "PLINK2/pgenlib_read.h"  // plink2 library for reading pgen files

// Global pgen reader object
static pgen::PgenReader* g_pgen_ptr = nullptr;
static arma::uvec g_sample_indices;
static size_t g_num_markers = 0;
static size_t g_num_samples = 0;

// [[Rcpp::export]]
void setgeno_pgen(std::string pgen_path,
                  std::string pvar_path, 
                  std::string psam_path,
                  Rcpp::StringVector sampleIDs,
                  double memoryChunk,
                  bool isDiagofKinSetAsOne = false)
{
  try {
    // Open pgen file
    g_pgen_ptr = new pgen::PgenReader(pgen_path);
    
    // Read psam file to get sample IDs
    std::vector<std::string> psam_ids;
    std::ifstream psam_file(psam_path);
    std::string line;
    
    // Skip header if present
    std::getline(psam_file, line);
    if(line.substr(0, 1) != "#") {
      // No header, rewind
      psam_file.clear();
      psam_file.seekg(0);
    }
    
    while(std::getline(psam_file, line)) {
      std::istringstream iss(line);
      std::string fid, iid;
      iss >> fid >> iid;
      psam_ids.push_back(iid);
    }
    psam_file.close();
    
    // Match requested sample IDs to psam file
    g_sample_indices.resize(sampleIDs.size());
    for(size_t i = 0; i < sampleIDs.size(); i++) {
      std::string target_id = Rcpp::as<std::string>(sampleIDs[i]);
      auto it = std::find(psam_ids.begin(), psam_ids.end(), target_id);
      if(it == psam_ids.end()) {
        Rcpp::stop("Sample ID not found in psam file: " + target_id);
      }
      g_sample_indices[i] = std::distance(psam_ids.begin(), it);
    }
    
    g_num_samples = sampleIDs.size();
    
    // Read pvar file to get number of variants
    std::ifstream pvar_file(pvar_path);
    g_num_markers = 0;
    while(std::getline(pvar_file, line)) {
      if(line.substr(0, 1) != "#") {
        g_num_markers++;
      }
    }
    pvar_file.close();
    
    Rcpp::Rcout << "Pgen file initialized: " << g_num_markers 
                << " variants, " << g_num_samples << " samples" << std::endl;
                
  } catch(std::exception& e) {
    Rcpp::stop(std::string("Error initializing pgen file: ") + e.what());
  }
}

// [[Rcpp::export]]
arma::vec Get_Geno_pgen(unsigned int markerIndex, 
                        double minMAF = 0.0,
                        double maxMissingRate = 1.0)
{
  if(g_pgen_ptr == nullptr) {
    Rcpp::stop("Pgen file not initialized. Call setgeno_pgen first.");
  }
  
  if(markerIndex >= g_num_markers) {
    Rcpp::stop("Marker index out of range");
  }
  
  try {
    // Read dosages for this variant
    std::vector<double> dosages(g_num_samples);
    g_pgen_ptr->Read(markerIndex, g_sample_indices, dosages);
    
    // Convert to Armadillo vector
    arma::vec geno(g_num_samples);
    size_t n_missing = 0;
    
    for(size_t i = 0; i < g_num_samples; i++) {
      if(std::isnan(dosages[i]) || dosages[i] < 0) {
        geno[i] = arma::datum::nan;
        n_missing++;
      } else {
        geno[i] = dosages[i];
      }
    }
    
    // Check missing rate
    double missing_rate = (double)n_missing / g_num_samples;
    if(missing_rate > maxMissingRate) {
      geno.fill(arma::datum::nan);
      return geno;
    }
    
    // Calculate MAF (excluding missing)
    arma::vec geno_nomiss = geno(arma::find_finite(geno));
    if(geno_nomiss.n_elem == 0) {
      geno.fill(arma::datum::nan);
      return geno;
    }
    
    double mean_dosage = arma::mean(geno_nomiss);
    double maf = mean_dosage / 2.0;
    if(maf > 0.5) maf = 1.0 - maf;
    
    // Check MAF threshold
    if(maf < minMAF) {
      geno.fill(arma::datum::nan);
      return geno;
    }
    
    // Mean impute missing values
    for(size_t i = 0; i < g_num_samples; i++) {
      if(!std::isfinite(geno[i])) {
        geno[i] = mean_dosage;
      }
    }
    
    return geno;
    
  } catch(std::exception& e) {
    Rcpp::stop(std::string("Error reading pgen variant: ") + e.what());
  }
}

// [[Rcpp::export]]
void closegeno_pgen()
{
  if(g_pgen_ptr != nullptr) {
    delete g_pgen_ptr;
    g_pgen_ptr = nullptr;
  }
  g_sample_indices.clear();
  g_num_markers = 0;
  g_num_samples = 0;
}

// [[Rcpp::export]]
unsigned int GetNumMarkers_pgen()
{
  return g_num_markers;
}

4. Modify Existing C++ Wrapper Functions

Update the existing wrapper functions to handle both file types:

// Modify existing setgeno to be a dispatcher
// [[Rcpp::export]]
void setgeno(std::string geno_path,
             Rcpp::StringVector sampleIDs,
             double memoryChunk,
             bool isDiagofKinSetAsOne = false)
{
  // This is the existing bed file reader - keep as is
  // ... existing bed implementation ...
}

// [[Rcpp::export]]
unsigned int GetNumMarkers()
{
  // Return markers from whichever format is initialized
  if(g_pgen_ptr != nullptr) {
    return GetNumMarkers_pgen();
  } else {
    // Return from bed format (existing code)
    return GetNumMarkers_bed();
  }
}

// [[Rcpp::export]]
void closegeno()
{
  // Close whichever format is open
  if(g_pgen_ptr != nullptr) {
    closegeno_pgen();
  } else {
    closegeno_bed();
  }
}

5. Namespace and Build Configuration

Update NAMESPACE:

# ... existing exports ...
export(setgeno_pgen)
export(Get_Geno_pgen)
export(closegeno_pgen)
export(GetNumMarkers_pgen)

Update src/Makevars:

PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lz
PKG_CXXFLAGS = -I../inst/include/PLINK2

Dependencies

Plink2 Library

SAIGE will need to link against the plink2 library (pgenlib). Options:

  1. Include pgenlib source in SAIGE (recommended):

    • Add inst/include/PLINK2/ directory
    • Include pgenlib headers and source
    • Modify during build

  2. System dependency:

    • Require users to install plink2 development libraries
    • Link against installed library

For step 2, SAIGE likely already has pgenlib included, so check:

  • inst/include/ directory
  • src/ for pgen-related source files

Testing

Test Script

#!/bin/bash

# Test with pgen input
Rscript step1_fitNULLGLMM.R \
  --pgenFile=./input/test_cohort \
  --phenoFile=./input/pheno.txt \
  --phenoCol=y_binary \
  --covarColList=age,sex \
  --sampleIDColinphenoFile=IID \
  --traitType=binary \
  --outputPrefix=./output/test_pgen_step1 \
  --nThreads=4 \
  --IsOverwriteVarianceRatioFile=TRUE

# Verify output files created
ls -lh ./output/test_pgen_step1.*

# Compare with bed format results
Rscript step1_fitNULLGLMM.R \
  --plinkFile=./input/test_cohort_bed \
  --phenoFile=./input/pheno.txt \
  --phenoCol=y_binary \
  --covarColList=age,sex \
  --sampleIDColinphenoFile=IID \
  --traitType=binary \
  --outputPrefix=./output/test_bed_step1 \
  --nThreads=4 \
  --IsOverwriteVarianceRatioFile=TRUE

# Compare variance ratios (should be very similar)
diff ./output/test_pgen_step1.varianceRatio.txt \
     ./output/test_bed_step1.varianceRatio.txt

Example Usage

After implementation, users can run:

# Step 0: Create sparse GRM (optional, for efficiency)
# Can use pgen or bed format here - your choice
Rscript createSparseGRM.R \
  --pgenFile=./data/wgs_common_variants \
  --phenoFile=./data/phenotypes.txt \
  --phenoCol=pheno1 \
  --sampleIDColinphenoFile=IID \
  --nThreads=16 \
  --outputPrefix=./output/sparseGRM \
  --relatednessCutoff=0.125

# Step 1: Fit null model using pgen format
Rscript step1_fitNULLGLMM.R \
  --pgenFile=./data/wgs_common_variants \
  --sparseGRMFile=./output/sparseGRM.mtx \
  --sparseGRMSampleIDFile=./output/sparseGRM.mtx.sampleIDs.txt \
  --useSparseGRMtoFitNULL=TRUE \
  --phenoFile=./data/phenotypes.txt \
  --phenoCol=pheno1 \
  --covarColList=age,sex,PC1,PC2,PC3,PC4,PC5 \
  --qCovarColList=sex \
  --sampleIDColinphenoFile=IID \
  --traitType=binary \
  --outputPrefix=./output/step1_result \
  --nThreads=16 \
  --IsOverwriteVarianceRatioFile=TRUE

# Step 2: Run association tests (already supports pgen!)
Rscript step2_SPAtests.R \
  --pgenPrefix=./data/wgs_all_variants \
  --SAIGEOutputFile=./output/associations_chr22.txt \
  --chrom=22 \
  --minMAF=0.0001 \
  --minMAC=1 \
  --GMMATmodelFile=./output/step1_result.rda \
  --varianceRatioFile=./output/step1_result.varianceRatio.txt \
  --sparseGRMFile=./output/sparseGRM.mtx \
  --sparseGRMSampleIDFile=./output/sparseGRM.mtx.sampleIDs.txt \
  --LOCO=TRUE \
  --is_Firth_beta=TRUE \
  --pCutoffforFirth=0.05

Benefits

  1. No duplicate storage: Use pgen format for both step 1 and step 2
  2. Better WGS support: Pgen handles high-density WGS data more efficiently
  3. Consistency: Same file format throughout pipeline
  4. Storage savings: Pgen is more space-efficient than bed for WGS data

Implementation Notes

  1. AlleleOrder: Pgen files use REF/ALT from pvar. The AlleleOrder parameter (if needed) should default to treating ALT as the effect allele (consistent with pvar format)

  2. Missing data: Pgen handles missing data natively. Ensure proper handling when calculating MAF and imputing missing genotypes

  3. Chromosome encoding: Pvar files can have different chromosome encodings (1-22, chr1-chr22, etc.). Handle flexibly

  4. Multi-allelic variants: Pgen can store multi-allelic sites. Decide on handling (skip or split)

  5. Memory efficiency: Pgen is already memory-efficient. The memoryChunk parameter may need adjustment for pgen vs bed

  6. Error handling: Add robust error checking for file existence, format validation, and sample matching

Validation Strategy

  1. Convert test dataset from bed to pgen using plink2
  2. Run step 1 with both formats on identical phenotype data
  3. Compare:
    • GRM matrices (should be identical)
    • Null model parameters (tau, coefficients)
    • Variance ratios (should match within numerical precision)
  4. Run step 2 with both step 1 outputs
  5. Compare association results (should be identical)

Additional Considerations

createSparseGRM.R

Also update createSparseGRM.R to accept pgen input:

option_list <- list(
  make_option("--plinkFile", type="character", default="",
    help="Path to plink file prefix"),
  make_option("--pgenFile", type="character", default="",
    help="Path to pgen file prefix"),  # NEW
  # ... other options ...
)

This ensures users can create sparse GRM from pgen files directly.

Summary of Changes

Files Modified:

  1. extdata/step1_fitNULLGLMM.R - Add --pgenFile option
  2. R/SAIGE_fitGLMM_fast.R - Add pgenFile parameter, handle both formats
  3. src/ (C++ files) - Add pgen reading functions
  4. NAMESPACE - Export new functions
  5. extdata/createSparseGRM.R - Add pgen support (optional but recommended)

New Dependencies:

  • pgenlib (likely already present from step 2 implementation)

Testing Requirements:

  • Validate identical results between bed and pgen formats
  • Test with various data types (binary/quantitative traits)
  • Test with sparse GRM options
  • Test LOCO functionality

This implementation will bring step 1 to feature parity with step 2 regarding pgen support, eliminating the need to maintain duplicate copies of large WGS datasets in different formats.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions