Whitelist Heme Driver Junctions

[kevinmhadi]

Description of feature

Keeping track of dev here.

Addy provided a set of whitelist regions in which junctions should be forced in to the JaBbA model, if any SVs overlap.

Shihab will write up a new process that will mark these junctions to force these in and that runs only in tumor-only mode.

[kevinmhadi]

@shihabdider Here is a template for the code I used to mark junctions, which I cached in an RDS and used as input to a dev run of JaBbA.

Each iteration of x is just a 1 row entry of a data.table that contains all of the input fields to JaBbA.

heme_gen = readRDS("/gpfs/data/imielinskilab/projects/Clinical_NYU/db/whitelist_genes_heme_gen.rds")

tiered_junctions = mclapply(
  mc.cores = 10,
  jabba_params_allbatches %>% split(1:NROW(.))
 ,
  function(x) {
    ## if (!(x$pair == "WG-24-18")) return(NULL)
    jpath = x$untiered_junctions
    jpath_tiered = glue::glue('{tools::file_path_sans_ext(jpath)}___tiered.rds')
    ra.all = readRDS(jpath)
    # Important part is below
    mcols_ra.all = mcols(ra.all)
    mcols_ra.all[[x$tfield]] = rep_len(2, NROW(mcols_ra.all))
    ix = unique((grl.unlist(ra.all) %&% heme_gen)$grl.ix)
    if (NROW(ix)) 
      mcols_ra.all[[x$tfield]][ix] = 1
    mcols(ra.all) = mcols_ra.all
    saveRDS(ra.all, jpath_tiered)
    # Important part is above
    ## return(mcols(ra.all)[[x$tfield]])
    jpath_tiered
  }
)