Example code for reading GATK4 GenomicsDBImport

[lima1]

Hi Nalini,

I was playing around with this package. Thanks a lot providing these R bindings! Much appreciated.

Do you have by any chance an example for iterating through a GATK4 GenomicsDBImport database?

I'm trying to build something similar to CreateSomaticPanelOfNormals which tries to fit beta-binomial distributions to all heterozygous variants. So simply extracting the AD field for all variants and all samples.

I'll figure it out if you don't have anything in place, but I thought I ask as it seems like a fairly standard task.

Thanks a lot,
Markus

[lima1]

The new query_variants_calls function is nice and does almost all I need. This is the code I came up with:

library(genomicsdb)
library(jsonlite)
reference <- "references/hg38/Reference/sequence.fa"

workspace <- "gatk4_pon_db/"

db <- genomicsdb::connect(workspace = workspace,
    vid_mapping_file = file.path(workspace, "vidmap.json"),
    callset_mapping_file=file.path(workspace, "callset.json"),
    reference_genome = reference,
    c("DP", "AD", "AF"))

j <- jsonlite::read_json(file.path(workspace, "callset.json"))

query <- query_variant_calls(db, array="1$3682058$243843487", column_ranges=list(c(3682058,243843487)), row_ranges=list(range(sapply(j$callsets, function(x) x$row_idx))))
head(query)
  ROW     COL        SAMPLE CHROM     POS     END REF  AD    AF ALT   DP
1   1 3682335 LIB-020359pc4     1 3682336 3682336   G 617 0.430 [A] 1105
2   4 3682335 LIB-020362pc4     1 3682336 3682336   G 626 0.445 [A] 1166
3   8 3682335 LIB-020366pc4     1 3682336 3682336   G 475 0.433 [A]  891
4   9 3682335 LIB-020367pc4     1 3682336 3682336   G 766 0.325 [A] 1191

One issue I have is that Mutect2 generates a DP field in both FORMAT and INFO, one is the unfiltered one the filtered. Looks like the unfiltered is picked. That's fine, but AD is a vector with the filtered REF and ALT reads, but only the REF reads make it into the data.frame. That's a bug, right? So I cannot get the filtered ALT count from the data.frame.

And quick question, is there an easy way to cycle through the available arrays + columns, say in batches of 100,000 columns? My database is small, so reading by chromosome is fine, but wondering if there is a cleaner way.

Best,
Markus

Update iterate through all arrays:

   workspace <- normalizePath(workspace, mustWork = TRUE)
   db <- genomicsdb::connect(workspace = workspace,
       vid_mapping_file = file.path(workspace, "vidmap.json"),
       callset_mapping_file=file.path(workspace, "callset.json"),
       reference_genome = reference,
       c("DP", "AD", "AF"))

   jcallset <- jsonlite::read_json(file.path(workspace, "callset.json"))
   jvidmap <- jsonlite::read_json(file.path(workspace, "vidmap.json"))

   # get all available arrays
   arrays <- sapply(dir(workspace, full.names=TRUE), file.path, "genomicsdb_meta_dir")
   arrays <- basename(names(arrays)[which(file.exists(arrays))])
   # get offsets and lengths
   contigs <- sapply(arrays, function(ary) strsplit(ary, "\\$")[[1]][1])
   contigs <- jvidmap$contigs[match(contigs, sapply(jvidmap$contigs, function(x) x$name))]
   idx <- order(rankSeqlevels(sapply(contigs, function(x) x$name)))

   all<- lapply(idx, function(i) {
       # avoid integer overflow
       c_offset <- as.numeric(contigs[[i]]$tiledb_column_offset)
       c_length <- as.numeric(contigs[[i]]$length)

       flog.info("Processing %s (offset %.0f, length %.0f)...",
           arrays[i], c_offset, c_length)
       query <- data.table(genomicsdb::query_variant_calls(db,
           array = arrays[i],
           column_ranges = list(c(c_offset, c_offset + c_length)),
           row_ranges = list(range(sapply(jcallset$callsets,
               function(x) x$row_idx)))))
   })

[nalinigans]

@lima1, sorry I was not watching this repository. Just wanted to let you know this is work in progress. Please let me know if you have solved your issue. Otherwise I will get back with answers soon.

And quick question, is there an easy way to cycle through the available arrays + columns, say in batches of 100,000 columns?
My database is small, so reading by chromosome is fine, but wondering if there is a cleaner way.

Are you open to using spark or mpi or any parallel processing paradigm?

[lima1]

Thanks @nalinigans . The code works well enough for my purpose, but would be curious to hear your suggestions once the code matured! Thanks for all your work, it’s really nice.