inst/doc/cpvSNP.R

In cpvSNP: Gene set analysis methods for SNP association p-values that lie in genes in given gene sets

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

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: loadLibrariesAndData
###################################################
library(cpvSNP)
data(geneSetAnalysis)
names(geneSetAnalysis)


###################################################
### code chunk number 3: createArrayData
###################################################
arrayDataGR <- createArrayData(geneSetAnalysis[["arrayData"]], 
    positionName="Position")
class(arrayDataGR)


###################################################
### code chunk number 4: geneSets
###################################################
geneSets <- geneSetAnalysis[["geneSets"]]
geneSets
length(geneSets)
head(geneIds(geneSets[[1]]))
details(geneSets[[1]])
head(geneIds(geneSets[[2]]))


###################################################
### code chunk number 5: geneToSNPList
###################################################
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
genesHg19 <- genes(txdb)

snpsGSC <- geneToSNPList(geneSets, arrayDataGR, genesHg19)
class(snpsGSC)


###################################################
### code chunk number 6: readyToRun
###################################################
arrayDataGR
snpsGSC


###################################################
### code chunk number 7: indep
###################################################
indep <- geneSetAnalysis[["indepSNPs"]]
head(indep)
dim(indep)


###################################################
### code chunk number 8: createIndepP
###################################################
pvals <- arrayDataGR$P[is.element(arrayDataGR$SNP, indep$V1)]
names(pvals) <- arrayDataGR$SNP[is.element(arrayDataGR$SNP, indep$V1)]
head(pvals)


###################################################
### code chunk number 9: glossi
###################################################
gRes <- glossi(pvals, snpsGSC)
gRes

gRes2 <- glossi(pvals, snpsGSC[[1]])
gRes2

pValue(gRes)
degreesOfFreedom(gRes)
statistic(gRes)


###################################################
### code chunk number 10: publishMultiple (eval = FALSE)
###################################################
## pvals <- p.adjust( unlist( pValue(gRes) ), method= "BH" )
## library(ReportingTools)
## report <- HTMLReport (shortName = "cpvSNP_glossiResult",
##         title = "GLOSSI Results", reportDirectory = "./reports")
## publish(geneSets, report, annotation.db = "org.Hs.eg",
##         setStats = unlist(statistic (gRes)), 
##         setPValues = pvals)
## finish(report)


###################################################
### code chunk number 11: hapmap (eval = FALSE)
###################################################
## download.file(
##     url="http://hapmap.ncbi.nlm.nih.gov/genotypes/hapmap3_r3/plink_format/hapmap3_r3_b36_fwd.consensus.qc.poly.ped.gz",
##     destfile="hapmap3_r3_b36_fwd.consensus.qc.poly.ped.gz")
## download.file(
##     url="http://hapmap.ncbi.nlm.nih.gov/genotypes/hapmap3_r3/plink_format/hapmap3_r3_b36_fwd.consensus.qc.poly.map.gz",
##     destfile="hapmap3_r3_b36_fwd.consensus.qc.poly.map.gz")
## system("gunzip hapmap3_r3_b36_fwd.consensus.qc.poly.ped.gz")
## system("gunzip hapmap3_r3_b36_fwd.consensus.qc.poly.map.gz")
## system("plink --file hapmap3_r3_b36_fwd.consensus.qc.poly --make-bed --chr 1")
## 
## genos <- read.plink(bed, bim, fam)
## genos$genotypes
## head(genos$map)
## x <- genos[,is.element(genos$map$snp.name,c(geneIds(snpsGSC[[2]])))]
## ldMat <- ld(x,y=x,stats="R.squared")


###################################################
### code chunk number 12: vegas (eval = FALSE)
###################################################
## ldMat <- geneSetAnalysis[["ldMat"]]
## vRes <- vegas(snpsGSC[1], arrayDataGR, ldMat)
## vRes
## summary(unlist(simulatedStats(vRes)))
## 
## pValue(vRes)
## degreesOfFreedom(vRes)
## statistic(vRes)


###################################################
### code chunk number 13: plotb
###################################################
plotPvals(gRes, main="GLOSSI P-values")


###################################################
### code chunk number 14: plot2
###################################################
pvals <- arrayDataGR$P
names(pvals) <- arrayDataGR$SNP
assocPvalBySetPlot(pvals, snpsGSC[[2]])


###################################################
### code chunk number 15: sessInfo
###################################################
toLatex(sessionInfo())


###################################################
### code chunk number 16: resetOptions
###################################################
options(prompt="> ", continue="+ ")