# NOT RUN {
# 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(HS = 4:5, HU = c(4,7), GR = c(1,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 = swapSex(linearPed(2), 5) # grandfather/granddaughter
G = addChildren(G, 1, 5, 1)
HSpat = swapSex(halfSibPed(), 5) # paternal half sibs
HSpat = addChildren(HSpat, 4, 5, 1)
HSmat = swapSex(HSpat, 1) # maternal half sibs
QHFC = quadHalfFirstCousins() # quad half first cousins
QHFC = addChildren(QHFC, 9, 10, nch = 1)
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