################################################### ### chunk number 1: wk1 ################################################### library(GGtools) if (!exists("hmceuB36.2021")) data(hmceuB36.2021) f1 = gwSnpTests(genesym("CPNE1")~male, hmceuB36.2021, chrnum(20)) topSnps(f1) ################################################### ### chunk number 2: f2 ################################################### plot(f1) ################################################### ### chunk number 3: lkgws ################################################### showMethods("gwSnpTests") ################################################### ### chunk number 4: lkf1 ################################################### f1 class(f1) getClass(class(f1)) ################################################### ### chunk number 5: lkhm ################################################### hmceuB36.2021 names(pData(hmceuB36.2021)) table(hmceuB36.2021$isFounder) ################################################### ### chunk number 6: lkg ################################################### s20 = smList(hmceuB36.2021)[["20"]] s20 as(s20, "matrix")[1:4,1:4] # show raw as(s20, "character")[1:4,1:4] as(s20, "numeric")[1:4,1:4] ################################################### ### chunk number 7: sol1 ################################################### #plot_EvG(genesym("CPNE1"), rsid("rs17093026"), hmceuB36.2021[chrnum(20),]) hmff = hmceuB36.2021[, which(hmceuB36.2021$isFounder)] hmfg = GTFfilter( hmff, .05 ) f1f = gwSnpTests(genesym("CPNE1")~male, hmfg, chrnum(20)) plot_EvG(genesym("CPNE1"), rsid("rs11698601"), hmfg[chrnum(20),]) ################################################### ### chunk number 8: dosrs ################################################### SS = snp.rhs.tests(exprs(hmceuB36.2021)["GI_23397697-A",]~male, data=pData(hmceuB36.2021), snp.data=s20, family="gaussian") -> SS class(SS) length(names(SS)) sum(is.na(p.value(SS))) SS[1:4,] SS["rs6060535",] ################################################### ### chunk number 9: lk eval=FALSE ################################################### ## table(as(smList(hmceuB36.2021)[[1]][,10], "character")) ################################################### ### chunk number 10: dot eval=FALSE ################################################### ## Mtests = sapply(1:100, function(x) topSnps(gwSnpTests(genesym("CPNE1")~male, ## sms=permEx(hmceuB36.2021), cnum=chrnum(20)))[[1]][1]) ################################################### ### chunk number 11: lk2 ################################################### gr1 = as(f1, "GRanges") gr1[1:4,] gr2 = as(f1f, "GRanges") ################################################### ### chunk number 12: getg ################################################### set.seed(1234) library(illuminaHumanv1.db) g20 = get("20", revmap(illuminaHumanv1CHR)) samp = sample(g20, size=50, replace=FALSE) hlit = hmceuB36.2021[ probeId(samp), ] ################################################### ### chunk number 13: domc ################################################### try(system("rm -rf dem50")) if (try(require(multicore))) { unix.time(e1 <- eqtlTests(hlit, ~male, geneApply = mclapply, targdir="dem50" )) } ################################################### ### chunk number 14: dom2 eval=FALSE ################################################### ## unix.time(e1 <- eqtlTests(hlit, ~male, geneApply = lapply, targdir="dem50" )) ################################################### ### chunk number 15: aaa ################################################### e1 class(e1) getClass(class(e1)) ################################################### ### chunk number 16: lkff ################################################### e1@fflist[[1]][1:4,1:4] e1@shortfac ################################################### ### chunk number 17: doi ################################################### e1[rsid("rs4814683"), probeId("GI_26051259-I")] ################################################### ### chunk number 18: lku ################################################### topFeats(probeid="GI_26051259-I", mgr=e1, ffind=1, anno="illuminaHumanv1.db") topFeats(rsid="rs708954", mgr=e1, ffind=1, anno="illuminaHumanv1.db") topFeats(sym="SRC", mgr=e1, ffind=1, anno="illuminaHumanv1.db") ################################################### ### chunk number 19: getg21 ################################################### g21 = get("21", revmap(illuminaHumanv1CHR)) samp21 = sample(g21, size=50, replace=FALSE) hlit21 = hmceuB36.2021[ probeId(samp21), ] ################################################### ### chunk number 20: dot ################################################### try(system("rm -rf dem50.21")) unix.time(e2 <- eqtlTests(hlit21, ~male, geneApply = mclapply, targdir="dem50.21" )) e2 ################################################### ### chunk number 21: lkdir ################################################### try(system("rm -rf d2021_probetab.ff")) try(system("rm -rf d2021_snptab.ff")) d1 = mkCisTransDirector(list(e1,e2), "d2021", "snptab", "probetab", "illuminaHumanv1.db" ) d1 d1[ "rs6090120", "GI_42734342-S" ] ################################################### ### chunk number 22: getmax ################################################### s = rep(NA,4) g = rep(NA,4) # we search the 2 chromosomes by 2 gene sets for highest score s[1] = which.max(apply(as.ram(e1@fflist[[1]]),1,max)) g[1] = which.max(apply(as.ram(e1@fflist[[1]]),2,max)) s[2] = which.max(apply(as.ram(e1@fflist[[2]]),1,max)) g[2] = which.max(apply(as.ram(e1@fflist[[2]]),2,max)) s[3] = which.max(apply(as.ram(e2@fflist[[1]]),1,max)) g[3] = which.max(apply(as.ram(e2@fflist[[1]]),2,max)) s[4] = which.max(apply(as.ram(e2@fflist[[2]]),1,max)) g[4] = which.max(apply(as.ram(e2@fflist[[2]]),2,max)) ind = which.max( c( e1@fflist[[1]][s[1],g[1]], e1@fflist[[2]][s[2],g[2]], e2@fflist[[1]][s[3],g[3]], e2@fflist[[2]][s[4],g[4]] )) if (ind %in% c(1,2)) { # max is in e1 pr = colnames(e1@fflist[[ind]])[g[ind]] rs = rownames(e1@fflist[[ind]])[s[ind]] } else { # max is in e2 pr = colnames(e2@fflist[[ind-2]])[g[ind]] rs = rownames(e2@fflist[[ind-2]])[s[ind]] } ################################################### ### chunk number 23: lkdd ################################################### if (ind %in% c(1,2)) plot_EvG(probeId(pr), rsid(rs), hlit) else plot_EvG(probeId(pr), rsid(rs), hlit21) if (ind %in% c(1,3)) topFeats(probeid=pr, mgr=d1, ffind=1, anno="illuminaHumanv1.db") else topFeats(probeid=pr, mgr=d1, ffind=2, anno="illuminaHumanv1.db") ################################################### ### chunk number 24: lkgg eval=FALSE ################################################### ## ggg = getGRanges(e2, 1, "GI_25777639-A", "chr20", getNamedLocs(chr="chr20")) ################################################### ### chunk number 25: rpup ################################################### library(GGtools) library(Rsamtools) pup238_6 = readPileup(system.file("pup/combn238_chr6.pup", package="GGtools"), variant="SNP") pup238_6[1:4,] ################################################### ### chunk number 26: getc ################################################### getPupCalls = function(pupfn) { tmp = strsplit(readLines(pupfn), "\t") sapply(tmp, "[", 9) } callp = getPupCalls( system.file("pup/combn238_chr6.pup", package="GGtools") ) elementMetadata(pup238_6)$callp = callp pup238_6[144:148,] ################################################### ### chunk number 27: getd ################################################### #hg18.txdb = loadFeatures(system.file("extdata/hg18.txdb.sqlite", # package="anno4NGS")) #if (.Platform$OS.type == "windows") # data(ex6) else library(GenomicFeatures) data(ex6) ex6[1:4,] #ex6 = exons(hg18.txdb, vals=list(exon_chrom="chr6")) exp_6 = pup238_6[ which( ranges(pup238_6) %in% ranges(ex6) ) , ] ################################################### ### chunk number 28: lkmsnp ################################################### sum( ranges(pup238_6) %in% ranges(ex6) ) isExonic = which( ranges(pup238_6) %in% ranges(ex6) ) pup238_6x = pup238_6[ isExonic, ] ################################################### ### chunk number 29: lkdb ################################################### library(SNPlocs.Hsapiens.dbSNP.20090506) s6 = getSNPlocs("chr6") s6[1:5,] knownLocs = IRanges(s6$loc, s6$loc) indbsnp = ranges(pup238_6x) %in% knownLocs sum(indbsnp) ################################################### ### chunk number 30: lkhz ################################################### p6xknown = pup238_6x[ which( indbsnp ), ] p6xnovel = pup238_6x[ -which( indbsnp ), ] nKnown = length(p6xknown) nhetKnown = sum( !(elementMetadata(p6xknown)$consensusBase %in% c("A", "C", "G", "T") ) ) nhetKnown/nKnown nNovel = length(p6xnovel) nhetNovel = sum( !(elementMetadata(p6xnovel)$consensusBase %in% c("A", "C", "G", "T") ) ) nhetNovel/nNovel ################################################### ### chunk number 31: lktab ################################################### data(degnerASE01) antab = with(degnerASE01, degnerASE01[ chr=="chr6" & indiv=="GM19238", 1:8] ) antab ################################################### ### chunk number 32: mkta ################################################### tabCalls = function (pup, ind) { ac = as.character empup = elementMetadata(pup) ref = ac(empup$referenceBase[ind]) maqcall = ac(empup$consensusBase[ind]) puptag = gsub("\\^.", "", empup$callp[ind]) list(ref = ref, maqcall = maqcall, calls = table(strsplit(puptag, ""))) } ################################################### ### chunk number 33: zzz ################################################### library(Biostrings) # has IUPAC_CODE_MAP decodeIU = function(x) { if(length(x)>1) warning("only handling first entry") strsplit(IUPAC_CODE_MAP[toupper(x)], "")[[1]] } ################################################### ### chunk number 34: lkra ################################################### ol = findOverlaps(IRanges(antab$loc, width=1), ranges(pup238_6x) ) ol ################################################### ### chunk number 35: lkcomp ################################################### for (i in 1:nrow(ol@matchMatrix)) { # ONEOFF cat("query", i, "\n") print(t1 <- antab[ ol@matchMatrix[i, "query"], ]) t2 <- tabCalls( pup238_6x, ol@matchMatrix[i, "subject"]) print(t2[["calls"]]) dpref = t1[1,2]/(t1[1,2]+t1[1,3]) cat("Degner refFreq =", round(dpref,4),"\n") dpuprefinds = which(toupper(names(t2[["calls"]])) %in% c(",", ".")) targcodes = decodeIU( t2[["maqcall"]] ) altcode = setdiff(targcodes, t2[["ref"]]) dpupaltinds = which(toupper(names(t2[["calls"]])) == altcode) dpuppref = sum(t2[["calls"]][dpuprefinds])/(sum(t2[["calls"]][dpuprefinds])+ sum(t2[["calls"]][dpupaltinds])) cat("Pileup refFreq =", round(dpuppref,4),"\n") cat("---\n") } ################################################### ### chunk number 36: lkd ################################################### library(hmyriB36) if (!exists("hmyriB36")) data(hmyriB36) library(illuminaHumanv1.db) pidVEGFA = get("VEGFA", revmap(illuminaHumanv1SYMBOL)) ve = exprs(hmyriB36)[ pidVEGFA, "NA19238"] ################################################### ### chunk number 37: lkgev ################################################### plot_EvG(genesym("VEGFA"), rsid("rs3025040"), hmyriB36) abline(h=ve) ################################################### ### chunk number 38: dos ################################################### toLatex(sessionInfo())