R/interpolate_cM.R

In jendelman/MapRtools: Tools for genetic mapping

#' Interpolate genetic distances
#' 
#' Creates monotone spline between physical (bp) and genetic (cM) distance
#' 
#' The input map can be generated by merging an existing genetic map with positions in bp and cM with additional markers with only bp information. Interpolation is based on minimizing the mean-squared error between the original and interpolated positions in cM. For the \code{df} and \code{max.extend} parameters, use either a single integer (if same for all chromosomes) or a vector of length equal to the number of chromosomes.
#'
#' @param map Data frame with four columns: marker, chrom, bp, cM 
#' @param df Degrees of freedom for the spline
#' @param max.extend Maximum distance (in cM) to extrapolate beyond the end of the input map
#'
#' @return List containing
#' \describe{
#' \item{map}{Map data frame with additional column named cM.spline}
#' \item{plot}{ggplot2 object}
#' }
#' 
#' @export
#' 
#' @import CVXR
#' @import ggplot2
#' @importFrom splines2 iSpline
#' 
interpolate_cM <- function(map,df=8,max.extend=5){
  stopifnot(ncol(map)==4)
  colnames(map) <- c("marker","chrom","bp","cM")
  map <- map[order(map$chrom,map$bp),]
  map$cM.spline <- numeric(nrow(map))
  chroms <- unique(map$chrom)
  n.chr <- length(chroms)
  if (length(df) >1){
    stopifnot(length(df)==n.chr)
  } else {
    df <- rep(df,n.chr)
  }
  if (length(max.extend) >1){
    stopifnot(length(df)==n.chr)
  } else {
    max.extend <- rep(max.extend,n.chr)
  }
  for (i in 1:n.chr) {
    ix <- which(map$chrom==chroms[i])
    x2 <- iSpline(x=map$bp[ix],intercept=T,df=df[i])
    q <- ncol(x2)
    x2 <- matrix(x2,ncol=q)
    beta <- Variable(q)
    y <- map$cM[ix]
    iq <- which(!is.na(y))
    m <- length(iq)
    max.y <- y[iq[m]] #length of input map
    objective <- Minimize(sum((y[iq]-x2[iq,] %*% beta)^2))
    prob <- Problem(objective,constraints=list(beta >= 0, x2 %*% beta <= max.y + max.extend[i]))
    result <- solve(prob)
    y.pred <- x2%*%result$getValue(beta)
    map$cM.spline[ix] <- round(y.pred,2)
  }
  plotme <- map
  plotme$type <- factor(ifelse(is.na(map$cM),1,0))
  plotme$Mb <- plotme$bp/1e6
  p <- ggplot(data=plotme) + geom_point(data=plotme[plotme$type==0,],aes(x=Mb,y=cM),colour="red") + geom_line(data=plotme[plotme$type==1,],aes(x=Mb,y=cM.spline),colour="black") + facet_wrap(~chrom) + theme_bw()
  return(list(map=map,plot=p))
}