create_gc_model.r

# Approximate way of creating a GC model from a different gcmodel
# This takes a while to run ...

library("data.table")

# We take Human 1M-Duev3 since it has the most SNPs
gcmodel.other <- data.table(read.table("run_ind_ex/lib/GC_hg19/Human1M-Duov3-0_hg19.gcmodel",h=T,sep="\t"))
setkey(gcmodel.other, "Chr", "Position")

# Read the map file 
lookup_nearest_gc <- function(chr, pos) {
	print(chr)
	print(pos)
	if (chr == "23")
		chr <- "X"
	else if (chr == "24") 
		chr <- "Y"

	gcmodel.other[list(as.character(chr), as.integer(pos)), roll="nearest"]$GC
}

map.core <- read.table("run_ind_ex/lib/MAP_hg19/HumanCoreExome-12v1-0_hg19.map", h=T, sep="\t")

gcmodel.core <- data.frame(
	Name = map.core$V2,
	Chr = map.core$V1,
	Position = map.core$V4,
        GC = mapply(lookup_nearest_gc, map.core$V1, map.core$V4) 
)

write.table(gcmodel.core, "run_ind_ex/lib/GC_hg19/HumanCoreExome-12v1-0_hg19.gcmodel", sep="\t", row.names=F, col.names=T, quot=F)