vpnagraj
commented
3 years ago @stephenturner some progress here. the helpers ive been working with seem to do the trick. long thread with all the code ...
first define the two helper functions (eventually to live in skater?):
## helper function to interpolate over segments
## inspired by py code http://faculty.washington.edu/sguy/ibd_relatedness.html
interpolate <- function(ibd_bp, chromgpos) {
tmp_set1 <-
chromgpos %>%
dplyr::filter(bp < ibd_bp)
tmp_set2 <-
chromgpos %>%
dplyr::filter(bp > ibd_bp)
if(nrow(tmp_set1) > 0 && nrow(tmp_set2) > 0) {
tmp_bp1 <- max(tmp_set1$bp)
tmp_bp2 <- min(tmp_set2$bp)
} else if (nrow(tmp_set1) > 0) {
tmp_bp2 <- max(tmp_set1$bp)
tmp_set1 <-
tmp_set1 %>%
dplyr::filter(bp != tmp_bp2)
tmp_bp1 <- max(tmp_set1$bp)
} else {
tmp_bp1 <- min(tmp_set2$bp)
tmp_set2 <-
tmp_set2 %>%
dplyr::filter(bp != tmp_bp1)
tmp_bp2 <- min(tmp_set2$bp)
}
if (tmp_bp1 == tmp_bp2) {
res <-
chromgpos %>%
filter(bp == tmp_bp1) %>%
pull(value)
return(res)
} else {
len1 <-
chromgpos %>%
filter(bp == tmp_bp1) %>%
pull(value)
len2 <-
chromgpos %>%
filter(bp == tmp_bp2) %>%
pull(value)
res <- len1+(len2-len1)*(ibd_bp - tmp_bp1)/(tmp_bp2-tmp_bp1)
return(res)
}
}
## function to loop over each chromosome and compute total length
## use 4 times that total length as denominator for pairwise shared length to get kinship coefficient
ibd2kin <- function(ibd_data, plinkmap) {
## get the min base pair for each chr
ibd_min_bp <-
ibd_data %>%
group_by(chr) %>%
summarise(bp = min(start), .groups = "drop")
## get the max base pair for each chromsome
ibd_max_bp <-
ibd_data %>%
group_by(chr) %>%
summarise(bp = max(end), .groups = "drop")
## loop over each chromosome to get length
## this will be summed up and serve as denominator for kinship coefficient calculations
totalchromlength <- vector()
for(i in 1:length(unique(ibd_min_bp$chr))) {
chrom <- unique(ibd_min_bp$chr)[i]
## these steps just get the min/max bp and map data for the given chromosome
tmp_ibd_min_bp <-
ibd_min_bp %>%
filter(chr == chrom)
tmp_ibd_max_bp <-
ibd_max_bp %>%
filter(chr == chrom)
tmp_chromgpos <-
plinkmap %>%
filter(chr == chrom)
## use interpolate to get centimorgan lengths
mincm <- interpolate(ibd_bp = tmp_ibd_min_bp$bp, chromgpos = tmp_chromgpos)
maxcm <- interpolate(ibd_bp = tmp_ibd_max_bp$bp, chromgpos = tmp_chromgpos)
## create chrom length and append to vector
tmp_totalchromlength <- maxcm - mincm
totalchromlength <- c(totalchromlength, tmp_totalchromlength)
}
## use ibd data and total chromlength as denominator above to get the kinship value
ibd_data %>%
group_by(id1,id2) %>%
summarise(totallength = sum(length), .groups = "drop") %>%
mutate(kinship = totallength/(4*sum(totalchromlength)))
}now using the same data you shared above calculate the kinship coefficient from the IBD segments and compare to theoretical kinship from "truth" data:
## read in fam
truth <- read_fam("5gens-everyone.fam.gz") %>%
fam2ped() %>%
mutate(x=map(ped, ped2kinpair)) %>%
select(x) %>%
unnest(x) %>%
mutate(d.truth=kin2degree(k, max_degree=8)) %>%
rename(k.truth=k) %>%
print()
## read in seg and convert to "ibd" output (i.e. similar to hap-ibd output)
seg <- read_tsv("5gens.seg.gz",
col_types="ccciicddd",
col_names=c("id1", "id2", "chr", "pstart", "pend", "type", "gstart", "gend", "glength"))
seglikeibd <-
seg %>%
filter(type == "IBD1" | type == "IBD2") %>%
select(id1, id2, chr, start = pstart, end = pend, length = glength)
## read in map file
plinkmap <-
read_delim("/data/projects/skate/data/hapmap/plink.allchr.GRCh37.map", delim = "\t", col_names = FALSE) %>%
select(chr = 1, value = 3, bp = 4)
## join up with truth data and see how we did
joined <-
ibd2kin(ibd_data = seglikeibd, plinkmap = plinkmap) %>%
mutate(degree = kin2degree(kinship)) %>%
left_join(truth) %>%
mutate(across(starts_with("d"), ~factor(replace_na(., "Unrelated"))))
## look at first 10 obs
joined %>%
head(10) %>%
knitr::kable()| id1 | id2 | totallength | kinship | degree | k.truth | d.truth |
|---|---|---|---|---|---|---|
| 5gens1_g1-b1-i1 | 5gens1_g2-b1-i1 | 3346.2989 | 0.2363204 | 1 | 0.2500 | 1 |
| 5gens1_g1-b1-i1 | 5gens1_g2-b1-i2 | 3346.2989 | 0.2363204 | 1 | 0.2500 | 1 |
| 5gens1_g1-b1-i1 | 5gens1_g2-b2-i1 | 3346.2989 | 0.2363204 | 1 | 0.2500 | 1 |
| 5gens1_g1-b1-i1 | 5gens1_g2-b2-i2 | 3346.2989 | 0.2363204 | 1 | 0.2500 | 1 |
| 5gens1_g1-b1-i1 | 5gens1_g3-b1-i1 | 1823.3637 | 0.1287685 | 2 | 0.1250 | 2 |
| 5gens1_g1-b1-i1 | 5gens1_g3-b1-i2 | 1846.1699 | 0.1303791 | 2 | 0.1250 | 2 |
| 5gens1_g1-b1-i1 | 5gens1_g3-b2-i1 | 1815.5866 | 0.1282193 | 2 | 0.1250 | 2 |
| 5gens1_g1-b1-i1 | 5gens1_g3-b2-i2 | 1744.5802 | 0.1232047 | 2 | 0.1250 | 2 |
| 5gens1_g1-b1-i1 | 5gens1_g4-b1-i1 | 887.1510 | 0.0626519 | 3 | 0.0625 | 3 |
| 5gens1_g1-b1-i1 | 5gens1_g4-b1-i2 | 876.6621 | 0.0619111 | 3 | 0.0625 | 3 |
questions
- what else do we need to do to confirm that this is "working"?
- is the "IBD1 or IBD2" logic appropriate to add up segments from the ped-sim .seg file (cc @genignored ) ? what about the HBD segments that are in the .seg file ... ignore those?
- when/how do we want to pull the helpers into
skater?
stephenturner
genignored
During the call A.W. brought up an idea that sounds reasonable -- comparing the calculated relatedness (kinship) against the actual simulated kinship based on the ibd segments ped-sim simulates, not just the expected kinship from the fam file. I think this can be done with the .seg file output by ped-sim. This is documented here: https://github.com/williamslab/ped-sim#output-ibd-segments-file.
I'm attaching a .fam file and a .seg file (compressed) from a simple 5 generation pedigree (not simulating any genetic data) using a pedigree definition that looks like this:
With a simple ped-sim simulation (that map and XO interference file are in /data/skate/skatesim):
This is the pedigree it produces:
You can make that plot and then count the number of degree relationships with a little skater code:
You can also read in that .seg file pretty easily as well. I started summing up the IBD1/2 segments, but not sure this makes sense.
@vpnagraj I'd like to assign this to you to at least take a look at. I think you did something very similar for getting k from hap-ibd, over at https://github.com/signaturescience/skatebench/blob/a6d6c6dc64c05639bfbd42fe9b9cadcccf529d3b/hap-ibd-test/hap-ibd-relatedness.R. If you have any ideas about how to get something working push up to a new branch, and maybe include the .seg and .fam in inst/extdata to use as examples. cc @genignored