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

###################################################
### code chunk number 1: getsam
###################################################
library(ind1KG)
if (!exists("n240"))
  data(n240)


###################################################
### code chunk number 2: lksc
###################################################
names(n240[[1]])


###################################################
### code chunk number 3: lkccc
###################################################
sapply(n240[[1]], class)


###################################################
### code chunk number 4: lkddd
###################################################
lapply(n240[[1]], head, 5)


###################################################
### code chunk number 5: lkgf
###################################################
library(GenomicFeatures)
library(TxDb.Hsapiens.UCSC.hg18.knownGene)
txdb = TxDb.Hsapiens.UCSC.hg18.knownGene
txdb


###################################################
### code chunk number 6: lktr
###################################################
tx6 <- transcripts(txdb, vals = list(tx_chrom = "chr6"))
chr6a <- head(unique(sort(ranges(tx6))), 50)
chr6a


###################################################
### code chunk number 7: getrs (eval = FALSE)
###################################################
## library(Rsamtools)
## p1 <- ScanBamParam(which=RangesList(`6`=chr6a))
## fl <- "/mnt/data/stvjc/1000GENOMES/NA19240.chrom6.SLX.maq.SRP000032.2009_07.bam"
## unix.time(cnt <- countBam(fl, param=p1))
## sum(cnt$records) # 2103439


###################################################
### code chunk number 8: dosc (eval = FALSE)
###################################################
## res <- scanBam(fl, param=p1)
## length(res)
## names(res[[1]])


###################################################
### code chunk number 9: lknd
###################################################
data(pup240_disc)
head(pup240_disc, 5)


###################################################
### code chunk number 10: om
###################################################
pup240_disc <- pup240_disc[ pup240_disc$ref != "*", ]
pup240_disc <- pup240_disc[ pup240_disc$ref %in% c("A", "C", "T", "G"), ]
table(pup240_disc$indiv)


###################################################
### code chunk number 11: lknd
###################################################
data(c6snp)
sum( !( pup240_disc$loc  %in% c6snp$chrPosFrom ) )


###################################################
### code chunk number 12: lkh
###################################################
nov <- pup240_disc[ !( pup240_disc$loc  %in% c6snp$chrPosFrom ), ]
table(nov$indiv)


###################################################
### code chunk number 13: lkd
###################################################
library(chopsticks)
data(yri240_6)
names(yri240_6)


###################################################
### code chunk number 14: lkhhh
###################################################
snps <- as(yri240_6[[1]], "character")
hets <- which(snps == "A/B") 
rshet <- colnames(snps)[hets]
smet <- yri240_6[[2]]
smethet <- smet[hets,]
head(smethet, 5)


###################################################
### code chunk number 15: lkpi
###################################################
data(pup240_500k)
head(pup240_500k, 2)


###################################################
### code chunk number 16: lkdu
###################################################
pup240_500ku <- pup240_500k[ !duplicated(pup240_500k[,1]),]


###################################################
### code chunk number 17: hp
###################################################
hpup <- pup240_500ku[ pup240_500ku[,1] %in% smethet[,"Position"], ]


###################################################
### code chunk number 18: lkhom
###################################################
hom <- hpup[ hpup[,2] == hpup[,3], ]
hom
smethet[ smethet[,"Position"] %in% hom[,1], ]


###################################################
### code chunk number 19: getd
###################################################
data(gw6c6.snp240)
head(gw6c6.snp240, 4)


###################################################
### code chunk number 20: lkloc
###################################################
hloc6 <- gw6c6.snp240[ gw6c6.snp240$call240 == 2, "physical_pos" ]


###################################################
### code chunk number 21: lkll
###################################################
inds <- which(pup240_500k[,1] %in% hloc6)
table(pup240_500k[inds, 3])


###################################################
### code chunk number 22: lko
###################################################
oloc6 <- gw6c6.snp240[ gw6c6.snp240$call240 %in% c(1,3), "physical_pos" ]
oinds <- which(pup240_500k[,1] %in% oloc6)
table(pup240_500k[oinds, 3])


###################################################
### code chunk number 23: dobr (eval = FALSE)
###################################################
## library(IRanges)
## nloc <- nov$loc
## nrd <- RangedData( IRanges(nloc, nloc) )
## names(nrd) <- "chr6"
## library(rtracklayer)
## br <- browserSession("UCSC")
## br[["novo"]] <- nrd
## v1 <- browserView(br, GenomicRanges(1, 1e7, "chr6"))


###################################################
### code chunk number 24: lrnam (eval = FALSE)
###################################################
## tn <- trackNames(br)
## grep("Genes", names(tn), value=TRUE) # many different gene sets
## tn["UCSC Genes"]   # resolve indirection


###################################################
### code chunk number 25: doal (eval = FALSE)
###################################################
## rsg <- track(br, "refGene")
## rsgdf <- as.data.frame(rsg)


###################################################
### code chunk number 26: lkrs
###################################################
data(rsgdf)
names(rsgdf)
rsgdf[1:3,1:7]


###################################################
### code chunk number 27: lkres
###################################################
library(org.Hs.eg.db)
rsgn <- as.character(rsgdf$name)
eid <- mget(rsgn, revmap(org.Hs.egREFSEQ), ifnotfound=NA)
eid <- na.omit(unlist(eid))
sym <- mget(eid, org.Hs.egSYMBOL, ifnotfound=NA)
head(unlist(sym), 10)


###################################################
### code chunk number 28: zzz
###################################################
nloc <- nov$loc  # this one is evaluated
nranges <- IRanges(nloc, nloc)
granges <- IRanges(rsgdf$start, rsgdf$end) # no guarantee of annotation
length(nranges)
length(granges)
sum(nranges %in% granges)
head(match(nranges,granges), 200)