magnusdv / pedtools Goto Github PK

View Code? Open in Web Editor NEW

21.0 21.0 3.0 1.07 MB

Tools for working with pedigrees in R

License: GNU General Public License v3.0

R 100.00%

pedtools's People

Contributors

Stargazers

Watchers

Forkers

luansheng asppagh

pedtools's Issues

distributeMarkers and NA

The NA-s and -s below appear strange to me:

library(pedtools)
x = distributeMarkers(nuclearPed(), n = 20, prefix = "M")
tail(getMap(x))
#>    CHROM MARKER          MB
#> 15    17    M15 100.2688989
#> 16    19    M16  18.4171132
#> 17    22    M17   0.9116074
#> 18    NA   <NA>  56.8220516
#> 19    NA   <NA>  40.7455858
#> 20    NA   <NA>  33.1612800

Created on 2022-03-29 by the reprex package (v2.0.1)

New errors breaking (marriage) loops

Currently, I am facing a new error when trying to break loops:
If the pedigree has marriage loops, I got the error:

Marriage loops detected, trying different selection method
Loop breakers: 001, 004
Error: Cannot restore to ped object: Disconnected input.

Breaking marriage loops worked for me the last time I used the package. At least I got an error message, when in a pedigree with 2 marriage loops, breaking them was not possible.

Furthermore, if the pedigree is split, and has _comp elements, I got the error-message "WARNING:numerical expression has 2 elements: only the first used". I am not sure, how this worked the last time,

Best wishes, Silke

`addMarker()` ped list fail

library(pedtools)

This works:

singleton(1) |> addMarker(alleles = 1:2)
#>  id fid mid sex <1>
#>   1   *   *   1 -/-

But not this:

list(singleton(1), singleton(2)) |> addMarker(alleles = 1:2)
#> Error: Genotypes must be named

^{Created on 2022-03-08 by the reprex package (v2.0.1)}

Adding children across components

library(pedtools)

list(singleton(1), singleton(2, sex = 2)) |> 
  addSon(parents = 1:2)
#> Error: Adding children across components is not implemented yet

^{Created on 2023-09-18 with reprex v2.0.2}

Unplottable pedigrees

The class of pedigrees allowed in pedtools is larger than for the kinship2 package, which does the plotting. Hence certain pedigrees cannot be plotted as of yet. In particular this concerns pedigrees with selfing:

x = pedtools::ped(1:2, fid=c(0,1), mid=c(0,1), sex=c(0,1))
x
#>  id fid mid sex
#>   1   *   *   0
#>   2   1   1   1
plot(x)
#> Warning in kinship2::pedigree(id = ped$id, dadid = ped$fid, momid =
#> ped$mid, : More than 25% of the gender values are 'unknown'
#> Error in kinship2::pedigree(id = ped$id, dadid = ped$fid, momid = ped$mid, : pedigree must have at least 2 founders

And a different error here:

library(pedtools)
y = nuclearPed(1, sex = 0)
y = addChildren(y, father = 3, mother = 3, nch = 1)
y
#>  id fid mid sex
#>   1   *   *   1
#>   2   *   *   2
#>   3   1   2   0
#>   4   3   3   1
plot(y)
#> Error in kinship2::pedigree(id = ped$id, dadid = ped$fid, momid = ped$mid, : Id not male, but is a father: 3

fill argument for `plot.ped()`

Disconnected random pedigrees

For example:

pedtools::randomPed(g = 2, founders = 3, seed = 123)
#> $`_comp1`
#>  id fid mid sex
#>   1   *   *   1
#>   3   *   *   2
#>   4   1   3   1
#>   5   4   3   1
#> 
#> $`_comp2`
#>  id fid mid sex
#>   2   *   *   1
#> 
#> attr(,"class")
#> [1] "pedList" "list"

randomPed() should perhaps have an argument connected, by default TRUE.

Some functions not exported?

The below documented functions appear not to be available

library(pedtools)
as.data.frame.ped
#> Error in eval(expr, envir, enclos): object 'as.data.frame.ped' not found
as.matrix.ped
#> Error in eval(expr, envir, enclos): object 'as.matrix.ped' not found
as.ped.data.frame
#> Error in eval(expr, envir, enclos): object 'as.ped.data.frame' not found
as.matrix.ped
#> Error in eval(expr, envir, enclos): object 'as.matrix.ped' not found

^{Created on 2019-04-24 by the reprex package (v0.2.1)}

plotPedList not run example

The 'not run' example of plotPedList reads
...
plotPedList(list(plot1, plot2, plot3, plot4), widths=c(2,3,1,2),
frames=list(1,2:3,4), available=T, newdev=T)
....
I f available=TRUE is removed, the example runs nicely without warnings.

Plot one marker

No big deal (there are other ways to do the below), but I seem to recall that the below used to work

library(pedtools)
x = nuclearPed()
m = marker(x,  "1" = 1)
plot(x, m)
#> Error in if (draw) {: the condition has length > 1

^{Created on 2023-02-28 with reprex v2.0.2}

Definition of generations

Dear Magnus,

I am using your pedtools-function generations() to determine the generations in a pedigree, and to calculate the generations in relation to a certain target. It seems, that all 'founders' are listed as generation 1 in the pedigree, even if there are married to a generation 2 (or even higher?) person. So a couple with children can consist of a parent from generation 1 and another parent from generation 2, the child is listed as generation 3.
When I calculate the difference between the child and their parents, the one parent is one (correctly) or the other one two generations (this seems incorrect to me) ahead.

I hope, I could make my point clear. I would find a generations-function very useful.

Best wishes, Silke

`setAfreq()` should update mutation models

library(pedtools)

x = nuclearPed() |> 
  addMarker(alleles = 1:2, afreq = c(0.5, 0.5)) |> 
  setMutmod(model = "prop", rate = 0.1)

mutmod(x, 1)
#> Unisex mutation matrix:
#>     1   2
#> 1 0.9 0.1
#> 2 0.1 0.9
#> 
#> Model: proportional 
#> Rate: 0.1 
#> Frequencies: 0.5, 0.5 
#> 

# Change allele freqs
y = setAfreq(x, marker = 1,  afreq = c("1" = 0.2, "2" = 0.8))

# Mutation model not updated
mutmod(y, 1)
#> Unisex mutation matrix:
#>     1   2
#> 1 0.9 0.1
#> 2 0.1 0.9
#> 
#> Model: proportional 
#> Rate: 0.1 
#> Frequencies: 0.5, 0.5 

# This works
z = setMutmod(y, update = TRUE)
mutmod(z, 1)
#> Unisex mutation matrix:
#>        1      2
#> 1 0.7500 0.2500
#> 2 0.0625 0.9375
#> 
#> Model: proportional 
#> Rate: 0.1 
#> Frequencies: 0.2, 0.8 
#>

^{Created on 2023-09-18 with reprex v2.0.2}

Failed to install 'pedtools' from GitHub

devtools::install_github("magnusdv/pedtools")
#> Error: Failed to install 'pedtools' from GitHub:
#>   cannot open URL 'https://api.github.com/repos/magnusdv/pedtools/commits/master'

^{Created on 2020-06-08 by the reprex package (v0.3.0)}

sorted output from getAlleles

In the first case below getAlleles() output is not sorted according to the names of individuals,
as expected. In the second case the output is sorted

# Case 1
library(pedtools)
x = nuclearPed(1, father = "a", mother ="b")
m = marker(x, a = 1, b = 2)
x = addMarkers(x, m)
getAlleles(x, c("b", "a"))
#>   na1.1 na1.2
#> b "2"   "2"  
#> a "1"   "1"

#Case 2
s1 = singleton(1)
s2 = singleton(2)
m1 = marker(s1, name =  "l1", "1" = 1)
s1 = addMarkers(s1,m1)
m2 = marker(s2, name = "l1", "2" = 2)
s2 = addMarkers(s2, m2)
getAlleles(list(s1,s2), ids = c("2", "1"))
#>   l1.1 l1.2
#> 1 "1"  "1" 
#> 2 "2"  "2"

^{Created on 2018-11-24 by the reprex package (v0.2.1)}

Changing mutation model

Below I use attr to change mutation rates. This leaves
the mutation matrix unchanged so it doesn't make sense
and should perhap not be allowed:

library(pedtools)
locus = list(list(name = "M", alleles = 1:2, mutmod = "equal", rate = 0.05))
x = setMarkers(nuclearPed(1), locus = locus)
attr(mutmod(x, marker = 1)[[1]], "rate" ) = 0.0 
attr(mutmod(x, marker = 1)[[2]], "rate" ) = 0.0
mutmod(x, marker = 1)
#> Unisex mutation matrix:
#>      1    2
#> 1 0.95 0.05
#> 2 0.05 0.95
#> 
#> Model: equal 
#> Rate: 0 
#> Frequencies: 0.5, 0.5 
#> 
#> Stationary: Yes 
#> Reversible: Yes 
#> Lumpable: Always

# The likelihood performs fine as it only uses the mutation matrix: 
amat = matrix(c(1,1,2,2), ncol = 2, byrow = T)
x = setAlleles(x, ids = c(1,3), markers = "M", alleles = amat) 
pedprobr::likelihood(x,1)
#> [1] 0.00625

^{Created on 2019-07-31 by the reprex package (v0.3.0)}

Markers with same names allowed

library(pedtools)
x = nuclearPed(mo = "Mo")
m1 = marker(x, name = "M1", Mo = "1/2")
m2 = marker(x, name = "M1", Mo = "1/1")
x = setMarkers(x, list(m1, m2))
genotype(x,  "M1", "Mo")
#> [1] "1" "2"

^{Created on 2020-08-02 by the reprex package (v0.3.0)}

Functions for manipulating the allele matrix

These would be useful to have:

getAlleles(x, ids, markers)
setAlleles(x, ids, markers, alleles)

For instance, to remove all genotypes of individuals "a" should be done by

x = setAlleles(x, ids = "a", alleles = 0)

Alternative function names: getAlleleMatrix, setAlleleMatrix.

`transferMarkers()` doesn't catch duplicated IDs

Moved here from magnusdv/forrel#11

When ids is explicitly given, it appears that transferMarkers() doesn't check for duplicated IDs.

library(pedtools)
x = singleton("a")
x = setMarkers(x, marker(x, a = 1:2, name = "M"))

# Error as expected
transferMarkers(list(x,x), x)
#> Error: ID label is not unique: a

# This should give the same error
transferMarkers(list(x,x), x, ids = "a")
#>  id fid mid sex   M
#>   a   *   *   1 1/2

^{Created on 2018-11-08 by the reprex package (v0.2.1)}

Release pedtools 0.9.0

Prepare for CRAN release:

Submit to CRAN:

usethis::use_version('minor')
Update cran-comments.md
devtools::submit_cran()
Approve email

Wait for CRAN...

`print.ped()` should return the pedigree (not the data frame) invisibly

x = pedtools::nuclearPed() |> print()
#>  id fid mid sex
#>   1   *   *   1
#>   2   *   *   2
#>   3   1   2   1

class(x)
#> [1] "data.frame"

^{Created on 2023-03-24 with reprex v2.0.2}

`swapSex()` messes up X-chromosomal markers

library(pedtools)

x = nuclearPed(1) |> 
  addMarker(chrom = "X", geno = c("a", "a/b", "a")) |> 
  swapSex(3)
x
#>  id fid mid sex <1>
#>   1   *   *   1   a
#>   2   *   *   2 a/b
#>   3   1   2   2   a

plot(x, marker = 1)

y = nuclearPed(1, sex = 2) |> 
  addMarker(chrom = "X", geno = c("a", "a/b", "a/b")) |> 
  swapSex(3)
y
#>  id fid mid sex <1> 
#>   1   *   *   1   a 
#>   2   *   *   2 a/b 
#>   3   1   2   1 a/b?

plot(y, marker = 1)

^{Created on 2023-01-29 with reprex v2.0.2}

getAlleles for X

getAlleles() and setAlleles() are very useful new functions!
If two different alleles are given for a male on X,
setmarkers() (and addMarkers()) silently (OK?!) chooses the first one
but getAlleles() returns both alleles:

library(pedtools)
x = nuclearPed(1)
m = marker(x, `1` = c(2,1), alleles = 1:2, name = "m",  chrom = "X")
x = setMarkers(x, m)
x
#>  id fid mid sex   m
#>   1   *   *   1   2
#>   2   *   *   2 -/-
#>   3   1   2   1   -
getAlleles(x)
#>   m.1 m.2
#> 1 "2" "1"
#> 2 NA  NA 
#> 3 NA  NA

^{Created on 2018-11-14 by the reprex package (v0.2.1)}

Feature requests for `transferMarkers`

We discussed some potential new features of the very useful function transferMarkers:
(i) add optional checking for sex consistency and (ii) optionally return first argument from with transferred
genotypes removed.

transferMarkers changing too much?

Below also the genotype of the father "FA" is changed for the second marker:

library(pedtools)
a = 1:2
v1 = singleton("V1")
mv1.1 = marker(v1, V1 = 2, alleles = a)
mv1.2 = marker(v1, V1 = 1, alleles = a)
v1 = setMarkers(v1, list(mv1.1, mv1.2))

x = nuclearPed(1, father = "FA", children = "V1")
mx.1 = marker(x, FA = 2, alleles = a)
mx.2 = marker(x, FA = 1, alleles = a)
x = setMarkers(x, list(mx.1,mx.2))
transferMarkers(v1 ,x , ids = c("V1"), erase = FALSE)
#>  id fid mid sex <1> <2>
#>  FA   *   *   1 2/2 2/2
#>   2   *   *   2 -/- -/-
#>  V1  FA   2   1 2/2 1/1

^{Created on 2018-11-03 by the reprex package (v0.2.1)}

posCm

Here is something I don't understand:

library(pedtools)
x = nuclearPed(1)
m = marker(x)
chrom(m) = 1 # OK, but
posCm(m) = 100
#> Error in posCm(m) = 100: could not find function "posCm<-"

Created on 2018-08-16 by the reprex package (v0.2.0).

validate_ped

library(pedtools)
x = nuclearPed(nch = 2, children = 3:4)
validate_ped(x) # OK, but
x = nuclearPed(nch = 2, children = 1:2)
validate_ped(x)
#> Error: Malformed pedigree.
#>  Duplicated ID label: 1
#>  Duplicated ID label: 2

Created on 2018-08-15 by the reprex package (v0.2.0).

Implement `as.data.frame()` and `writePed()` for lists

`randomPed()` error

randomPed() plots a pedigree when it stops on error:

library(pedtools) # 2.1.1

randomPed(n = 5, maxGenerationGap = 0, seed = 1)
#> Error in sample.int(length(x), ...): invalid first argument

# # Other examples
# randomPed(n = 10, maxGenerationGap = 0, seed = 1)
# randomPed(n = 5, maxGenerationGap = 1, seed = 58)
# randomPed(n = 10, maxGenerationGap = 1, seed = 2)
# randomPed(n = 10, maxGenerationGap = 2, seed = 63)
# randomPed(n = 10, maxGenerationGap = 3, seed = 202)

^{Created on 2023-02-28 with reprex v2.0.2}

`readPed()` fails when ID labels have spaces

library(pedtools)
x = singleton("A B C")
tmp <- writePed(x, tempfile(), what = "ped")

readPed(tmp)
#> Error in read.table(pedfile, header = header, colClasses = "character", : more columns than column names
readPed(tmp, sep = "\t")
#> Error in read.table(pedfile, header = header, colClasses = "character", : more columns than column names

^{Created on 2021-08-27 by the reprex package (v2.0.1)}

The correct solution readPed(..., sep = "\t") fails because the sep argument is already used for something else (separating alleles).

Rename to alleleSep?

consolidate function - feature request

A function, say consolidate ensuring that ped objects are defined for
the same markers would be nice. Example:

library(forrel)
#> Loading required package: pedtools
x1 = singleton("s1")
x1 = setMarkers(x1, locusAttributes = NorwegianFrequencies[1:2])
x2 = singleton("s2")
x2 = setMarkers(x2, locusAttributes = NorwegianFrequencies[2:3])

# The function should ensure that  x1 and x2 are defined for the samme markers

miss1 = setdiff(name(x2), name(x1))
miss2 = setdiff(name(x1), name(x2))

x1 = x1 |> addMarkers(locusAttributes = NorwegianFrequencies[miss1])
x2 = x2 |> addMarkers(locusAttributes = NorwegianFrequencies[miss2])
list(x1, x2)
#> [[1]]
#>  id fid mid sex D3S1358 TH01 D21S11
#>  s1   *   *   1     -/-  -/-    -/-
#> 
#> [[2]]
#>  id fid mid sex TH01 D21S11 D3S1358
#>  s2   *   *   1  -/-    -/-     -/-

^{Created on 2022-04-27 by the reprex package (v2.0.1)}

`addParents()` should add parents before the child

... not append them at the end.

library(pedtools)
addParents(singleton(1), 1)
#> Father: Creating new individual with ID = 2
#> Mother: Creating new individual with ID = 3
#>  id fid mid sex
#>   1   2   3   1
#>   2   *   *   1
#>   3   *   *   2

Duplicated individuals

Apparenty nuclearPed and likelihood(https://github.com/magnusdv/pedprobr) tolerate some duplication, not plot :

library(pedtools)
library(pedprobr)
x = nuclearPed(children = c("sib1", "sib1"))
m = marker(x, sib1 = 1:2, sib1 = 1, alleles = 1:2)
likelihood(x,m)
#> [1] 0.25
plot(x,m)
#> Error in kinship2::pedigree(id = ped$id, dadid = ped$fid, momid = ped$mid, : Duplicate subject id: sib1

Created on 2018-08-13 by the reprex package (v0.2.0).

Malformed random pedigrees with selfing

library(pedtools)

x = randomPed(2, 2, selfing = TRUE, seed = 2)
x
#>  id fid mid sex
#>   1   *   *   1
#>   2   *   *   2
#>   3   1   2   1
#>   4   1   1   0

validatePed(x)
#> Error: Malformed pedigree.
#>  Individual 1 is male, but appear as the mother of 4

^{Created on 2022-03-08 by the reprex package (v2.0.1)}

`setFrequencyDatabase` and marker name

I did not manage to get the variable names right (below files
are in https://github.com/thoree/exclusion/tree/master/magnus
and I used reprex_addin() for the below):

library(pedtools)
x = readPed("pedfile.ped", sep ="/")
x = setFrequencyDatabase(x, database = "dbfile.txt")
#> M1 M2
#> 1 0.5 NA
#> 2 0.5 0.8
#> 3 NA 0.2
#> Error: Unknown marker name: M1, M2
Created on 2019-08-09 by the reprex package (v0.3.0)

`distributeMarkers` and `dist` argument

Return an error or only one marker if dist exceeds the the total chromosome length?:

library(pedtools)
x = distributeMarkers(nuclearPed(), dist = 3000, prefix = "M")
getMap(x)
#> CHROM MARKER MB
#> 1 1 M1 0
#> 2 2 M2 0
#> 3 3 M3 0
#> 4 4 M4 0
#> 5 5 M5 0
#> 6 6 M6 0
#> 7 7 M7 0
#> 8 8 M8 0
#> 9 9 M9 0
#> 10 10 M10 0
#> 11 11 M11 0
#> 12 12 M12 0
#> 13 13 M13 0
#> 14 14 M14 0
#> 15 15 M15 0
#> 16 16 M16 0
#> 17 17 M17 0
#> 18 18 M18 0
#> 19 19 M19 0
#> 20 20 M20 0
#> 21 21 M21 0
#> 22 22 M22 0
Created on 2022-03-29 by the reprex package (v2.0.1)

linearPed

I heard rumour of a new ped, linearPed (?), but can't find it ...

`randomPed()` error

pedtools::randomPed(n = 5, f = 2, seed = 58)
#> Error in sample.int(length(x), ...): invalid first argument

^{Created on 2023-04-20 with reprex v2.0.2}

proportional model in `mutationMatrix()`

So far only the "equal" model is supported. The "proportional" model should follow; use the draft (and tests) suggested by @thoree in magnusdv/forrel#1 (comment).

componentwise `nMarkers()`

nMarkers() should have an option to report number of markers for each component.

library(pedtools)

x = list(singleton("A") |> addMarker(),
         singleton("B"))

nMarkers(x)
#> Error: The pedigree components have different number of markers attached

^{Created on 2023-09-06 with reprex v2.0.2}

Problem installing with `build_vignettes = TRUE`

Installation failed just now iff I asked for the vignette:

devtools::install_github("magnusdv/pedtools",
build_vignettes = TRUE, force = TRUE)
#> Downloading GitHub repo magnusdv/pedtools@master
#>
#>

checking for file 'C:\Users\theg\AppData\Local\Temp\RtmpkVIyn5\remotes2fe447ce73ab\magnusdv-pedtools-93add0e/DESCRIPTION' ...

v checking for file 'C:\Users\theg\AppData\Local\Temp\RtmpkVIyn5\remotes2fe447ce73ab\magnusdv-pedtools-93add0e/DESCRIPTION'
#>

preparing 'pedtools':
#> checking DESCRIPTION meta-information ...

checking DESCRIPTION meta-information ...

v checking DESCRIPTION meta-information
#>

installing the package to build vignettes
#>

creating vignettes ...

creating vignettes ...

E creating vignettes (8.4s)
#> pandoc.exe: <>
#> Error: processing vignette 'pedtools.Rmd' failed with diagnostics:
#> pandoc document conversion failed with error 1
#> Execution halted
#>
#> Error: Failed to install 'pedtools' from GitHub:
#> System command error, exit status: 1, stdout + stderr:
#> E> * checking for file 'C:\Users\theg\AppData\Local\Temp\RtmpkVIyn5\remotes2fe447ce73ab\magnusdv-pedtools-93add0e/DESCRIPTION' ... OK
#> E> * preparing 'pedtools':
#> E> * checking DESCRIPTION meta-information ... OK
#> E> * installing the package to build vignettes
#> E> * creating vignettes ... ERROR
#> E> pandoc.exe: <>
#> E> Error: processing vignette 'pedtools.Rmd' failed with diagnostics:
#> E> pandoc document conversion failed with error 1
#> E> Execution halted
Created on 2019-08-11 by the reprex package (v0.3.0)

request for getSex for list of pedigrees and singletons

It would be nice with a version of getSex working
for a list of pedigrees and singletons:

library(pedtools)
v = list(singleton(1),singleton(2))
getSex(v,"1")
#> Error: Input is not a `ped` object

^{Created on 2018-12-05 by the reprex package (v0.2.1)}

Check for duplicated allele names?

library(pedprobr)
#> Loading required package: pedtools
x = nuclearPed()
m1 = marker(x,"1" = "a/a",   afreq = c(a = 0.9, a = 0.1))
likelihood(x,m1)
#> [1] 0.81

^{Created on 2020-08-02 by the reprex package (v0.3.0)}

`getAlleles()` output and marker names

The output from getAlleles() below is not as I expected.
The name given to the second marker, m2, appears to be ignored

#Case 2
library(pedtools)
s1 = singleton(1)
s2 = singleton(2)
m1 = marker(s1, name=  "locus", "1" = 1)
s1 = addMarkers(s1,m1)
m2 = marker(s2, name = "other locus", "2" = 2)
s2 = addMarkers(s2,m2)
getAlleles(list(s1,s2), ids =c("2", "1"))
#>   locus.1 locus.2
#> 1 "1"     "1"    
#> 2 "2"     "2"

^{Created on 2018-11-24 by the reprex package (v0.2.1)}

`custom` mutation model

Is it possible to also set the custom model. I unsuccessfully tried:

library(pedtools)
library(pedmut, quietly = TRUE)
mutmat = matrix(c(0.99, 0.1, 0.01, 0.9), ncol = 2, dimnames = list(1:2, 1:2)) 
locus = list(list(name = "M", alleles = 1:2, mutmod = "custom", matrix = mutmat))
x = setMarkers(nuclearPed(1), locus = locus)
#> Error: The "custom" model requires the argument `matrix` to be non-NULL

^{Created on 2019-07-31 by the reprex package (v0.3.0)}

implementation of missing parents

Hi Magnus, I was wondering if you were planning to extend your package to allow for cases where only one parent is missing? Many thanks!