bugNumerical Library for Haploid Genetics - Bugs: bug #29904, rec_gen_table produces impossible...

 
 

bug #29904: rec_gen_table produces impossible offspring

Submitter:  Joel J. Adamson <trashbird1240>
Submitted:  Tue 18 May 2010 03:44:17 PM UTC
   
 
Category:  Library Severity:  5 - Blocker
Item Group:  None Status:  Fixed
Privacy:  Public Assigned to:  trashbird1240
Open/Closed:  Closed
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 10 Jun 2010 07:14:58 PM UTC, comment #13: 

This is now fixed thanks to a discussion with Reinhard Burger and a more aggressive (albeit simplistic) classification of cases.

Joel J. Adamson <trashbird1240>
Group administrator
Thu 03 Jun 2010 05:49:16 PM UTC, comment #12: 

One pattern is that the totally recombinant offspring are all between zero and one.  Perhaps it is only the cases that include non-recombinant offspring that needs "adjustment."


Offspring 1:
0:0.00000000 1.50000000 0.00000000 0.15625000 0.00000000 0.25000000 0.00000000 0.12500000
1:1.50000000 1.00000000 1.34375000 1.50000000 1.25000000 1.50000000 1.21875000 1.34375000
2:0.00000000 1.34375000 0.00000000 0.00000000 0.00000000 0.12500000 0.00000000 0.00000000
3:0.15625000 1.50000000 0.00000000 0.00000000 0.03125000 0.15625000 0.00000000 0.00000000
4:0.00000000 1.25000000 0.00000000 0.03125000 0.00000000 0.00000000 0.00000000 0.00000000
5:0.25000000 1.50000000 0.12500000 0.15625000 0.00000000 0.00000000 0.00000000 0.00000000
6:0.00000000 1.21875000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
7:0.12500000 1.34375000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000


As we can see, the values over 1 are in the row and column that produces non-recombinant targets.  This pattern occurs for every target offspring.

Here are the exact words from Burger:

"...where $R_I (j_Ij_J, k_Ik_J -> i_Ii_J)$ is the probability that a randomly chosen gamete produced by a jk individual is i and each of the sets I and J is passed to the next generation without recombination."  In other words, there is only recombination between some locus in I and some locus in J.

Joel J. Adamson <trashbird1240>
Group administrator
Thu 03 Jun 2010 05:16:46 PM UTC, comment #11: 

One solution is to view the R(jk->i) as a relative frequency and divide by the total summed across possible offspring.  Then we'd get a probability.  We could only apply this to values that have been calculated (e.g. the diagonal intersection value 1.0, is not calculated).  I predict this would be computationally expensive, but it would be worth it.  There may be a way to hold the values in cache and execute the division at the end instead of looping twice.

The test for this idea is this: does it still yield valid results under the tlta model?

Joel J. Adamson <trashbird1240>
Group administrator
Wed 02 Jun 2010 07:21:24 PM UTC, comment #10: 

As I suspected:

Offspring 0:
0:1.00000000 1.50000000 1.50000000 1.34375000 1.50000000 1.25000000 1.34375000 1.21875000
1:1.50000000 0.00000000 0.15625000 0.00000000 0.25000000 0.00000000 0.12500000 0.00000000
2:1.50000000 0.15625000 0.00000000 0.00000000 0.15625000 0.03125000 0.00000000 0.00000000
3:1.34375000 0.00000000 0.00000000 0.00000000 0.12500000 0.00000000 0.00000000 0.00000000
4:1.50000000 0.25000000 0.15625000 0.12500000 0.00000000 0.00000000 0.00000000 0.00000000
5:1.25000000 0.00000000 0.03125000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
6:1.34375000 0.12500000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
7:1.21875000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000



The formula is wrong somehow (or my interpretation of it).  Unless the delta's are calculated differently from how I think they are, the formula is wrong.

Joel J. Adamson <trashbird1240>
Group administrator
Wed 02 Jun 2010 06:18:48 PM UTC, comment #9: 

Here's the output from the new algorithm:


Offspring 0:
0:1.00000000 0.00000000 0.00000000 0.50000000
1:0.00000000 0.00000000 0.00000000 0.00000000
2:0.00000000 0.00000000 0.00000000 0.00000000
3:0.50000000 0.00000000 0.00000000 0.00000000

Offspring 1:
0:0.00000000 0.00000000 0.00000000 0.00000000
1:0.00000000 1.00000000 0.50000000 0.00000000
2:0.00000000 0.50000000 0.00000000 0.00000000
3:0.00000000 0.00000000 0.00000000 0.00000000

Offspring 2:
0:0.00000000 0.00000000 0.00000000 0.00000000
1:0.00000000 0.00000000 0.50000000 0.00000000
2:0.00000000 0.50000000 1.00000000 0.00000000
3:0.00000000 0.00000000 0.00000000 0.00000000

Offspring 3:
0:0.00000000 0.00000000 0.00000000 0.50000000
1:0.00000000 0.00000000 0.00000000 0.00000000
2:0.00000000 0.00000000 0.00000000 0.00000000
3:0.50000000 0.00000000 0.00000000 1.00000000
x[0] = 0.25000000
x[1] = 0.25000000
x[2] = 0.25000000
x[3] = 0.25000000
Call to rec_mating () took 0.00000000 sec
PASS: rec_test
rec_test3: src/rec.c:119: rec_total: Assertion `total <= 1.0' failed.
/bin/sh: line 5:  3289 Aborted                 (core dumped) ${dir}$tst
FAIL: rec_test3
Value by hand: 0.170442
Value by ld_from_geno: 0.170442


promising, but I had already uncovered by hand that Bürger's formula appears to be wrong (or I completely misunderstand it).

Joel J. Adamson <trashbird1240>
Group administrator
Tue 01 Jun 2010 08:35:30 PM UTC, comment #8: 

Partitions are easily achieved by iterating over odd natural numbers from 1 (0001) to

(((1<<geno) - 1) & ~2) & ((1 << geno) - 1) (e.g. 1101).

It's so simple!  No need for a separate function; no need for an extra data structure.

Joel J. Adamson <trashbird1240>
Group administrator
Sat 29 May 2010 03:22:12 PM UTC, comment #7: 

The algorithm is not guaranteed to work in any way.  What we need is to follow an established algorithm instead of making wild guesses. 

Now the task is to find an algorithm that can generate the partitions needed to use the formula in Reinhard Burger's book: R_{ij \rightarrow k} = \sum_I R_I

the algorithm we need is to generate the sets I, which are all the non-empty disjoint pairwise partitions.  In other words, we need all the possible ways to split up the genome.  Once we find that we can generate the deltas easily (differences between the genomes of the parents and the offspring).

Joel J. Adamson <trashbird1240>
Group administrator
Fri 28 May 2010 09:11:25 PM UTC, comment #6: 

I've made a few more changes to rec_total(), meaning we need to carefully check the three-locus result.  The two-locus result is fine.

Joel J. Adamson <trashbird1240>
Group administrator
Thu 20 May 2010 07:16:47 PM UTC, comment #5: 

New algorithm identifies cases that need to call rec_total() completely; values appear correct for three loci.  We cannot close this bug until we have a systematic check (such as an externally-generated table that agrees).

Joel J. Adamson <trashbird1240>
Group administrator
Thu 20 May 2010 03:03:17 AM UTC, comment #4: 

A new algorithm is in place; however we cannot officially close this until we produce some way to check the three-locus case.  The two-locus case is easy to check by hand.  I will attempt to check it by hand tomorrow.

Joel J. Adamson <trashbird1240>
Group administrator
Wed 19 May 2010 08:42:21 PM UTC, comment #3: 

We need to verify the output of rec_test3; or we need to prove this rather than through calculation (e.g. Reiersol algebra).

Joel J. Adamson <trashbird1240>
Group administrator
Wed 19 May 2010 07:38:30 PM UTC, comment #2: 

rec_total:  TLTA problem solved; we now need to reduce larger genome problems to the two locus case.  That is, we need to loop over pairs of loci while properly multiplying the total probability

Joel J. Adamson <trashbird1240>
Group administrator
Wed 19 May 2010 02:48:09 PM UTC, comment #1: 

Fixed logic so that rec_gen_table no longer produces impossible offspring; however, inspection of table for offspring 0 from rec_test3 shows that probabilities are incorrect

Joel J. Adamson <trashbird1240>
Group administrator
Tue 18 May 2010 03:44:17 PM UTC, original submission:  

Definition: impossible offspring -- offspring that two parents cannot
produce through recombination

While developing my age dependent sexual selection simulation, I have
found that rec_gen_table produces non-zero offspring probabilities for
parents that cannot produce that particular offspring.  For example,
for offspring 010 (2):


Offspring 2:
0.00000000 0.25000000 0.43750000 0.25000000 0.25000000 0.00000000 0.25000000 0.25000000
0.25000000 0.00000000 0.31250000 0.12500000 0.12500000 0.00000000 0.12500000 0.12500000
0.43750000 0.31250000 1.00000000 0.43750000 0.43750000 0.18750000 0.31250000 0.43750000
0.25000000 0.12500000 0.43750000 0.00000000 0.25000000 0.00000000 0.25000000 0.25000000
0.25000000 0.12500000 0.43750000 0.25000000 0.00000000 0.00000000 0.25000000 0.25000000
0.00000000 0.00000000 0.18750000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.25000000 0.12500000 0.31250000 0.25000000 0.25000000 0.00000000 0.00000000 0.12500000
0.25000000 0.12500000 0.43750000 0.25000000 0.25000000 0.00000000 0.12500000 0.00000000


The entry at position (0,1) (and its transpose) is incorrect; it
should be zero.  I discovered this from seeing that rec_mating was
producing offspring frequencies for offspring genotypes that should
not be possible: I entered allele frequencies of zero, and was getting
out non-zero allele frequencies (impossble, again).

The tabular output was produced by copying sparse.h into my source
directory and using the same procedure as in tests/prtable.c of the
haploid source tree.

Joel J. Adamson <trashbird1240>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by trashbird1240 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-06-10 trashbird1240 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2010-05-29 trashbird1240 StatusFixed In Progress
    2010-05-29 trashbird1240 Severity4 - Important 5 - Blocker
    2010-05-20 trashbird1240 Severity5 - Blocker 4 - Important
    2010-05-20 trashbird1240 StatusIn Progress Fixed
    2010-05-19 trashbird1240 CategoryNone Library
        StatusNone In Progress

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code