This page gives a compact mental model for misha. Use it as the first
quick read before the full Manual vignette.
Most analyses follow the same pattern:
In misha this is usually one call to gextract,
gscreen, or gsummary.
You are not limited to raw track names. You can pass full
expressions, for example log(dense_track + 1),
dense_track / (chip.sum + 1e-6), or
pmin(dense_track, 2).
All examples below assume the bundled examples database:
A track is genomic signal organized over coordinates.
dense_track in the examples DB).Useful starter commands:
gtrack.ls() # list tracks in the examples DB
#> [1] "array_track" "dense_track" "rects_track"
#> [4] "sparse_track" "subdir.dense_track2"
gtrack.info("dense_track") # inspect type/metadata
#> $type
#> [1] "dense"
#>
#> $dimensions
#> [1] 1
#>
#> $size.in.bytes
#> [1] 80012
#>
#> $format
#> [1] "per-chromosome"
#>
#> $bin.size
#> [1] 50
gtrack.info("sparse_track")
#> $type
#> [1] "sparse"
#>
#> $dimensions
#> [1] 1
#>
#> $size.in.bytes
#> [1] 32012
#>
#> $format
#> [1] "per-chromosome"For intuition, you can think of dense_track as a
ChIP-seq-like coverage signal.
An interval set defines genomic regions
(chrom, start, end) where you
want to work.
The iterator is the stepping policy inside the scope.
iterator = 100 -> fixed 100 bp binsiterator = "some_sparse_track" -> iterate over that
track’s intervalsiterator = some_intervals_df -> iterate over
explicit regionsiterator = "my_intervals_set" -> iterate directly
over an intervals setThink of it as: scope says where, iterator says in what chunks.
out <- gextract("dense_track", regions, iterator = 100)
head(out)
#> chrom start end dense_track intervalID
#> 1 chr1 0 100 0.1688889 1
#> 2 chr1 100 200 0.1700000 1
#> 3 chr1 200 300 0.1800000 1
#> 4 chr1 300 400 0.1600000 1
#> 5 chr1 400 500 0.1100000 1
#> 6 chr1 500 600 0.0400000 1
log_out <- gextract("log(dense_track + 1)", regions, iterator = 100)
head(log_out)
#> chrom start end log(dense_track + 1) intervalID
#> 1 chr1 0 100 0.15605364 1
#> 2 chr1 100 200 0.15700375 1
#> 3 chr1 200 300 0.16551444 1
#> 4 chr1 300 400 0.14842000 1
#> 5 chr1 400 500 0.10436001 1
#> 6 chr1 500 600 0.03922071 1Create and use an intervals set as an iterator:
A virtual track is a named on-the-fly transformation, not stored as a physical track file.
Examples:
gvtrack.create("chip.sum", "dense_track", "sum")
out <- gextract("chip.sum", regions, iterator = 200)
head(out)
#> chrom start end chip.sum intervalID
#> 1 chr1 0 200 0.67777777 1
#> 2 chr1 200 400 0.68000001 1
#> 3 chr1 400 600 0.29999998 1
#> 4 chr1 600 800 0.04000000 1
#> 5 chr1 800 1000 0.09999999 1
#> 6 chr1 1000 1200 0.12000000 1You can also shift the iterator window used by the virtual track:
gvtrack.create("chip.shifted", "dense_track", "sum")
gvtrack.iterator("chip.shifted", sshift = -100, eshift = 100)
out <- gextract("chip.shifted", regions, iterator = 200)
head(out)
#> chrom start end chip.shifted intervalID
#> 1 chr1 0 200 1.037778 1
#> 2 chr1 200 400 1.240000 1
#> 3 chr1 400 600 0.660000 1
#> 4 chr1 600 800 0.140000 1
#> 5 chr1 800 1000 0.200000 1
#> 6 chr1 1000 1200 0.220000 1Here, each iterator interval is expanded by 100 bp on both sides
before evaluating dense_track.
Virtual tracks are session objects (easy to list with
gvtrack.ls and delete with gvtrack.rm).
library(misha)
gdb.init_examples()
# 1) pick scope
regions <- gintervals(1, 0, 50000)
# 2) inspect available tracks
gtrack.ls()
#> [1] "array_track" "dense_track" "rects_track"
#> [4] "sparse_track" "subdir.dense_track2"
# 3) extract signal with a chosen iterator
chip <- gextract("dense_track", regions, iterator = 100)
head(chip)
#> chrom start end dense_track intervalID
#> 1 chr1 0 100 0.1688889 1
#> 2 chr1 100 200 0.1700000 1
#> 3 chr1 200 300 0.1800000 1
#> 4 chr1 300 400 0.1600000 1
#> 5 chr1 400 500 0.1100000 1
#> 6 chr1 500 600 0.0400000 1
# 4) screen high-signal bins (as a simple peak-like filter).
# Pick the threshold from the data: dense_track only reaches 0.34 in this
# scope, so a threshold above that would silently screen nothing.
hi_chip <- gscreen("dense_track > 0.2", regions, iterator = 100)
head(hi_chip)
#> chrom start end
#> 1 chr1 17200 17300
#> 2 chr1 20000 20100
#> 3 chr1 23300 23400
#> 4 chr1 26200 26300
#> 5 chr1 32600 32800
#> 6 chr1 32900 33000
# 5) summarize distribution/coverage
gsummary("dense_track", regions, iterator = 100)
#> Total intervals NaN intervals Min Max Sum
#> 500.00000000 0.00000000 0.00000000 0.34000000 43.68888862
#> Mean Std dev
#> 0.08737778 0.06208607A PWM/PSSM is a motif model over A/C/G/T. In misha, a common pattern is:
regions <- gintervals(1, c(1000, 2000), c(1020, 2020))
seqs <- gseq.extract(regions)
seqs
#> [1] "cctcagtaatccgaaaagcc" "CTGCATGTAACTTAATACCA"
pssm <- matrix(c(
0.80, 0.05, 0.10, 0.05,
0.10, 0.10, 0.70, 0.10,
0.05, 0.80, 0.05, 0.10,
0.10, 0.10, 0.10, 0.70
), ncol = 4, byrow = TRUE)
colnames(pssm) <- c("A", "C", "G", "T")
scores <- gseq.pwm(seqs, pssm, mode = "lse")
scores
#> [1] -2.198000 -2.119781If your database has motif files under pssms/, you can
create a genome-wide PWM-energy track with
gtrack.create_pwm_energy(...).