### R code from vignette source 'vignettes/CGEN/inst/doc/vignette.Rnw' ################################################### ### code chunk number 1: start ################################################### library(CGEN) ################################################### ### code chunk number 2: load data ################################################### data(Xdata, package="CGEN") Xdata[1:5, ] ################################################### ### code chunk number 3: dummy vars ################################################### for (a in unique(Xdata[, "age.group"])) { dummyVar <- paste("age.group_", a, sep="") Xdata[, dummyVar] <- 0 temp <- Xdata[, "age.group"] == a if (any(temp)) Xdata[temp, dummyVar] <- 1 } ################################################### ### code chunk number 4: table ################################################### table(Xdata[, "case.control"], Xdata[, "age.group"], exclude=NULL) ################################################### ### code chunk number 5: snp.logistic ################################################### mainVars <- c("oral.years", "n.children", "age.group_1", "age.group_2", "age.group_3", "age.group_5") fit <- snp.logistic(Xdata, "case.control", "BRCA.status", main.vars=mainVars, int.vars=c("oral.years", "n.children"), strata.var="ethnic.group") ################################################### ### code chunk number 6: summary table ################################################### getSummary(fit) ################################################### ### code chunk number 7: Wald test ################################################### getWaldTest(fit, c("BRCA.status", "BRCA.status:oral.years", "BRCA.status:n.children")) ################################################### ### code chunk number 8: snp.list ################################################### g <- system.file("sampleData", "SNPdata.rda", package="CGEN") snp.list <- list(file=g, file.type=1, delimiter="|", in.miss="NA") ################################################### ### code chunk number 9: pheno.list ################################################### f <- system.file("sampleData", "Xdata.txt", package="CGEN") pheno.list <- list(file=f, id.var="id", file.type=3, delimiter="\t", response.var="case.control", strata.var="ethnic.group", main.vars=c("oral.years", "n.children"), int.vars="oral.years") ################################################### ### code chunk number 10: options list ################################################### out.file <- "snp.scan.logistic.output.file_h5jb47j.txt" op <- list(out.file=out.file) ################################################### ### code chunk number 11: snp.scan.logistic ################################################### ret <- snp.scan.logistic(snp.list, pheno.list, op=op) ################################################### ### code chunk number 12: read data ################################################### x <- read.table(out.file, sep="\t", header=TRUE) x[1:5, ] ################################################### ### code chunk number 13: qq plot ################################################### ret <- QQ.plot(x[, "UML.omnibus.pvalue"]) ################################################### ### code chunk number 14: system.file ################################################### map <- system.file("sampleData", "LocusMapData.txt", package="CGEN") ################################################### ### code chunk number 15: locus map list ################################################### read.table(map, sep="\t", header=1, nrows=5) lmap.list <- list(file=map, file.type=3, header=1, delimiter="\t", snp.var="SNP", chrm.var="CHROMOSOME", loc.var="LOCATION") ################################################### ### code chunk number 16: plot setup ################################################### plot.vars <- c("UML.omnibus.pvalue", "CML.omnibus.pvalue", "EB.omnibus.pvalue") op <- list(pch=c(21, 22, 23), alt.colors=1) #plot.new() ################################################### ### code chunk number 17: chromosome.plot ################################################### ret <- chromosome.plot(out.file, plot.vars, lmap.list, op=op) ################################################### ### code chunk number 18: top hits ################################################### temp <- sort(x[, "UML.omnibus.pvalue"], index.return=TRUE)$ix topHits <- x[temp, ][1:10,] topHits ################################################### ### code chunk number 19: print table ################################################### table(Xdata$case.control, Xdata$ethnic.group) ################################################### ### code chunk number 20: getMatchedSets ################################################### library("cluster") size <- ifelse(Xdata$ethnic.group == 3, 3, 4) d <- daisy(Xdata[,c("age.group_1","gynSurgery.history","BRCA.history")]) mx <- getMatchedSets(d, CC=TRUE, NN=TRUE, ccs.var = Xdata$case.control, strata.var = Xdata$ethnic.group, size = size, fixed = FALSE) ################################################### ### code chunk number 21: Xdata ################################################### mx$CC[1:10] mx$tblCC Xdata <- cbind(Xdata, CCStrat = mx$CC, NNStrat = mx$NN) Xdata[1:5,] ################################################### ### code chunk number 22: snp.matched ################################################### intVars <- ~ oral.years + n.children snpVars <- ~ BRCA.status fit <- snp.matched(Xdata, "case.control", snp.vars=snpVars, main.vars=intVars, int.vars=intVars, cc.var="CCStrat", nn.var="NNStrat") ################################################### ### code chunk number 23: summary table 2 ################################################### getSummary(fit, method = c("CLR", "CCL")) ################################################### ### code chunk number 24: Wald test 2 ################################################### getWaldTest(fit$HCL, c( "BRCA.status" , "BRCA.status:oral.years" , "BRCA.status:n.children")) ################################################### ### code chunk number 25: sessionInfo ################################################### sessionInfo()