gsynth.sample: Sample a synthetic genome from a trained Markov model

In misha: Toolkit for Analysis of Genomic Data

Sample a synthetic genome from a trained Markov model

Description

Generates a synthetic genome by sampling from a trained stratified Markov-5 model. The generated genome preserves the k-mer statistics of the original genome within each stratification bin.

Usage

gsynth.sample(
  model,
  output_path = NULL,
  output_format = c("misha", "fasta", "vector"),
  mask_copy = NULL,
  seed = NULL,
  intervals = NULL,
  n_samples = 1,
  bin_merge = NULL
)

Arguments

model	A gsynth.model object from gsynth.train
output_path	Path to the output file (ignored when output_format = "vector")
output_format	Output format: "misha": .seq binary format (default) "fasta": FASTA text format "vector": Return sequences as a character vector (does not write to file)
mask_copy	Optional intervals to copy from the original genome instead of sampling. Use this to preserve specific regions (e.g., repeats, regulatory elements) exactly as they appear in the reference. Regions not in mask_copy will be sampled using the Markov model. Note: mask_copy intervals should be non-overlapping and sorted by start position within each chromosome. Overlapping intervals may result in only the first overlapping region being copied, with subsequent overlaps skipped due to cursor advancement during sequential processing.
seed	Random seed for reproducibility. If NULL, uses current random state.
intervals	Genomic intervals to sample. If NULL, uses all chromosomes.
n_samples	Number of samples to generate per interval. Default is 1. When n_samples > 1 and output_format = "fasta", headers include "_sampleN". When output_format = "vector", returns n_samples * n_intervals sequences.
bin_merge	Optional list of bin merge specifications to apply during sampling, one per dimension (length must equal model$n_dims). Each element should be: A list of merge specifications (same format as in gsynth.train: each spec is list(from = ..., to = ...)) Or NULL to use the bin mapping from training for that dimension This allows merging sparse bins at sampling time without re-training. Example for a 2D model: bin_merge = list( # Dimension 1: merge bins with values >= 0.8 to bin [0.7, 0.8) list(list(from = c(0.8, Inf), to = c(0.7, 0.8))), # Dimension 2: use training-time bin_map (no override) NULL )

Details

N bases during sampling: When the sampler needs to initialize the first 5-mer context and encounters regions with only N bases, it falls back to uniform random base selection until a valid context is established. Similarly, if a bin has no learned statistics (sparse bin with NA CDF), uniform sampling is used for that position.

Sparse bins: If the model has sparse bins (from min_obs during training), a warning is issued when sampling regions that fall into these bins. Consider using bin_merge to redirect sparse bins to well-populated ones.

Value

When output_format is "misha" or "fasta", returns invisible NULL and writes the synthetic genome to output_path. When output_format is "vector", returns a character vector of sequences (length = n_intervals * n_samples).

See Also

gsynth.train, gsynth.save

Examples

gdb.init_examples()

# Create virtual tracks
gvtrack.create("g_frac", NULL, "kmer.frac", kmer = "G")
gvtrack.create("c_frac", NULL, "kmer.frac", kmer = "C")
gvtrack.create("cg_frac", NULL, "kmer.frac", kmer = "CG")
gvtrack.create("masked_frac", NULL, "masked.frac")

# Define repeat mask (regions to preserve from original)
repeats <- gscreen("masked_frac > 0.5",
    intervals = gintervals.all(),
    iterator = 100
)

# Train model (excluding repeats from training)
model <- gsynth.train(
    list(expr = "g_frac + c_frac", breaks = seq(0, 1, 0.025)),
    list(expr = "cg_frac", breaks = c(0, 0.01, 0.02, 0.03, 0.04, 0.2)),
    mask = repeats,
    iterator = 200,
    min_obs = 1000
)

# Sample with mask_copy to preserve repeats from original genome
temp_dir <- tempdir()
synthetic_genome_file <- file.path(temp_dir, "synthetic_genome.fa")
gsynth.sample(model, synthetic_genome_file,
    output_format = "fasta",
    mask_copy = repeats,
    seed = 60427,
    bin_merge = list(
        list(list(from = 0.7, to = c(0.675, 0.7))),
        list(list(from = 0.04, to = c(0.03, 0.04)))
    )
)

misha documentation built on Jan. 28, 2026, 1:07 a.m.