# map<-read.table("http://orpheus/Mapview/master_map.dat", h=T, na=".")

map<-read.table("newmap20060618.dat", h=T, na=".")
pseudo <- read.table("http://orpheus/Mapview/pseudoautosomal_map.dat", h=T, na=".")


#
# Predict Rutgers cM
#

interp.cM <- function(chrom) {
  require(locfit)
  chromdata<-map[map$chr==chrom & !is.na(map$I.bp),]
  if (nrow(chromdata)==0) {
    stop(paste("No usable physical map positions chrom",chrom))
  }
  plot(m<-locfit(R.cM ~ I.bp, data=chromdata, alpha=0.15), get.data=TRUE,
       main=paste("Chromosome",chrom))
  p<-predict(m, newdata=data.frame(I.bp=chromdata$I.bp))
  points(chromdata$I.bp, p, pch=18)
#
# chromdata$I.cM<-p-min(p, na.rm=TRUE) <= no good for acrocentric chroms
#
  chromdata$I.cM<-p
  chromdata
}
interp.cM.X <- function() {
  require(locfit)
  chromdata<-map[map$chr==23 & !is.na(map$I.bp),]
  if (nrow(chromdata)==0) {
    stop(paste("No usable physical map positions chrom",chrom))
  }
  plot(m<-locfit(R.cM.F ~ I.bp, data=chromdata, alpha=0.05), get.data=TRUE,
       main=paste("Chromosome X"))
  p<-predict(m, newdata=data.frame(I.bp=chromdata$I.bp))
  points(chromdata$I.bp, p, pch=18)
  chromdata$I.cM<-p-min(p, na.rm=TRUE)
  chromdata
}

interp.all23 <- function() {
  newmap<-data.frame(NULL)
  for(i in 1:22) {
    newmap<-rbind(newmap, interp.cM(i))
  } 
  newmap<-rbind(newmap, interp.cM.X())
  new<-rep(NA,nrow(map))
  new[as.integer(rownames(newmap))]<-newmap$I.cM
  new[new<0]<-0
  map$I.cM<-new
  idx<-is.na(map$I.cM) | (!is.na(map$R.cM) & is.na(map$build34))
  map$I.cM[idx]<-map$R.cM[idx]
  map<-map[order(map$chr, map$I.cM),]
  map$I.bp<-round(map$I.bp)
  map$I.cM<-round(map$I.cM,3)
  write.table(map,"newmap.table", quote=FALSE, na=".", row.names=FALSE)
  map
}
#
# Predict male PAR cM (Lien et al 2000)
#

interp.cM.Male <- function() {
  require(locfit)
  chromdata<-pseudo[!is.na(pseudo$I.bp),]
  if (nrow(chromdata)==0) {
    stop(paste("No usable physical pseudo positions chrom",chrom))
  }
  plot(m<-locfit(lien.cM ~ I.bp, data=chromdata, subset=(pseudo$I.bp<2750000)), get.data=TRUE,
       main="Pseudoautosomal Region 1-2")
  p<-predict(m, newdata=data.frame(I.bp=chromdata$I.bp))
  pab <- which(chromdata$name1=="PABX")
  max.par1 <- p[pab]
  p[p<0]<-0
  p[seq(pab+1,length(p))]<-max.par1
  points(chromdata$I.bp, p, pch=18)
  chromdata$I.cM.M<-p
  chromdata$I.cM.M[chromdata$name=="DXS1107"]<-chromdata$lien.cM[chromdata$name=="DXS1107"]
  chromdata$I.cM.M<-round(chromdata$I.cM.M,3)
  write.table(chromdata,"newpar.table", quote=FALSE, na=".", row.names=FALSE)
  chromdata
}


map2 <- interp.all23()
pseudo$I.cM.F <- map2$I.cM[match(pseudo$name1, map2$name1)]
pseudo2 <- interp.cM.Male()