Clean up package structure and optimize script

[jeff-mandell]

Hi,
Thanks for making this package. I develop integrative methods to quantify cancer evolution (see https://townsend-lab-yale.github.io/cancereffectsizeR/index.html), and your armSCNA peak-finding algorithm and fitness measure will likely be useful for my work.

In going through your code, I saw that the GEV calculations could be optimized to make the algorithm run a lot faster. After making the change, I'm able to run the pan-cancer analysis in a few minutes on a single core, compared to several hours on multiple cores. The codebase could also be improved by restructuring the file organization to match a standard R package (which would enable simpler installation), and it would also be nice if a lot of hard-coded options and file directories were changed so that default outputs don't go into the package directory itself.

If I make a version with these changes, would you be willing to review a pull request?

Thanks,
Jeff

[shahab-sarmashghi]

Hi Jeff,

I'd love to see such improvements in running time and codebase/UI, and I'm sure the whole community will benefit from that. That being said, let me first discuss this internally before you put more efforts into this. This might take a week or two but I will surely get back to you.

Best,
Shahab

[shahab-sarmashghi]

@jeff-mandell we welcome your contribution. I'll be more than happy to review your PR and accept it after some testing to make sure the algorithm works as expected.

[jeff-mandell]

Great, I will try to have something over by next week.

[jeff-mandell]

I created a pull request; see #6.

[jeff-mandell]

Hi @shahab-sarmashghi, is there anything I could do that would be helpful in facilitating review of the pull request? There's no rush on my end, although eventually it would be nice to point people to the original package rather than my fork

[shahab-sarmashghi]

Hi @jeff-mandell, sorry I've just been very busy. Will try to get through these in the next week or two.

[jeff-mandell]

Hi @shahab-sarmashghi, just checking if there is any chance of the pull request getting reviewed at some point? If there is any way I could help streamline that process, let me know.

[shahab-sarmashghi]

Hi @jeff-mandell I'm sorry that I forgot about this. Did you try to run the original version and yours on the same exact input that we use in the paper (pan-cancer mode)? How did the peaks compare? Are they identical?

[jeff-mandell]

I compared the results reported in your paper (Table S3A) with the output produced by my BISCUT branch using PANCAN_ISAR.seg.txt as input, acquired from your tutorial. I'm not sure if this is exactly the input used in your paper; if not, I can run again if you are able to provide the input.

Due to the stochasticity of finding peaks and their endpoints, results do not exactly match. One way to compare between the two sets of results is by comparing peak_ids (which contain chromsome arm, amp/del, tel/cent, the chain of selection directions, and the iteration number; e.g., 9q_amp_tel_n-rn_2). By peak ID, 83% of peaks found in Table S3A were re-identified in the new run. When considering just first iteration peaks, 97.9% of the original peaks were re-identified.

Below is the code used for comparison. Let me know if you want to run this, and I'll make a new branch that has the files in place for this to be quickly reproduced.

Code

# Run from BISCUT dev directory
library(devtools)
load_all() # load dev version of package
library(data.table)

# Get pan-TCGA peaks from Shih et al., table S3A
original = as.data.table(readxl::read_excel('for_compare/41586_2023_6266_MOESM5_ESM.xlsx', 
                                            sheet = 1, skip = 17))
original = original[, peak_id := paste(arm, direction, telcent, code, iter, sep = "_")]


# Compare to peaks running on altered version of BISCUT-py3
# Downloaded from the link given in the BISCUT README: 
# https://drive.google.com/drive/folders/1-VZ_A0uodEs04Jg-Gphkl3HU5VdEyICW?usp=share_link
segment_file = 'for_compare/PANCAN_ISAR.seg.txt.gz'

 # You may need to adjust cores depending on your system
make_breakpoint_files(segment_file = segment_file, output_dir = 'for_compare/breakpoint_for_compare', cores = 8)
results = do_biscut(breakpoint_file_dir = 'for_compare/breakpoint_for_compare/', 
                    results_dir = 'for_compare/biscut_output', cores = 8, seed = 999)


original_peaks = unique(original[, .(peak_id, code, iter, direction, telcent, negpos, combined_sig)])
new_peaks = results$peaks[, .(peak_id, code, iter, direction, telcent, negpos, combined_sig)]



mean(original_peaks$peak_id %in% new_peaks$peak_id)
# [1] 0.8341969

all_peak_id = unique(c(original_peaks$peak_id, new_peaks$peak_id))
peak_compare = data.table(peak_id = all_peak_id,
                          in_original = all_peak_id %in% original_peaks$peak_id,
                          in_new = all_peak_id %in% new_peaks$peak_id)


# 83.4% of the peaks in the original results match peaks found in the new run.
peak_compare[in_original == T, mean(in_new)]
# [1] 0.8341969

table(peak_compare[, .(in_original, in_new)])
#             in_new
# in_original FALSE TRUE
# FALSE       0   70
# TRUE        32  161

# Considering only 1st iteration peaks, 97.9% of original peaks are reproduced in the new results.
peak_compare[, iter := sub('.*_', '', peak_id)]
first_iter_compare = peak_compare[iter == 1]
first_iter_compare[in_original == T, mean(in_new)]
# [1] 0.9789474

table(first_iter_compare[, .(in_original, in_new)])
#             in_new
# in_original FALSE TRUE
# FALSE       0   20
# TRUE        2   93

[jeff-mandell]

Hi again,
I made a new branch on my fork show the comparison between original and refactored outputs. It's the same codebase as my pull request except there is a for_compare folder with a script that demonstrates the comparisons given in my last comment here. If you install my forked version as described on the pull request, you can run through the code in compare_results.R if you like.

I added a check for how many gene-level ampdel/telcent peaks agree in direction of selection. Of the 4,921 gene/ampdel/telcent records that are shared between the original peak results (reported in Shih et al.) and the version that I generated, 97% agree on negative/positive selection according to the negpos field.

I'm continuing to work on my own projects related to SCNA selection, which use orthogonal methods and are not focused on anchored events. If I include comparisons with gene-level BISCUT results in write-ups, it would be nice to be able to say that I used BISCUT, rather than that I used a modified version of BISCUT. I'm also interested in making it easy for users of my software to also run BISCUT, and I would prefer to advise them to install your package rather than a fork of it.

[jeff-mandell]

Hi @shahab-sarmashghi, are you still the maintainer?

[jeff-mandell]

Hi @rberouk, is there anyone who might be able to review this issue/pull request? Sorry that it's a bit of a handful, but it would be nice if BISCUT could incorporate the improvements that I made to facilitate its future use. See #6 for details.

[rberouk]

Hi Jeff, thanks so much for working on this. Shahab is indeed still the maintainer and we just discussed including your improvements last week. Hopefully this can be done soon! And we may be working towards similar goals on improved methods to your goals. So it might be useful to keep in contact as we move towards papers in case we can support each other.

[shahab-sarmashghi]

Hi @jeff-mandell, thank you so much for all the great work on the codebase and performance improvements. I was able to successfully test the new version. We're happy to merge the pull request but there are a few points I wanted to address while doing that. Here are the main points, I'll leave the more detailed ones under the pull request:

• The results are mostly overlapping but it seems in some cases there are differences. The measure of overlap between two sets of peaks (say A and B) that we prefer the most is the number of more significant peaks in set A that are overlapped by any peaks in set B in the same direction (amp-pos, amp-neg, del-pos, or del-neg). Out of 170 paper peaks with log10ksby>2, 160 are overlapped by peaks from a different seeding of the current version (by at least 50kb), whereas that number is 150, 152, and 154 for three different seeds of the new version (999, 9, 99). This is not a dealbreaker but it's worthwhile to understand the reason and mention it in the README and the new version's release note. I tried to look into some of the differences to understand this better. One example is an amp-tel-pos peak on 5p15.33 (chr5:1020354-1287540) which is missing in the new version and instead there is an amp-tel-neg peak detected at that locus (chr5:1006073-1373493). When you look at the sequence of peaks detected on that arm, more differences show up, starting with a very different peak boundaries around the peak detected in the first iteration. I have attached the plots here for comparison. Could that be the reason? My impression of the changes to the peak boundary finding was that it's only refactored for speedup, but could it be that the behavior has changed beyond random seeding?

Current version's first iteration 5p amp tel pos peak

New version's first iteration 5p amp tel pos peak

• The current version makes some outputs and plots that have changed in the new version, which are mostly fine. However, there were a groups of outputs that were particularly useful where each peak was a column with varying number of rows depending on the number genes it included. I think they are generated by make_column_results () and extract_cols() functions in combine_BISCUT_results.py. Is it possible to add those?

[jeff-mandell]

Hi @shahab-sarmashghi,
Sorry for the slow response. I'm in a transition period and may not be very available for the next month or two, but I'll do what I can.

I don't know why the two versions produce different peaks at 5p. You're correct that the boundary-finding methods are intended to be conceptually the same. Actually, for my results comparisons, I've been using Table S3A from Shih et al. (named 41586_2023_6266_MOESM5_ESM.xlsx), and this table does not report the peak that you're describing (chr5:1,020,354-1,287,540). It instead reports an amp-tel-pos-iter1 peak at 28,581,733-31,773,999 might be able to investigate more if you can provide the exact inputs that can be used to produce the peak in your plot.

As for wide output, here's some code to transform the existing output in something like the output of make_column_results(). I can make this a function in the package (what should it be called?) or have it automatically included in do_biscut() output (again, what would you want to call it?), if this is in fact the format you're looking for.

# Starting with BISCUT output
# results = do_biscut(...)
library(data.table)
peaks = copy(results$peaks)
peaks[results$genes_by_peak, cytoband := Cytoband, on = 'peak_id']
peaks = peaks[, .(peak_id, cytoband, peak_location = peak_interval,
                  combined_sig, log10_ksby, ks_stat,
                  n_events, direction, telcent, selection = negpos,
                  code, `TS or onco-like` = type_of_selection)]

column_names = c(names(peaks), 'genes')
peaks_to_gene = results$genes_by_peak[, .(gene = list(Gene)), by = 'peak_id']
tmp2 = apply(peaks, 1, as.list)
tmp3 = lapply(tmp2, function(x) unlist(unname(c(x, as.list(peaks_to_gene[peak_id == x$peak_id, gene])))))
req_length = max(sapply(tmp3, length))
tmp4 = lapply(tmp3, function(x) c(x, rep(NA, req_length - length(x))))
tmp5 = as.data.table(tmp4)
tmp6 = setnames(tmp5[2:.N, -"V1"], new = unlist(tmp5[1, -"V1"]))
final_column_names = c(setdiff(column_names, 'peak_id'), rep(NA, req_length - length(column_names)))

output = cbind(final_column_names, tmp6)
setnames(output, 'final_column_names', 'peak_id')

# For example
output[1:20, 1:3]