## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----setup-------------------------------------------------------------------- library(misha) gdb.init_examples() ## ----------------------------------------------------------------------------- # The Dataset API needs more than one database. For this vignette the bundled # examples database is the working database, and a second, small database under # tempdir() plays the role of a shared, read-only annotation source. A dataset # only has to share the working database's chrom_sizes.txt. my_project <- .misha$GROOT shared_annotations <- file.path(tempdir(), "shared_annotations") unlink(shared_annotations, recursive = TRUE) dir.create(file.path(shared_annotations, "tracks"), recursive = TRUE) invisible(file.copy(file.path(my_project, "chrom_sizes.txt"), shared_annotations)) gtrack.copy("dense_track", "annotation_track", db = shared_annotations) # Set your working database gsetroot(my_project) # Load an additional dataset gdataset.load(shared_annotations) # List all sources (working db + loaded datasets) gdataset.ls() # Get detailed information gdataset.ls(dataframe = TRUE) ## ----error = TRUE------------------------------------------------------------- try({ # A second dataset that also carries a track named "dense_track": other_dataset <- file.path(tempdir(), "other_dataset") unlink(other_dataset, recursive = TRUE) dir.create(file.path(other_dataset, "tracks"), recursive = TRUE) invisible(file.copy(file.path(my_project, "chrom_sizes.txt"), other_dataset)) gtrack.copy("dense_track", db = other_dataset) gdataset.load(other_dataset) }) ## ----------------------------------------------------------------------------- # Working database always wins over datasets # (for dataset-to-dataset collisions, later-loaded wins) gdataset.load(other_dataset, force = TRUE) # Check which source provides a track gtrack.dataset("dense_track") # See all sources where a track exists (for debugging) gtrack.dbs("dense_track") gdataset.unload(other_dataset) ## ----------------------------------------------------------------------------- gsetroot(my_project) gdataset.load(shared_annotations) # Extract tracks from different sources in a single call result <- gextract(c("dense_track", "annotation_track"), gintervals(1, 0, 10000), iterator = 100) head(result) # Use track expressions across sources normalized <- gextract("dense_track - annotation_track", gintervals(1, 0, 10000), iterator = 100) head(normalized) # Access track attributes from any source gtrack.attr.get("dense_track", "description") gtrack.attr.get("annotation_track", "description") ## ----------------------------------------------------------------------------- # Create virtual tracks from different sources gvtrack.create("vt_signal", "dense_track", "avg") gvtrack.create("vt_annotation", "annotation_track", "avg") # Combine them in expressions head(gextract("vt_signal / vt_annotation", gintervals(1, 0, 10000))) ## ----------------------------------------------------------------------------- # Single track - returns source path gtrack.dataset("annotation_track") # All sources containing a track (useful for debugging shadowed tracks) gtrack.dbs("dense_track") # Multiple tracks (vectorized) gtrack.dataset(c("dense_track", "sparse_track", "annotation_track")) # Intervals gintervals.dataset("annotations") gintervals.dbs("annotations") ## ----------------------------------------------------------------------------- # List tracks from a specific source only gtrack.ls(db = shared_annotations) # List intervals from working database gintervals.ls(db = my_project) ## ----------------------------------------------------------------------------- gsetroot(my_project) # Create a dataset with selected tracks and intervals chipseq_dataset <- file.path(tempdir(), "my_chipseq_dataset") unlink(chipseq_dataset, recursive = TRUE) gdataset.save( path = chipseq_dataset, description = "ChIP-seq-like tracks", tracks = gtrack.ls("dense"), # Pattern matching intervals = "annotations" ) # Options: # - symlinks = TRUE: Create symlinks instead of copying (saves space) # - copy_seq = TRUE: Copy seq/ directory instead of symlinking ## ----------------------------------------------------------------------------- gdataset.info(chipseq_dataset) # Returns: description, author, created date, track/interval counts, genome hash ## ----------------------------------------------------------------------------- # Create linked database with symlinks to parent's seq and chrom_sizes my_db <- file.path(tempdir(), "my_db") unlink(my_db, recursive = TRUE) gdb.create_linked(my_db, parent = my_project) # Use as your working database gsetroot(my_db) # Load datasets from the parent gdataset.load(my_project) # Create your own tracks gtrack.create("my_analysis", "Analysis results", "dense_track * 2") head(gextract("my_analysis", gintervals(1, 0, 500))) ## ----error = TRUE------------------------------------------------------------- try({ # Unload a dataset (tracks/intervals become unavailable) gdataset.unload(my_project) # Safe to call even if not loaded (no error by default) gdataset.unload("/nonexistent/path") # Error if validate=TRUE and not loaded gdataset.unload("/nonexistent/path", validate = TRUE) }) ## ----------------------------------------------------------------------------- gsetroot(my_project) gtrack.copy("dense_track", "old_name") # Rename a track gtrack.mv("old_name", "new_name") # Move to a different namespace (directory) gtrack.mv("new_name", "results.new_name") gtrack.ls("new_name") gtrack.rm("results.new_name", force = TRUE) ## ----------------------------------------------------------------------------- # Copy a track within the same database gtrack.copy("dense_track", "copy_track") gtrack.exists("copy_track") gtrack.rm("copy_track", force = TRUE) # Copy from a loaded dataset to the working database gsetroot(my_project) gdataset.load(shared_annotations) gtrack.copy("annotation_track", "my_local_copy") # Copy to working db gtrack.dataset("my_local_copy") gtrack.rm("my_local_copy", force = TRUE) ## ----------------------------------------------------------------------------- gsetroot(my_project) gdataset.load(shared_annotations) # Create virtual track referencing track from working db gvtrack.create("vt1", "dense_track", "avg") # Create virtual track referencing track from dataset gvtrack.create("vt2", "annotation_track", "max") # Use both in same expression head(gextract("vt1 + vt2", gintervals(1, 0, 10000))) ## ----------------------------------------------------------------------------- gsetroot(my_project) # Works unchanged gdb.init(my_project) # Equivalent, also works gdataset.ls() ## ----------------------------------------------------------------------------- gtrack.copy("dense_track", "my_track") # Convert a track to indexed format gtrack.convert_to_indexed("my_track") # Check track format info <- gtrack.info("my_track") print(info$format) # "indexed" or "per-chromosome" gtrack.rm("my_track", force = TRUE) ## ----------------------------------------------------------------------------- # Only big interval sets are stored per chromosome, so lower the threshold to # get a big set out of the small example database (the default is 1,000,000). options(gbig.intervals.size = 10) gintervals.save("my_intervals", gscreen("dense_track > 0.3")) gintervals.save("my_2d_intervals", gextract("rects_track", gintervals.2d.all())[, 1:6]) options(gbig.intervals.size = 1e6) # Convert 1D interval set to indexed format gintervals.convert_to_indexed("my_intervals") # Convert 2D interval set to indexed format gintervals.2d.convert_to_indexed("my_2d_intervals") # Convert and remove old per-chromosome files gintervals.convert_to_indexed("my_intervals", remove.old = TRUE, force = TRUE) gintervals.rm("my_intervals", force = TRUE) gintervals.rm("my_2d_intervals", force = TRUE) ## ----------------------------------------------------------------------------- # 'annotations' is an intervals set saved in Genomic Database gintervals.intersect("annotations", gintervals(2)) ## ----------------------------------------------------------------------------- gvtrack.create("myvtrack", "dense_track") ## ----------------------------------------------------------------------------- gvtrack.create("myvtrack", "dense_track", "global.percentile") ## ----------------------------------------------------------------------------- gvtrack.create("myvtrack", "array_track", "sum") gvtrack.array.slice("myvtrack", c("col2", "col5"), "max") ## ----------------------------------------------------------------------------- gvtrack.iterator("myvtrack", sshift = -100, eshift = 200) ## ----------------------------------------------------------------------------- gvtrack.create("myvtrack", "dense_track") gvtrack.iterator("myvtrack", dim = 2) ## ----------------------------------------------------------------------------- gvtrack.create("myvtrack", "annotations", "distance") intervs <- gscreen("dense_track > 0.45") head(gextract("myvtrack", .misha$ALLGENOME, iterator = intervs)) ## ----------------------------------------------------------------------------- # Create a data frame with intervals and numeric values intervals_with_values <- data.frame( chrom = "chr1", start = c(100, 300, 500), end = c(200, 400, 600), score = c(10, 20, 30) ) # Use as value-based sparse track gvtrack.create("myvtrack", intervals_with_values, "avg") gvtrack.create("myvtrack_max", intervals_with_values, "max") ## ----------------------------------------------------------------------------- options(gbuf.size = 1) getOption("gbuf.size") options(gbuf.size = 1000) # back to the default for the rest of this vignette ## ----------------------------------------------------------------------------- gextract("dense_track", gintervals(2, 340, 520)) ## ----------------------------------------------------------------------------- intervs <- gintervals.2d(1, 200, 800, 1, 100, 1000) intervs <- rbind(intervs, gintervals.2d(1, 900, 950, 1, 0, 200)) intervs <- rbind(intervs, gintervals.2d(1, 0, 100, 1, 0, 400)) intervs <- rbind(intervs, gintervals.2d(1, 900, 950, 2, 0, 200)) intervs gintervals.2d.band_intersect(intervs, band = c(500, 1000)) ## ----------------------------------------------------------------------------- intervs <- gintervals.2d(1, c(100, 400), c(300, 490), 1, c(120, 180), c(200, 500)) gtrack.2d.create("test2d", "test 2D track", intervs, c(10, 20)) gextract("test2d", .misha$ALLGENOME) gextract("test2d", .misha$ALLGENOME, iterator = gintervals.2d(1, 0, 1000, 1, 0, 1000)) gintervals.2d.band_intersect(intervs, band = c(150, 1000)) gextract("test2d", .misha$ALLGENOME, iterator = gintervals.2d(1, 0, 1000, 1, 0, 1000), band = c(150, 1000)) gtrack.rm("test2d", force = TRUE) ## ----------------------------------------------------------------------------- set.seed(60427) r1 <- gsample("dense_track", 10) r2 <- gsample("dense_track", 10) # r2 differs from r1 set.seed(60427) r3 <- gsample("dense_track", 10) # r3 == r1 identical(r1, r2) identical(r1, r3)