# GEMMA output differs from R/qtl2

# Tags

* assigned: pjotrp, davea
* priority: high
* type: bug, enhancement
* status: unclear
* keywords: database, gemma, reaper, rqtl2

# Description

When running trait BXD_21526 results differ significantly.

=> https://genenetwork.org/show_trait?trait_id=21526&dataset=BXDPublish
=> https://genenetwork.org/show_trait?trait_id=21529&dataset=BXDPublish

So I confirm I am getting the same results as Dave in GN for GEMMA.

# Tasks

I run GEMMA for precompute on the command line and that I confirmed to
be the same as what we see in the browser. This suggests either data
or method is different with Dave's approach.

I confirmed that gemma in GN matches Dave's results. It is interesting
to see that running without LOCO has some impact, but not as bad as
the R/qtl2 difference. First we should check the genotype files to see
if they match. I checked that the phenotypes match.

Our inputs are different if I count genotypes (first yours, the other on production):

```
     1  2184941 B
     2  2132744 D
     3   628980 H
     1  2195662 B
     2  2142959 D
     3   650168 H
```

The number of rows/markers is the same. So we probably added some
genometypes, but if we miss one that would matter. Dave you can find
the file in /home/wrk/BXD.geno on tux02 if you want to look.

Dave, I notice that you don't use H in the R/qtl2 control file. That
might make a difference though it probably won't explain what we see
now. BTW I also correlated the LOD scores from GEMMA and R/qtl2 in
your spreadsheet and at 0.7 that is too low.  So it is probably not
just a magnitude problem. The results differ a lot in your
spreadsheet.

Next step is that I need to run R/qtl2 using the script in your
dropbox and see what Karl's code does. The exercise does not hurt
because it will help us bring R/qtl2 to GN.