[gbrsales]

Summary

Speed up the GSVA scoring stage by evaluating the KS random walk only at the ranks of the genes in the set under test.

Changes Made

• Rewrote the dense gsva_rnd_walk to accumulate over the k ranked genes belonging to the set under exam. The sparse and NA-handling variants are unchanged.
• Moved rank fetching and random-walk scoring into separate compilation units (ranks.c, rnd_walk.c).
• Rank fetching uses a per-call context with a pre-allocated, reused rank buffer instead of allocating one per column.
• Golden test output (test_gsva_newapi.R) is now rounded to 5 significant digits: the rewritten walk sums in a different order, so results differ from previous versions in trailing digits only while remaining numerically equivalent.
• Added a benchmark suite (benchmarks), excluded from the package build via .Rbuildignore (could be dropped from the PR if not wanted outside development).

Speedup results

Benchmark of the scoring stage (gsvaColScores) on PBMC68K profiles with different MSigDB collections. Medians of 2–5 bench::mark iterations. baseline = devel, current = this branch.

Tier	Gene sets	Cells	Baseline (ms)	Current (ms)	Speedup
T1 Hallmark	50	66,722	135,531.8	25,052.1	5.41×
T2 C7 (top 1000)	984	66,722	2,613,065.2	564,877.0	4.63×
T3 C2 (top 1000)	917	66,722	2,415,864.5	383,911.2	6.29×
T4 C5:BP (top 2000)	1,472	66,722	3,857,668.7	599,567.2	6.43×
T5 Hallmark+C5 (top 2000)	1,472	66,722	3,871,121.2	605,780.6	6.39×

[rcastelo]

Thanks @gbrsales !! I'll review it as soon as possible and get back to you.

[rcastelo]

Hi @gbrsales thanks again, please wait to push changes until I'm done with the review, otherwise makes the review harder for me.

[rcastelo]

Hi @gbrsales thanks again for this PR with this first optimization of the GSVA random walk for dense expression data. While reviewing this PR, I found a few bugs in the original R and C code, which I have fixed in devel (and in RELEASE_3_23), and which they are important to fix before merging this PR:

• Undefined random walks. A random walk score is undefined when the gene set contains all genes, or when all its genes have a zero rank statistic (e.g., a single gene at the middle rank), because they result in a 0 value in the denominator of the running CDF. These now return NA explicitly, are counted, and trigger a warning.
• Duplicated genes in gene sets. These are now removed, with a warning, when the parameter object is built, and also when gene sets are mapped to features.
• Memory leak. An early return in gsva_rnd_walk_nas() skipped freeing two R_Calloc() buffers.
• Integer overflow. nr*j in the dense column function fetchers, overflows rank matrices with more than 2^31 elements.
• Walks with missing values. These ran over all p positions instead of the p - nnas non-missing ones, so missing genes acted as phantom out-of-set genes at the bottom of the ranking. This biased the use="na.rm" and use="everything" scores.

Several of these touch code that this PR moves into ranks.c and rnd_walk.c, which is how I discovered them, so the PR needs to be brought up to date with devel. Details below. Please follow the following steps.

1. Bring the branch up to date with devel

I tried this locally. Rebasing your 7 commits one by one runs into conflicts starting at the second commit, Move functions for ranks to a separate compilation unit. They come back in later commits, because the fixes touch code the PR first moves and then rewrites. The simplest way is to squash the PR into a single commit first, so you resolve the conflicts only once:

git remote add upstream https://github.com/rcastelo/GSVA.git   # if not there yet
git fetch upstream
git checkout <your-PR-branch>
git branch pr285-backup                      # safety net
git reset --soft $(git merge-base HEAD upstream/devel)
git commit -m "Faster dense GSVA random walk"  # one commit with all your changes
git rebase upstream/devel                    # conflicts only in src/ks_test.c

If you'd rather not rewrite history, git merge upstream/devel gives exactly the same single conflict, and you can push without forcing.

Resolving src/ks_test.c

There are three conflicting code chunks:

1. extern declarations at the top. Keep your side and add the new symbol:

   extern SEXP GSVA_attrNAsSym,
               GSVA_attrUndefWalksSym;

2. Local variables in gsva_score_genesets_R(). Keep your side and add:

   double   nundefwalks = 0; /* double to avoid overflowing m x n */

3. The three random-walk calls. Keep your signatures, but take the return value and pass p - nnas to the NA walk:

   undefwalk = gsva_rnd_walk_nas(gsetidx, k, decordstat_col,
                                 symrnkstat_col, p - nnas, tau, nause,
                                 minsize, &walkstatpos, &walkstatneg,
                                 &wna);
   ...
   undefwalk = gsva_sparse_rnd_walk(gsetidx, k, decordstat_col,
                                    symrnkstat_col, p, tau, nzeros_col,
                                    zerosymrnkstat, &walkstatpos,
                                    &walkstatneg);
   ...
   undefwalk = gsva_rnd_walk(walkctx, gsetidx, k, decordstat_col,
                             symrnkstat_col, p, tau, &walkstatpos,
                             &walkstatneg);

   The rest of gsva_score_genesets_R() (nnas = 0;, undefwalk, nundefwalks, the attrUndefWalks attribute, and the all.obs check) merges automatically.

Fixes that git will not carry over for you

The devel side of the first code chunk contains the fixed versions of the functions this PR moves out of ks_test.c. When you keep your side, those fixes are silently dropped. Please port them by hand:

• src/rnd_walk.c and src/rnd_walk.h. The three walk functions now return int: 1 if the walk is undefined, 0 otherwise. Change the prototypes to match.
  • gsva_rnd_walk(): when k == 0 || k == n, or when total_pos_sum == 0.0, set *walkstatpos = *walkstatneg = NA_REAL and return 1, instead of returning 0.0 scores. Return 0 at the end. A score of 0 claims "no enrichment" for a statistic that is actually undefined. It would also disagree with the sparse and missing-value paths, which return NA for the same gene sets.
  • gsva_rnd_walk_nas(): free gsetrnk and gsetidx_wonas before the early return taken when na_use < 3, and return 0 there. Set undefined = !(stepcdfingeneset[n-1] > 0 && stepcdfoutgeneset[n-1] > 0) and return it at the end.
  • gsva_sparse_rnd_walk(): compute and return undefined in the same way.
  • The comments above these functions in devel's ks_test.c describe this contract. Please move them along with the code.
• src/ranks.c. In fetch_intcol_intmatrix() and fetch_intcol_dblmatrix(), compute the column offset as R_xlen_t offset = (R_xlen_t)ctx->p * j; and use X + offset / X[offset + i].

Your new gsva_rnd_walk() makes the reasonable assumption that gene indices in a set are unique, which means that duplicates would be double-counted. Since devel now removes duplicates before they reach C, it would useful for us and any other future developer to have a one-line comment stating that assumption, and that is correct to make it because duplicates are removed upstream.

2. Revert the test changes

After applying the fixes described in step 1, devel's original regression test output tests/test_gsva_newapi.Rout.save should match unchanged, and so should all the other test outputs. I checked this with my own resolution, and the scores agreed with devel to within 3.3e-15, including the NA and edge cases. So, in this PR, there is no need to change any tests/*.Rout.save. Because differences are around 1e-16, so neither test change is needed:

• inst/unitTests/test_gsvaColRanks.R: please drop tolerance = 1e-6. The C score and the R gsvaEnrichment() score differ by about 1e-16, well within the default tolerance of checkEqualsNumeric().
• tests/test_gsva_newapi.R and tests/test_gsva_newapi.Rout.save: please revert the rounding to 5 digits. The unrounded output reproduces exactly with your implementation, and rounding only makes the regression test less sensitive.

Please also update the PR description: the results are not "rounded to 5 significant digits", they are exact up to floating-point summation order.

3. Smaller C issues

• ranks.c, ranks_ctx_create(): eval(lang2(install("class"), XR), R_BaseEnv) evaluates an unprotected call, and its result is also unprotected (rchk would flag this). It also duplicates get_class_name(), which still lives in ks_test.c and is used by average.c. Please call get_class_name(XR) instead, which is also faster.
• ks_test.c: <assert.h> and <R_ext/Rdynload.h> are included but not used.
• ranks.h: please add include guards, as in rnd_walk.h.
• ranks_ctx_t: the type field is set but never read. Either drop it or use it.

4. Benchmarks and DESCRIPTION

Keeping benchmarks/ in the repository is fine, since it's in .Rbuildignore. For the same reason, nothing in the built package uses bench or msigdbr, so please remove them from Suggests. Otherwise every user and builder pulls them in for no benefit, and R CMD check / BiocCheck may complain about unused suggested packages. A comment at the top of benchmarks/setup.R listing the packages it needs would be enough.

Thanks again. Once the branch is rebased and these points are addressed, I'll be happy to merge.

[rcastelo]

A quick heads-up: devel has moved on by one commit since my review (926e1e5, which replaces bpiterate() with bplapply() for the DoparParam back-end). It doesn't touch any of the code in this PR, but it adds foreach and doParallel to the end of Suggests in DESCRIPTION, the same line where this PR adds bench and msigdbr.

So when you rebase onto upstream/devel (or merge it), you'll get one extra conflict in DESCRIPTION, besides the one in src/ks_test.c described above in my review:

<<<<<<< upstream/devel
          batchtools, foreach, doParallel
=======
          batchtools, bench, msigdbr
>>>>>>> (your commit)

If you follow the suggestion to drop bench and msigdbr from Suggests, just keep the devel line, batchtools, foreach, doParallel. The ks_test.c conflict is exactly the same as before, so the instructions in my review still apply unchanged.