karyoSim: karyoSim

In mskilab/skitools: Various mskilab R utilties

karyoSim

Description

Simulate (random) evolution of rearrangements according to input junctions, which are provided as a GRangesList, and grouped into "events" by events list (list of numeric vectors or of lists of numeric vectors indexing "junctions")

Goal of the simulation is to instantiate a collection of junctions (+/- approach some copy number profile) through a sequence of rearrangements and whole-chromosome copy changes

Junctions are sequences of signed reference intervals that are contiguous on the on the tumor genome (usually pairs)

Each event consists of either (1) a "quasi reciprocal sequence" (QRS) of junctions, implemented during a single "cell cycle", and are specified by vectors of junction indices, (2) a set of sets of junctions, specified as a list of list of junction indices, again without repetition, corresponding to complex events spanning multiple QRS's or "cell cycles" e.g. a BFB, which involve a replication step in between each QRS. The subsequent QRS's (attempted to be) instantiated in cis to the last item in the previous QRS (3) a GRanges object specifying a reference locus and meta data field $type = "loss" or "gain" specifying one or more pieces of reference genome that should be lost or gained at a given step.

- Events are interpreted as strings of one or more "quasi reciprocal sequence" (QRS) of junctions which may be closed / cyclic (if they begin and end with the same junction index) or open. in which case they will result in at least some interval loss. We restrict QRS's to contain at most one repeated junction, and this has to be the first and the last item in the sequence. "Quasi" reciprocal means that we allow some sequence to be lost or gained in between breaks. - Every QRS is instantiated in the current genome, by mapping junctions, which are specified in haploid reference coordinates to intervals on the current genome. By default, instantiation is chosen so that the source interval of every junction in the QRS is on the same chromosome as the target interval on the previous junction in the QRS. If this is the case, then we say that the current junction "follows" the previous one in this WQRS instantiation. Quasi-reciprocality is then applied by possibly adding intervals at the site of a break (i.e. if the target interval of the previous junction is upstream of the source interval of the next junction). In situations where an instance of a subsequent junction cannot be found to followthe current junction, then the chain is either (1) interpreted as "broken", i.e. equiv to an unbalanced rearrangement (if strict.chain = F or (2) the event is discarded (if strict.chain = T) - Junctions / events can have many possible instantiations at a given round of evolution. This is because a given haploid interval on the reference can be associated with many loci on the tumor chromosome (in the simplest case, two homologues of the same chromosome) By default, the following preferences are exercised for choosing junction instantiations: (1) if a chromosome strand can be found that contains all the intervals in the junction (2) a chromosome whose both strands contain all the intervals in the junction (3) a set of chromosome that instantiates the event as a chain of junctions these prefereences can be over-ridden by specifying instant.local and instant.chain flags - After every cycle we do a "clean up" which involves (1) rejoining any pairs of broken ends that were partners at the previous iteration (2) removing any fragments that lack a telomere (if req.tel = T) or lack other req.gr (3) replacing reverse complements of chromosomes from previous iteration that were rearranged in the previous iteration with the reverse complements of their alteration products in the current iteration. - Every junction is implemented <exactly once> during the evolutionary history, i.e. lightning does not strike twice, infinite sites model

p.chrom = prob of chrom event at each simulation step p.chromloss = probability of chromosomal loss | chrom event (default 0.5) p.chromgain = probability of chromosomal gain | chrom event (default 0.5) lambda.chrom = poisson lambda of number of different chromosomes gained or lost at a chromosomal event lambda.chromgain = poisson lambda of number of chromosomes gained at each "gain" event (default lambda.chrom) labmda.chromloss = poisson lambda of number of chromosomes lostd at each "loss" event (default lambda.chrom)

p.wgd = prob of whole genome doubling at each simulation step

Optionally can provide a copy profile cn.profile (GRanges tiling genome with $cn meta data field) and heuristic will be applied to attempt to "evolve" the simulation towards the observed copy profile (to be implemented)

Output is provided as - (if full = F) list with fields $chroms = Named GRangesList of final tumor chromosomes in terms of reference intervals $gChain = gChain mapping reference to tumor genome $cn = gRanges in reference genome of copy counts of reference intervals $events = data frame of event indices with field $id (for event id), $desc (see below for description) - (if full = T) List of lists, each item k corresponding to each stage k of evolution and containing the following fields: $chroms = Named GRangesList of tumor chromosomes at that stage in terms of reference intervals $gChain = gChain mapping reference to current genome k $gChain.last = gChain mapping last evolution step to current (from reference in first item of history) $cn = gRanges in reference genome of copy counts of reference intervals $event = list with $id, $desc that gives the id and description of event for chromosomal loss / gain $id = 'chromgain', or 'chromloss', $desc = indices chromosomes for ra event, $id event id, desc = junctions involved

done.events = rowSums(event2junction)==0

Usage

karyoSim(
  junctions,
  events = NULL,
  p.chrom = 0,
  p.wgd = 0,
  p.chromgain = 0.5,
  cn = NULL,
  req.gr = NULL,
  req.tel = TRUE,
  neo.tel = NULL,
  haploid = T,
  local.junction = T,
  local.qrs = T,
  force.event = T,
  lambda.chrom = 0,
  lambda.chromgain = lambda.chrom,
  lambda.chromloss = lambda.chrom,
  full = F,
  random.event = T,
  precedence = NULL,
  dist = 1000,
  verbose = T,
  ...
)

mskilab/skitools documentation built on Aug. 31, 2023, 1:13 p.m.