R/aSPUs.R

Defines functions aSPUs

Documented in aSPUs

#' Sum of Powered Score (SPUs) tests and adaptive SPU (aSPUs) test for single trait - SNP set association with GWAS summary statistics.
#'
#' It gives p-values of the SPUs tests and aSPUs test with GWAS summary statistics.
#'
#' @param Zs Z-scores for each SNPs. It could be P-values if the Ps option is TRUE. 
#'
#' @param corSNP Correlation matirx of the SNPs to be tested; estimated from a
#' reference panel (based on the same set of the reference alleles as
#' used in calculating Z-scores).
#'
#' @param pow power used in SPU test. A vector of the powers.
#'
#' @param n.perm number of permutations or bootstraps.
#'
#' @param Ps TRUE if input is p-value, FALSE if input is Z-scores. The default is FALSE.
#'
#' @param prune if it is TRUE, do pruing before the test using pruneSNP function. 
#'
#' @return A list object, Ts : test statistics for the SPU tests (in the order of the specified pow) and finally for the aSPU test.
#'         pvs : p-values for the SPUs and aSPUs tests.
#'
#' @author Il-Youp Kwak and Wei Pan
#'
#' @references
#'
#' Il-Youp Kwak, Wei Pan (2015)
#' Adaptive Gene- and Pathway-Trait Association Testing with GWAS Summary Statistics,
#' Bioinformatics, 32(8), 1178-1184.
#'
#' @examples
#'
#' data(kegg9)
#' ## example analysis using aSPUM test.
#' g <- kegg9$gene.info[1,1]  #  SOAT1 
#' ## Take snps mapped on gene "SOAT1" from the information of gene.info and snp.info. 
#' snps <- which( ( kegg9$snp.info[,2] == kegg9$gene.info[kegg9$gene.info[,1] == g, 2] ) &
#'                  (kegg9$snp.info[,3] > kegg9$gene.info[kegg9$gene.info[,1] == g, 3] ) &
#'                  (kegg9$snp.info[,3] < kegg9$gene.info[kegg9$gene.info[,1] == g, 4] )  )
#' ## Take subsets
#' newP <- kegg9$nP[snps] ;
#' ldsub <- kegg9$ldmatrix[snps, snps];
#'
#' ## Get p-value for gene SOAT1. Read vignette for details.
#' out <- aSPUs(newP, corSNP=ldsub , pow=c(1:8, Inf), n.perm=100, Ps=TRUE)
#'
#' out$Ts
#' # This is a vector of Test Statistics for SPUM and aSPUM tests.
#' # SPUs1 to SPUsInf corresponds with the option pow=c(1:8, Inf)
#' # They are SPUs test statistics.
#' # The last element aSPUs is minimum of them, aSPUs statistic.
#'
#' out$pvs
#' # This is a vector of p-values for SPUs and aSPUs tests.
#' # SPUs1 to SPUsInf corresponds with the option pow=c(1:8, Inf)
#' # They are p-values for corresponding SPUs tests.
#' # The last element is p-value of aSPUs test.
#'
#' @seealso \code{\link{aSPUw}} \code{\link{aSPU}} \code{\link{aSPUsPath}}
#'
#' @export

aSPUs <- function(Zs, corSNP, pow = c(1:8, Inf), n.perm = 1000, Ps = FALSE, prune = TRUE)
{

    if(prune== TRUE) {
        pr <- pruneSNP(corSNP)
        if( length(pr$to.erase) > 0 ) {
            Zs <- Zs[-pr$to.erase]
            corSNP <- corSNP[-pr$to.erase, -pr$to.erase]
        }
    }
    
    if( length(Zs) <= 1 ) {
        stop("less than 1 SNP.")
    }

    if ( !( length(Zs) == dim(corSNP)[1] & dim(corSNP)[2] == dim(corSNP)[1] ) ) {
        stop("dimension do not match. Check dimension of Zs and corSNP.")
    }
    
    n <- length(Zs)
    k <- n
    if(Ps == TRUE)
      Zs <- qnorm(1 - Zs/2)

    U <- Zs
    CovS <- corSNP

    eS <- eigen(CovS, symmetric = TRUE)
    ev <- eS$values

    CovSsqrt <- eS$vectors %*% diag(sqrt(pmax(ev, 0)), k)

#    svd.CovS <- svd(CovS)
#    CovSsqrt <- svd.CovS$u %*% diag(sqrt(svd.CovS$d))

    Ts = rep(NA, length(pow))
    for (j in 1:length(pow)) {
        if (pow[j] < Inf)
            Ts[j] = sum(U^pow[j])
        else Ts[j] = max(abs(U))
    }
    pPerm0 = rep(NA, length(pow))
    T0s = numeric(n.perm)
    s <- sample(1:10^5, 1)
    for (j in 1:length(pow)) {
        set.seed(s)
        for (b in 1:n.perm) {
            U00 <- rnorm(k, 0, 1)
            U0 <- CovSsqrt %*% U00
            if (Ps == TRUE)
              U0 <- abs(U0)
            if (pow[j] < Inf) {
	       T0s[b] = round(sum(U0^pow[j]), digits = 8)
            }
            if (pow[j] == Inf) {
                T0s[b] = round(max(abs(U0)), digits = 8)
            }
        }
        pPerm0[j] = sum(abs(Ts[j]) <= abs(T0s))/n.perm
        P0s = ((n.perm - rank(abs(T0s))) + 1)/(n.perm)
        if (j == 1)
            minp0 = P0s
        else minp0[which(minp0 > P0s)] = P0s[which(minp0 > P0s)]
    }
    Paspu <- (sum(minp0 <= min(pPerm0)) + 1)/(n.perm + 1)
    pvs <- c(pPerm0, Paspu)
    Ts <- c(Ts, min(pPerm0))
    names(Ts) <- c(paste("SPUs", pow, sep = ""), "aSPUs")
    names(pvs) = names(Ts)
    list(Ts = Ts, pvs = pvs)
}

Try the aSPU package in your browser

Any scripts or data that you put into this service are public.

aSPU documentation built on June 29, 2021, 1:06 a.m.