This document will show the relationship between the Haldane and
Kosambi mapping functions with respect to recombination rate and genetic
distance. It will then show how each function is represented using the
gamma model for recombination. Finally, it will show how recombination
parameters in AlphaSimR.
Haldane function vs Kosambi function
# Haldane's function to convert genetic distance (Morgans) to
# recombination rate.
calcRate_Haldane = function(d){
(1 - exp(-2*d)) / 2
}
# Haldane's function to convert recombination rate to
# genetic distance (Morgans).
calcDist_Haldane = function(r){
-log(1 - 2*r) / 2
}
# Kosambi's function to convert genetic distance (Morgans) to
# recombination rate.
calcRate_Kosambi = function(d){
(1 - exp(-4*d)) / (2 * (1 + exp(-4*d)))
}
# Kosambi's function to convert recombination rate to
# genetic distance (Morgans).
calcDist_Kosambi = function(r){
log( (1 + 2*r) / (1 - 2*r) ) / 4
}
distance = seq(0, 1, length.out=100)
rateHaldane = calcRate_Haldane(distance)
rateKosambi = calcRate_Kosambi(distance)
plot(distance, rateKosambi, type="l",
xlab="Genetic Distance (Morgans)",
ylab="Recombination Rate",
col="red")
lines(distance, rateHaldane)

rate = seq(0, 0.5, length.out=100)
distHaldane = calcDist_Haldane(rate)
distKosambi = calcDist_Kosambi(rate)
plot(rate, distHaldane, type="l",
xlab="Recombination Rate",
ylab="Genetic Distance (Morgans)"
)
lines(rate, distKosambi, col="red")

Gamma model for recombination
# Function to sample recombination positions for a
# chromosome with 1 Morgan genetic length using the
# gamma model with a variable interference parameter (v).
sampleRecomb = function(v){
# Sample deviates from a gamma distribution
x = rgamma(n=50, shape=v, rate=2*v)
# Convert deviates to map positions
# Subtracting by 10 to produce a burn-in to avoid too
# few recombinations at the beginning of the chromosome
pos = cumsum(x) - 10
# Trimming positions to only those within the chromosome
# Best practice would be to check that at least one position
# is greater than 1 and sample more deviates if not.
pos = pos[pos>0 & pos<1]
# Thin out positions by one half
# This accounts for recombinations between non-target chromatids
keep = rbinom(n=length(pos), size=1, prob=0.5)
pos = pos[as.logical(keep)]
# Return the positions
return(pos)
}
# Look at distribution of number of recombinations under
# Haldane and approximate Kosambi assumptions.
nHaldaneCO = nKosambiCO = numeric(10000)
for(i in 1:10000){
nHaldaneCO[i] = length( sampleRecomb(v=1) )
nKosambiCO[i] = length( sampleRecomb(v=2.6) )
}
mean(nHaldaneCO)
[1] 0.9914
mean(nKosambiCO)
[1] 0.9957
hist(nHaldaneCO, breaks = c(-1,0:max(nHaldaneCO))+0.5)

hist(nKosambiCO, breaks = c(-1,0:max(nKosambiCO))+0.5)

# Show distance between crossovers
distHaldaneCO = distKosambiCO = vector("list", length=10000)
for(i in 1:10000){
distHaldaneCO[[i]] = diff( sampleRecomb(v=1) )
distKosambiCO[[i]] = diff( sampleRecomb(v=2.6) )
}
distHaldaneCO = do.call("c", distHaldaneCO)
distKosambiCO = do.call("c", distKosambiCO)
hist(distHaldaneCO)

hist(distKosambiCO)

AlphaSimR
library(AlphaSimR)
Loading required package: R6
founderPop = quickHaplo(nInd=1, nChr=1, segSites=11, genLen=1)
SP = SimParam$new(founderPop)
# Interference parameter
SP$v
[1] 2.6
# Proportion of non-interferring chiasmata
SP$p
[1] 0
# Genetic map
getGenMap()
# Turn on recombination tracking option
SP$setTrackRec(TRUE)
Parent = newPop(founderPop)
A = self(Parent, nProgeny=5)
# Pulling IBD haplotypes
hap = pullIbdHaplo(A)
hap
1_1 1_2 1_3 1_4 1_5 1_6 1_7 1_8 1_9 1_10 1_11
2_1 2 1 1 1 2 2 2 2 2 2 2
2_2 2 1 1 1 2 2 2 2 2 2 2
3_1 2 2 2 2 2 2 2 1 1 1 1
3_2 2 2 1 1 1 1 1 2 2 2 2
4_1 1 2 2 1 1 1 1 1 1 1 1
4_2 2 2 2 2 2 1 1 1 1 1 2
5_1 1 1 1 1 1 1 1 2 2 2 2
5_2 1 1 1 1 2 2 2 2 2 2 2
6_1 1 1 1 1 1 2 2 1 1 1 1
6_2 1 1 1 2 2 2 2 2 2 2 2
# Function to transpose and flip a matrix for visualizing with
# the `image` function. See ?image
tf = function(m) t(m)[, nrow(m):1]
image(tf(hap),
col = rainbow(2),
useRaster = TRUE,
axes = FALSE,
xlab = "Locus",
ylab = "Haplotype")

LS0tDQp0aXRsZTogIkdlbmV0aWMgTWFwcyBhbmQgUmVjb21iaW5hdGlvbiINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNClRoaXMgZG9jdW1lbnQgd2lsbCBzaG93IHRoZSByZWxhdGlvbnNoaXAgYmV0d2VlbiB0aGUgSGFsZGFuZSBhbmQgS29zYW1iaSBtYXBwaW5nIGZ1bmN0aW9ucyB3aXRoIHJlc3BlY3QgdG8gcmVjb21iaW5hdGlvbiByYXRlIGFuZCBnZW5ldGljIGRpc3RhbmNlLiBJdCB3aWxsIHRoZW4gc2hvdyBob3cgZWFjaCBmdW5jdGlvbiBpcyByZXByZXNlbnRlZCB1c2luZyB0aGUgZ2FtbWEgbW9kZWwgZm9yIHJlY29tYmluYXRpb24uIEZpbmFsbHksIGl0IHdpbGwgc2hvdyBob3cgcmVjb21iaW5hdGlvbiBwYXJhbWV0ZXJzIGluIEFscGhhU2ltUi4NCg0KIyMgSGFsZGFuZSBmdW5jdGlvbiB2cyBLb3NhbWJpIGZ1bmN0aW9uDQoNCmBgYHtyfQ0KIyBIYWxkYW5lJ3MgZnVuY3Rpb24gdG8gY29udmVydCBnZW5ldGljIGRpc3RhbmNlIChNb3JnYW5zKSB0byANCiMgcmVjb21iaW5hdGlvbiByYXRlLg0KY2FsY1JhdGVfSGFsZGFuZSA9IGZ1bmN0aW9uKGQpew0KICAoMSAtIGV4cCgtMipkKSkgLyAyDQp9DQoNCiMgSGFsZGFuZSdzIGZ1bmN0aW9uIHRvIGNvbnZlcnQgcmVjb21iaW5hdGlvbiByYXRlIHRvIA0KIyBnZW5ldGljIGRpc3RhbmNlIChNb3JnYW5zKS4NCmNhbGNEaXN0X0hhbGRhbmUgPSBmdW5jdGlvbihyKXsNCiAgLWxvZygxIC0gMipyKSAvIDINCn0NCg0KIyBLb3NhbWJpJ3MgZnVuY3Rpb24gdG8gY29udmVydCBnZW5ldGljIGRpc3RhbmNlIChNb3JnYW5zKSB0byANCiMgcmVjb21iaW5hdGlvbiByYXRlLg0KY2FsY1JhdGVfS29zYW1iaSA9IGZ1bmN0aW9uKGQpew0KICAoMSAtIGV4cCgtNCpkKSkgLyAoMiAqICgxICsgZXhwKC00KmQpKSkNCn0NCg0KIyBLb3NhbWJpJ3MgZnVuY3Rpb24gdG8gY29udmVydCByZWNvbWJpbmF0aW9uIHJhdGUgdG8gDQojIGdlbmV0aWMgZGlzdGFuY2UgKE1vcmdhbnMpLg0KY2FsY0Rpc3RfS29zYW1iaSA9IGZ1bmN0aW9uKHIpew0KICBsb2coICgxICsgMipyKSAvICgxIC0gMipyKSApIC8gNA0KfQ0KDQoNCmRpc3RhbmNlID0gc2VxKDAsIDEsIGxlbmd0aC5vdXQ9MTAwKQ0KcmF0ZUhhbGRhbmUgPSBjYWxjUmF0ZV9IYWxkYW5lKGRpc3RhbmNlKQ0KcmF0ZUtvc2FtYmkgPSBjYWxjUmF0ZV9Lb3NhbWJpKGRpc3RhbmNlKQ0KDQpwbG90KGRpc3RhbmNlLCByYXRlS29zYW1iaSwgdHlwZT0ibCIsDQogICAgIHhsYWI9IkdlbmV0aWMgRGlzdGFuY2UgKE1vcmdhbnMpIiwNCiAgICAgeWxhYj0iUmVjb21iaW5hdGlvbiBSYXRlIiwgDQogICAgIGNvbD0icmVkIikNCmxpbmVzKGRpc3RhbmNlLCByYXRlSGFsZGFuZSkNCg0KDQpyYXRlID0gc2VxKDAsIDAuNSwgbGVuZ3RoLm91dD0xMDApDQpkaXN0SGFsZGFuZSA9IGNhbGNEaXN0X0hhbGRhbmUocmF0ZSkNCmRpc3RLb3NhbWJpID0gY2FsY0Rpc3RfS29zYW1iaShyYXRlKQ0KDQpwbG90KHJhdGUsIGRpc3RIYWxkYW5lLCB0eXBlPSJsIiwNCiAgICAgeGxhYj0iUmVjb21iaW5hdGlvbiBSYXRlIiwNCiAgICAgeWxhYj0iR2VuZXRpYyBEaXN0YW5jZSAoTW9yZ2FucykiDQogICAgICkNCmxpbmVzKHJhdGUsIGRpc3RLb3NhbWJpLCBjb2w9InJlZCIpDQoNCmBgYA0KDQojIyBHYW1tYSBtb2RlbCBmb3IgcmVjb21iaW5hdGlvbg0KDQpgYGB7cn0NCiMgRnVuY3Rpb24gdG8gc2FtcGxlIHJlY29tYmluYXRpb24gcG9zaXRpb25zIGZvciBhIA0KIyBjaHJvbW9zb21lIHdpdGggMSBNb3JnYW4gZ2VuZXRpYyBsZW5ndGggdXNpbmcgdGhlIA0KIyBnYW1tYSBtb2RlbCB3aXRoIGEgdmFyaWFibGUgaW50ZXJmZXJlbmNlIHBhcmFtZXRlciAodikuDQpzYW1wbGVSZWNvbWIgPSBmdW5jdGlvbih2KXsNCiAgIyBTYW1wbGUgZGV2aWF0ZXMgZnJvbSBhIGdhbW1hIGRpc3RyaWJ1dGlvbg0KICB4ID0gcmdhbW1hKG49NTAsIHNoYXBlPXYsIHJhdGU9Mip2KQ0KICANCiAgIyBDb252ZXJ0IGRldmlhdGVzIHRvIG1hcCBwb3NpdGlvbnMNCiAgIyBTdWJ0cmFjdGluZyBieSAxMCB0byBwcm9kdWNlIGEgYnVybi1pbiB0byBhdm9pZCB0b28gDQogICMgZmV3IHJlY29tYmluYXRpb25zIGF0IHRoZSBiZWdpbm5pbmcgb2YgdGhlIGNocm9tb3NvbWUNCiAgcG9zID0gY3Vtc3VtKHgpIC0gMTANCiAgDQogICMgVHJpbW1pbmcgcG9zaXRpb25zIHRvIG9ubHkgdGhvc2Ugd2l0aGluIHRoZSBjaHJvbW9zb21lDQogICMgQmVzdCBwcmFjdGljZSB3b3VsZCBiZSB0byBjaGVjayB0aGF0IGF0IGxlYXN0IG9uZSBwb3NpdGlvbg0KICAjIGlzIGdyZWF0ZXIgdGhhbiAxIGFuZCBzYW1wbGUgbW9yZSBkZXZpYXRlcyBpZiBub3QuDQogIHBvcyA9IHBvc1twb3M+MCAmIHBvczwxXQ0KICANCiAgIyBUaGluIG91dCBwb3NpdGlvbnMgYnkgb25lIGhhbGYNCiAgIyBUaGlzIGFjY291bnRzIGZvciByZWNvbWJpbmF0aW9ucyBiZXR3ZWVuIG5vbi10YXJnZXQgY2hyb21hdGlkcw0KICBrZWVwID0gcmJpbm9tKG49bGVuZ3RoKHBvcyksIHNpemU9MSwgcHJvYj0wLjUpDQogIHBvcyA9IHBvc1thcy5sb2dpY2FsKGtlZXApXQ0KICANCiAgIyBSZXR1cm4gdGhlIHBvc2l0aW9ucw0KICByZXR1cm4ocG9zKQ0KfQ0KDQojIExvb2sgYXQgZGlzdHJpYnV0aW9uIG9mIG51bWJlciBvZiByZWNvbWJpbmF0aW9ucyB1bmRlciANCiMgSGFsZGFuZSBhbmQgYXBwcm94aW1hdGUgS29zYW1iaSBhc3N1bXB0aW9ucy4NCm5IYWxkYW5lQ08gPSBuS29zYW1iaUNPID0gbnVtZXJpYygxMDAwMCkNCmZvcihpIGluIDE6MTAwMDApew0KICBuSGFsZGFuZUNPW2ldID0gbGVuZ3RoKCBzYW1wbGVSZWNvbWIodj0xKSApDQogIG5Lb3NhbWJpQ09baV0gPSBsZW5ndGgoIHNhbXBsZVJlY29tYih2PTIuNikgKQ0KfQ0KDQptZWFuKG5IYWxkYW5lQ08pDQptZWFuKG5Lb3NhbWJpQ08pDQoNCmhpc3QobkhhbGRhbmVDTywgYnJlYWtzID0gYygtMSwwOm1heChuSGFsZGFuZUNPKSkrMC41KQ0KaGlzdChuS29zYW1iaUNPLCBicmVha3MgPSBjKC0xLDA6bWF4KG5Lb3NhbWJpQ08pKSswLjUpDQoNCiMgU2hvdyBkaXN0YW5jZSBiZXR3ZWVuIGNyb3Nzb3ZlcnMNCmRpc3RIYWxkYW5lQ08gPSBkaXN0S29zYW1iaUNPID0gdmVjdG9yKCJsaXN0IiwgbGVuZ3RoPTEwMDAwKQ0KZm9yKGkgaW4gMToxMDAwMCl7DQogIGRpc3RIYWxkYW5lQ09bW2ldXSA9IGRpZmYoIHNhbXBsZVJlY29tYih2PTEpICkNCiAgZGlzdEtvc2FtYmlDT1tbaV1dID0gZGlmZiggc2FtcGxlUmVjb21iKHY9Mi42KSApDQp9DQoNCmRpc3RIYWxkYW5lQ08gPSBkby5jYWxsKCJjIiwgZGlzdEhhbGRhbmVDTykNCmRpc3RLb3NhbWJpQ08gPSBkby5jYWxsKCJjIiwgZGlzdEtvc2FtYmlDTykNCg0KaGlzdChkaXN0SGFsZGFuZUNPKQ0KaGlzdChkaXN0S29zYW1iaUNPKQ0KDQpgYGANCg0KIyMgQWxwaGFTaW1SDQoNCmBgYHtyfQ0KbGlicmFyeShBbHBoYVNpbVIpDQoNCmZvdW5kZXJQb3AgPSBxdWlja0hhcGxvKG5JbmQ9MSwgbkNocj0xLCBzZWdTaXRlcz0xMSwgZ2VuTGVuPTEpDQoNClNQID0gU2ltUGFyYW0kbmV3KGZvdW5kZXJQb3ApDQoNCiMgSW50ZXJmZXJlbmNlIHBhcmFtZXRlcg0KU1Akdg0KDQojIFByb3BvcnRpb24gb2Ygbm9uLWludGVyZmVycmluZyBjaGlhc21hdGENClNQJHANCg0KIyBHZW5ldGljIG1hcA0KZ2V0R2VuTWFwKCkNCg0KIyBUdXJuIG9uIHJlY29tYmluYXRpb24gdHJhY2tpbmcgb3B0aW9uDQpTUCRzZXRUcmFja1JlYyhUUlVFKQ0KDQpQYXJlbnQgPSBuZXdQb3AoZm91bmRlclBvcCkNCg0KQSA9IHNlbGYoUGFyZW50LCBuUHJvZ2VueT01KQ0KDQojIFB1bGxpbmcgSUJEIGhhcGxvdHlwZXMNCmhhcCA9IHB1bGxJYmRIYXBsbyhBKQ0KaGFwDQoNCiMgRnVuY3Rpb24gdG8gdHJhbnNwb3NlIGFuZCBmbGlwIGEgbWF0cml4IGZvciB2aXN1YWxpemluZyB3aXRoDQojIHRoZSBgaW1hZ2VgIGZ1bmN0aW9uLiBTZWUgP2ltYWdlDQp0ZiA9IGZ1bmN0aW9uKG0pIHQobSlbLCBucm93KG0pOjFdDQppbWFnZSh0ZihoYXApLCANCiAgICAgIGNvbCA9IHJhaW5ib3coMiksIA0KICAgICAgdXNlUmFzdGVyID0gVFJVRSwNCiAgICAgIGF4ZXMgPSBGQUxTRSwNCiAgICAgIHhsYWIgPSAiTG9jdXMiLA0KICAgICAgeWxhYiA9ICJIYXBsb3R5cGUiKQ0KDQpgYGANCg0K