inst/doc/example4.R

In SEAGLE: Scalable Exact Algorithm for Large-Scale Set-Based Gene-Environment Interaction Tests

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(SEAGLE)

## -----------------------------------------------------------------------------
# Read in gene list
dir <- "../inst/extdata/"
genelist <- read.delim(paste0(dir, "glist-hg19_chr22_example"), sep=" ", header=FALSE)

# Identify number of genes in genelist
num_genes <- dim(genelist)[1]

# Generate synthetic phenotype and covariate data
n <- 1092 # number of study participants

set.seed(1)
y <- 2 * rnorm(n)

set.seed(2)
X <- as.matrix(rnorm(n))

set.seed(3)
E <- as.matrix(rnorm(n))

## -----------------------------------------------------------------------------
# Initialize output containers for T and p-value for each gene
T.list <- numeric(length=num_genes)
pv.list <- numeric(length=num_genes)

# Run SEAGLE on each gene in genelist
for (i in 1:num_genes) {
  
  # Identify current gene
  gene_name <- genelist[i,4]
  
  # Obtain G
  Gtemp <- read.delim(paste0(dir, gene_name, ".raw"), sep=" ")
  G <- as.matrix(Gtemp[,-c(1:6)])
  L <- dim(G)[2]                ## number of SNPs
  
  # Make weights
  avg_newsnp <- colMeans(G)
  fA  = avg_newsnp/2            ## freq of allele "A"
  fa  = 1-fA                    ## freq of allele "a"
  maf = fA; maf[fA>0.5]= fa[fA>0.5]## maf should be b/w 0 and 0.5
  G = G[ ,maf>0]                ## only keep those SNPs with MAF>0
  maf <- maf[maf > 0]           ## Keep only MAF > 0)
  wt   = sqrt(maf^(-3/4))       ## Note we take the square root here
  if (length(wt) > 1) {
    G_final    = G %*% diag(wt) ## Input this G
  } else {
    T.list[i] <- NA
    pv.list[i] <- NA
    next
  }

  # Run SEAGLE
  objSEAGLE <- prep.SEAGLE(y=y, X=X, intercept=0, 
                           E=E, G=G_final)
  res <- SEAGLE(objSEAGLE, init.tau=0.5, init.sigma=0.5)
  
  # Save T and p-value into output lists
  T.list[i] <- res$T
  pv.list[i] <- res$pv
}

## -----------------------------------------------------------------------------
resMat <- cbind(T.list, pv.list)
colnames(resMat) <- c("T", "p-value")
resMat