sim_do_kinship.R

# simulate DO data, to look at kinship coefficient for siblings vs others


library(simcross) # install with devtools::install_github("kbroman/simcross")
library(qtl2)


# example DO data with megaMUGA genotypes, used for map + founder genotypes
file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/master/DO_Svenson291/svenson.zip")
do_mm <- read_cross2(file)

# drop markers where there are missing founder genotypes
fg <- do.call("cbind", do_mm$founder_geno)
do_mm <- drop_markers(do_mm, colnames(fg)[colSums(fg==0) > 0])

# simulate DO genotypes, without X chromosome
do_ped <- sim_do_pedigree(nkids_per=2)
do_gen <- sim_from_pedigree(do_ped, sapply(do_mm$gmap, max)[-20])
do_gen2 <- collapse_do_alleles(do_gen)
do_gen2 <- lapply(do_gen2, function(g) g[do_ped$do & do_ped$gen==12])

# create SNP genotypes
gen <- lapply(seq_along(do_gen2), function(chr) convert2geno(do_gen2[[chr]], do_mm$gmap[[chr]],
                                                             founder_geno=do_mm$founder_geno[[chr]]))
names(gen) <- names(do_gen2)

# paste into cross object
do_mm$geno <- gen
do_mm$founder_geno <- do_mm$founder_geno[-20]
do_mm$covar <- NULL
do_mm$is_female <- setNames(rep(TRUE, 288), rownames(do_mm$geno[[1]]))
do_mm$cross_info <- as.matrix(rep(12L, 288))
rownames(do_mm$cross_info) <- rownames(do_mm$geno[[1]])
do_mm$gmap <- do_mm$gmap[-20]
do_mm$pmap <- do_mm$pmap[-20]
do_mm$is_x_chr <- setNames(rep(FALSE, 19), names(do_mm$gmap))

# now look at kinship matrix
pr <- calc_genoprob(do_mm, cores=0) # <- really slow
k <- calc_kinship(pr, cores=0)

# heat map in gray scale
image(k, col=rev(gray((0:256)/256)))

# histogram
hist(k, breaks=50)

# kinship coefficients for first couple of sib-pairs
k[1,2]
k[3,4]