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

###################################################
### code chunk number 1: getlib
###################################################
library(gwascat)


###################################################
### code chunk number 2: lk1
###################################################
if (length(grep("gwascat", search()))>0) detach("package:gwascat")


###################################################
### code chunk number 3: lk2
###################################################
library(gwascat)
objects("package:gwascat")


###################################################
### code chunk number 4: lkgr
###################################################
data(ebicat37)


###################################################
### code chunk number 5: lktr
###################################################
topTraits(ebicat37)


###################################################
### code chunk number 6: lklocs
###################################################
subsetByTraits(ebicat37, tr="LDL cholesterol")[1:3]


###################################################
### code chunk number 7: lkm
###################################################
gwtrunc = ebicat37
requireNamespace("S4Vectors")
mcols = S4Vectors::mcols
mlpv = mcols(ebicat37)$PVALUE_MLOG
mlpv = ifelse(mlpv > 25, 25, mlpv)
S4Vectors::mcols(gwtrunc)$PVALUE_MLOG = mlpv
library(GenomeInfoDb)
seqlevelsStyle(gwtrunc) = "UCSC"
gwlit = gwtrunc[ which(as.character(seqnames(gwtrunc)) %in% c("chr4", "chr5", "chr6")) ]
library(ggbio)
mlpv = mcols(gwlit)$PVALUE_MLOG
mlpv = ifelse(mlpv > 25, 25, mlpv)
S4Vectors::mcols(gwlit)$PVALUE_MLOG = mlpv


###################################################
### code chunk number 8: dodie (eval = FALSE)
###################################################
## pdf(file="litman.pdf", width=8, height=2)


###################################################
### code chunk number 9: dorunap (eval = FALSE)
###################################################
## methods:::bind_activation(FALSE)
## autoplot(gwlit, geom="point", aes(y=PVALUE_MLOG), xlab="chr4-chr6")


###################################################
### code chunk number 10: dodie (eval = FALSE)
###################################################
## dev.off()


###################################################
### code chunk number 11: lktrm
###################################################
args(traitsManh)


###################################################
### code chunk number 12: lkn (eval = FALSE)
###################################################
## pdf(file="annman.pdf", width=8, height=2)


###################################################
### code chunk number 13: runit (eval = FALSE)
###################################################
## traitsManh(gwtrunc)


###################################################
### code chunk number 14: dodie (eval = FALSE)
###################################################
## dev.off()


###################################################
### code chunk number 15: getgd
###################################################
gffpath = system.file("gff3/chr17_43000000_45000000.gff3", package="gwascat")
library(rtracklayer)
c17tg = import(gffpath)


###################################################
### code chunk number 16: lkgvt
###################################################
c17td = c17tg[ which(S4Vectors::mcols(c17tg)$type == "Degner_dsQTL") ]
library(Gviz)
dsqs = DataTrack( c17td, chrom="chr17", genome="hg19", data="score",
  name="dsQTL")


###################################################
### code chunk number 17: dost
###################################################
g2 = GRanges(seqnames="chr17", IRanges(start=4.3e7, width=2e6))

seqlevelsStyle(ebicat37) = "UCSC"
basic = gwcex2gviz(basegr = ebicat37, contextGR=g2, plot.it=FALSE) 


###################################################
### code chunk number 18: getstr
###################################################
c17ts = c17tg[ which(S4Vectors::mcols(c17tg)$type == "Stranger_eqtl") ]
eqloc = AnnotationTrack(c17ts,  chrom="chr17", genome="hg19", name="Str eQTL")
displayPars(eqloc)$col = "black"
displayPars(dsqs)$col = "red"
integ = list(basic[[1]], eqloc, dsqs, basic[[2]], basic[[3]])


###################################################
### code chunk number 19: dogr
###################################################
pdf(file="integ.pdf", width=9, height=4)


###################################################
### code chunk number 20: doplot
###################################################
plotTracks(integ)


###################################################
### code chunk number 21: dono
###################################################
dev.off()


###################################################
### code chunk number 22: dobig (eval = FALSE)
###################################################
## library(pd.genomewidesnp.6)
## con = pd.genomewidesnp.6@getdb()
## locon6 = dbGetQuery(con, 
##    "select dbsnp_rs_id, chrom, physical_pos from featureSet limit 10000")


###################################################
### code chunk number 23: doloc
###################################################
data(locon6)
rson6 = as.character(locon6[[1]])
rson6[1:5]


###################################################
### code chunk number 24: lkdtab
###################################################
intr = ebicat37[ intersect(getRsids(ebicat37), rson6) ]
sort(table(getTraits(intr)), decreasing=TRUE)[1:10]


###################################################
### code chunk number 25: lkexp
###################################################
gr6.0 = GRanges(seqnames=ifelse(is.na(locon6$chrom),0,locon6$chrom), 
       IRanges(ifelse(is.na(locon6$phys),1,locon6$phys), width=1))
S4Vectors::mcols(gr6.0)$rsid = as.character(locon6$dbsnp_rs_id)
seqlevels(gr6.0) = paste("chr", seqlevels(gr6.0), sep="")


###################################################
### code chunk number 26: dosub
###################################################
ag = function(x) as(x, "GRanges")
ovraw = suppressWarnings(subsetByOverlaps(ag(ebicat37), gr6.0))
length(ovraw)
ovaug = suppressWarnings(subsetByOverlaps(ag(ebicat37+500), gr6.0))
length(ovaug)


###################################################
### code chunk number 27: dosub2
###################################################
rawrs = mcols(ovraw)$SNPS
augrs = mcols(ovaug)$SNPS
ebicat37[augrs]


###################################################
### code chunk number 28: lkrelax
###################################################
setdiff( getTraits(ebicat37[augrs]), getTraits(ebicat37[rawrs]) )


###################################################
### code chunk number 29: lkcout
###################################################
data(gg17N) # translated from GGdata chr 17 calls using ABmat2nuc
gg17N[1:5,1:5]


###################################################
### code chunk number 30: dorun
###################################################
h17 = riskyAlleleCount(gg17N, matIsAB=FALSE, chr="ch17",
 gwwl = ebicat37)
h17[1:5,1:5]
table(as.numeric(h17))


###################################################
### code chunk number 31: domo
###################################################
gwr = ebicat37
gwr = gwr[colnames(h17),]
S4Vectors::mcols(gwr) = cbind(mcols(gwr), DataFrame(t(h17)))
sn = rownames(h17)
gwr[,c("DISEASE/TRAIT", sn[1:4])]


###################################################
### code chunk number 32: getbase
###################################################
data(low17)
low17


###################################################
### code chunk number 33: lkggd
###################################################
data(g17SM)
g17SM


###################################################
### code chunk number 34: dog
###################################################
data(gw6.rs_17)
g17SM = g17SM[, intersect(colnames(g17SM), gw6.rs_17)]
dim(g17SM)


###################################################
### code chunk number 35: lkrul
###################################################
if (!exists("rules_6.0_1kg_17")) data(rules_6.0_1kg_17)
rules_6.0_1kg_17[1:5,]


###################################################
### code chunk number 36: lksum
###################################################
summary(rules_6.0_1kg_17)


###################################################
### code chunk number 37: lkov
###################################################
length(intersect(colnames(g17SM), mcols(ebicat37)$SNPS))


###################################################
### code chunk number 38: doimp
###################################################
exg17 = impute.snps(rules_6.0_1kg_17, g17SM)


###################################################
### code chunk number 39: lkl
###################################################
length(intersect(colnames(exg17), mcols(ebicat37)$SNPS))


###################################################
### code chunk number 40: getdo
###################################################
library(DO.db)
DO()


###################################################
### code chunk number 41: getallt
###################################################
alltob = unlist(mget(mappedkeys(DOTERM), DOTERM))
allt = sapply(alltob, Term)
allt[1:5]


###################################################
### code chunk number 42: dohit
###################################################
cattra = mcols(ebicat37)$Disease.Trait
mat = match(tolower(cattra), tolower(allt))
catDO = names(allt)[mat]
na.omit(catDO)[1:50]
mean(is.na(catDO))


###################################################
### code chunk number 43: dogr
###################################################
unique(cattra[is.na(catDO)])[1:20]
nomatch = cattra[is.na(catDO)]
unique(nomatch)[1:5]


###################################################
### code chunk number 44: dopar
###################################################
hpobo = gzfile(dir(system.file("obo", package="gwascat"), pattern="hpo", full=TRUE))
HPOgraph = obo2graphNEL(hpobo)
close(hpobo)


###################################################
### code chunk number 45: dohpot
###################################################
requireNamespace("graph")
hpoterms = unlist(graph::nodeData(HPOgraph, graph::nodes(HPOgraph), "name"))
hpoterms[1:10]


###################################################
### code chunk number 46: lkint
###################################################
intersect(tolower(nomatch), tolower(hpoterms))


###################################################
### code chunk number 47: chkcadd (eval = FALSE)
###################################################
## g3 = as(ebicat37, "GRanges")
## bg3 = bindcadd_snv( g3[which(seqnames(g3)=="chr3")][1:20] )
## inds = ncol(mcols(bg3))
## bg3[, (inds-3):inds]


###################################################
### code chunk number 48: getrs
###################################################
if ("SNPlocs.Hsapiens.dbSNP144.GRCh37" %in% installed.packages()[,1]) {
 library(SNPlocs.Hsapiens.dbSNP144.GRCh37)
 chklocs("20", ebicat37)
}