# \dontshow{
options(gmax.processes = 2)
# }
gdb.init_examples()
gvtrack.create("vtrack1", "dense_track", "max")
gvtrack.create("vtrack2", "dense_track", "quantile", 0.5)
gextract("dense_track", "vtrack1", "vtrack2",
gintervals(1, 0, 10000),
iterator = 1000
)
gvtrack.create("vtrack3", "dense_track", "global.percentile")
gvtrack.create("vtrack4", "annotations", "distance")
gdist(
"vtrack3", seq(0, 1, l = 10), "vtrack4",
seq(-500, 500, 200)
)
gvtrack.create("cov", "annotations", "coverage")
gextract("cov", gintervals(1, 0, 1000), iterator = 100)
pssm <- matrix(
c(
0.7, 0.1, 0.1, 0.1, # Example PSSM
0.1, 0.7, 0.1, 0.1,
0.1, 0.1, 0.7, 0.1,
0.1, 0.1, 0.7, 0.1,
0.1, 0.1, 0.7, 0.1,
0.1, 0.1, 0.7, 0.1
),
ncol = 4, byrow = TRUE
)
colnames(pssm) <- c("A", "C", "G", "T")
gvtrack.create(
"motif_score", NULL, "pwm",
list(pssm = pssm, bidirect = TRUE, prior = 0.01)
)
gvtrack.create("max_motif_score", NULL, "pwm.max",
pssm = pssm, bidirect = TRUE, prior = 0.01
)
gvtrack.create("max_motif_pos", NULL, "pwm.max.pos",
pssm = pssm
)
gextract(
c(
"dense_track", "motif_score", "max_motif_score",
"max_motif_pos"
),
gintervals(1, 0, 10000),
iterator = 500
)
# Kmer counting examples
gvtrack.create("cg_count", NULL, "kmer.count", kmer = "CG", strand = 1)
gvtrack.create("cg_frac", NULL, "kmer.frac", kmer = "CG", strand = 1)
gextract(c("cg_count", "cg_frac"), gintervals(1, 0, 10000), iterator = 1000)
gvtrack.create("at_pos", NULL, "kmer.count", kmer = "AT", strand = 1)
gvtrack.create("at_neg", NULL, "kmer.count", kmer = "AT", strand = -1)
gvtrack.create("at_both", NULL, "kmer.count", kmer = "AT", strand = 0)
gextract(c("at_pos", "at_neg", "at_both"), gintervals(1, 0, 10000), iterator = 1000)
# GC content
gvtrack.create("g_frac", NULL, "kmer.frac", kmer = "G")
gvtrack.create("c_frac", NULL, "kmer.frac", kmer = "C")
gextract("g_frac + c_frac", gintervals(1, 0, 10000),
iterator = 1000,
colnames = "gc_content"
)
# Masked base pair counting
gvtrack.create("masked_count", NULL, "masked.count")
gvtrack.create("masked_frac", NULL, "masked.frac")
gextract(c("masked_count", "masked_frac"), gintervals(1, 0, 10000), iterator = 1000)
# Combined with GC content (unmasked regions only)
gvtrack.create("gc", NULL, "kmer.frac", kmer = "G")
gextract("gc * (1 - masked_frac)",
gintervals(1, 0, 10000),
iterator = 1000,
colnames = "gc_unmasked"
)
# Value-based track examples
# 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 (functions like sparse track)
gvtrack.create("value_track", intervals_with_values, "avg")
gvtrack.create("value_track_max", intervals_with_values, "max")
gextract(c("value_track", "value_track_max"),
gintervals(1, 0, 10000),
iterator = 1000
)
# Spatial PWM examples
# Create a PWM with higher weight in the center of intervals
pssm <- matrix(
c(
0.7, 0.1, 0.1, 0.1,
0.1, 0.7, 0.1, 0.1,
0.1, 0.1, 0.7, 0.1,
0.1, 0.1, 0.1, 0.7
),
ncol = 4, byrow = TRUE
)
colnames(pssm) <- c("A", "C", "G", "T")
# Spatial factors: low weight at edges, high in center
# For 200bp intervals with 40bp bins: bins 0, 40, 80, 120, 160
spatial_weights <- c(0.5, 1.0, 2.0, 1.0, 0.5)
gvtrack.create(
"spatial_pwm", NULL, "pwm",
list(
pssm = pssm,
bidirect = TRUE,
spat_factor = spatial_weights,
spat_bin = 40L
)
)
# Compare with non-spatial PWM
gvtrack.create(
"regular_pwm", NULL, "pwm",
list(pssm = pssm, bidirect = TRUE)
)
gextract(c("spatial_pwm", "regular_pwm"),
gintervals(1, 0, 10000),
iterator = 200
)
# Using spatial parameters with iterator shifts
gvtrack.create(
"spatial_extended", NULL, "pwm.max",
pssm = pssm,
spat_factor = c(0.5, 1.0, 2.0, 2.5, 2.0, 1.0, 0.5),
spat_bin = 40L
)
# Scan window will be 280bp (100bp + 2*90bp)
gvtrack.iterator("spatial_extended", sshift = -90, eshift = 90)
gextract("spatial_extended", gintervals(1, 0, 10000), iterator = 100)
# Using spat_min/spat_max to restrict scanning to a window
# For 500bp intervals, scan only positions 30-470 (440bp window)
gvtrack.create(
"window_pwm", NULL, "pwm",
pssm = pssm,
bidirect = TRUE,
spat_min = 30, # 1-based position
spat_max = 470 # 1-based position
)
gextract("window_pwm", gintervals(1, 0, 10000), iterator = 500)
# Combining spatial weighting with window restriction
# Scan positions 50-450 with spatial weights favoring the center
gvtrack.create(
"window_spatial_pwm", NULL, "pwm",
pssm = pssm,
bidirect = TRUE,
spat_factor = c(0.5, 1.0, 2.0, 2.5, 2.0, 1.0, 0.5, 1.0, 0.5, 0.5),
spat_bin = 40L,
spat_min = 50,
spat_max = 450
)
gextract("window_spatial_pwm", gintervals(1, 0, 10000), iterator = 500)
# --- Edit distance examples ---
# (Using the same 4-position PSSM defined above)
# First, check the PWM score distribution to pick sensible thresholds
gvtrack.create("pwm_score", NULL, "pwm.max", pssm = pssm)
head(gextract("pwm_score", gintervals(1, 500, 520), iterator = 1))
# Scores range from ~-8.3 (poor) to ~-0.8 (strong match)
# --- direction = "above": how many edits to CREATE a motif match? ---
# "How many substitutions until this window scores >= -5?"
gvtrack.create("edist_above", NULL, "pwm.edit_distance",
pssm = pssm, score.thresh = -5
)
# At 1bp: each position gets its own edit distance (0 if already matching)
head(gextract(c("pwm_score", "edist_above"),
gintervals(1, 500, 520),
iterator = 1
))
# --- direction = "below": how many edits to DISRUPT a motif match? ---
#
# score.thresh is the target: "push the score below this value"
# score.min is a pre-filter: "only consider windows currently scoring above this"
# These are independent — score.min controls WHICH windows to evaluate,
# score.thresh controls WHAT score to target.
#
# Example: among windows scoring >= -5 (strong matches),
# how many edits to push below -7 (disrupt the match)?
gvtrack.create("edist_below", NULL, "pwm.edit_distance",
pssm = pssm, score.thresh = -7, score.min = -5,
direction = "below"
)
head(gextract(c("pwm_score", "edist_above", "edist_below"),
gintervals(1, 500, 520),
iterator = 1
))
# score >= -5 (e.g. -2.7): gets edit distance (2 edits to go below -7)
# score < -5 (e.g. -6.4): NA (filtered out by score.min)
# Without score.min, ALL windows are evaluated — windows already
# below -7 return 0 edits:
gvtrack.create("edist_below_all", NULL, "pwm.edit_distance",
pssm = pssm, score.thresh = -7, direction = "below"
)
head(gextract(c("pwm_score", "edist_below", "edist_below_all"),
gintervals(1, 500, 520),
iterator = 1
))
# score >= -5: same result as edist_below (2 edits)
# score between -7 and -5 (e.g. -6.4): edist_below=NA, edist_below_all=0
# score < -7: both return 0 (already below threshold)
# --- max_edits: cap the search depth ---
# Return NA if more than 2 substitutions are needed
gvtrack.create("edist_max2", NULL, "pwm.edit_distance",
pssm = pssm, score.thresh = -5, max_edits = 2
)
# Positions needing 3+ edits now return NA instead of 3
head(gextract(c("pwm_score", "edist_above", "edist_max2"),
gintervals(1, 500, 520),
iterator = 1
))
# --- Aggregation over intervals ---
# With a larger iterator, the minimum edit distance across the
# interval is returned.
gextract(c("pwm_score", "edist_above", "edist_below_all"),
gintervals(1, 0, 2000),
iterator = 200
)
Run the code above in your browser using DataLab