path: root/topics/systems/mariadb/precompute-publishdata.gmi

# Precompute PublishData

Based on the QTL_Reaper_cal_lrs.py aka QTL_Reaper_v8_PublishXRef.py. This script simply updates PublishXRef table with a highest hit as computed by qtlreaper.

In a first attempt to update the database we are going to do just that using GEMMA.

For the new script we will pass in the genotype file as well as the phenotype file, so gemma-wrapper can process it. I wrote quite a few scripts already

=> https://github.com/genetics-statistics/gemma-wrapper/tree/master/bin

So we can convert a .geno file to BIMBAM. I need to extract GN traits to a R/qtl2 or lmdb trait format file and use that as input.

* [X] Visit use of PublishXRef
* [X] geno -> BIMBAM (BXD first)
* [X] Get PublishData trait(s) and convert to gemma, R/qtl2 or lmdb
* - [X] see scripts/lmdb-publishdata-export.scm
* - [X] see scripts for ProbeSetData
* - [X] Make sure the BXDs are mappable
* [X] Run gemma-wrapper
* [X] We should map by trait-id, data id is not intuitive: curl http://127.0.0.1:8091/dataset/bxd-publish/values/8967044.json > 10002-pheno.json
* [X] Check why Zach/GN JSON file lists different mappable BXDs
* [X] Update DB on run-server
* [X] Add batch run and some metadata so we can link back from results
* [ ] Check lmdb duplicate key warning
* [ ] Create a DB/table containing hits and old reaper values
* [ ] Update PublishXRef and store old reaper value(?)
* [ ] Correctly Handle escalating errors
* [ ] Make sure the trait fetcher handles authorization or runs localhost only
* [ ] gemma-wrapper --force does not work for GRM and re-check GRM does not change on phenotype

For the last we should probably add a few columns. Initially we'll only store the maximum hit.

After

* [ ] provide distributed storage of files using https

# Visit use of PublishXRef

In GN2 this table is used in search, auth, and router. For search it is to look for trait hits (logically). For the router it is to fetch train info as well as dataset info.

In GN3 this table is used for partial correlations. Also to fetch API trait info and to build the search index.

In GN1 usage is similar.

# geno -> BIMBAM

We can use the script in gemma-wrapper

=> https://github.com/genetics-statistics/gemma-wrapper/blob/master/bin/gn-geno-to-gemma.py

there is probably something similar in GN2. And I have another version somewhere.

To identify the geno file the reaper script uses

```python
cursor.execute('select Id, Name from InbredSet')
results = cursor.fetchall()
InbredSets = {}
for item in results:
	InbredSets[item[0]] = genotypeDir+str(item[1])+'.geno'
```

which assumes one single geno file for the BXD that is indexed by the InbredSetID (a number). Note it ignores the many genotype files we have per inbredset (today). Also there is a funny hardcoded

```python
	if InbredSetId==3:
		InbredSetId=1
```

(no comment).

Later we'll output to lmdb when GEMMA supports it.

There are about 100 InbredSets. Genotype files can be found on production in
/export/guix-containers/genenetwork/var/genenetwork/genotype-files/genotype. For the BXD alone there are

```
BXD.2.geno               BXD-Heart-Metals_old.geno   BXD-Micturition.6.geno
BXD.4.geno               BXD-JAX-AD.4.geno           BXD-Micturition.8.geno
BXD.5.geno               BXD-JAX-AD.8.geno           BXD-Micturition.geno
BXD.6.geno               BXD-JAX-AD.geno             BXD-Micturition_old.4.geno
BXD.7.geno               BXD-JAX-AD_old.geno         BXD-Micturition_old.6.geno
BXD.8.geno               BXD-JAX-OFS.geno            BXD-Micturition_old.geno
BXD-AE.4.geno            BXD-Longevity.4.geno        BXD_mm8.geno
BXD-AE.8.geno            BXD-Longevity.8.geno        BXD-NIA-AD.4.geno
BXD-AE.geno              BXD-Longevity.9.geno        BXD-NIA-AD.8.geno
BXD-AE_old.geno          BXD-Longevity.array.geno    BXD-NIA-AD.geno
BXD-Bone.geno            BXD-Longevity.classic.geno  BXD-NIA-AD_old2.geno
BXD-Bone_orig.geno       BXD-Longevity.geno          BXD-NIA-AD_old.geno
BXD.geno                 BXD-Longevity_old.4.geno    BXD_Nov_23_2010_before_polish_101_102_103.geno
BXD-Harvested.geno       BXD-Longevity_old.8.geno    BXD_Nov_24_2010_before_polish_55_81.geno
BXD-Heart-Metals.4.geno  BXD-Longevity_old.geno      BXD_old.geno
BXD-Heart-Metals.8.geno  BXD-MBD-UTHSC.geno          BXD_unsure.geno
BXD-Heart-Metals.geno    BXD-Micturition.4.geno      BXD_UT-SJ.geno
```

Not really reflected in the DB:

```
MariaDB [db_webqtl]> select Id, Name from InbredSet where name like '%BXD%';
+----+------------------+
| Id | Name             |
+----+------------------+
|  1 | BXD              |
| 58 | BXD-Bone         |
| 64 | BXD-Longevity    |
| 68 | BXD_Dev          |
| 76 | DOD-BXD-GWI      |
| 84 | BXD-Heart-Metals |
| 86 | BXD-AE           |
| 91 | BXD-Micturition  |
| 92 | BXD-JAX-AD       |
| 93 | BXD-NIA-AD       |
| 94 | CCBXD-TM         |
| 96 | BXD-JAX-OFS      |
| 97 | BXD-MBD-UTHSC    |
+----+------------------+
```

Bit of a mess. Looks like some files are discarded. Let's see what the reaper script does.

We should also look into distributed storage. One option is webdav.

# Get PublishData trait(s) and convert to R/qtl2 or lmdb

Let's see how the scripts do it. Note that we already did that for the probeset script in

=> precompute-mapping-input-data

The code is reflected in

=> https://git.genenetwork.org/gn-guile/tree/scripts/precompute/list-traits-to-compute.scm

Now I need to do the exact same thing, but for PublishData.

Let's connect to a remote GN DB:

```
ssh -L 3306:127.0.0.1:3306 -f -N tux02.genenetwork.org
```

and follow

=> https://github.com/genenetwork/genenetwork2/blob/testing/scripts/maintenance/QTL_Reaper_v8_PublishXRef.py

the script takes a number of values 'PublishFreezeIds'. Alternatively it picks it up by SpeciesId (hard effing coded, of course).

=> https://github.com/genenetwork/genenetwork2/blob/fcde38b0f37f12508a01b16b7820029aa951bded/scripts/maintenance/QTL_Reaper_v8_PublishXRef.py#L62

Next it picks the geno file from the InbredSetID with

```
select InbredSetId  from PublishFreeze  where PublishFreeze.Id = 1;
```

Here we are initially going to focus on BXD=1 datasets only.

```
MariaDB [db_webqtl]> select Id,InbredSetId  from PublishFreeze  where InbredSetId = 1;
+----+-------------+
| Id | InbredSetId |
+----+-------------+
|  1 |           1 |
+----+-------------+
```

(we are half way the script now). Next we capture some metadata

```
MariaDB [db_webqtl]> select PhenotypeId, Locus, DataId, Phenotype.Post_publication_description from PublishXRef, Phenotype where PublishXRef.PhenotypeId = Phenotype.Id and InbredSetId=1 limit 5;
+-------------+----------------+---------+----------------------------------------------------------------------------------------------------------------------------+
| PhenotypeId | Locus          | DataId  | Post_publication_description                                                                                               |
+-------------+----------------+---------+----------------------------------------------------------------------------------------------------------------------------+
|           4 | rs48756159     | 8967043 | Central nervous system, morphology: Cerebellum weight, whole, bilateral in adults of both sexes [mg]                       |
|          10 | rsm10000005699 | 8967044 | Central nervous system, morphology: Cerebellum weight after adjustment for covariance with brain size [mg]                 |
|          15 | rsm10000013713 | 8967045 | Central nervous system, morphology: Brain weight, male and female adult average, unadjusted for body weight, age, sex [mg] |
|          20 | rs48756159     | 8967046 | Central nervous system, morphology: Cerebellum volume [mm3]                                                                |
|          25 | rsm10000005699 | 8967047 | Central nervous system, morphology: Cerebellum volume, adjusted for covariance with brain size [mm3]                       |
+-------------+----------------+---------+----------------------------------------------------------------------------------------------------------------------------+
```

it captures LRS

```
MariaDB [db_webqtl]> select LRS from PublishXRef where PhenotypeId=4 and InbredSetId=1;
+--------------------+
| LRS                |
+--------------------+
| 13.497491147108706 |
+--------------------+
```

and finally the trait values that are used for mapping

```
select Strain.Name, PublishData.value from Strain, PublishData where Strain.Id = PublishData.StrainId and PublishData.Id = 8967043;
+-------+-----------+
| Name  | value     |
+-------+-----------+
| BXD1  | 61.400002 |
| BXD2  | 49.000000 |
| BXD5  | 62.500000 |
| BXD6  | 53.099998 |
| BXD8  | 59.099998 |
| BXD9  | 53.900002 |
| BXD11 | 53.099998 |
| BXD12 | 45.900002 |
| BXD13 | 48.400002 |
| BXD14 | 49.400002 |
| BXD15 | 47.400002 |
| BXD16 | 56.299999 |
| BXD18 | 53.599998 |
| BXD19 | 50.099998 |
| BXD20 | 48.200001 |
| BXD21 | 50.599998 |
| BXD22 | 53.799999 |
| BXD23 | 48.599998 |
| BXD24 | 54.900002 |
| BXD25 | 49.599998 |
| BXD27 | 47.400002 |
| BXD28 | 51.500000 |
| BXD29 | 50.200001 |
| BXD30 | 53.599998 |
| BXD31 | 49.700001 |
| BXD32 | 56.000000 |
| BXD33 | 52.099998 |
| BXD34 | 53.700001 |
| BXD35 | 49.700001 |
| BXD36 | 44.500000 |
| BXD38 | 51.099998 |
| BXD39 | 54.900002 |
| BXD40 | 49.900002 |
| BXD42 | 59.400002 |
+-------+-----------+
```

Note that we need to filter out the parents - the original reaper script does not do that! My gn-guile code does handle that:

```
SELECT StrainId,Strain.Name FROM Strain, StrainXRef WHERE StrainXRef.StrainId = Strain.Id AND StrainXRef.InbredSetId =1 AND Used_for_mapping<>'Y' limit 5;
+----------+----------+
| StrainId | Name     |
+----------+----------+
|        1 | B6D2F1   |
|        2 | C57BL/6J |
|        3 | DBA/2J   |
|      150 | A/J      |
|      151 | AXB1     |
+----------+----------+
etc.
```

Also Bonz' script

=> https://git.genenetwork.org/gn-guile/tree/scripts/lmdb-publishdata-export.scm

has an interesting query:

```
MariaDB [db_webqtl]>
SELECT DISTINCT PublishFreeze.Name, PublishXRef.Id FROM PublishData
  INNER JOIN Strain ON PublishData.StrainId = Strain.Id
  INNER JOIN PublishXRef ON PublishData.Id = PublishXRef.DataId
  INNER JOIN PublishFreeze ON PublishXRef.InbredSetId = PublishFreeze.InbredSetId
  LEFT JOIN PublishSE ON PublishSE.DataId = PublishData.Id AND PublishSE.StrainId = PublishData.StrainId
  LEFT JOIN NStrain ON NStrain.DataId = PublishData.Id AND NStrain.StrainId = PublishData.StrainId
  WHERE PublishFreeze.public > 0 AND PublishFreeze.confidentiality < 1
  ORDER BY PublishFreeze.Id, PublishXRef.Id limit 5;
+------------+-------+
| Name       | Id    |
+------------+-------+
| BXDPublish | 10001 |
| BXDPublish | 10002 |
| BXDPublish | 10003 |
| BXDPublish | 10004 |
| BXDPublish | 10005 |
+------------+-------+
5 rows in set (0.239 sec)
```

that shows we have 13689 BXDPublish datasets. It also has

```
SELECT
JSON_ARRAYAGG(JSON_ARRAY(Strain.Name, PublishData.Value)) AS data,
 MD5(JSON_ARRAY(Strain.Name, PublishData.Value)) as md5hash
FROM
    PublishData
    INNER JOIN Strain ON PublishData.StrainId = Strain.Id
    INNER JOIN PublishXRef ON PublishData.Id = PublishXRef.DataId
    INNER JOIN PublishFreeze ON PublishXRef.InbredSetId = PublishFreeze.InbredSetId
LEFT JOIN PublishSE ON
    PublishSE.DataId = PublishData.Id AND
    PublishSE.StrainId = PublishData.StrainId
LEFT JOIN NStrain ON
    NStrain.DataId = PublishData.Id AND
    NStrain.StrainId = PublishData.StrainId
WHERE
    PublishFreeze.Name = "BXDPublish" AND
    PublishFreeze.public > 0 AND
    PublishData.value IS NOT NULL AND
    PublishFreeze.confidentiality < 1
ORDER BY
    LENGTH(Strain.Name), Strain.Name LIMIT 5;
```

best to pipe that to a file. It outputs JSON and an MD5SUM straight from mariadb. Interesting.

Finally, let's have a look at the existing GN API

```
SELECT
                            Strain.Name, Strain.Name2, PublishData.value, PublishData.Id, PublishSE.error, NStrain.count
                        FROM
                            (PublishData, Strain, PublishXRef, PublishFreeze)
                        LEFT JOIN PublishSE ON
                            (PublishSE.DataId = PublishData.Id AND PublishSE.StrainId = PublishData.StrainId)
                        LEFT JOIN NStrain ON
                            (NStrain.DataId = PublishData.Id AND
                            NStrain.StrainId = PublishData.StrainId)
                        WHERE
                            PublishXRef.InbredSetId = 1 AND
                            PublishXRef.PhenotypeId = 4 AND
                            PublishData.Id = PublishXRef.DataId AND
                            PublishData.StrainId = Strain.Id AND
                            PublishXRef.InbredSetId = PublishFreeze.InbredSetId AND
                            PublishFreeze.public > 0 AND
                            PublishFreeze.confidentiality < 1
                        ORDER BY
                            Strain.Name;
 +-------+-------+-----------+---------+-------+-------+
| Name  | Name2 | value     | Id      | error | count |
+-------+-------+-----------+---------+-------+-------+
| BXD1  | BXD1  | 61.400002 | 8967043 |  2.38 | NULL  |
| BXD11 | BXD11 | 53.099998 | 8967043 |   1.1 | NULL  |
| BXD12 | BXD12 | 45.900002 | 8967043 |  1.09 | NULL  |
| BXD13 | BXD13 | 48.400002 | 8967043 |  1.63 | NULL  |
...
```

which actually blocks non-public sets and shows std err, as well as counts when available(?) It does not exclude the parents for mapping (btw). That probably happens on the mapping page itself.

Probably the most elegant query is in GN3 API:

```
SELECT st.Name, ifnull(pd.value, 'x'), ifnull(ps.error, 'x'), ifnull(ns.count, 'x')
    FROM PublishFreeze pf JOIN PublishXRef px ON px.InbredSetId = pf.InbredSetId
        JOIN PublishData pd ON pd.Id = px.DataId JOIN Strain st ON pd.StrainId = st.Id
        LEFT JOIN PublishSE ps ON ps.DataId = pd.Id AND ps.StrainId = pd.StrainId
        LEFT JOIN NStrain ns ON ns.DataId = pd.Id AND ns.StrainId = pd.StrainId
    WHERE px.PhenotypeId = 4 limit 5;
+------+-----------------------+-----------------------+-----------------------+
| Name | ifnull(pd.value, 'x') | ifnull(ps.error, 'x') | ifnull(ns.count, 'x') |
+------+-----------------------+-----------------------+-----------------------+
| BXD1 | 61.400002             | 2.38                  | x                     |
| BXD2 | 49.000000             | 1.25                  | x                     |
| BXD5 | 62.500000             | 2.32                  | x                     |
| BXD6 | 53.099998             | 1.22                  | x                     |
| BXD8 | 59.099998             | 2.07                  | x                     |
+------+-----------------------+-----------------------+-----------------------+
```

written by Zach and Bonface. See

=> https://github.com/genenetwork/genenetwork3/blame/main/gn3/db/sample_data.py



We can get a list of the 13689 BXD datasets we can use. Note that we start with public data because we'll feed it to AI and all privacy will be gone after. We'll design an second API that makes use of Fred's authentication/authorization later.
Let's start with the SQL statement listed on:


We can run mysql through an ssh tunnel with

```
ssh -L 3306:127.0.0.1:3306 -f -N tux02.genenetwork.org
mysql -A -h 127.0.0.1 -uwebqtlout -pwebqtlout db_webqtl
```

and test the query, i.e.

```
MariaDB [db_webqtl]> SELECT DISTINCT PublishFreeze.Name, PublishXRef.Id FROM PublishData
    ->   INNER JOIN Strain ON PublishData.StrainId = Strain.Id
    ->   INNER JOIN PublishXRef ON PublishData.Id = PublishXRef.DataId
    ->   INNER JOIN PublishFreeze ON PublishXRef.InbredSetId = PublishFreeze.InbredSetId
    ->   LEFT JOIN PublishSE ON PublishSE.DataId = PublishData.Id AND PublishSE.StrainId = PublishData.StrainId
    ->   LEFT JOIN NStrain ON NStrain.DataId = PublishData.Id AND NStrain.StrainId = PublishData.StrainId
    ->   WHERE PublishFreeze.public > 0 AND PublishFreeze.confidentiality < 1
    ->   ORDER BY PublishFreeze.Id, PublishXRef.Id limit 5;
+------------+-------+
| Name       | Id    |
+------------+-------+
| BXDPublish | 10001 |
| BXDPublish | 10002 |
| BXDPublish | 10003 |
| BXDPublish | 10004 |
| BXDPublish | 10005 |
```

Let's take this apart a little. First of all PublishFreeze has only one record for BXDPublish where ID=1. PublishData may be used to check valid fields, but the real information is in PublishXRef. A simple

```
 select count(*) from PublishXRef WHERE InbredSetId=1;
+----------+
| count(*) |
+----------+
|    13711 |
+----------+
```

counts a few extra datasets (it was 13689). It may mean that PublishXRef contains some records that are still not public? Anyway,
let's go for the full dataset for precompute right now. We'll add an API endpoint to gn-guile so it can be used later.

Note GN2 on the menu search

=> https://genenetwork.org/search?species=mouse&group=BXD&type=Phenotypes&dataset=BXDPublish&search_terms_or=*&search_terms_and=&accession_id=None&FormID=searchResult

gives 13,729 entries, including recent BXD_51094. That is because that production database is newer. If we look at our highest records:

```
select * from PublishXRef WHERE InbredSetId=1 ORDER BY ID DESC limit 3;
+-------+-------------+-------------+---------------+----------+-------------------+----------------+--------------------+--------------------+----------+----------+
| Id    | InbredSetId | PhenotypeId | PublicationId | DataId   | mean              | Locus          | LRS                | additive           | Sequence | comments |
+-------+-------------+-------------+---------------+----------+-------------------+----------------+--------------------+--------------------+----------+----------+
| 51060 |           1 |       45821 |         39794 | 41022015 |              NULL | rsm10000000968 | 13.263934206457122 | 2.1741201177177185 |        1 |          |
| 51049 |           1 |       45810 |         39783 | 41022004 | 8.092333210508029 | rsm10000014174 |   16.8291804498215 | 18.143229769230775 |        1 |          |
| 51048 |           1 |       45809 |         39782 | 41022003 | 6.082199917286634 | rsm10000009222 | 14.462661474938166 |  4.582111488461538 |        1 |          |
+-------+-------------+-------------+---------------+----------+-------------------+----------------+--------------------+--------------------+----------+----------+
```

You can see they match that list (51060 got updated on production). The ID matches record BXD_51060 on the production search table.
We can look at the DataId with

```
select Id,PhenotypeId,DataId from PublishXRef WHERE InbredSetId=1 ORDER BY ID DESC limit 3;
+-------+-------------+----------+
| Id    | PhenotypeId | DataId   |
+-------+-------------+----------+
| 51060 |       45821 | 41022015 |
| 51049 |       45810 | 41022004 |
| 51048 |       45809 | 41022003 |
+-------+-------------+----------+
```

And get the actual values with

```
select * from PublishData WHERE Id=41022003;
+----------+----------+-----------+
| Id       | StrainId | value     |
+----------+----------+-----------+
| 41022003 |        2 |  9.136000 |
| 41022003 |        3 |  4.401000 |
| 41022003 |        9 |  4.360000 |
| 41022003 |       29 | 15.745000 |
| 41022003 |       98 |  4.073000 |
| 41022003 |       99 | -0.580000 |
```

which match the values on

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

The phenotypeid is useful for some metadata:


```
select * from Phenotype WHERE ID=45809;
| 45809 | Central nervous system, metabolism, nutrition, toxicology: Difference score for Iron (Fe) concentration in cortex (CTX) between 20 to 120-day-old and 300 to 918-day-old males mice fed Envigo diet 7912 containing 240, 93, and 63 ppm Fe, Cu and Zn, respectively [µg/g wet weight]  | Central nervous system, metabolism, nutrition, toxicology: Difference score for Iron (Fe) concentration in cortex (CTX) between 20 to 120-day-old and 300 to 918-day-old males mice fed Envigo diet 7912 containing 240, 93, and 63 ppm Fe, Cu and Zn, respectively [µg/g wet weight]  | Central nervous system, metabolism, nutrition, toxicology: Difference score for Iron (Fe) concentration in cortex (CTX) between 20 to 120-day-old and 300 to 918-day-old males mice fed Envigo diet 7912 containing 240, 93, and 63 ppm Fe, Cu and Zn, respectively [µg/g wet weight]  | [ug/mg wet weight] | Fe300-120CTXMale             | Fe300-120CTXMale              | NULL     | acenteno  | Jones B | joneslab         |
```

Since I am going for the simpler query I'll add an API endpoint named
datasets/bxd-publish/list (so others can use that too).  We'll return
tuples for each entry so we can extend it later. First we need the
DataID so we can point into PublishData. We expect the endpoint to
return something like

```
+-------+-------------+----------+
| Id    | PhenotypeId | DataId   |
+-------+-------------+----------+
| 51060 |       45821 | 41022015 |
| 51049 |       45810 | 41022004 |
| 51048 |       45809 | 41022003 |
...
```

Alright, let's write some code. The following patch returns on the endpoint:

```
[
  {
    "Id": 10001,
    "PhenotypeId": 4,
    "DataId": 8967043
  },
  {
    "Id": 10002,
    "PhenotypeId": 10,
    "DataId": 8967044
  },
  {
    "Id": 10003,
    "PhenotypeId": 15,
    "DataId": 8967045
  },
...
```

in about 3 seconds. It will run a lot faster on a local network. But for our purpose it is fine. The code I wrote is here:

=> https://git.genenetwork.org/gn-guile/commit/?id=1590be15f85e30d7db879c19d2d3b4bed201556a

Note the simple SQL query (compared to the first one).
Next step is to fetch the trait values we can feed to GEMMA. The full query using the PhenotypeId and DataId in GN is:

```
SELECT Strain.Name, Strain.Name2, PublishData.value, PublishData.Id, PublishSE.error, NStrain.count
  FROM
      (PublishData, Strain, PublishXRef, PublishFreeze)
  LEFT JOIN PublishSE ON
      (PublishSE.DataId = PublishData.Id AND PublishSE.StrainId = PublishData.StrainId)
  LEFT JOIN NStrain ON
      (NStrain.DataId = PublishData.Id AND
      NStrain.StrainId = PublishData.StrainId)
  WHERE
      PublishXRef.InbredSetId = 1 AND
      PublishXRef.PhenotypeId = 4 AND
      PublishData.Id = PublishXRef.DataId AND
      PublishData.StrainId = Strain.Id AND
      PublishXRef.InbredSetId = PublishFreeze.InbredSetId AND
      PublishFreeze.public > 0 AND
      PublishFreeze.confidentiality < 1;
+-------+-------+-----------+---------+-------+-------+
| Name  | Name2 | value     | Id      | error | count |
+-------+-------+-----------+---------+-------+-------+
| BXD1  | BXD1  | 61.400002 | 8967043 |  2.38 | NULL  |
| BXD2  | BXD2  | 49.000000 | 8967043 |  1.25 | NULL  |
| BXD5  | BXD5  | 62.500000 | 8967043 |  2.32 | NULL  |
| BXD6  | BXD6  | 53.099998 | 8967043 |  1.22 | NULL  |
...
```

(result includes parents). We can simplify this for GEMMA because it only wants the name and (mean) value.

The short version when you have the data ID is:

```
SELECT Strain.Name, PublishData.value FROM Strain, PublishData WHERE PublishData.Id=41022003 and Strain.Id=StrainID;
+----------+-----------+
| Name     | value     |
+----------+-----------+
| C57BL/6J |  9.136000 |
| DBA/2J   |  4.401000 |
| BXD9     |  4.360000 |
| BXD32    | 15.745000 |
| BXD43    |  4.073000 |
| BXD44    | -0.580000 |
| BXD48    | -1.810000 |
| BXD51    |  4.294000 |
| BXD60    | -0.208000 |
| BXD62    | -0.013000 |
| BXD63    |  3.221000 |
| BXD66    |  2.472000 |
| BXD69    | 12.886000 |
| BXD70    | -1.973000 |
| BXD78    | 19.511999 |
| BXD79    |  7.845000 |
| BXD73a   |  3.201000 |
| BXD87    | -3.054000 |
| BXD48a   | 11.585000 |
| BXD100   |  7.088000 |
| BXD102   |  8.485000 |
| BXD124   | 13.442000 |
| BXD170   | -1.274000 |
| BXD172   | 18.587000 |
| BXD186   | 10.634000 |
+----------+-----------+
```

which matches GN perfectly (some individuals where added). Alright, let's add an endpoint for this named
'dataset/bxd-publish/values/dataid/41022003'. Note we only deal with public data (so far). Later we may come up with more generic
end points and authorization. At this point the API is either on the local network (this one is) or public.

The first version returns this data from the endpoint:

```
time curl http://127.0.0.1:8091/dataset/bxd-publish/values/41022003
[{"Name":"C57BL/6J","value":9.136},{"Name":"DBA/2J","value":4.401},{"Name":"BXD9","value":4.36},{"Name":"BXD32","value":15.745},{"Name":"BXD43","value":4.073},{"Name":"BXD44","value":-0.58},{"Name":"BXD48","value":-1.81},{"Name":"BXD51","value":4.294},{"Name":"BXD60","value":-0.208},{"Name":"BXD62","value":-0.013},{"Name":"BXD63","value":3.221},{"Name":"BXD66","value":2.472},{"Name":"BXD69","value":12.886},{"Name":"BXD70","value":-1.973},{"Name":"BXD78","value":19.511999},{"Name":"BXD79","value":7.845},{"Name":"BXD73a","value":3.201},{"Name":"BXD87","value":-3.054},{"Name":"BXD48a","value":11.585},{"Name":"BXD100","value":7.088},{"Name":"BXD102","value":8.485},{"Name":"BXD124","value":13.442},{"Name":"BXD170","value":-1.274},{"Name":"BXD172","value":18.587},{"Name":"BXD186","value":10.634}]
real    0m0.537s
user    0m0.002s
sys     0m0.005s
```

Note it includes the parents. We should drop them. In this case we can simple check for (string-contains name "BXD"). The database records allow for a filter, so we get

```
curl http://127.0.0.1:8091/dataset/bxd-publish/mapping/values/41022003
[{"Name":"BXD9","value":4.36},{"Name":"BXD32","value":15.745},{"Name":"BXD43","value":4.073},{"Name":"BXD44","value":-0.58},{"Name":"BXD48","value":-1.81},{"Name":"BXD51","value":4.294},{"Name":"BXD60","value":-0.208},{"Name":"BXD62","value":-0.013},{"Name":"BXD63","value":3.221},{"Name":"BXD66","value":2.472},{"Name":"BXD69","value":12.886},{"Name":"BXD70","value":-1.973},{"Name":"BXD78","value":19.511999},{"Name":"BXD79","value":7.845},{"Name":"BXD73a","value":3.201},{"Name":"BXD87","value":-3.054},{"Name":"BXD48a","value":11.585},{"Name":"BXD100","value":7.088},{"Name":"BXD102","value":8.485},{"Name":"BXD124","value":13.442},{"Name":"BXD170","value":-1.274},{"Name":"BXD172","value":18.587},{"Name":"BXD186","value":10.634}]
```

That code went in as

=> https://git.genenetwork.org/gn-guile/commit/?id=9ad0793eb477611c700f4a5b02f60ac793bfae96

It took a bit longer than I wanted because I made a mistake converting the results to a hash table. It broke the JSON conversion and the error was not so helpful.

To write a CSV it turns out I have written

=> https://git.genenetwork.org/gn-guile/tree/gn/runner/gemma.scm?id=9ad0793eb477611c700f4a5b02f60ac793bfae96#n18

which takes the GN BXD.json file and our trait file. BXD.json captures the genotype information GN has:

```
{
        "mat": "C57BL/6J",
        "pat": "DBA/2J",
        "f1s": ["B6D2F1", "D2B6F1"],
        "genofile" : [{
                "title" : "WGS-based (Mar2022)",
                "location" : "BXD.8.geno",
                "sample_list" : ["BXD1", "BXD2", "BXD5", "BXD6", "BXD8", "BXD9", "BXD11", "BXD12", "BXD13", "BXD14", "BXD15", "BXD16", "BXD18", "BXD19", "BXD20", "BXD21", "BXD22", "BXD23", "BXD24", "BXD24a", "BXD25", "BXD27", "BXD28", "BXD29", "BXD30", "BXD31", "BXD32", "BXD33", "BXD34", "BXD35", "BXD36", "BXD37", "BXD38", "BXD39", "BXD40", "BXD41", "BXD42", "BXD43", "BXD44", "BXD45", "BXD48", "BXD48a", "BXD49", "BXD50", "BXD51", "BXD52", "BXD53", "BXD54", "BXD55", "BXD56", "BXD59", "BXD60", "BXD61",
(...)
"BXD065xBXD077F1", "BXD069xBXD090F1", "BXD071xBXD061F1", "BXD073bxBXD065F1", "BXD073bxBXD077F1", "BXD073xBXD034F1", "BXD073xBXD065F1", "BXD073xBXD077F1", "BXD074xBXD055F1", "BXD077xBXD062F1", "BXD083xBXD045F1", "BXD087xBXD100F1", "BXD065bxBXD055F1", "BXD102xBXD077F1", "BXD102xBXD73bF1", "BXD170xBXD172F1", "BXD172xBXD197F1", "BXD197xBXD009F1", "BXD197xBXD170F1"]
```

The code maps the traits values I generated against these columns to see what inviduals overlap which corrects for unmappable individuals (anyway).

The function 'write-pheno-file', listed above, does not work however because of the format of the endpoint. Remember it generates

```
[{"Name":"BXD9","value":4.36},{"Name":"BXD32","value":15.745}...]
```

While this function expects the shorter

```
{"BXD9":4.36,"BXD23":15.745...}
```

Now, for endpoints there is no real standard. We have written ideas up here:

=> https://git.genenetwork.org/gn-docs/tree/api

and, most recently

=> https://git.genenetwork.org/gn-docs/tree/api/GN-REST-API-v2.md

Where I make a case for having the metadata as a separate endpoint that can be reasoned on by people and machines (and AI).
That means I should default to the short version of the data and describe that layout using metadata. This we can do later.

I modified the endpoint to return the shorter hash:

```
time curl http://127.0.0.1:8091/dataset/bxd-publish/values/41022003
{"BXD9":4.36,"BXD23":15.745...}
```

Next, to align with

=> https://github.com/genenetwork/gn-docs/blob/master/api/GN-REST-API-v2.md

I gave the API the json extension, so we have http://127.0.0.1:8091/dataset/bxd-publish/values/41022003.json

This allows writing a special handler for GEMMA output (.gemma extension) downloading the pheno file with

```
curl http://127.0.0.1:8091/dataset/bxd-publish/values/41022003.gemma
NA
NA
NA
NA
NA
4.36NA
NA
NA
NA
(...)
``

that GEMMA can use directly and matches the order of the individuals in the BXD.8.geno file and the founders/parents are not included. Note that all of this now only works for the BXD (on PublishData) and I am using BXD.json as described in

=> https://issues.genenetwork.org/topics/systems/mariadb/precompute-mapping-input-data

I.e., it is Zach's listed stopgap solution. Code is here:

=> https://git.genenetwork.org/gn-guile/log/

Next step run gemma as we are on par with my earlier work on ProbeSetData. I wrote a gemma runner for that too at

=> https://git.genenetwork.org/gn-guile/tree/gn/runner/gemma.scm#n79

Now here I use guile to essentially script running GEMMA. There is no real advantage for that, so I will simply tell gemma-wrapper to use the output of above .gemma endpoint to fetch the trait values. Basically gemma-wrapper can specify the standard gemma -p switch, or pass in --phenotypes, that are used for permutations.

Now the new method we want to introduce is that the trait values are read from a REST API, instead of a file. The dirty way is to provide that functionality directly to gemma-wrapper, but we plan to get rid of that code (useful as it is -- it duplicates what Arun's ravanan does and ravanan has the advantage that it can be run on a cluster).

So we simply download the data and write it to a file with a small script. To run:

```
curl http://127.0.0.1:8091/dataset/bxd-publish/values/41022003.gemma > 41022003-pheno.txt
```

Next we create a container for gemma-wrapper (and includes the gemma that GN uses):

```
. .guix-deploy
env TMPDIR=tmp ruby ./bin/gemma-wrapper --force --json \
        --loco -- \
        -g BXD.8_geno.txt.gz \
        -p 41022003-pheno.txt \
        -a BXD.8_snps.txt \
        -gk > K.json
```

this bailed out with

Executing: parallel --results /tmp/test --joblog /tmp/test/5f3849a9e61b70e3d562b20c5eade5a699923c68-parallel.log < /tmp/test/parallel-commands.txt
Command exited with non-zero status 20

When running an individual chromosome (from the parallel log) we get two warnings and an error:

```
**** WARNING: The maximum genotype value is not 2.0 - this is not the BIMBAM standard and will skew l_lme and effect sizes
**** WARNING: Columns in geno file do not match # individuals in phenotypes
ERROR: Enforce failed for not enough genotype fields for marker in src/gemma_io.cpp at line 1470 in BimbamKin
```

Looks familiar!
The first warning we'll ignore for now, as we just want the hits initially. The second warning relates to the error that there is a mismatch in number of inds.

This topic I have covered in the past, particularly trying to debug Dave's conflicting results:

=> https://issues.genenetwork.org/topics/lmms/gemma/permutations

It makes somewhat depressive reading, though we have a solution.

Note the correct conversion we only have to do once (basically the code I wrote earlier
to fetch BXD traits needs to work with the latest BXD genotypes).
The real problem is that gemma itself does not compare individual names (at all), so any corrections need to be done beforehand. In this case our pheno file contains 212 inds from the earlier BXD.json file.

```
wc -l 41022003-pheno.txt
212 41022003-pheno.txt
```

And that is off. Let's try the tool I wrote during that exercise. It can create a different json file after parsing BXD.geno
that has in the header:

> # Date Modified: April 23, 2024 by Arthur Centeno, Suheeta Roy. March 22, 2022 by Rob Williams, David Ashbrook, and Danny Arends to remove excessive cross-over events in strains BXD42 (Chr9), BXD81 (Chrs1, 5, 10), BXD99 (Chr1), and BXD100 (Chrs2 and 6); and to add Taar1 maker on Chr 10 for T. Phillips-Richards.   Jan 19, 2017: Danny Arends computed BXD cM values and recombinations between markers. Rob W. Williams fixed errors on most chromosomes and added Affy eQTL markers. BXD223 now has been added based on David Ashbrook's spreadsheet genotype information.

```
md5sum BXD.geno:
  a78aa312b51ac15dd8ece911409c5b98  BXD.geno
gemma-wrapper$ ./bin/gn-geno-to-gemma.py BXD.geno > BXD.geno.txt
```

creates a .json file (that is different from Zach/GN's) and a bimbam file GEMMA can use. Now in the next step I need to adapt above code to use this format. What I *should* have done, instead of writing gemma phenotypes directly, is write the R/qtl2 format that includes the ind names (so we can compare and validate those) and *then* parse that data against our new JSON file created by gn-geno-to-gemma.py using the rqtl2-pheno-to-gemma.py script. Both Python scripts are already part of gemma-wrapper:

=> https://github.com/genetics-statistics/gemma-wrapper/blob/master/bin/gn-geno-to-gemma.py
=> https://github.com/genetics-statistics/gemma-wrapper/blob/master/bin/rqtl2-pheno-to-gemma.py

The idea was to create the rqtl2 API endpoint, or I'll adapt the 2nd script to take the endpoint as input and then correct for GEMMA's requirements.

OK, updated the endpoints and the code for rqtl2-pheno-to-gemma.py so it accepts a URL instead of a file. So the idea is
to run

```
./bin/rqtl2-pheno-to-gemma.py BXD_pheno_Dave.csv --json BXD.geno.json > BXD_pheno_matched.txt
```

A line in BXD_pheno_Dave.csv is:

```
BXD113,24.52,205.429001,3.643,2203.312012,3685.907959,1.199,2.019,29.347143,0.642857,205.428574,24.520409,3.642857,2203
.312012,3685.908203,1.198643,2.018643,0.642857,33.785709,1.625,2,1.625,1,22.75
```

Now if I read the Rqtl2 docs it says:

> We split the numeric phenotypes from the mixed-mode covariates, as two separate CSV files. Each file forms a matrix of individuals × phenotypes (or covariates), with the first column being individual IDs and the first row being phenotype or covariate names. Sex and line IDs (if needed) can be columns in the covariate data.

This differs from the BXD Dave layout (it is transposed). Karl added in the docs:

> All of these CSV files may be transposed relative to the form described below. You just need to include, in the control file, a line like: "geno_transposed: true". So, OK, we can use the transposed form. First we make it possible to parse json:

```
curl http://127.0.0.1:8091/dataset/bxd-publish/values/41022003.json > 41022003-pheno.json
jq < 41022003-pheno.json
{
  "C57BL/6J": 9.136,
  "DBA/2J": 4.401,
  "BXD9": 4.36,
  "BXD32": 15.745,
(...)
```

note it includes the parents. Feed it to

```
./bin/rqtl2-pheno-to-gemma.py 41022003-pheno.json --json BXD.geno.json
```

where BXD.geno.json is not the Zach/GN json file, but the actual BXDs in GEMMA's bimbam file.

One question is why Zach's JSON file gives a different number of mappable BXDs. I made of note of that to check.

I wrote a new script and we had our first GEMMA run with lmdb output:

```
wrk@napoli /export/local/home/wrk/iwrk/opensource/code/genetics/gemma-wrapper [env]$ tar tvf /tmp/3fddda2374509c7b346>
-rw-r--r-- wrk/users    294912 2025-08-06 05:49 3fddda2374509c7b346b7819ae358ed23be9cb46-gemma-GWA.mdb
```

The script is just 10 lines of code (after the command line handler)

=> https://github.com/genetics-statistics/gemma-wrapper/blob/master/bin/gn-pheno-to-gemma.rb

Excellent, now we can run gemma and the next step is to look at the largest hit.

So the trait we try to run is 41022003 = https://genenetwork.org/show_trait?trait_id=51048&dataset=BXDPublish. The inputs match up. When we run GEMMA in GN it has a 4.0 score on chr 12 and 3.9 on chr 19.

Running gemma-wrapper we get

```
LOCO K computation with caching and JSON output

gemma-wrapper --json --force --loco -- -g BXD.geno.txt -p BXD_pheno.txt -a BXD.8_snps.txt -n 2 -gk -debug > K.json

LMM's using the K's captured in K.json using the --input switch

gemma-wrapper --json --force --lmdb --loco --input K.json -- -g BXD.geno.txt -p BXD_pheno.txt -a BXD.8_snps.txt -lmm 9 -maf 0.1 -n 2 -debug > GWA.json
```

We can view the lmdb file with something like:

```
./bin/view-gemma-mdb --sort /tmp/66b8c19be87e9566358ce904682a56250eb05748-gemma-GWA.tar.xz --anno BXD.8_snps.txt > test.out
/tmp/3fddda2374509c7b346b7819ae358ed23be9cb46-gemma-GWA.tar.xz
chr,pos,marker,af,beta,se,l_mle,l_lrt,-logP
7,67950073,rsm10000004928,0.543,1.5226,1.3331,100000.0,0.0002,3.79
7,68061665,rs32453663,0.543,1.5226,1.3331,100000.0,0.0002,3.79
7,68111284,rs32227186,0.543,1.5226,1.3331,100000.0,0.0002,3.79
19,30665443,rsm10000014129,0.522,2.2128,1.0486,100000.0,0.0002,3.77
19,30671753,rs31207057,0.522,2.2128,1.0486,100000.0,0.0002,3.77
12,40785621,rsm10000009222,0.565,2.8541,1.3576,100000.0,0.0002,3.75
12,40786657,rs29124638,0.565,2.8541,1.3576,100000.0,0.0002,3.75
12,40842857,rs13481410,0.565,2.8541,1.3576,100000.0,0.0002,3.75
12,40887762,rsm10000009223,0.565,2.8541,1.3576,100000.0,0.0002,3.75
12,40887894,rsm10000009224,0.565,2.8541,1.3576,100000.0,0.0002,3.75
12,40900825,rs50979658,0.565,2.8541,1.3576,100000.0,0.0002,3.75
12,41054766,rs46705481,0.565,2.8541,1.3576,100000.0,0.0002,3.75
```

Interestingly the hits are very similar to what is on production now, though not the same! That points out that I am not using the production database on this recent dataset. Let's try an older one. BXD_10002 has data id 8967044

```
curl http://127.0.0.1:8091/dataset/bxd-publish/values/8967044.json > 10002-pheno.json
./bin/gn-pheno-to-gemma.rb -p 10002-pheno.json --geno-json BXD.geno.json > 10002-pheno.txt
gemma-wrapper --json --force --loco -- -g BXD.geno.txt -p 10002-pheno.txt -a BXD.8_snps.txt -n 2 -gk -debug > K.json
gemma-wrapper --json --force --lmdb --loco --input K.json -- -g BXD.geno.txt -p 10002-pheno.txt -a BXD.8_snps.txt -lmm 9 -maf 0.1 -n 2 -debug > GWA.json
./bin/view-gemma-mdb --sort /tmp/c4ffedf358698814c6e29a54a2a51cb6c66328d0-gemma-GWA.tar.xz --anno BXD.8_snps.txt > test.out
```

Luckily this is a perfect match:

```
1,179861787,rsm10000000444,0.559,0.8837,0.3555,100000.0,0.0,4.99
1,179862838,rs30712622,0.559,0.8837,0.3555,100000.0,0.0,4.99
1,179915631,rsm10000000787,0.559,0.8837,0.3555,100000.0,0.0,4.99
1,179919811,rsm10000000788,0.559,0.8837,0.3555,100000.0,0.0,4.99
(...)
8,94479237,rs32095272,0.441,1.0456,0.4362,100000.0,0.0,4.75
8,94765445,rsm10000005684,0.441,1.0456,0.4362,100000.0,0.0,4.75
8,94785223,rsm10000005685,0.441,1.0456,0.4362,100000.0,0.0,4.75
8,94840921,rsm10000005686,0.441,1.0456,0.4362,100000.0,0.0,4.75
```

The lmdb file contains the full vector and compresses to 100K. For 13K traits that equals about 1Gb.

First I wanted to check how Zach's list of mappable inds compares to mine. A simple REPL exercise shows:

```
zach = JSON.parse(File.read('BXD.json'))
pj = JSON.parse(File.read('BXD.geno.json'))
s1 = zach["genofile"][0]["sample_list"]
=> ["BXD1", "BXD2", "BXD5", "BXD6", "BXD8", "BXD9", "BXD11", "BXD12", "BXD13", "BXD14", "BXD15", "BXD16", "BXD18",...
s2 = pj["samples"]
=> ["BXD1", "BXD2", "BXD5", "BXD6", "BXD8", "BXD9", "BXD11", "BXD12", "BXD13", "BXD14", "BXD15", "BXD16", "BXD18",...
s1.size()
=> 235
s2.size()
=> 237
 s2-s1
=> ["BXD077xBXD065F1", "BXD065xBXD102F1"]
```

So it turns out the newer geno file contains these two new inds that are *also* in the .geno file and confuses the hell out of my scripts ;). The GN2 webserver probably uses the header of the geno file to fetch the correct number. The trait page also lists these inds, so (I guess) the BXD.json file ought to be updated.

Now that is explained and we are good.

## Running at scale

In the next step we need to batch run GEMMA. Initially we'll run on one server. gemma-wrapper takes care of running only once, so we can restart the pipeline at any point (we'll move to ravanan after to run on the cluster). At this point the API uses the dataid to return the trait values. I think that is not so intuitive, so I modified the endpoint to give the same results for:

```
curl http://127.0.0.1:8091/dataset/bxd-publish/values/10002.json > 10002-pheno.json
curl http://127.0.0.1:8091/dataset/bxd-publish/dataid/values/8967044.json > 10002-pheno.json
```

Now that works we can get a list of all BXDPublish datasets that I wrote earlier:

```
curl http://127.0.0.1:8091/dataset/bxd-publish/list > bxd-publish.json
[
  {
    "Id": 10001,
    "PhenotypeId": 4,
    "DataId": 8967043
  },
  {
    "Id": 10002,
    "PhenotypeId": 10,
    "DataId": 8967044
  },
  {
    "Id": 10003,
    "PhenotypeId": 15,
    "DataId": 8967045
  },
```

so we can use this to create our batch list. There are 13711 datasets listed on this DB. We can use jq to extract all Ids

```
jq ".[] | .Id" < bxd-publish.json > ids.txt
```

All set to run our first batch! Now we replicate our guix-wrapper environment, start the gn-guile server and fire up a batch script that pulls the data from the database and runs gemma for every step.


To get precompute going we need a server set up with a recent database. I don't want to use the production server. The fastest other server we have is balg01, and it is not busy right now, so let's use that. First we recover a DB from our backup, as described in

=> topics/systems/mariadb/precompute-mapping-input-data

(btw that examples show we started on precompute since November 2023 - 1.5 years ago). On that server mariadb is running as
/usr/local/guix-profiles/gn-latest/bin/mariadbd --datadir=/export/mariadb/tux01. We can simply overwrite that database as it
is an installation of Feb 18 2024. We extract:

```
borg extract --progress /export/backup/bacchus/drop/tux04/genenetwork::borg-tux04-sql-20250807-04:16-Thu
```

After extracting the backup we need to update permissions and point mariadb to the new dir: balg01:/export/mariadb/tux04/latest/.
Restarting the DB and it all appears to work.

Before I move the code across we need to make sure metadata on the traits get added to the lmdb mapping data. I actually wrote the code for that here. This adds the metadata to lmdb:

=> https://github.com/genetics-statistics/gemma-wrapper/blob/a0eb8ed829072cb539b32affe135a7930989ca30/bin/gemma2lmdb.py#L99

gemma-wrapper writes data like this:

```
  "meta": {
    "type": "gemma-wrapper",
    "version": "0.99.7-pre1",
    "population": "BXD",
    "name": "HC_U_0304_R",
    "trait": "101500_at",
    "url": "https://genenetwork.org/show_trait?trait_id=101500_at&dataset=HC_U_0304_R",
    "archive_GRM": "46bfba373fe8c19e68be6156cad3750120280e2e-gemma-cXX.tar.xz",
    "archive_GWA": "779a54a59e4cd03608178db4068791db4ca44ab3-gemma-GWA.tar.xz",
    "dataid": 75629,
    "probesetid": 1097,
    "probesetfreezeid": 7
    }
```

This was done for probesetdata and needs to be adapted for our BXD PublishData exercise. Also I want the archive_GWA file name to include the trait name/ID so we can find it quickly on the storage (without having to parse/query all lmdb files).

From the gemma-wrapper invocation you can see I added a few switches to pass in this information:

=> https://git.genenetwork.org/gn-guile/tree/gn/runner/gemma.scm#n97

```
        --meta NAME                  Pass in metadata as JSON file
        --population NAME            Add population identifier to metadata
        --name NAME                  Add dataset identifier to metadata
        --id ID                      Add identifier to metadata
        --trait TRAIT                Add trait identifier to metadata
```

We can add BXD as population and BXDPublish as a dataset identifier. Set id with dataid, and trait id with PublishXRefID and point it back to GN, so we can click

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

Another thing I want to add are the existing qtlreaper hit values.  That way we can assess where the biggest impact was of using gemma over qtlreaper. To achieve this we will create a new API endpoint that can serve that data. Remember we get the trait values with:

=> http://127.0.0.1:8091/dataset/bxd-publish/values/10002.json

so we can add an endpoint that lists the mapping results

=> http://127.0.0.1:8091/dataset/bxd-publish/trait-hits/10002.json

we also will have

=> http://127.0.0.1:8091/dataset/bxd-publish/trait-info/10002.json

That will return more metadata and point into our RDF store. Note that this is now all very specific to bxd-publish. Later we'll have to think how to generalise these endpoints. We are just moving forward to do the BXD precompute run.

Interestingly GN2 shows this information (well, only the highest hit) on the search page, but not on the trait page. As we can get hits from multiple sources we should (eventually) account for that with something like:

```
=> http://127.0.0.1:8091/dataset/bxd-publish/trait-hits/10002.json
{ "qtlreaper-hk":
  {
    [
      { "name":..., "chr": ..., "pos":..., "LRS":..., "additive":..., }
    ]
  }
  "gemma-loco":
  {
    [
      { "name":..., "chr": ..., "pos":..., "LRS":..., "additive":..., }
      { "name":..., "chr": ..., "pos":..., "LRS":..., "additive":..., }
      { "name":..., "chr": ..., "pos":..., "LRS":..., "additive":..., }
    ]
  }
}
```

Eventually we may list gemma, Rqtl2 hits with and without LOCO and with and without covariates. Once we build this support we can adapt our search tools.

Obviously this won't fit the current PublishXRef format, so -- for now -- we will just mirror its contents:

```
{ "qtlreaper-hk":
  {
    [
      { "name":..., "chr": ..., "pos":..., "LRS":..., "additive":..., }
    ]
  }
}
```

To get compute going I am going to skip above because we can update the lmdb files later.
The first fix is to add the trait name to the file names and the following record to lmdb:

  "meta": {
    "type": "gemma-wrapper",
    "version": "0.99.7-pre1",
    "population": "BXD",
    "name": "BXDPublish",
    "table": "PublishData",
    "traitid": 10002, // aka PublishXrefId
    "url": "https://genenetwork.org/show_trait?trait_id=51048&dataset=BXDPublish,
    "archive_GRM": "46bfba373fe8c19e68be6156cad3750120280e2e-gemma-cXX.tar.xz",
    "archive_GWA": "779a54a59e4cd03608178db4068791db4ca44ab3-BXDPublish-10002-gemma-GWA.tar.xz",
    "dataid": 8967044,
    }

This required modifications to gemma-wrapper.

Running:

```
gemma-wrapper --json --force --loco -- -g BXD.geno.txt -p BXD_pheno.txt -a BXD.8_snps.txt -n 2 -gk -debug > K.json
gemma-wrapper --json --force --lmdb --population BXD --name BXDPublish --trait 10002 --loco --input K.json -- -g BXD.geno.txt -p BXD_pheno.txt -a BXD.8_snps.txt -lmm 9 -maf 0.1 -n 2 -debug > GWA.json
```

begets '66b8c19be87e9566358ce904682a56250eb05748-BXDPublish-10002-gemma-GWA.tar.xz'. When I check the meta data in the lmdb file it is set to

```
"meta": {"type": "gemma-wrapper", "version": "1.00-pre1", "population": "BXD", "name": "BXDPublish", "trait": "10002", "geno_filename": "BXD.geno.txt", "geno_hash": "3b65ed252fa47270a3ea867409b0bdc5700ad6f6", "loco": true, "url": "https://genenetwork.org/show_trait?trait_id=10002&dataset=BXDPublish", "archive_GRM": "185eb08dc3897c7db5d7ea987170898035768f93-gemma-cXX.tar.xz", "archive_GWA":"66b8c19be87e9566358ce904682a56250eb05748-BXDPublish-10002-gemma-GWA.tar.xz", "table": "PublishData", "traitid": 10002, "dataid": 0}
```

which is good enough (for now). I may still add the dataid, but it requires a SQL call. Code is here:

=> https://github.com/genetics-statistics/gemma-wrapper/commit/49587523fc93bdcf0265da9da97f8d6d2a9e1008

I should note that up to this point I would have had no advantage from AI programming. I know there are topics that I'll work on where I may benefit, but this type of architecturing, with very little code writing, does not really help. I certainly have the intention of using AI! Next steps, unfortunately, there is still little to be gained. Where we'll probably gain is:

- Using the RDF data store and documenting the endpoint(s)
- Refactoring some of GN2's code to introduce lmdb\
- Deduplicating GN2/GN3 SQL code
- Improving the REST API and writing documentation and tests
- Analysing existing code bases, such as GEMMA itself

Next step is getting the data churn going! After that we'll list all the hits which requires processing the lmdb output.

Precompute of 13K traits has its first test run on balg01.

It is going at 30 gemma runs per minute, so perhaps 8 hours for the full run if it keeps going. But I am hitting errors.

Afther that will be to digest hits from the precomputed vectors in lmdb.

## Yesterday's tux02 crash

All servers work on tux02 except for BNW.

I tried to restart BNW, but it is giving an error, including the mystifying shepherd error (that I have as a sticker on my laptop):

> 2025-08-11 01:13:41 error in finalization thread: Success

It is on our end, so no need to ping Yan. I'll fix it when I have time (I did below).

## Precompute

To get precompute up and running I need to create the environment on balg01. The DB I updated a few days ago, so that should be fine.

First we check out the guile webserver:

```
git clone tux02.genenetwork.org:/home/git/public/gn-guile gn-guile-8092
```

Now gn-guile is already running serving aliases, so we want to run this as an internal endpoint right now with something like

```
unset GUIX_PROFILE
. /usr/local/guix-profiles/guix-pull/etc/profile
guix shell -L ~/guix-bioinformatics --container --network --file=guix.scm -- guile -L . --fresh-auto-compile -e main web/webserver.scm 8092
```

so, this renders

```
curl http://127.0.0.1:8092/dataset/bxd-publish/values/10002.json
{"BXD1":54.099998,"BXD2":50.099998,"BXD5":53.299999,"BXD6":55.099998
```

Next step is to set up gemma-wrapper. Now this failed because guix was not happy. We have been updating things these last weeks. Rather than trying to align with recent changes I could have rolled back to the version I am using on my desktop. But I decided not to let those bits rot and updated guix from

guix describe Thu Mar 14 21:33:55 2024

to

guix describe Sun Aug 10 18:18:20 2025

Should use a newer version first! Let's try

```
guix pull --url=https://codeberg.org/guix/guix  -p ~/opt/guix-pull
```

(that took a while, so I took the opportunity to fix BNW -- turns out someone disabled BNW in shepherd by creating a systemd version that did not start properly).

After the pull there were quite a few problems with gemma dependencies that needed fixing. First problem

```
guix package: warning: failed to load '(gn packages gemma)':
In procedure abi-check: #<record-type <git-reference>>: record ABI mismatch; recompilation needed
```

required

```
find ~/.cache/guile -name "*.go" -delete
```

I also had to point guix-past to the new codeberg record! And now, magically, things started working.

So, now I have an identical setup on my desktop and on the balg server. Next is to write a script that will batch run gemma-wrapper for every BXD PublishData ID. We created that list with jq earlier.

```
curl http://127.0.0.1:8092/dataset/bxd-publish/list > bxd-publish.json
jq ".[] | .Id" < bxd-publish.json > ids.txt
```

For every ID in that list we are going to fetch the trait values with

```
#! /bin/env sh

export TMPDIR=./tmp
curl http://127.0.0.1:8092/dataset/bxd-publish/list > bxd-publish.json
jq ".[] | .Id" < bxd-publish.json > ids.txt
./bin/gemma-wrapper --force --json --loco -- -g BXD.geno.txt -p BXD_pheno.txt -a BXD.8_snps.txt -n 2 -gk > K.json

for id in `cat ids.txt` ; do
  echo Precomputing $id
  curl http://127.0.0.1:8092/dataset/bxd-publish/values/$id.json > pheno.json
  ./bin/gn-pheno-to-gemma.rb --phenotypes pheno.json --geno-json BXD.geno.json > BXD_pheno.txt
  ./bin/gemma-wrapper --json --lmdb --population BXD --name BXDPublish --trait $id --loco --input K.json -- -g BXD.geno.txt -p BXD_pheno.txt -a BXD.8_snps.txt -lmm 9 -maf 0.1 -n 2 -debug > GWA.json
done
```

I hard copied the following files

```
BXD.geno.json
BXD.geno.txt
BXD.8_snps.txt
```

One thing I need to check is that the GRM is actually a constant. I forgot what GEMMA does.

We hit an error

```
/gnu/store/vvl1g1l0j19w39kry2xcsawvlhbyb87j-ruby-3.4.4/lib/ruby/3.4.0/json/common.rb:221:in 'JSON::Ext::Parser.parse':
unexpected token at '' (JSON::ParserError)
FATAL ERROR: gemma-wrapper bailed out with pid 340588 exit 20
./bin/gemma-wrapper:494:in 'block (2 levels) in <main>'
./bin/gemma-wrapper:479:in 'IO.open'
./bin/gemma-wrapper:479:in 'block in <main>'
./bin/gemma-wrapper:832:in '<main>'Precomputing 10137
```

The JSON file is empty 10136. Hmmm.

I also see

```
WARNING: failed to update lmdb record with key b'\r\x02n\x7f\x10' -- probably a duplicate 13:40795920 (b'\r':40795920)
```

For the first the webserver actually stopped on `In procedure accept: Too many open files`. The problem looks similar to

=> https://issues.guix.gnu.org/60226

and Arun's patch

=> https://cgit.git.savannah.gnu.org/cgit/guix/mumi.git/commit/?id=897967a84d3f51da2b1cc8c3ee942fd14f4c669b

I raised ulimit, but may need to restart the webserver several time. We are computing though:

```
-rw-r--r-- 1 wrk wrk   82968 Aug 11 05:16 ab51d69f79601cfa7399feebca619ea1a71c1270-BXDPublish-10146-gemma-GWA.tar.xz
-rw-r--r-- 1 wrk wrk   82772 Aug 11 05:16 e6739ace8ca4931fc51baa1844b3b5ceac592104-BXDPublish-10147-gemma-GWA.tar.xz
-rw-r--r-- 1 wrk wrk   81848 Aug 11 05:16 60880fc7e8c86dffb17f28664e478204ea26f827-BXDPublish-10148-gemma-GWA.tar.xz
-rw-r--r-- 1 wrk wrk   79336 Aug 11 05:16 c914d6221b004dec98d60e08c0fdf8791c09cb41-BXDPublish-10149-gemma-GWA.tar.xz
-rw-r--r-- 1 wrk wrk   83536 Aug 11 05:16 3d72b19730edab29bdc593cb6a1a86dd789d351f-BXDPublish-10150-gemma-GWA.tar.xz
-rw-r--r-- 1 wrk wrk   69060 Aug 11 05:16 0e965f1778425071a5497d0fe69f2dc2e534ef60-BXDPublish-10151-gemma-GWA.tar.xz
-rw-r--r-- 1 wrk wrk   69072 Aug 11 05:16 4de26e62a75727bc7edd6b266dfcd7753d185f1a-BXDPublish-10152-gemma-GWA.tar.xz
(...)
```

There are some scarily small datasets:

```
GET /dataset/bxd-publish/values/10198.json
;;; ("8967240")

;;; ((("C57BL/6J" . 1.62) ("BXD1" . 2.37) ("BXD5" . 2.73) ("BXD9" . 3.52) ("BXD11" . 0.18) ("BXD12" . 3.69) ("BXD16" . 0.29) ("BXD21" . 2.34) ("BXD27" . 3.38) ("BXD32" . 0.24)))
```

i.e. https://genenetwork.org/show_trait?trait_id=10198&dataset=BXDPublish

Not sure we should be running GEMMA on those!