magnusdv / pedtools

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

Changing mutation model #19

Closed thoree closed 5 years ago

thoree commented 5 years ago

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)

magnusdv commented 5 years ago

First of all, I don't think there is a simple/useful way to stop you from tinkering with the attributes of an R object.

Secondly, as you already noticed, it is the mutation matrix that matters in computations; the rate parameter is mainly used to create this matrix. Changing the rate means changing the whole model; hence it makes sense (to me, at least) that the way to do this is to replace the whole model. (See code example below.)

Thirdly, in the particular case r=0 you could of course just remove the mutation model entirely.

library(pedtools)

# Attach marker locus with mutation model 
locus = list(alleles = 1:2, mutmod = "equal", rate = 0.05)
x = setMarkers(nuclearPed(1), locus = list(locus))

# Setting a new rate
mutmod(x, marker = 1) = pedmut::mutationModel("equal", alleles = alleles(x, 1), rate = 0)

# Inspect it:
mutmod(x, marker = 1)
#> Unisex mutation matrix:
#>   1 2
#> 1 1 0
#> 2 0 1
#> 
#> Model: equal 
#> Rate: 0 
#> 
#> Lumpable: Always

# Or simply remove the mutation model
mutmod(x, marker = 1) = NULL

Created on 2019-08-01 by the reprex package (v0.3.0)

thoree commented 5 years ago

Thanks. One could also go "trivial":

library(pedtools)

# Attach marker locus with mutation model 
locus = list(alleles = 1:2, mutmod = "equal", rate = 0.05)
x = setMarkers(nuclearPed(1), locus = list(locus))

# Setting a new rate
mutmod(x, marker = 1) = "trivial"

# Inspect it:
mutmod(x, marker = 1)
#> Unisex mutation matrix:
#>   1 2
#> 1 1 0
#> 2 0 1
#> 
#> Model: trivial 
#> Rate: NA 
#> 
#> Lumpable: Always

Created on 2019-08-02 by the reprex package (v0.3.0)

magnusdv commented 5 years ago

Ah, yes! I had completely forgotten about the "trivial" model. :)