Skip to content

Instantly share code, notes, and snippets.

@kbroman
Created March 21, 2023 14:38
Show Gist options
  • Save kbroman/378445ee60d4435413e7ca82d356a14e to your computer and use it in GitHub Desktop.
Save kbroman/378445ee60d4435413e7ca82d356a14e to your computer and use it in GitHub Desktop.
smooth genetic map by mixing it with a bit of constant recombination (separate rec rate for each chromosome)
# smooth genetic map by mixing it with a bit of constant recombination (separate rec rate for each chromosome)
# (the input gmap and pmap are lists of vectors of marker positions)
smooth_gmap <-
function(gmap, pmap, alpha=0.02)
{
stopifnot(is.list(gmap), is.list(pmap))
stopifnot(length(gmap)==length(pmap), all(names(gmap) == names(pmap)))
stopifnot(length(alpha)==1, alpha >= 0, alpha <= 1)
for(chr in seq(along=gmap)) {
g <- gmap[[chr]]
p <- pmap[[chr]]
stopifnot(length(g) == length(p), all(names(g) == names(p)) )
gd <- diff(g)
pd <- diff(p)
rr <- diff(range(g)) / diff(range(p))
gd <- gd*(1-alpha) + alpha*pd*rr
gmap[[chr]] <- setNames(cumsum(c(g[1], gd)), names(g))
}
gmap
}
library(qtl2convert)
cox <- read.csv("https://raw.githubusercontent.com/kbroman/CoxMapV3/main/cox_v3_map.csv")
cox$cM_coxV3_ave[cox$chr=="X"] <- cox$cM_coxV3_female[cox$chr=="X"]
cox$Mbp <- cox$bp_grcm39/1e6
cox_gmap <- map_df_to_list(cox, pos_column="cM_coxV3_ave")
cox_pmap <- map_df_to_list(cox, pos_column="Mbp")
cox_gmap_adj <- smooth_gmap(cox_gmap, cox_pmap)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment