inst/doc/vignette.R

In ASSET: An R package for subset-based association analysis of heterogeneous traits and subtypes

### R code from vignette source 'vignette.Rnw'

###################################################
### code chunk number 1: start
###################################################
library(ASSET)


###################################################
### code chunk number 2: data file
###################################################
datafile <- system.file("sampleData", "vdata.rda", package="ASSET")


###################################################
### code chunk number 3: load data
###################################################
load(datafile)
data1[1:5, ]

SNPs    <- paste("SNP", 1:3, sep="")
nSNP    <- length(SNPs)
studies <- paste("STUDY", 1:4, sep="") 
nStudy  <- length(studies)


###################################################
### code chunk number 4: ncase ncontrol
###################################################
case    <- matrix(data=NA, nrow=nSNP, ncol=nStudy) 
control <- matrix(data=NA, nrow=nSNP, ncol=nStudy) 

for (i in 1:nStudy) {
  data <- eval(parse(text=paste("data", i, sep="")))
  caseVec <- data[, "CC"] == 1
  controlVec <- !caseVec
  for (j in 1:nSNP) {
    temp <- !is.na(data[, SNPs[j]])
    case[j, i] <- sum(caseVec & temp, na.rm=TRUE)
    control[j, i] <- sum(controlVec & temp, na.rm=TRUE)
  }
}
case
control


###################################################
### code chunk number 5: log reg
###################################################
beta  <- matrix(data=NA, nrow=nSNP, ncol=nStudy) 
sigma <- matrix(data=NA, nrow=nSNP, ncol=nStudy) 
for (i in 1:nStudy) {
  data <- eval(parse(text=paste("data", i, sep="")))
  for (j in 1:nSNP) {
    data[, "SNP"] <- data[, SNPs[j]]
    fit <- glm(CC ~ AGE + SNP, data=data, family=binomial())
    coef <- summary(fit)$coefficients
    beta[j, i] <- coef["SNP", 1]
    sigma[j, i] <- coef["SNP", 2]
  }
}
beta
sigma



###################################################
### code chunk number 6: h.traits
###################################################
res <- h.traits(SNPs, studies, beta, sigma, case, control, meta=TRUE)


###################################################
### code chunk number 7: summary table
###################################################
h.summary(res)


###################################################
### code chunk number 8: subset function
###################################################
sub.def <- function(logicalVec) {
  sum <- sum(logicalVec)  
  ret <- all(logicalVec[1:sum])
  ret
}


###################################################
### code chunk number 9: h.traits 2
###################################################
res <- h.traits(SNPs, studies, beta, sigma, case, control, meta=TRUE, 
         zmax.args=list(sub.def=sub.def), pval.args=list(sub.def=sub.def))


###################################################
### code chunk number 10: summary table 2
###################################################
h.summary(res)


###################################################
### code chunk number 11: combine data
###################################################
data <- NULL
for (i in 1:nStudy) {
  temp <- eval(parse(text=paste("data", i, sep="")))
  temp[, "STUDY"]  <- i
  data <- rbind(data, temp)
}


###################################################
### code chunk number 12: dummy vars
###################################################
for (i in 1:nStudy) {
  dvar <- paste("STUDY_", i, sep="")
  data[, dvar] <- as.numeric(data[, "STUDY"] %in% i)
}


###################################################
### code chunk number 13: variables
###################################################
snp.vars <- paste("SNP", 1:3, sep="")
adj.vars <- c("AGE", "STUDY_1", "STUDY_2", "STUDY_3")
types.lab <- paste("TYPE_", 1:3, sep="")


###################################################
### code chunk number 14: h.types
###################################################
ret <- h.types(data, "TYPE", snp.vars, adj.vars, types.lab, "CONTROL",
               logit=TRUE)


###################################################
### code chunk number 15: h.types summary
###################################################
h.summary(ret)


###################################################
### code chunk number 16: sessionInfo
###################################################
sessionInfo()