Bound populateKinship.r peak memory and remove the kinship2 dependency

[ankurjuneja]

Rationale

Bound the memory the nightly kinship calculation charges against the web server's JVM budget, and take kinship2 out of the coefficient calculation.

The R process runs as a child of the web server, so its peak counts against the same budget as the JVM and has been driving instance sizing. Computing one family at a time and appending its rows to disk bounds the peak by the largest family rather than the whole colony. Running the recursion in the script removes what remains: the duplicate matrix kinship2 allocates for every family, and the placeholder parent rows it invents for animals with only one known parent, which at one colony are 26,529 of 66,848 rows in a single species.

Measured against five production pedigrees between 2,773 and 62,757 animals, peak memory falls 2.5x to 3.8x. The largest colony cannot be completed by the current script on a 48 GB machine at all; it now finishes in 16 seconds at 11 GB.

Output is unchanged. On every pedigree the current script can process, the rows imported into ehr.kinship are identical.

Removing kinship2 from the calculation also closes a source of cross-server inconsistency: install.r pins no versions, and the same script run against two different kinship2 builds produces coefficients that differ in their last bits, which was the original diagnosis behind Issue 47002.

Related Pull Requests

None.

Changes

• Compute kinship one family at a time and append each family's rows to disk, instead of accumulating the whole colony in memory and sorting at the end.
• Run the kinship recursion in the script rather than calling kinship2, avoiding both the duplicate allocation of the largest matrix and the placeholder parent rows the library invents.
• Write to a temporary file and rename it on success, so a failed run cannot leave a partial file for the importer to load over good data, and remove that temporary file when a run fails.
• Stop writing self-pairs and the species column, neither of which the importer uses.
• Validate the pedigree before computing so a contradictory record names the animal at fault.
• Warn when a parent's recorded gender contradicts the role it is used in, which the previous implementation silently corrected and never reported.
• Stop rejecting a species whose animals all have unknown gender, which autosomal kinship does not use.

[bbimber]

@ankurjuneja and @labkey-martyp: have you confirmed this produces identical results using the real data from production NPRCs?

[ankurjuneja]

@ankurjuneja and @labkey-martyp: have you confirmed this produces identical results using the real data from production NPRCs?

Yes, tested with data from three centers and output is identical.

[bbimber]

@ankurjuneja and @labkey-martyp: have you confirmed this produces identical results using the real data from production NPRCs?

Yes, tested with data from three centers and output is identical.

OK, thanks for confirming. All of this is quite old code and it would not surprise me if there was big room for efficiency both here, and at import time.

[labkey-martyp]

@ankurjuneja and @labkey-martyp: have you confirmed this produces identical results using the real data from production NPRCs?

Yes, tested with data from three centers and output is identical.

@bbimber if you want to send us study.Pedigree export and your kinship.txt, we can test yours as well.

[bbimber]

@ankurjuneja and @labkey-martyp: have you confirmed this produces identical results using the real data from production NPRCs?

Yes, tested with data from three centers and output is identical.

@bbimber if you want to send us study.Pedigree export and your kinship.txt, we can test yours as well.

@ankurjuneja wrote that this was tested on three centers - was ONPRC not one of them?

[bbimber]

If you're touching this code, the first message I got from R was that kinship2 is deprecated in favor of this, which seems like it might be a drop-in replacement: https://louislenezet.github.io/Pedixplorer/

[bbimber]

@ankurjuneja: I ran this on our data. Some high level observations first:

• Some of the refactors, like changing filtering style, make a lot of sense (again, this was written >15 years ago).
• One of the ways this brings efficiency is using the makefam() function to split families, rather than process the entire species at once. On our production data, in practice this basically means one massive family/species; however, it drops some barely related subjects, which reduces the max matrix size a little. The result is that the output omits 20-30K coefficient=0 records. It might be fine to break that expectation, but if we're going to do that, we should get the maximum value from it. More below.
• When this was originally written, makefam() didnt seem to make a major difference on runtime, and there is the possibility of it splitting families incorrectly, which could result in missing values. This may not matter in practice.
• Piping tidyr::pivot_longer() and dplyr::filter() and using sparse matrices might accomplish the same matrix melting operation in an efficient manner with less custom code, but that's nothing wrong with the change per se.

Specific suggestions:

• If we want to bring efficiency to both R script memory and import time, we should have a conversation on how to handle rows where coefficient=0. In the current form, the script writes, and the SQL table stores, one row for each pair within a species, by design. In most species, the majority of these rows are zeros. It would probably help R memory and definitely reduce the SQL import time to only report rows where coefficient>0. Consumers of these data would need to make the assumption that lack of row equals zero kinship. This is a change, but I cant currently think of a reason we couldnt make that change in behavior. This script has a MIN_COEFFICIENT variable that doesnt appear to be used, which might have been intended to go that direction.

• If we are going to use makefam in R, I also wonder if it would be useful to write a TSV mapping SubjectId->FamilyId with each run in R. This would at least preserve an artifact of what IDs were calculated together if a question arises. We could consider storing the family ID as a new column in the kinship table; however, that's probably more than warranted here. If we did store family ID somewhere, we could omit zeros in ehr.kinship and accurately 're-hydrate' the kinship table. My current thinking is that it would be OK in practice to write non-zero rows and let callers make the assumption that no row = zero, without needing to check how families were calculated in R.

[ankurjuneja]

Thanks for running this on your data.

Agreed on the big one. "Basically one massive family per species" matches what I measured: the largest family is 76%, 98% and 100% of the largest species at the three centers I tested. That's exactly why this PR only gets 1.2x at the largest site and the peak is set by one dense matrix, and splitting families doesn't shrink it. The split is kept because it lets each matrix be written and released before the next is allocated, not because it lowers the largest allocation.

One correction on the dropped rows. I think the 20-30K records are self-pairs, not zeros.
I ran both scripts and counted:

| | total rows | coefficient == 0 |

| current script | 1,681,615 | 0 |
| this PR | 1,678,444 | 0 |
Neither emits a zero. The current script drops them too - as(temp.kin, "dgCMatrix") discards zeros during the coercion. The 3,171-row delta is entirely the matrix diagonal (an animal against itself), and those coefficients are 0.5 or higher.

Those rows were already dead weight: GeneticCalculationsImportTask has if (fields[0].equalsIgnoreCase(fields[1])) continue; //dont import self-kinship. They were being written and then thrown away on import, so nothing reaches ehr.kinship differently. Old vs new output is byte-identical after normalising for that, at all three colonies.

Could you check coefficient == 0 on your current output? If you're seeing real zeros, your
pedigree hits something mine don't and I'd like to look at it.

MIN_COEFFICIENT is wired up — populateKinship.r:94-95. It's set to 0, which disables it.

On thresholding, I think that's the right conversation and it's the real remaining win, but
there's a catch worth settling first: the colony-wide kinship average queries divide by an
independent population count rather than by the rows present, so dropping distant pairs silently
biases those averages downward. Worth confirming which centers depend on that before we enable it. @labkey-martyp

On pivot_longer/filter: that materialises the whole-colony triplet frame, which is the
specific allocation causing the nightly spike and it's what the current script does via bind_rows.
Also worth noting dplyr was never listed in install.r, so it was an undeclared dependency this
PR removes.

The SubjectId -> FamilyId TSV is a good idea and cheap, happy to add it. For what it's worth, if makefamid were splitting incorrectly I'd expect missing pairs, and output is identical to the current script at all three centers, so it doesn't appear to be happening in practice.

[bbimber]

@ankurjuneja: is that that AI-written or guided?

To be clear: I am not saying anything is omitting coefficient = 0. I am saying that they different in the number of unrelated (coefficient==0) rows they output. If your script changes the composition of the matrices passed to the kinship function, it will change the set of ID pairs in output. This isnt an NA coercion question.

I am also saying that the original script had an implicit covenant to report all within-species pairs that have been tested, which allows the downstream consumer to differentiate 'not compared' with 'not related'. That is an important difference in this PR from prior behavior. If this PR is partially walking back that covenant (it is), I am saying that we might consider going all the way. If we omit all unrelated rows (which implies that downstream code can assume no-data = unrelated), that massively reduces a lot of steps in this process.

[ankurjuneja]

yes, the response is claude assisted. I'm not an R expert, so I'm using Claude for the analysis and verification. In my understanding on the "Zeros" neither the changed script or existing script writes a row for an unrelated pair (the old one throws them and mine never generates them).

I think in the existing script, every animal appeared at least once paired with itself but they never reached the database because the Java importer throws self-pairs away on import - ehr/src/org/labkey/ehr/pipeline/GeneticCalculationsImportTask.java - L273-274

if (fields[0].equalsIgnoreCase(fields[1])) continue; //dont import self-kinship
So self-pairs have been written by R and discarded by Java on every run, long before this PR.

The change in my PR for the animals with no relatives at all, whose only row was a self-pair, no longer appear in the file.

[bbimber]

@ankurjuneja:

I think in the existing script, every animal appeared at least once paired with itself but they never reached the database because the Java importer throws self-pairs away on import - ehr/src/org/labkey/ehr/pipeline/GeneticCalculationsImportTask.java - L273-274

The row count of the output changed with this PR. there's explainable reasons for this, but it has implications. This has nothing to do with IDs paired to themselves.

The reason is changed is that different combinations of Id1/Id2 are written. the reason for this is clearly because rather than comparing the entire species at once (therefore creating a matrix where each dimension includes every animal that ever existed for that the species), it is now partitioning this into families and only calling kinship() per family. That means that total number of pairs in the output is quite different. If makefam() is accurate, all of those missing records are unrelated (coefficient=0).

Historically, the output of this script ensured that every pair within a species was calculated. This means that within a species, every pair of Ids should be represented in the output. Calculating kinship() in R across the full species ensures this. A sizable number of animals do have no kinship (i.e. coefficient = 0). Nonetheless, omitting a record for a given pair and reporting a pair where coefficient is zero are not the same thing. It's a reasonable question to ask whether this matters. This PR quasi-arbitrarily stops reporting some of those pairs, because it splits up IDs by makefam(). That is arguably a reasonable thing to do; however, if you're going to do this, you might as well go all the way and get a substantial benefit on import, rather than a marginal one. If consumers can no longer guarantee that lack-of-data for a given pair means that kinship was calculated and is zero, then I dont think there is a lot of value for storing millions of coefficient=zero records.

if (fields[0].equalsIgnoreCase(fields[1])) continue; //dont import self-kinship So self-pairs have been written by R and discarded by Java on every run, long before this PR.

This is not relevant.

[ankurjuneja]

can you confirm whether you ran both scripts against the same pedigree.txt, or compared against an existing production kinship.txt?

If you did run both on the same input, could you post three numbers from your data?

• row count from the existing script
• row count from this PR
• count of rows in the existing output where Id and Id2 are the same animal

[bbimber]

Hi @ankurjuneja: My apologies. I made a mistake when I first looked at these outputs. My comments above are wrong. You are correct that the only difference in output is same-ID rows, which would get dropped by the java import anyway.

I havent yet wrapped my head around how splitting families doesnt result in failure to report some pairs (i.e., two animals in the same species but not listed in the same family), but in practice maybe all living animals are lumped into one family. Anyway, this is seems good. Again, sorry for the mistake above.

These might be good changes to explore if you go any further, but the script is largely working for our purposes and I dont I feel especially strongly about this: 1) switching from the deprecated kinhsip2 package to Pedixplorer, 2) using sparse matrices, 3) only write and store non-zero coefficients (big win for import time), 4) write a TSV with the ID->family assignments for debugging/QC.

[labkey-martyp]

@ankurjuneja @bbimber I've done another refactor here to take the parts of kinship2 needed for our computation and put them directly in our script with some memory optimizations and removed kinship2 dependency. That removes the unpinned dependency and lets us make optimizations that match our implementation. I verified identical results on four different pedigrees.

@bbimber try it out on your data to see if you see any issues. Should produce identical results.

Specific script changes on the copied kinship2 implementation (from claude):

• Same algorithm, same complexity. Identical recursion, still O(n²) — output is bit-identical (identical(), max diff exactly 0) on families up to 24,186 animals.
• No trailing matrix copy. kinship2 allocates (n+1)² then copies to n², holding 8.9 GB live to produce a 4.5 GB result; the in-script version allocates n × n once and never copies.
• Bounded generation temporaries. kinship2 updates a whole generation in one expression — about 1 GB per temporary at center B's 5,506-animal generation — where chunking to 512 rows caps them near 100 MB regardless of generation size.
• No placeholder rows. pedigree()'s "both parents or neither" rule forces a phantom founder for every one-parent animal; dropping them shrinks n itself, so the saving is quadratic — at the largest colony, 53,223 → 30,359 animals and 21.1 GB → 6.9 GB.
• Measured decomposition. Isolating kinship() on B's largest family, peak falls from 3.2× the matrix size to 1.56× from the first two items, then a further 16% once the placeholders go.

[bbimber]

This matrix is going to remain sparse, right? Would:

kmat <- Matrix::.sparseDiagonal(x = 0.5, n = n)

be more efficient?

[labkey-martyp]

So claude did analyze this and determined that a sparse matrix would give a very small benefit to the amount of memory used for the matrix but due to the way it's accessing the matrix it would cause a very significant slow down. So maybe slightly more efficient in terms of memory but the trade off in performance would be large.

[bbimber]

I know this isnt used, but I think it would make more sense to set this to -1, and change the logic below to " > MIN_COEFFICIENT". I suspect a logical use-case would be to report only non-zero pairs, so the existing GTE logic becomes cumbersome.

[labkey-martyp]

So just to clarify, the script already outputs only non-zero pairs regardless of this coefficient (filtered around line 216). I think this is more for filtering out very small coefficients if further optimization is needed.

[bbimber]

See comment above about allowing MIN_COEFFICIENT=0, which if the simplest way to let the user report all non-zero coefficients. Since most are zero, I dont see a lot of value in picking a really, really small number instead.

I suggest line 241 should use GTE. Line 242 should use GT logic.

[labkey-martyp]

As mentioned above the zero pairs should already be filtered out. I'm not sure if there's really a good use case for filtering out small numbers now that this PR is a pretty significant reduction in memory usage. I could be convinced to remove this filter entirely.

[bbimber]

What's the reasoning for not just failing? I dont quite understand the comment above this method, but even if kinship2's existing code ignored sex discrepancies during the kinship calculation, they clearly indicate something is incorrect.

[labkey-martyp]

This basically does the same behavior as before, the kinship is calculated by the role (dam/sire) not the sex. The sex is used only for validation. If there's a female sire or male dam, the sex is ignored and the animal is added to the pedigree based on their role, dam or sire. This actually adds a warning that was not there before to inform users of the inconsistency. I chose to make it a warning because there are currently about 60 animals across the five centers I tested that have the wrong sex for their role. This will make the users aware without breaking their nightly calculations.

[bbimber]

The resulting text file is quite large. Would it make sense to just write a gzip output?

con <- gzfile(TEMP_FILE, open = 'w')

[labkey-martyp]

Yeah that's a good idea. That would require some changes on the java side. I can investigate for a future PR.