magnusdv / forrel

Forensic pedigree analysis and relatedness inference
GNU General Public License v2.0
10 stars 0 forks source link

Error when plotting pedigrees with `exclusionPower()` #17

Closed knifecake closed 5 years ago

knifecake commented 5 years ago

Hi!

I ran into this problem when trying to run exclusionPower. Below is the output from reprex.

library(pedtools)
library(paramlink)
#> 
#> Attaching package: 'paramlink'
#> The following objects are masked from 'package:pedtools':
#> 
#>     addDaughter, addParents, addSon, ancestors, branch,
#>     breakLoops, connectedComponents, cousinPed, cousins,
#>     descendants, doubleCousins, doubleFirstCousins,
#>     findLoopBreakers, findLoopBreakers2, fullSibMating,
#>     getMarkers, grandparents, halfCousinPed, halfSibStack,
#>     is.singleton, leaves, marker, mendelianCheck, mergePed,
#>     nephews_nieces, nuclearPed, offspring, parents, plotPedList,
#>     quadHalfFirstCousins, randomPed, relabel, removeIndividuals,
#>     removeMarkers, setMarkers, siblings, singleton, spouses,
#>     swapSex, tieLoops, unrelated
library(forrel, quietly = T)
#> 
#> Attaching package: 'forrel'
#> The following objects are masked from 'package:paramlink':
#> 
#>     IBDestimate, IBDtriangle, LR, exclusionPower, markerSim,
#>     readFamiliasLoci, showInTriangle, simpleSim

############################################
### A standard case paternity case:
### Compute the power of exclusion when the claimed father is in fact
### unrelated to the child.
############################################

# Claim: Individual 1 is the father of indiv 3
claim = nuclearPed(nch = 1, sex = 2)
#> Error in nuclearPed(nch = 1, sex = 2): unused argument (nch = 1)

# Truth: 1 and 3 are unrelated
true = list(singleton(id = 1), singleton(id = 3, sex = 2))

# Indivs 1 and 3 are available for genotyping
ids = c(1, 3)

# An equifrequent SNP
als = 2
afreq = c(0.5, 0.5)

# The exclusion power assuming no known genotypes
PE1 = exclusionPower(claim, true, ids = ids, alleles = als, afreq = afreq)
#> Error in is.ped(ped_claim): object 'claim' not found

Created on 2019-06-05 by the reprex package (v0.3.0)

However, if I don't load the pedtools and paramlink packages I get a different error:

library(forrel, quietly = T)

############################################
### A standard case paternity case:
### Compute the power of exclusion when the claimed father is in fact
### unrelated to the child.
############################################

# Claim: Individual 1 is the father of indiv 3
claim = nuclearPed(nch = 1, sex = 2)

# Truth: 1 and 3 are unrelated
true = list(singleton(id = 1), singleton(id = 3, sex = 2))

# Indivs 1 and 3 are available for genotyping
ids = c(1, 3)

# An equifrequent SNP
als = 2
afreq = c(0.5, 0.5)

# The exclusion power assuming no known genotypes
PE1 = exclusionPower(claim, true, ids = ids, alleles = als, afreq = afreq)
#> Error in kinship2::pedigree(id = ped$id, dadid = ped$fid, momid = ped$mid, : Wrong length for affected

Created on 2019-06-05 by the reprex package (v0.3.0)

I was using R version 3.5.0 (2018-04-23) -- "Joy in Playing". As for the packages, I was using these versions:

> packageVersion('paramlink')
[1] '1.1.2'
> packageVersion('pedtools')
[1] '0.1.0'
> packageVersion('forrel')
[1] '0.1.0'

Thanks for your time!

magnusdv commented 5 years ago

Hi, Many thanks for supplying reprexes!

1) The reason for the first error is very clear: pedtools is not compatible with paramlink.

The story is as follows: The paramlink package was split (and re-implemented/improved) into several new packages collectively called the "ped suite". The main package here is pedtools, but you rarely have to load this explicitly, since it is automatically loaded by all the others (including forrel).

Paramlink still exists, but is no longer developed, and I of course strongly advice switching to the ped suite.

2) As for your second error I suspect it may be caused by breaking changes in kinship2. This was recently updated, and I haven't had time to test it yet.

In the meantime, could you try installing the old version 1.6.4 of kinship2, and try then?

magnusdv commented 5 years ago

Ok, this should be fixed now, and it should work with any version of kinship2.

magnusdv commented 5 years ago

(To be clear, you need to install the latest github version of pedtools.)