R/staar2aistaar_nullmodel.R

In xihaoli/STAARpipeline: STAARpipeline for Analyzing Whole-Genome/Whole-Exome Sequencing Data

#' Transforming the null model object fitted for STAAR to the null model object to be used for the ancestry-informed (AI) framework. 
#'
#' The \code{staar2aistaar_nullmodel} function takes in the object from fitting the null model for STAAR analyses
#' and transforms it to the object from fitting the null model to be used for AI framework.
#' @param obj_nullmodel_staar an object from fitting the null model, which is the
#' output from \code{fit_nullmodel} function in the STAAR package.
#' @param pop.groups a vector of defined ancestries for all individuals. 
#' @param B a positive numerical value for the number of base tests for ancestry-informed ensemble testing. 
#' @param seed a numerical value to set the initial seed for generating ensemble weights. 
#' @return An object from fitting the null model for related samples to be used for the AI framework,
#' which is an option for output from \code{\link{fit_nullmodel}} function.
#' @references Li, X., Li, Z., et al. (2020). Dynamic incorporation of multiple
#' in silico functional annotations empowers rare variant association analysis of
#' large whole-genome sequencing studies at scale. \emph{Nature Genetics}, \emph{52}(9), 969-983.
#' (\href{https://doi.org/10.1038/s41588-020-0676-4}{pub})
#' @references Li, Z., Li, X., et al. (2022). A framework for detecting
#' noncoding rare-variant associations of large-scale whole-genome sequencing
#' studies. \emph{Nature Methods}, \emph{19}(12), 1599-1611.
#' (\href{https://doi.org/10.1038/s41592-022-01640-x}{pub})
#' @export

staar2aistaar_nullmodel <- function(obj_nullmodel_staar, pop.groups = NULL, B = NULL, seed = 7590){

    obj_nullmodel_aistaar <- NULL

    obj_nullmodel_aistaar <- obj_nullmodel_staar

    K <- length(unique(pop.groups))
    if(K >= 2 && !is.null(B))
    {
        set.seed(seed)
        obj_nullmodel_aistaar$pop_weights_1_1 <- cbind(rep(1, K), matrix(abs(rnorm(K * B)), nrow = K))
        obj_nullmodel_aistaar$pop_weights_1_25 <- cbind(rep(1, K), matrix(abs(rnorm(K * B)), nrow = K))
        obj_nullmodel_aistaar$pop.groups <- pop.groups
    }

    return(obj_nullmodel_aistaar)
}