# bug #29904: rec_gen_table produces impossible offspring

 Submitted by: Joel J. Adamson Submitted on: Tue 18 May 2010 03:44:17 PM UTC Category: Library Severity: 5 - Blocker Item Group: None Status: Fixed Privacy: Public Assigned to: Joel J. Adamson Open/Closed: Closed

You are not logged in

## Discussion

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.

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.

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?

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.

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).

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.

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).

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.

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).

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.

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).

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

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

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.

## Attached Files

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

Attach Files:

Comment:

No files currently attached

## Dependencies

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)
•

Do you think this task is very important?
This task has 0 encouragements so far.

Only logged-in users can vote.