Skip to content

Commit 9f61013

Browse files
smasongarrisonsmasongarrison
committed
Support optional mtdna and mit_val=NULL in slicing
Add mtdna parameter to sliceFamilies (default TRUE) and allow .write_bin_data to accept mit_val = NULL to produce bin files that are not split by mitRel. Update sliceFamilies to call .write_bin_data with mit_val=NULL when mtdna is FALSE and add verbose messages. Adjust .write_bin_data to build file paths for mit_val present or NULL and only write when rows exist. Add tests to verify mitRel=0/1 and mit_val=NULL output files and their contents. Also micro-optimizations: replace repeated ifelse(u %in% v$i, ...) patterns in makeLinks with match-based indexing to avoid %in% overhead and NA issues, and speed up readGedcom's pattern counting by using vapply over file$X1 with fixed = TRUE. Update NEWS.md with a note about the optimized gedcom reader and sliceFamilies change.
1 parent e70790e commit 9f61013

7 files changed

Lines changed: 99 additions & 16 deletions

File tree

‎NEWS.md‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,9 @@
22

33
# Development version:
44
## BGmisc 1.7.0.1.1
5+
* Optimized gedcom reader, com2links for speed and memory usage, with a focus on large pedigrees
56
* Fixed bug in gedcom reader that resulted in document records being added to the final person in the pedigree
6-
* Optimized sliceFamilies to be more abstract
7+
* Optimized sliceFamilies to be more abstract, and no longer require mtdna
78
* Created `.require_openmx()` to make it easier to use OpenMx functions without making OpenMx a dependency
89
* Smarter string ID handling for ped2id
910
* Fixed how different-sized matrices are handled by `com2links()`
@@ -22,7 +23,7 @@
2223
* Allow confidence intervals for pedigree mx wrappers
2324

2425
# BGmisc 1.6.0.1
25-
## CRAN submission
26+
* CRAN submission
2627
* Add OpenMx pedigree model builders and docs
2728
* Added vignette for OpenMx pedigree model builders
2829
* Add option for MZ twins in the additive genetic matrix

‎R/helpSliceFamilies.R‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -62,15 +62,26 @@
6262
# @param data_directory Output directory path
6363
# @param verbose Print file names if TRUE
6464
# @keywords internal
65-
.write_bin_data <- function(data, range_min, range_max, mit_val, data_directory, verbose = FALSE) {
65+
.write_bin_data <- function(data, range_min, range_max,
66+
mit_val=NULL,
67+
data_directory, verbose = FALSE) {
68+
if (!is.null(mit_val)) {
6669
range_data <- data[
6770
base::round(data$addRel, 6) >= range_min &
6871
base::round(data$addRel, 6) < range_max &
6972
data$mitRel == mit_val,
7073
]
74+
file_path <- file.path(data_directory, paste0("df_mt", mit_val, "_r", range_min, "-r", range_max, ".csv"))
75+
} else {
76+
range_data <- data[
77+
base::round(data$addRel, 6) >= range_min &
78+
base::round(data$addRel, 6) < range_max
79+
]
80+
file_path <- file.path(data_directory, paste0("df_r", range_min, "-r", range_max, ".csv"))
81+
}
7182

7283
if (base::nrow(range_data) > 0) {
73-
file_name <- file.path(data_directory, paste0("df_mt", mit_val, "_r", range_min, "-r", range_max, ".csv"))
84+
file_name <- file_path
7485
if (verbose) {
7586
message(file_name)
7687
}

‎R/makeLinks.R‎

Lines changed: 11 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -434,9 +434,12 @@ process_all_three <- function(
434434
if (length(u) > 0) {
435435
ID1 <- ids[u]
436436
tds <- data.frame(ID1 = ID1, ID2 = ID2)
437-
tds[[name1]] <- if (!is.null(v1)) ifelse(u %in% v1$i, v1$x[match(u, v1$i)], 0) else 0
438-
tds[[name2]] <- if (!is.null(v2)) ifelse(u %in% v2$i, v2$x[match(u, v2$i)], 0) else 0
439-
tds[[name3]] <- if (!is.null(v3)) ifelse(u %in% v3$i, v3$x[match(u, v3$i)], 0) else 0
437+
tds[[name1]] <- if (!is.null(v1)) { idx <- match(u, v1$i); ifelse(is.na(idx), 0, v1$x[idx]) } else 0
438+
tds[[name2]] <- if (!is.null(v2)) { idx <- match(u, v2$i); ifelse(is.na(idx), 0, v2$x[idx]) } else 0
439+
tds[[name3]] <- if (!is.null(v3)) { idx <- match(u, v3$i); ifelse(is.na(idx), 0, v3$x[idx]) } else 0
440+
# tds[[name1]] <- if (!is.null(v1)) ifelse(u %in% v1$i, v1$x[match(u, v1$i)], 0) else 0
441+
# tds[[name2]] <- if (!is.null(v2)) ifelse(u %in% v2$i, v2$x[match(u, v2$i)], 0) else 0
442+
# tds[[name3]] <- if (!is.null(v3)) ifelse(u %in% v3$i, v3$x[match(u, v3$i)], 0) else 0
440443

441444
if (drop_upper_triangular) {
442445
tds <- tds[tds$ID1 <= tds$ID2, ]
@@ -446,7 +449,7 @@ process_all_three <- function(
446449
if (writetodisk) {
447450
write_buffer[[length(write_buffer) + 1L]] <- tds
448451
if (length(write_buffer) >= write_buffer_size) {
449-
utils::write.table(do.call(rbind, write_buffer),
452+
utils::write.table(data.table::rbindlist(write_buffer),
450453
file = rel_pairs_file, row.names = FALSE,
451454
col.names = FALSE, append = TRUE, sep = ","
452455
)
@@ -531,8 +534,10 @@ process_two <- function(
531534
if (length(u) > 0) {
532535
ID1 <- ids[u]
533536
tds <- data.frame(ID1 = ID1, ID2 = ID2)
534-
tds[[name1]] <- if (!is.null(v1)) ifelse(u %in% v1$i, v1$x[match(u, v1$i)], 0) else 0
535-
tds[[name2]] <- if (!is.null(v2)) ifelse(u %in% v2$i, v2$x[match(u, v2$i)], 0) else 0
537+
tds[[name1]] <- if (!is.null(v1)) { idx <- match(u, v1$i); ifelse(is.na(idx), 0, v1$x[idx]) } else 0
538+
tds[[name2]] <- if (!is.null(v2)) { idx <- match(u, v2$i); ifelse(is.na(idx), 0, v2$x[idx]) } else 0
539+
# tds[[name1]] <- if (!is.null(v1)) ifelse(u %in% v1$i, v1$x[match(u, v1$i)], 0) else 0
540+
# tds[[name2]] <- if (!is.null(v2)) ifelse(u %in% v2$i, v2$x[match(u, v2$i)], 0) else 0
536541

537542
if (drop_upper_triangular) {
538543
tds <- tds[tds$ID1 <= tds$ID2, ]

‎R/readGedcom.R‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -137,9 +137,12 @@ readGedcom <- function(file_path,
137137
birth = c("birth_date", "birth_lat", "birth_long", "birth_place"),
138138
death = c("death_caus", "death_date", "death_lat", "death_long", "death_place"),
139139
attributes = c(
140-
"attribute_caste", "attribute_children", "attribute_description", "attribute_education",
141-
"attribute_idnumber", "attribute_marriages", "attribute_nationality", "attribute_occupation",
142-
"attribute_property", "attribute_religion", "attribute_residence", "attribute_ssn",
140+
"attribute_caste", "attribute_children",
141+
"attribute_description", "attribute_education",
142+
"attribute_idnumber", "attribute_marriages",
143+
"attribute_nationality", "attribute_occupation",
144+
"attribute_property", "attribute_religion",
145+
"attribute_residence", "attribute_ssn",
143146
"attribute_title"
144147
),
145148
relationships = c("FAMC", "FAMS")
@@ -449,14 +452,16 @@ extract_info <- function(line, type) {
449452
#' @param file A data frame with a column \code{X1} containing GEDCOM lines.
450453
#' @return A list with counts of specific GEDCOM tag occurrences.
451454
countPatternRows <- function(file) {
452-
pattern_counts <- sapply(
455+
x <- file$X1
456+
pattern_counts <- vapply(
453457
c(
454458
"@ INDI", " NAME", " GIVN", " NPFX", " NICK", " SURN", " NSFX", " _MARNM",
455459
" BIRT", " DEAT", " SEX", " CAST", " DSCR", " EDUC", " IDNO", " NATI",
456460
" NCHI", " NMR", " OCCU", " PROP", " RELI", " RESI", " SSN", " TITL",
457461
" FAMC", " FAMS", " PLAC", " LATI", " LONG", " DATE", " CAUS"
458462
),
459-
function(pat) sum(grepl(pat, file$X1))
463+
function(pat) sum(grepl(pat, x, fixed = TRUE)),
464+
integer(1L)
460465
)
461466
num_rows <- list(
462467
num_indi_rows = pattern_counts["@ INDI"],

‎R/sliceFamilies.R‎

Lines changed: 28 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
#' @param data_directory Directory where output files will be saved. If NULL, it is constructed based on `outcome_name` and `folder_prefix`.
1919
#' @param verbose Logical; whether to print progress messages (default FALSE)
2020
#' @param addRel_ceiling Numeric. Maximum relatedness value to bin to. Default is 1.5
21+
#' @param mtdna Logical. Whether to separate bins by mitochondrial relatedness (mitRel) value. Default is TRUE
2122
#' @param error_handling Logical. Should more aggressive error handling be attempted? Default is FALSE
2223
#' @param max_retries Integer. Number of retry attempts with halved chunk size when error_handling is TRUE. Default is 2
2324
#' @return NULL. Writes CSV files to disk and updates progress logs.
@@ -31,6 +32,7 @@ sliceFamilies <- function(
3132
chunk_size = 2e7,
3233
max_lines = 1e13,
3334
addRel_ceiling = 1.5,
35+
mtdna= TRUE,
3436
input_file = NULL,
3537
folder_prefix = "data",
3638
progress_csv = "progress.csv",
@@ -43,6 +45,11 @@ sliceFamilies <- function(
4345
) {
4446
bin_width_string <- as.character(bin_width * 100)
4547

48+
# if(mtdna ==FALSE & "mitRel" %in% file_column_names) {
49+
# file_column_names <- file_column_names[!file_column_names %in% "mitRel"]
50+
# message("mtdna is set to FALSE, so 'mitRel' column will be ignored in processing.")
51+
# }
52+
4653
if (is.null(data_directory)) {
4754
# Set the data directory based on the outcome name and folder prefix
4855
link_suffix <- if (biggest == TRUE) "links_" else "links_allbut_"
@@ -167,6 +174,10 @@ sliceFamilies <- function(
167174
range_max <- addRel_maxs[i]
168175
range_min <- addRel_mins[i]
169176

177+
if(mtdna == TRUE) {
178+
if (verbose == TRUE) {
179+
message("Processing bin: ", range_min, " to ", range_max, " with mitRel = 1")
180+
}
170181
# filter the data for the current bin
171182
.write_bin_data(
172183
data = dataRelatedPair_merge,
@@ -184,7 +195,22 @@ sliceFamilies <- function(
184195
data_directory = data_directory,
185196
verbose = verbose
186197
)
198+
} else {
199+
if (verbose == TRUE) {
200+
message("Processing bin: ", range_min, " to ", range_max)
201+
}
202+
.write_bin_data(
203+
data = dataRelatedPair_merge,
204+
range_min = range_min,
205+
range_max = range_max,
206+
mit_val = NULL,
207+
data_directory = data_directory,
208+
verbose = verbose
209+
)
210+
}
187211
}
212+
213+
188214
base::message(start_line)
189215
df_nrows <- base::nrow(dataRelatedPair_merge)
190216
if (verbose == TRUE) {
@@ -243,4 +269,5 @@ sliceFamilies <- function(
243269

244270
# Close the progress file
245271
base::close(progress_status_conn)
246-
}
272+
}
273+

‎man/sliceFamilies.Rd‎

Lines changed: 3 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎tests/testthat/test-sliceFamilies.R‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -288,6 +288,37 @@ test_that(".write_bin_data creates file only when matching data exists", {
288288

289289
files_after <- list.files("bin_test", pattern = "\\.csv$")
290290
expect_equal(length(files_after), 1) # still only 1 file
291+
292+
# Write mitRel=0 bin that matches addRel ~0.5 - should create a new file
293+
BGmisc:::.write_bin_data(test_dt,
294+
range_min = 0.45, range_max = 0.55, mit_val = 0,
295+
data_directory = "bin_test", verbose = FALSE
296+
)
297+
files_final <- list.files("bin_test", pattern = "\\.csv$")
298+
expect_equal(length(files_final), 2) # now we should have 2 files
299+
expect_all_true(files_final %in% c("df_mt0_r0.45-r0.55.csv","df_mt1_r0.45-r0.55.csv"))
300+
301+
BGmisc:::.write_bin_data(test_dt,
302+
range_min = 0.45, range_max = 0.55, mit_val = NULL,
303+
data_directory = "bin_test", verbose = FALSE
304+
)
305+
files_final_final <- list.files("bin_test", pattern = "\\.csv$")
306+
expect_equal(length(files_final_final), 3) # now we should have 3 files
307+
expect_all_true(files_final_final %in% c( "df_r0.45-r0.55.csv","df_mt0_r0.45-r0.55.csv","df_mt1_r0.45-r0.55.csv"))
308+
309+
written1 <- data.table::fread(file.path("bin_test", "df_r0.45-r0.55.csv"))
310+
written2 <- data.table::fread(file.path("bin_test", "df_mt0_r0.45-r0.55.csv"))
311+
written3 <- data.table::fread(file.path("bin_test", "df_mt1_r0.45-r0.55.csv"))
312+
expect_equal(nrow(written1), 2)
313+
expect_equal(nrow(written2), 1)
314+
expect_equal(nrow(written3), 1)
315+
316+
expect_all_true(written1$V1 %in% c(written2$V1, written3$V1))
317+
expect_all_true(written2$V1 %in% written1$V1)
318+
expect_all_true(written3$V1 %in% written1$V1)
319+
expect_false(any(written2$V1 %in% written3$V1)) # mitRel=0 vs mitRel=1 should not overlap
320+
expect_false(any(written3$V1 %in% written2$V1))
321+
expect_false(any(written1$V1 %in% c(4,5,6))) # ID1=2,3 should not be in the 0.45-0.55 bin
291322
})
292323

293324
test_that("sliceFamilies uses file.path correctly for output paths (no trailing slash needed)", {

0 commit comments

Comments
 (0)