### R code from vignette source 'vignettes/ind1KG/inst/doc/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)