R/shaman.R
shaman_score_hic_track.Rdshaman_score_hic_track
shaman_score_hic_track(
track_db,
work_dir,
score_track_nm,
obs_track_nms,
exp_track_nms = paste0(obs_track_nms, "_shuffle"),
points_track_nms = obs_track_nms,
near_cis = 5000000,
expand = 2000000,
k = 100,
max_jobs = 100
)Directory of the misha database.
Centralized directory to store temporary files.
Score track that will be created.
Names of observed 2D genomic tracks for the hic data. Pooling of multiple observed tracks is supported.
Names of expected (shuffled) 2D genomic tracks. Pooling of multiple expected tracks is supported.
Names of 2D genomic tracks that contain points on which to compute normalized score. Pooling points from multiple tracks is supported.
Size of matrix in grid.
Size of expansion, points to include outside the matrix for accurate computing of the score. Note that for each observed point, its k-nearest neighbors must be included in the expanded matrix.
The number of neighbor distances used for the score. For higher resolution maps, increase k. For lower resolution maps, decrease k.
Maximal number of qsub jobs.
No return value, called for side effects: creates the 2D score track score_track_nm in track_db.
This function generates a 2D score track based on observed and expected hic data.
The score is computed by generating a grid of small matrices spanning all chromosomes
and computing the score of each matrix independently.
The model for computing the score relies on the KS D statistic computed for each observed point,
over the distances of the k-nearest neighbors in the observed compared to the expected.
High scores represent contact enrichment while low scores depict insulation.
Note that this function requires either sge (qsub) or multicore to compute in a timely manner.
Parameters can be set via shaman.sge_support or shaman.mc_support in shaman.conf file.
Sun Grid Engine mode, options(shaman.sge_support = 1), is preferred; work_dir must then be
accessible by all jobs. Multi-core mode is options(shaman.mc_support = 1).
options(shaman.score.threads = N) computes the kNN distances and scores of each matrix on N
threads (default 1); in multi-core mode each of the max_jobs processes uses N threads, and in SGE mode
each job does, so shaman.sge_flags should ask for N slots.
Each extraction of a 2D track can keep about the whole chromosome's track file in memory;
misha's options(gtrack.num.chunks = 1000) bounds that (2-3x lower peak memory on dense
matrices in our runs) without changing the scores.
Each step creates temporary files of the matrix scores which are then joined to a track. Temporary files are deleted upon track creation. Matrices that still have no score after 3 rounds of jobs stop the run with an error that names them. The finished score files stay in work_dir, so a rerun with the same work_dir computes only the missing ones.
# The example below runs on the test misha db provided with shaman.
# Note that this is a toy db sampled from K562 ela data -
# scoring based on the observed and expected tracks will not produce the score track,
# as most of the genome is missing.
library(misha)
track_db <- shaman_get_test_track_db()
gsetroot(track_db)
old_opts <- options(shaman.mc_support = 1) # configuring multi-core mode
if (gtrack.exists("hic_score_new")) {
gtrack.rm("hic_score_new", force = TRUE)
gdb.reload()
}
ret <- shaman_score_hic_track(track_db,
# work_dir can be tempdir() only in multi-core mode.
# For sge mode, work_dir must be accessible by all jobs.
work_dir = tempdir(),
score_track_nm = "hic_score_new",
obs_track_nms = "hic_obs",
exp_track_nms = "hic_exp",
near_cis = 1e09, # this test db contains very little data, can increase the size of each job
k = 20, # default is set to 100
max_jobs = 2
) # increase number of jobs for optimal runtime when running in sge mode
#> missing 1 score files
#> kk norm on 90300 points
#> obs = hic_obs
#> n_obs = 97618, n_exp = 200292, k_exp = 40
#> missing 0 score files
gdb.reload()
gtrack.ls("hic_score_new") # new score track that was created
#> [1] "hic_score_new"
options(old_opts)