The pedigree routines came out of a simple need – to quickly draw a Pedigree structure on the screen, within R, that was “good enough” to help with debugging the actual routines of interest, which were those for fitting mixed effecs Cox models to large family data. As such the routine had compactness and automation as primary goals; complete annotation (monozygous twins, multiple types of affected status) and most certainly elegance were not on the list. Other software could do that much better.
It therefore came as a major surprise when these routines proved useful to others. Through their constant feedback, application to more complex pedigrees, and ongoing requests for one more feature, the routine has become what it is today. This routine is still not suitable for really large pedigrees, nor for heavily inbred ones such as in animal studies, and will likely not evolve in that way. The authors fondest hope is that others will pick up the project.
The Pedigree function is the first step, creating an object of class
Pedigree.
It accepts the following input
Note that a factor variable is not listed as one of the choices for the subject identifier. This is on purpose. Factors were designed to accomodate character strings whose values came from a limited class – things like race or gender, and are not appropriate for a subject identifier. All of their special properties as compared to a character variable turn out to be backwards for this case, in particular a memory of the original level set when subscripting is done.
However, due to the awful decision early on in S to automatically turn every character into a factor — unless you stood at the door with a club to head the package off — most users have become ingrained to the idea of using them for every character variable.
(I encourage you to set the global option
stringsAsFactors = FALSE
to turn off autoconversion – it
will measurably improve your R experience).
Therefore, to avoid unnecessary hassle for our users the code will
accept a factor as input for the id variables, but the final structure
does not retain it.
Gender and relation do become factors. Status follows the pattern of the
survival routines and remains an integer.
Based on the dataframe given for ped_df and rel_df and their corresponding named list, the columns are renamed for them to be used correctly. The renaming is done as follow
rel_df <- data.frame(
indId1 = c("110", "204"),
indId2 = c("112", "205"),
code = c(1, 2),
family = c("1", "2")
)
cols_ren_rel <- list(
id1 = "indId1",
id2 = "indId2",
famid = "family"
)
## Rename columns rel
old_cols <- as.vector(unlist(cols_ren_rel))
new_cols <- names(cols_ren_rel)
cols_to_ren <- match(old_cols, names(rel_df))
names(rel_df)[cols_to_ren[!is.na(cols_to_ren)]] <-
new_cols[!is.na(cols_to_ren)]
print(rel_df)
## id1 id2 code famid
## 1 110 112 1 1
## 2 204 205 2 2
If the normalisation process is selected
normalize = TRUE
, then both dataframe will be checked by
their dedicated normalization function. It will ensure that all
modalities are written correctly and set up the right way. If a famid
column is present in the dataframe, then it will be aggregated to the id
of each individual and separated by an ’’_’’ to ensure the uniqueness of
the individuals identifiers.
## sex id avail
## Min. :1.000 Length:55 Min. :0.0000
## 1st Qu.:1.000 Class :character 1st Qu.:0.0000
## Median :1.000 Mode :character Median :0.0000
## Mean :1.491 Mean :0.4364
## 3rd Qu.:2.000 3rd Qu.:1.0000
## Max. :2.000 Max. :1.0000
## sex id avail
## male :28 Length:55 Mode :logical
## female :27 Class :character FALSE:31
## unknown : 0 Mode :character TRUE :24
## terminated: 0
If any error is detected after the normalisation process, then the normalised dataframe is gave back to the user with errors column added describing the encountered problems.
## Warning in .local(obj, ...): The relationship informations are not valid. Here is the normalised
## relationship informations with the identified problems
## id1 id2 code famid error
## 1 1_110 1_112 MZ twin 1 <NA>
## 2 2_204 2_205 <NA> 2 CodeNotRecognise
Now that the data for the Pedigree object creation are ready, they are given to a new Pedigree object, trigerring the validation process.
This validation step will check up for many errors such as:
After validation an S4 object is generated. This new concept make it possible to easily setup methods for this new type of object. The controls of the parameters is also more precise.
The Pedigree object contains 4 slots, each of them contains a different S4 object containing a specific type of information used for the Pedigree construction.
For more information on each object:
help(Ped)
help(Rel)
help(Scales)
help(Hints)
As the Pedigree object is now an S4 class, we have made available a number of accessors. Most of them can be used as a getter or as a setter to modify a value in the correponding slot of the object
The mcols() accessors is the one you should use to add more informations to your individuals.
## DataFrame with 55 rows and 5 columns
## error sterilisation vitalStatus affection_mods avail_mods
## <character> <logical> <logical> <numeric> <numeric>
## 1 NA NA NA 0 0
## 2 NA NA NA 1 0
## 3 NA NA NA 1 0
## 4 NA NA NA 0 0
## 5 NA NA NA NA 0
## ... ... ... ... ... ...
## 51 NA NA NA 0 0
## 52 NA NA NA 0 1
## 53 NA NA NA 0 1
## 54 NA NA NA 0 0
## 55 NA NA NA 1 1
## Add new columns as a threshold if identifiers of individuals superior
## to a given threshold for example
mcols(ped)$idth <- ifelse(as.numeric(mcols(ped)$indId) < 200, "A", "B")
mcols(ped)$idth
## [1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A"
## [25] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B" "B" "B"
## [49] "B" "B" "B" "B" "B" "B" "B"
With this new S4 object comes multiple methods to ease the use of it:
[
## We can change the family name based on an other column
ped <- upd_famid_id(ped, mcols(ped)$idth)
## We can substract a given family
pedA <- ped[famid(ped) == "A"]
## Plot it
plot(pedA, cex = 0.5)
## Pedigree object with
## [1] "Ped object with 41 individuals and 13 metadata columns"
## [1] "Rel object with 0 relationshipswith 0 MZ twin, 0 DZ twin, 0 UZ twin, 0 Spouse"
## $id
## [1] "A_101" "A_102" "A_103" "A_104" "A_105" "A_106" "A_107" "A_108" "A_109" "A_110" "A_111" "A_112"
## [13] "A_113" "A_114" "A_115" "A_116" "A_117" "A_118" "A_119" "A_120" "A_121" "A_122" "A_123" "A_124"
## [25] "A_125" "A_126" "A_127" "A_128" "A_129" "A_130" "A_131" "A_132" "A_133" "A_134" "A_135" "A_136"
## [37] "A_137" "A_138" "A_139" "A_140" "A_141"
##
## $dadid
## [1] NA NA "A_135" NA NA NA NA NA "A_101" "A_103" "A_103" "A_103"
## [13] NA "A_103" "A_105" "A_105" NA "A_105" "A_105" "A_107" "A_110" "A_110" "A_110" "A_110"
## [25] "A_112" "A_112" "A_114" "A_114" "A_117" "A_119" "A_119" "A_119" "A_119" "A_119" NA NA
## [37] NA "A_135" "A_137" "A_137" "A_137"
##
## $momid
## [1] NA NA "A_136" NA NA NA NA NA "A_102" "A_104" "A_104" "A_104"
## [13] NA "A_104" "A_106" "A_106" NA "A_106" "A_106" "A_108" "A_109" "A_109" "A_109" "A_109"
## [25] "A_118" "A_118" "A_115" "A_115" "A_116" "A_120" "A_120" "A_120" "A_120" "A_120" NA NA
## [37] NA "A_136" "A_138" "A_138" "A_138"
## Shrink it to keep only the necessary information
lst1_s <- shrink(pedA, max_bits = 10)
plot(lst1_s$pedObj, cex = 0.5)
## 10 x 10 sparse Matrix of class "dsCMatrix"
## [[ suppressing 10 column names 'A_101', 'A_102', 'A_103' ... ]]
##
## A_101 0.50 . . . . . . . 0.25 .
## A_102 . 0.50 . . . . . . 0.25 .
## A_103 . . 0.50 . . . . . . 0.25
## A_104 . . . 0.50 . . . . . 0.25
## A_105 . . . . 0.5 . . . . .
## A_106 . . . . . 0.5 . . . .
## A_107 . . . . . . 0.5 . . .
## A_108 . . . . . . . 0.5 . .
## A_109 0.25 0.25 . . . . . . 0.50 .
## A_110 . . 0.25 0.25 . . . . . 0.50
## Get the useful individuals
pedA <- useful_inds(pedA, informative = "AvAf")
as.data.frame(ped(pedA))["useful"][1:10,]
## [1] FALSE FALSE TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Pedixplorer_1.1.0 BiocStyle_2.33.1
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.9 utf8_1.2.4 generics_0.1.3 tidyr_1.3.1
## [5] stringi_1.8.4 lattice_0.22-6 digest_0.6.36 magrittr_2.0.3
## [9] evaluate_0.24.0 grid_4.4.1 fastmap_1.2.0 plyr_1.8.9
## [13] jsonlite_1.8.8 Matrix_1.7-0 brio_1.1.5 BiocManager_1.30.23
## [17] purrr_1.0.2 fansi_1.0.6 scales_1.3.0 jquerylib_0.1.4
## [21] cli_3.6.3 rlang_1.1.4 munsell_0.5.1 withr_3.0.0
## [25] cachem_1.1.0 yaml_2.3.10 tools_4.4.1 dplyr_1.1.4
## [29] colorspace_2.1-1 ggplot2_3.5.1 BiocGenerics_0.51.0 buildtools_1.0.0
## [33] vctrs_0.6.5 R6_2.5.1 stats4_4.4.1 lifecycle_1.0.4
## [37] stringr_1.5.1 S4Vectors_0.43.2 pkgconfig_2.0.3 pillar_1.9.0
## [41] bslib_0.8.0 gtable_0.3.5 glue_1.7.0 Rcpp_1.0.13
## [45] highr_0.11 xfun_0.46 tibble_3.2.1 tidyselect_1.2.1
## [49] sys_3.4.2 knitr_1.48 htmltools_0.5.8.1 rmarkdown_2.27
## [53] maketools_1.3.0 testthat_3.2.1.1 compiler_4.4.1 quadprog_1.5-8