Skip to content

Commit 6e0a956

Browse files
add option to remove alternative isChild
1 parent 92deafe commit 6e0a956

2 files changed

Lines changed: 56 additions & 3 deletions

File tree

R/convertPedigree.R

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ ped2com <- function(ped, component,
3131
standardize.colnames = TRUE,
3232
transpose_method = "tcrossprod",
3333
adjacency_method = "indexed",
34+
isChild_method = "classic",
3435
saveable = FALSE,
3536
resume = FALSE,
3637
save_rate = 5,
@@ -222,14 +223,22 @@ ped2com <- function(ped, component,
222223
isChild <- readRDS(checkpoint_files$isChild)
223224
} else {
224225
# isChild is the 'S' matrix from RAM
226+
227+
if(isChild_method == "partialparent"){
225228
isChild <- apply(ped[, c("momID", "dadID")], 1, function(x) {
226229
.5 + .25*sum(is.na(x)) # 2 parents -> .5, 1 parent -> .75, 0 parents -> 1
227230
})
231+
}else{
232+
isChild <- apply(ped[, c("momID", "dadID")], 1, function(x) {
233+
2^(-!all(is.na(x)))
234+
})
235+
236+
}
228237
if (saveable) {
229238
saveRDS(isChild, file = checkpoint_files$isChild)
230239
}
231-
}
232240

241+
}
233242
# --- Step 2: Compute Relatedness Matrix ---
234243
if (resume && file.exists(checkpoint_files$r_checkpoint) && file.exists(checkpoint_files$gen_checkpoint) && file.exists(checkpoint_files$mtSum_checkpoint) && file.exists(checkpoint_files$newIsPar_checkpoint) &&
235244
file.exists(checkpoint_files$count_checkpoint)
@@ -659,9 +668,10 @@ compute_parent_adjacency <- function(ped, component,
659668
update_rate = update_rate,
660669
parList = parList,
661670
lens=lens,
662-
save_rate_parlist=save_rate_parlist,
671+
save_rate_parlist =save_rate_parlist,
663672
...)
664673
}
674+
665675
} else if (adjacency_method == "indexed") { # Garrison version
666676
if (lastComputed < nr) {
667677
list_of_adjacency <- .adjIndexed(ped=ped,
@@ -702,7 +712,7 @@ compute_parent_adjacency <- function(ped, component,
702712
if (saveable) {
703713
saveRDS(parList, file = checkpoint_files$parList)
704714
saveRDS(lens, file = checkpoint_files$lens)
705-
if (verbose) cat("Final checkpoint saved for adjacency matrix.\n")
715+
if (verbose) {cat("Final checkpoint saved for adjacency matrix.\n")}
706716
}
707717
return(list_of_adjacency)
708718
}

tests/testthat/test-convertPedigree.R

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -320,3 +320,46 @@ test_that("adjacency_method 'indexed', 'loop', and direct produce the same resu
320320
expect_equal(ped_gen_loop, ped_gen_direct, tolerance = tolerance)
321321
expect_equal(ped_gen_indexed, ped_gen_direct, tolerance = tolerance)
322322
})
323+
324+
test_that("isChild_method product the same results for add matrix", {
325+
data(inbreeding)
326+
df <- inbreeding
327+
df$momID[df$ID == 6] <- NA
328+
tolerance <- 1e-10
329+
# add
330+
ped_add_partial <- ped2com(df, isChild_method= "partialparent",
331+
component = "additive",
332+
adjacency_method = "direct")
333+
ped_add_classic <- ped2com(df, isChild_method= "classic",
334+
component = "additive", adjacency_method = "direct")
335+
336+
expect_equal(ped_add_partial[6,6], 1, tolerance = tolerance)
337+
expect_equal(ped_add_classic[6,6], .75, tolerance = tolerance)
338+
difference <- ped_add_partial - ped_add_classic
339+
# expect_equal(ped_add_partial, ped_add_classic, tolerance = tolerance)
340+
341+
difference <- ped_add_partial - ped_add_classic
342+
343+
expect_gt(sum(abs(difference)),0)
344+
})
345+
346+
347+
348+
test_that("isChild_method product the same results for mtdna matrix", {
349+
data(hazard)
350+
df <- hazard
351+
df$momID[df$ID == 4] <- NA
352+
tolerance <- 1e-10
353+
# maternal
354+
ped_mit_partial <- ped2com(df, isChild_method= "partialparent",
355+
component = "mitochondrial",
356+
adjacency_method = "indexed")
357+
ped_mit_classic <- ped2com(df, isChild_method= "classic",
358+
component = "mitochondrial", adjacency_method = "indexed")
359+
# should be the same
360+
expect_equal(ped_mit_partial, ped_mit_classic, tolerance = tolerance)
361+
362+
})
363+
364+
365+

0 commit comments

Comments
 (0)