## bug #29904: rec_gen_table produces impossible offspring

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

Add a New Comment (Rich Markup):

## You are not logged in

Please log in, so followups can be emailed to you.

### document.write('<a onclick="document.getElementById(\'hidsubpartcontentdiscussion\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkhidediscussion\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkshowdiscussion\').style.display=\'block\';" id="hidsubpartlinkhidediscussion" style="display: block" href="#discussion"><span class="minusorplus">(-)</span> Discussion</a>'); document.write('<a onclick="document.getElementById(\'hidsubpartcontentdiscussion\').style.display=\'inline\'; document.getElementById(\'hidsubpartlinkhidediscussion\').style.display=\'block\'; document.getElementById(\'hidsubpartlinkshowdiscussion\').style.display=\'none\';" id="hidsubpartlinkshowdiscussion" style="display: none" href="#discussion"><span class="minusorplus">(+)</span> Discussion</a>'); Discussion

Thu Jun 10 19:14:58 2010, 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>
Thu Jun 3 17:49:16 2010, 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."

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>
Thu Jun 3 17:16:46 2010, 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>
Wed Jun 2 19:21:24 2010, comment #10:

As I suspected:

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>
Wed Jun 2 18:18:48 2010, comment #9:

Here's the output from the new algorithm:

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>
Tue Jun 1 20:35:30 2010, 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>
Sat May 29 15:22:12 2010, 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>
Fri May 28 21:11:25 2010, 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>
Thu May 20 19:16:47 2010, 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>
Thu May 20 03:03:17 2010, 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>
Wed May 19 20:42:21 2010, 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>
Wed May 19 19:38:30 2010, 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>
Wed May 19 14:48:09 2010, 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>
Tue May 18 15:44:17 2010, 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):

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>

### document.write('<a onclick="document.getElementById(\'hidsubpartcontentattached\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkhideattached\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkshowattached\').style.display=\'block\';" id="hidsubpartlinkhideattached" style="display: block" href="#attached"><span class="minusorplus">(-)</span> Attached Files</a>'); document.write('<a onclick="document.getElementById(\'hidsubpartcontentattached\').style.display=\'inline\'; document.getElementById(\'hidsubpartlinkhideattached\').style.display=\'block\'; document.getElementById(\'hidsubpartlinkshowattached\').style.display=\'none\';" id="hidsubpartlinkshowattached" style="display: none" href="#attached"><span class="minusorplus">(+)</span> Attached Files</a>'); Attached Files

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

Attach File(s):

Comment:

No files currently attached

### document.write('<a onclick="document.getElementById(\'hidsubpartcontentdependencies\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkhidedependencies\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkshowdependencies\').style.display=\'block\';" id="hidsubpartlinkhidedependencies" style="display: block" href="#dependencies"><span class="minusorplus">(-)</span> Dependencies</a>'); document.write('<a onclick="document.getElementById(\'hidsubpartcontentdependencies\').style.display=\'inline\'; document.getElementById(\'hidsubpartlinkhidedependencies\').style.display=\'block\'; document.getElementById(\'hidsubpartlinkshowdependencies\').style.display=\'none\';" id="hidsubpartlinkshowdependencies" style="display: none" href="#dependencies"><span class="minusorplus">(+)</span> Dependencies</a>'); Dependencies

Depends on the following items: None found

Items that depend on this one: None found

### document.write('<a onclick="document.getElementById(\'hidsubpartcontentcc\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkhidecc\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkshowcc\').style.display=\'block\';" id="hidsubpartlinkhidecc" style="display: none" href="#cc"><span class="minusorplus">(-)</span> Mail Notification Carbon-Copy List</a>'); document.write('<a onclick="document.getElementById(\'hidsubpartcontentcc\').style.display=\'inline\'; document.getElementById(\'hidsubpartlinkhidecc\').style.display=\'block\'; document.getElementById(\'hidsubpartlinkshowcc\').style.display=\'none\';" id="hidsubpartlinkshowcc" style="display: block" href="#cc"><span class="minusorplus">(+)</span> Mail Notification Carbon-Copy List</a>'); Mail Notification Carbon-Copy List

Carbon-Copy List
• -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.

Please enter the title of George Orwell's famous dystopian book (it's a date):

### document.write('<a onclick="document.getElementById(\'hidsubpartcontenthistory\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkhidehistory\').style.display=\'none\'; document.getElementById(\'hidsubpartlinkshowhistory\').style.display=\'block\';" id="hidsubpartlinkhidehistory" style="display: none" href="#history"><span class="minusorplus">(-)</span> History</a>'); document.write('<a onclick="document.getElementById(\'hidsubpartcontenthistory\').style.display=\'inline\'; document.getElementById(\'hidsubpartlinkhidehistory\').style.display=\'block\'; document.getElementById(\'hidsubpartlinkshowhistory\').style.display=\'none\';" id="hidsubpartlinkshowhistory" style="display: block" href="#history"><span class="minusorplus">(+)</span> History</a>'); History

Follow 8 latest changes.

Date Changed By Updated Field Previous Value => Replaced By
Thu Jun 10 19:14:58 2010trashbird1240StatusIn ProgressFixed
Open/ClosedOpenClosed
Sat May 29 15:25:29 2010trashbird1240StatusFixedIn Progress
Sat May 29 15:22:12 2010trashbird1240Severity4 - Important5 - Blocker
Thu May 20 19:16:47 2010trashbird1240Severity5 - Blocker4 - Important
Thu May 20 03:03:17 2010trashbird1240StatusIn ProgressFixed
Wed May 19 14:47:27 2010trashbird1240CategoryNoneLibrary
StatusNoneIn Progress