R/shaman.R
shaman_shuffle_hic_mat_for_track.Rdshaman_shuffle_hic_mat_for_track
shaman_shuffle_hic_mat_for_track(
track_db,
track,
work_dir,
chrom,
start1,
end1,
start2,
end2,
min_dist = 1024,
max_dist = max(gintervals.all()$end),
dist_resolution = NA,
decay_smooth = NA,
proposal_iterations = 10000000,
shuffle = 80,
hic_mcmc_max_resolution = 400,
raw_ext = "raw",
shuffled_ext = "shuffled",
grid_small = 500000,
grid_high = 1000000,
grid_increase = 500000,
grid_step_iter = 40,
sort_uniq = FALSE,
seed = NULL
)Directory of the misha database.
Name of observed 2D genomic track for the hic data.
Centralized directory to store temporary files.
The chromosome of the matrix.
The start coordinate of the first dimension.
The end coordinate of the first dimension.
The start coordinate of the second dimension.
The end coordinate of the second dimension.
The minimum distance between contact end points.
The maximum distance between contact end points.
Number of bins in each log2 distance unit. If NA, value is determined based on observed data (recommended).
Number of bins to use for smoothing the MCMC target function: the decay curve. If NA, value is determined based on observed data (recommended).
Number of MCMC sampling iterations between proposal corrections.
Number of shuffling rounds for each observed point.
Maximum number of bins for each log2 unit
File extension of the observed data.
File extension of the shuffled data.
Initial size of maximum distance between contact pairs considered for switching
Final size of maximum distance between contact pairs considered for switching
Grid increase size
Number of iterations in each grid size
Binary flag, indicating whether the shuffled matrix file should be sorted and contacts combined. This is required prior to importing the track to misha, and should be applied to full chromosomes only.
Seed for the shuffler: NULL (default) seeds it from the current time (in seconds), or a whole number between 0 and 2^31 - 1. The shuffler prints the seed it uses.
The seed the shuffler used, or NA if the matrix was not shuffled here (no or too few contacts, or the shuffled file already existed).
This function generates an expected 2D hic matrix from observed hic data. Should not be called externally, but rather as jobs when creating a 2D shuffled track. The observed data is a combination of observed contacts in scope plus already shuffled near-cis contacts (stored in work_dir) which we sample from to maintain the decay probability curve.
Each step creates temporary files of the shuffled matrices which are then joined to a track. Temporary files are deleted upon track creation.
# The example below runs on the test misha db provided with shaman.
# Called by shaman_shuffle_hic_track() for each chromosome; here it shuffles the contacts
# of chr2:176.5e06-177e06.
library(misha)
track_db <- shaman_get_test_track_db()
work_dir <- tempfile("shaman_shuffle_")
dir.create(work_dir)
seed <- shaman_shuffle_hic_mat_for_track(track_db, "hic_obs", work_dir,
chrom = "chr2", start1 = 176.5e06, end1 = 177e06, start2 = 176.5e06, end2 = 177e06,
shuffle = 1, # default is set to 80
grid_step_iter = 1, # default is set to 40
seed = 1, # the same seed gives the same shuffle
sort_uniq = TRUE
)
#> writing 48810 raw data, dist resolution= 13
seed
#> [1] 1
# the shuffled contacts (in both orientations) and their counts
shuffled <- read.delim(list.files(work_dir, pattern = "uniq$", full.names = TRUE))
head(shuffled)
#> chrom1 start1 end1 chrom2 start2 end2 obs
#> 1 chr2 176500017 176500018 chr2 176501976 176501977 1
#> 2 chr2 176500017 176500018 chr2 176503835 176503836 1
#> 3 chr2 176500017 176500018 chr2 176504088 176504089 1
#> 4 chr2 176500017 176500018 chr2 176506177 176506178 1
#> 5 chr2 176500017 176500018 chr2 176507402 176507403 2
#> 6 chr2 176500017 176500018 chr2 176507863 176507864 1
unlink(work_dir, recursive = TRUE)