magnusdv / pedtools

Tools for working with pedigrees in R
GNU General Public License v3.0
23 stars 3 forks source link

transferMarkers changing too much? #9

Closed thoree closed 5 years ago

thoree commented 5 years ago

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)

magnusdv commented 5 years ago

Thanks. The problem here was in fact wrong input that wasn't properly caught. By default, when erase = FALSE, markers are matched on name. In your example the markers were unnamed, which should have (and now does) raise an error.

If the markers are unnamed, but you are sure they match and are in the same order, then you can avoid the name matching by setting matchNames = FALSE.

magnusdv commented 5 years ago

An example similar to yours but a bit differently coded:

library(pedtools)
x = nuclearPed(father = "FA", children = "CH")
m1 = marker(x, FA = 2, alleles = 1:2, name = 'locus1')
m2 = marker(x, FA = 'a', alleles = c('a', 'b'), name = 'locus2')
x = setMarkers(x, list(m1, m2))
x
#>  id fid mid sex locus1 locus2
#>  FA   *   *   1    2/2    a/a
#>   2   *   *   2    -/-    -/-
#>  CH  FA   2   1    -/-    -/-

# Extract child as singleton (including locus annotations)
ch = subset(x, "CH")

# Add some genotypes and transfer back to the trio
genotype(ch, "locus1", "CH") = 1
genotype(ch, "locus2", "CH") = 'b'
transferMarkers(ch, x, erase = FALSE)
#>  id fid mid sex locus1 locus2
#>  FA   *   *   1    2/2    a/a
#>   2   *   *   2    -/-    -/-
#>  CH  FA   2   1    1/1    b/b
thoree commented 5 years ago

Thanks, I understand where I went wrong. I'll continue using transferMarkers for my application.

thoree commented 5 years ago

transferMarkers works nicely. Here's something that puzzled me:

library(pedtools)
x1 = nuclearPed(father = "FA1", children = "CH1")
m1 = marker(x1, FA1 = 2, alleles = 1:2, name = 'locus1')
m2 = marker(x1, FA1 = 'a', alleles = c('a', 'b'), name = 'locus2')
x1 = setMarkers(x1, list(m1, m2))

x2 = nuclearPed(father = "FA2", children = "CH2")
m1 = marker(x2, FA2 = 1, alleles = 1:2, name = 'locus1')
m2 = marker(x2, FA2 = 'b', alleles = c('a', 'b'), name = 'locus2')
x2 = setMarkers(x2, list(m1, m2))

ch1 = subset(x1, "CH1")
genotype(ch1, "locus1", "CH1") = 1
genotype(ch1, "locus2", "CH1") = 'b'

ch2 = subset(x2, "CH2")
genotype(ch2, "locus1", "CH2") = 2
genotype(ch2, "locus2", "CH2") = 'a'
transferMarkers(list(ch1, ch2), list(x1, x2), erase = FALSE) #Fine:
#> [[1]]
#>   id fid mid sex locus1 locus2
#>  FA1   *   *   1    2/2    a/a
#>    2   *   *   2    -/-    -/-
#>  CH1 FA1   2   1    1/1    b/b
#> 
#> [[2]]
#>   id fid mid sex locus1 locus2
#>  FA2   *   *   1    1/1    b/b
#>    2   *   *   2    -/-    -/-
#>  CH2 FA2   2   1    2/2    a/a
# However:
transferMarkers(list(ch1, ch2), list(x1, x2), ids = c("CH1", "CH2"), erase = FALSE)
#> Error: Unknown ID label: CH2

Created on 2018-11-04 by the reprex package (v0.2.1)

magnusdv commented 5 years ago

Ah, more unit tests are needed for this, clearly.

Fixed now.