Adding difference and division

[eotamm]

This pull request adds now two new transformations to the transformAssay() function: difference (computes all pairwise differences (x - y) between features (rows) across samples) and division (computes all pairwise ratios (x / y) between features)

[antagomir]

See the suggestions.

Unit tests should be added as well.

Examples section could also show how to do logx-logy by using the combination of log and difference transformations.

[eotamm]

Added rowname check, zero value warning, and switched to vapply in difference and division functions.

[eotamm]

I updated transformCounts unit tests and documentation for division and difference.

[TuomasBorman]

I started to wonder if these current methods are meaningful. As far as I know, usually, some reference is used instead of calculating all pairwise differences/ratios (reference for instance mean of all features).

This also causes technical problem as number of pairs increases quickly and the resulting matrix is huge.

[antagomir]

There are uses for full log-ratios, for instance this one: https://academic.oup.com/nargab/article/6/2/lqae038/7658053 and we would like to add comparisons to it. And one might have similar needs in settings where number of features is more limited.

[TuomasBorman]

There are uses for full log-ratios, for instance this one: https://academic.oup.com/nargab/article/6/2/lqae038/7658053 and we would like to add comparisons to it. And one might have similar needs in settings where number of features is more limited.

OK. I think we should set limit for number of features, and give error if it exceeds it. Otherwise, users just keep on crashing their sessions.

[eotamm]

What do you think would be a reasonable maximum number of features?

[antagomir]

This depends on the computing capacity so it is tricky to set a definite hard limit.

How about throwing an informative warning instead?

Also: make sure not to double calculate the ratios (x/y vs. y/x),

[eotamm]

I added warning if the matrix has too many rows, supported "-" and "/" aliases for method selection, and switched to using feature indices.

[TuomasBorman]

Can you still explore with SparseAssays?

Problem also is usage of combn() as it creates a huge data.frame which is loaded into memory. Can it create SparseAssay?

If not, you could first initialize SparseArray and then create a nested loop that goes through the feature pairs. It is not the nicest solution, but should preserve memory which is the problem here. Could you check if that works?

Another option is C++...

[eotamm]

Thanks for the suggestions! I will try to check if using a SparseAssay is possible, either directly or with a nested loop.

[eotamm]

I replaced combn() with a nested loop implementation. As far as I understand, SparseAssay doesn’t support direct use with combn(). Updated unit tests too.

[TuomasBorman]

Can you make comparisons: this new approach vs the old one. Test these methods with different number of samples, let's say 100, 200, 500, 1000, 2000, 5000, 10000. Record the time it takes to run foe each. (looping is also suboptimal solution, but might be safer and more practical than combn)

[eotamm]

I compared the combn()-based and nested loop implementations using 100 features and varying the number of samples. Below are the runtimes (in seconds):

100 samples: combn = 0.03 s, loop = 10.63 s
200 samples: combn = 0.05 s, loop = 21.61 s
500 samples: combn = 0.06 s, loop = 66.34 s
1000 samples: combn = 0.14 s, loop = 170.16 s
2000 samples: combn = 0.23 s, loop = 827.42 s
5000 samples: combn = 1.58 s, loop = 1787.61 s
10000 samples: combn = 1.70 s, loop = 3397.47 s

Would it make sense to add a sparse = TRUE argument or something like that (defaulting to FALSE, i.e., the faster combn-based version), so that users can explicitly choose the sparse loop-based implementation when needed? Or would it be better to always use just one of the methods?

[TuomasBorman]

I compared the combn()-based and nested loop implementations using 100 features and varying the number of samples. Below are the runtimes (in seconds):

100 samples: combn = 0.03 s, loop = 10.63 s 200 samples: combn = 0.05 s, loop = 21.61 s 500 samples: combn = 0.06 s, loop = 66.34 s 1000 samples: combn = 0.14 s, loop = 170.16 s 2000 samples: combn = 0.23 s, loop = 827.42 s 5000 samples: combn = 1.58 s, loop = 1787.61 s 10000 samples: combn = 1.70 s, loop = 3397.47 s

Would it make sense to add a sparse = TRUE argument or something like that (defaulting to FALSE, i.e., the faster combn-based version), so that users can explicitly choose the sparse loop-based implementation when needed? Or would it be better to always use just one of the methods?

The loop seems to scale very badly, however, N features (instead of samples) might be more interesting as that is the bottleneck in the current implementation (too big data.frame).

Anyways, loop seems not to be feasible. Also, the current combn() is problematic. Can you still check options? For instance, DelayedArray, C++, or datatable (https://stackoverflow.com/questions/26828301/faster-version-of-combn)

[eotamm]

The data.table-based approach seems slightly faster than combn(), but since it returns a dense matrix, I'm not entirely sure about its memory efficiency in practice. The C++ implementation, which directly constructs a sparse matrix, appears to be similarly fast and is likely the most memory-efficient option of those tested. I translated the sparse matrix logic from R to C++, but I don’t have previous experience with C++, so I’m not fully confident in the implementation details.

[antagomir]

If the tests produce correct results then it should be fine?

[TuomasBorman]

The data.table-based approach seems slightly faster than combn(), but since it returns a dense matrix, I'm not entirely sure about its memory efficiency in practice. The C++ implementation, which directly constructs a sparse matrix, appears to be similarly fast and is likely the most memory-efficient option of those tested. I translated the sparse matrix logic from R to C++, but I don’t have previous experience with C++, so I’m not fully confident in the implementation details.

Could you check faith index implementation (along with its PR)? @jepasan implemented that some time ago, if you propose initial draft I could help you from that

#697

[eotamm]

I just pushed my own changes there as well at the same time. I'll check that everything follows the same style and let you know once it's all good. At least all the tests passed successfully.

[eotamm]

I added the corresponding entries to RcppExports, following the same approach as in the Faith index implementation. I'm not entirely sure whether these generated entries in RcppExports.R and RcppExports.cpp should be committed to Git, since they are automatically created (at least locally), but I included them just in case. It should now work.

[TuomasBorman]

Awesome!!!

For the sake of maintenance, can you move functionality that can be done in R to R functions?

C++ functions should only contain the calculation and they should output sparse matrix. R can do the rest.

[TuomasBorman]

Also add your name to DESCRIPTION file

[eotamm]

I moved the warnings and rowname generation logic to R, only the calculations are now handled in C++.

[TuomasBorman]

Great, thanks!!