magnusdv / paramlink2

Parametric Linkage Analysis
3 stars 0 forks source link

incorrect lod score #4

Closed liuhankui closed 1 year ago

liuhankui commented 1 year ago

Assume we have a family that contain a father, a mother and five kids The father and the first two kids are healthy The mother and the last three kids have a same disease.

Some times we did not observe the genotypes of parents so we set their genotypes as missing (0/0)

modAD = diseaseModel("AD",penetrances = c(0,1,1)) x = nuclearPed(5) |> addMarker(geno = c('1/1','1/2','1/1','1/1','1/2','1/2','1/2')) |> #segregation marker1 addMarker(geno = c('1/2','2/2','1/2','1/2','2/2','2/2','2/2')) |> #segregation marker2 but not fellow AD model addMarker(geno = c('0/0','0/0','1/1','1/1','1/2','1/2','1/2')) |> #marker1 without parental genotypes addMarker(geno = c('0/0','0/0','1/2','1/2','2/2','2/2','2/2')) #marker2 without parental genotypes

aff<-c(1,2,1,1,2,2,2) lod(x,aff,modAD)

CHROM MARKER MB LOD 1 1 NA 1.2041 2 2 NA 0.0000 3 3 NA 0.8573 4 4 NA 0.9031

see the lod scores of the marker1 and marker2 are correct, but the marker4 have a higher lod score than marker3

we know that the kids genotypes of marker4 is as same as marker2, the lod score of marker4 was supposed to be 0 or lower than marker3.

magnusdv commented 1 year ago

Hi,

I don't see any obvious errors here. Marker2 has lod=0 because the mother is homozygous and therefore uninformative for linkage. For marker4 she is informative because she has nonzero probability of being 1/2.

When you say "the lod score of marker4 was supposed to be ...", which authority are you referring to?

liuhankui commented 1 year ago

I mean marker3 provide co-segregated genotypes of five kids and it got a 0.85 lod score but marker4 only provide probability of mother being 1/2 but it got a 0.93 lod score

marker3 is more like a associated locus compared with marker4

magnusdv commented 1 year ago

I find it hard to follow hand-waving arguments like that when it comes to lod scores... I'm afraid you have to dig into the math. :)

The lod scores are correct, and the difference between m3 and m4 are in the denominator, i.e., assuming no linkage. The following calculation shows that genotype probabilities are different for the two markers:

library(pedsuite, quietly = T)

x = nuclearPed(5) |>
  addMarker(geno = c('0/0','0/0','1/1','1/1','1/2','1/2','1/2')) |>
  addMarker(geno = c('0/0','0/0','1/2','1/2','2/2','2/2','2/2'))

likelihood(x)
#> [1] 0.009765625 0.008789062

Created on 2023-04-04 with reprex v2.0.2

liuhankui commented 1 year ago

I find it hard to follow hand-waving arguments like that when it comes to lod scores... I'm afraid you have to dig into the math. :)

The lod scores are correct, and the difference between m3 and m4 are in the denominator, i.e., assuming no linkage. The following calculation shows that genotype probabilities are different for the two markers:

library(pedsuite, quietly = T)

x = nuclearPed(5) |>
  addMarker(geno = c('0/0','0/0','1/1','1/1','1/2','1/2','1/2')) |>
  addMarker(geno = c('0/0','0/0','1/2','1/2','2/2','2/2','2/2'))

likelihood(x)
#> [1] 0.009765625 0.008789062

Created on 2023-04-04 with reprex v2.0.2

yes, i should calculate the lod(0.5) of null hypothesis. thanks for giving me the likelihood(x), this function work.