
library(ceu1kg)
data(ceu1kgMeta_20)
geneLimits = function( anno="org.Hs.eg.db", chr="20" ) {
 require(anno, character.only=TRUE)
 clnanno = gsub(".db", "", anno)
 gna = get(chr, revmap( get(paste(clnanno, "CHR", sep="")) ) )
 gs = mget(gna, get(paste(clnanno, "CHRLOC", sep="")) ) 
 ge = mget(gna, get(paste(clnanno, "CHRLOCEND", sep="")) ) 
 gs = abs(sapply(gs, "[", 1))
 ge = abs(sapply(ge, "[", 1))
 bad = which(is.na(gs) | is.na(ge))
 ans = GRanges( IRanges(gs[-bad], ge[-bad]), seqnames=paste("chr", chr, sep=""))
 names(ans) = gna[-bad]
 ans
}
gl = geneLimits( anno="illuminaHumanv1.db" )

extendGR = function( gr, siz=1e6 )
  punion(shift(gr, -siz), shift(gr, siz), fill.gap=TRUE)

egl = extendGR(gl)
save(egl, file="egl.rda")

containedSnps = function(geneRanges, snpRanges) {
 #
 # returns a list with one element per gene in geneRanges
 # giving the names of the SNP that lie within the gene's range
 #
  mm = findOverlaps( geneRanges, snpRanges )@matchMatrix
  sinds = split( mm[,2], mm[,1] )  # split snp indexes by gene index
  ugi = unique(mm[,1])
  gn = names(geneRanges)[ugi]
  ans = lapply( 1:length(sinds), function(x) names(snpRanges)[sinds[[x]] ] )
  names(ans) = gn
  ans
}

cisScores2 = function( mgr, geneRanges, snpRanges, applier=lapply ) {
 snm = unlist(lapply(mgr@fflist, rownames))
 snpRanges = snpRanges[ intersect(snm, names(snpRanges)) ]
 snps = containedSnps( geneRanges, snpRanges )
 gn = names(geneRanges)
 applier( 1:length(snps), function(x) mgr[ rsid(snps[[x]]), probeId(gn[x]) ])
}
 
library(GGtools)
data(hmceuB36.2021)
gl20 = names(egl[ which(seqnames(egl) == "chr20") ] )[1:20]
p2 = intersect( gl20, featureNames(hmceuB36.2021) )
system("rm -rf foo")
ee = eqtlTests( hmceuB36.2021[ probeId(p2), ], ~1 )
save(ee, file="ee.rda")
cisScores2( ee, egl[1:10], ceu1kgMeta_20 )
