# Simulation parameters used in the below examples.
map = uniformMap(M = 10) # recombination map
N = 5 # number of sims
# For more realistic results, replace with e.g.:
# map = loadMap("decode19")
# N = 1000
#################################################################
# EXAMPLE 1
# Comparison of IBD segment distributions
# between paternal and maternal half siblings.
#################################################################
# Define the pedigrees
xPat = halfSibPed()
xMat = swapSex(xPat, 1)
simPat = ibdsim(xPat, N = N, map = map)
simMat = ibdsim(xMat, N = N, map = map)
# By default, the IBD segments of the "leaves" are computed and plotted
plotSegmentDistribution(simPat, simMat, type = "ibd1", ids = 4:5,
labels = c("HSpat", "HSmat"))
#################################################################
# EXAMPLE 2
# Half siblings vs half uncle vs grandparent/grandchild
#################################################################
# Only one pedigree needed here
x = addSon(halfSibPed(), 5)
s = ibdsim(x, N = N, map = map)
# Indicate the pairs explicitly this time.
ids = list(GR = c(2,7), HS = 4:5, HU = c(4,7))
# List names are used as labels in the plot
plotSegmentDistribution(s, type = "ibd1", ids = ids, shape = 1:3)
#################################################################
# EXAMPLE 3
# Comparison of autozygosity distributions in various individuals
# with the same expected inbreeding coefficient (f = 1/8)
#################################################################
G = linearPed(2) |> swapSex(5) |> addSon(c(1,5)) # grandfath/granddaughter
HSpat = halfSibPed(sex2 = 2) |> addSon(4:5) # paternal half sibs
HSmat = swapSex(HSpat, 2) # maternal half sibs
QHFC = quadHalfFirstCousins() # quad half first cousins
QHFC = addSon(QHFC, 9:10)
peds = list(G = G, HSpat = HSpat, HSmat = HSmat, QHFC = QHFC)
plotPedList(peds, newdev = TRUE)
dev.off()
# Simulations
s = lapply(peds, function(p)
ibdsim(p, N = N, ids = leaves(p), verbose = FALSE, map = map))
# Plot distributions
plotSegmentDistribution(s, type = "autoz", title = "Autozygous segments")
Run the code above in your browser using DataLab