Drycleaned fragCounter input for JaBbA causes drop of SVs during edge selection

[pblaney]

Hello @mskilab,

I see on the JaBbA repo that a 'best practices' are in the works so this might have come up in already.
However, I figured it'd best to detail my thoughts and get an opinion on if this is reasonable.

The behavior I am experiencing occurs when using newly drycleaned fragCounter coverage with JaBbA, the algorithm excludes numerous (but not all) junctions due to "both breakpoints are in NA coverage regions". From what I see in the code, this has a direct impact on the edges selected in from the graph and should to be included (those lost include both inter- and intra- chromosomal breakpoints).

My quick assessment was the input SVs were properly read as input and then dropped as the drycleaned fragCounter coverage contained gaps in the coverage where this would lead to a NA value for the foreground. Considering the original fragCounter input (pre-drycleaning) is binned across the entire genome, I suspected there should be no gaps especially after drycleaning (I believed this is one of the reasons for this step).

To see what occurred between the drycleaning and original fragCounter coverage, I ran start_wash_cycle changing only the germline filter to FALSE:

# Original Run
start_wash_cycle(cov = prewashed_tumor_cov,
		               mc.cores = 1,
		               detergent.pon.path = path_to_detergent_pon,
		               verbose = TRUE,
		               whole_genome = TRUE,
		               germline.filter = TRUE)

# Second Run, without germline.filter
start_wash_cycle(cov = prewashed_tumor_cov,
		              mc.cores = 1,
		              detergent.pon.path = path_to_detergent_pon,
		              verbose = TRUE,
		              whole_genome = TRUE,
		              germline.filter = FALSE)

Interestingly, this revealed a drycleaned fragCounter coverage profile with the same number of bins as the input (i.e. there weren't any gaps in the GRanges). Given this change, I searched the code base and found lines 784-806 in dryclean.R:

setnames(cov, "signal", "input.read.counts")
cov = cbind(decomposed[[2]], cov)
colnames(cov)[1] = 'foreground.log'
cov[is.na(input.read.counts), foreground.log := NA]
cov[, foreground := exp(foreground.log)]
cov[input.read.counts == 0, foreground := 0]
cov[is.na(input.read.counts), foreground := NA]
cov = cbind(decomposed[[1]], cov)
colnames(cov)[1] = 'background.log'
cov[is.na(input.read.counts), background.log := NA]
cov[, background := exp(background.log) ]
cov[input.read.counts == 0, background := 0]
cov[is.na(input.read.counts), background := NA]
cov[, log.reads := log(input.read.counts)]
cov[is.infinite(log.reads), log.reads := NA]

if (germline.filter){
    germ.file = rpca.1$inf_germ
    cov$germline.status = germ.file$germline.status
    cov[germline.status == TRUE, foreground := NA]
    cov[germline.status == TRUE, foreground.log := NA]
    cov = na.omit(cov)
}

The problem is all the NA records in the cov object. Under either condition for germline.filter, there coverage regions that are evaluated to NA that are used for input in JaBbA which causes issues downstream.

My thought process in solving this two-sided issue was:

• If we are confident to set foreground == 0 when input.read.counts == 0, then a reasonable assumption would be that in all scenarios where is.na(input.read.counts) == TRUE we can should set foreground to some small non-zero value. I chose exp(-5) as this was the lowest I saw for the range of foreground.log in some of my fragCounter cov files.

• The main issue is the NA for segments where germline.status == TRUE, while this will achieve the desired effect of not having these FPs appear as enriched SCNAs but something like GISTIC. However, the NAs can cover regions where a somatic SVs breakpoint resides, leading to its exclusion. I reasoned an alternative would be that in all scenarios where germline.status == TRUE we can should set foreground to median.chr. I feel would achieve the same effect of reducing the potential FP enrichment of SCNAs for these regions in a patient series (with potential exceptions in extreme aneuploidy of single chromosomes perhaps).

After some troubleshooting runs with dryclean and JaBbA, I have gotten stable results in a variety of JaBbA execution parameters. All of these runs have recovered the nearly all of the previous excluded SVs in the model.

I will continue to test but would interested to hear thoughts on this.

Best,
Patrick

[ZWael]

Hello, @pblaney
where can i found the start_wash_cycle function ?