###################################################
### chunk number 1: getsam
###################################################
#line 118 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
library(ind1KG)
if (!exists("n240"))
  data(n240)


###################################################
### chunk number 2: lksc
###################################################
#line 125 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
names(n240[[1]])


###################################################
### chunk number 3: lkccc
###################################################
#line 129 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
sapply(n240[[1]], class)


###################################################
### chunk number 4: lkddd
###################################################
#line 133 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
lapply(n240[[1]], head, 5)


###################################################
### chunk number 5: lkgf
###################################################
#line 150 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
library(GenomicFeatures)
if (file.exists("hg18.txdb.sqlite"))
txdb = loadFeatures("hg18.txdb.sqlite") else {
txdb <-
  makeTranscriptDbFromUCSC(genome="hg18",
                           tablename="knownGene")
}
txdb


###################################################
### chunk number 6: lktr
###################################################
#line 164 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
tx6 <- transcripts(txdb, vals = list(tx_chrom = "chr6"))
chr6a <- head(unique(sort(ranges(tx6))), 50)
chr6a


###################################################
### chunk number 7: getrs eval=FALSE
###################################################
## #line 179 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
## 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


###################################################
### chunk number 8: dosc eval=FALSE
###################################################
## #line 189 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
## res <- scanBam(fl, param=p1)
## length(res)
## names(res[[1]])


###################################################
### chunk number 9: lknd
###################################################
#line 220 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
data(pup240_disc)
head(pup240_disc, 5)


###################################################
### chunk number 10: om
###################################################
#line 228 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
pup240_disc <- pup240_disc[ pup240_disc$ref != "*", ]
pup240_disc <- pup240_disc[ pup240_disc$ref %in% c("A", "C", "T", "G"), ]
table(pup240_disc$indiv)


###################################################
### chunk number 11: lknd
###################################################
#line 236 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
data(c6snp)
sum( !( pup240_disc$loc  %in% c6snp$chrPosFrom ) )


###################################################
### chunk number 12: lkh
###################################################
#line 242 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
nov <- pup240_disc[ !( pup240_disc$loc  %in% c6snp$chrPosFrom ), ]
table(nov$indiv)


###################################################
### chunk number 13: lkd
###################################################
#line 254 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
library(chopsticks)
data(yri240_6)
names(yri240_6)


###################################################
### chunk number 14: lkhhh
###################################################
#line 262 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
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)


###################################################
### chunk number 15: lkpi
###################################################
#line 273 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
data(pup240_500k)
head(pup240_500k, 2)


###################################################
### chunk number 16: lkdu
###################################################
#line 278 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
pup240_500ku <- pup240_500k[ !duplicated(pup240_500k[,1]),]


###################################################
### chunk number 17: hp
###################################################
#line 282 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
hpup <- pup240_500ku[ pup240_500ku[,1] %in% smethet[,"Position"], ]


###################################################
### chunk number 18: lkhom
###################################################
#line 287 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
hom <- hpup[ hpup[,2] == hpup[,3], ]
hom
smethet[ smethet[,"Position"] %in% hom[,1], ]


###################################################
### chunk number 19: getd
###################################################
#line 298 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
data(gw6c6.snp240)
head(gw6c6.snp240, 4)


###################################################
### chunk number 20: lkloc
###################################################
#line 303 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
hloc6 <- gw6c6.snp240[ gw6c6.snp240$call240 == 2, "physical_pos" ]


###################################################
### chunk number 21: lkll
###################################################
#line 308 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
inds <- which(pup240_500k[,1] %in% hloc6)
table(pup240_500k[inds, 3])


###################################################
### chunk number 22: lko
###################################################
#line 313 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
oloc6 <- gw6c6.snp240[ gw6c6.snp240$call240 %in% c(1,3), "physical_pos" ]
oinds <- which(pup240_500k[,1] %in% oloc6)
table(pup240_500k[oinds, 3])


###################################################
### chunk number 23: dobr eval=FALSE
###################################################
## #line 331 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
## 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"))


###################################################
### chunk number 24: lrnam eval=FALSE
###################################################
## #line 367 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
## tn <- trackNames(br)
## grep("Genes", names(tn), value=TRUE) # many different gene sets
## tn["UCSC Genes"]   # resolve indirection


###################################################
### chunk number 25: doal eval=FALSE
###################################################
## #line 375 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
## rsg <- track(br, "refGene")
## rsgdf <- as.data.frame(rsg)


###################################################
### chunk number 26: lkrs
###################################################
#line 380 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
data(rsgdf)
names(rsgdf)
rsgdf[1:3,1:7]


###################################################
### chunk number 27: lkres
###################################################
#line 388 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
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)


###################################################
### chunk number 28: zzz
###################################################
#line 401 "vignettes/ind1KG/inst/doc/ind1KG.Rnw"
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)