Package 'jackalope'

Description

Random chromosomes are generated to create a new ref_genome object. Note that this function will never generate empty chromosomes.

Usage

create_genome(
  n_chroms,
  len_mean,
  len_sd = 0,
  pi_tcag = rep(0.25, 4),
  n_threads = 1
)

Arguments

Value

A ref_genome object.

Examples

genome <- create_genome(10, 100e3, 100, pi_tcag = c(0.1, 0.2, 0.3, 0.4))

Description

Uses one of multiple methods to create variant haplotypes from a reference genome. See haps_functions for the methods available.

Usage

create_haplotypes(
  reference,
  haps_info,
  sub = NULL,
  ins = NULL,
  del = NULL,
  epsilon = 0.03,
  n_threads = 1,
  show_progress = FALSE
)

Arguments

Value

A haplotypes object.

References

Cao, Y., D. T. Gillespie, and L. R. Petzold. 2006. Efficient step size selection for the tau-leaping simulation method. The Journal of Chemical Physics 124(4): 044109.

Doob, J. L. 1942. Topics in the theory of markoff chains. Transactions of the American Mathematical Society 52(1): 37–64.

Gillespie, D. T. 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics 22(4): 403–434.

Wieder, N., R. H. Fink, and F. von Wegner. 2011. Exact and approximate stochastic simulation of intracellular calcium dynamics. Journal of Biomedicine and Biotechnology 2011: 572492.

Examples

r <- create_genome(10, 1000)
v_phylo <- create_haplotypes(r, haps_phylo(ape::rcoal(5)), sub_JC69(0.1))
v_theta <- create_haplotypes(r, haps_theta(0.001, 5), sub_K80(0.1, 0.2))

Description

From Table 1 in Sung et al. (2016).

Usage

evo_rates

Format

A data frame with 15 rows and 8 variables:

domain

Either Bacteria or Eukarya for what type of organism the species is.

species

Species name.

Ge

Effective genome size using only coding DNA.

Gc_Gnc

Effective genome size using coding DNA and non-coding DNA that is under purifying selection.

indels

Rate of insertions and deletions ($10^{-10}$ events per site per generation).

subs

Base-substitution mutation rate ($10^{-10}$ events per site per generation).

Ne

Effective population size ($\times 10^{6}$).

theta_s

Population mutation rate estimated using $\theta_s$.

pi_s

Population mutation rate estimated using $\pi_s$.

Source

doi:10.1534/g3.116.030890

References

Sung, W., M. S. Ackerman, M. M. Dillon, T. G. Platt, C. Fuqua, V. S. Cooper, and M. Lynch. 2016. Evolution of the insertion-deletion mutation rate across the tree of life. G3: Genes | Genomes | Genetics 6:2583–2591.

Description

Interactive wrapper for a pointer to a C++ object that stores information about variant haplotypes from a single reference genome.

Details

This class should NEVER be created using haplotypes$new. Only use create_haplotypes. Because this class wraps a pointer to a C++ object, there are no fields to manipulate directly. All manipulations are done through this class's methods.

Connections to ref_genome objects

Regarding the ref_genome object you use to create a haplotypes object, you should note the following:

• This point is the most important. Both the ref_genome and haplotypes objects use the same underlying C++ object to store reference genome information. Thus, if you make any changes to the ref_genome object, those changes will also show up in the haplotypes object. For example, if you make a haplotypes object named V based on an existing ref_genome object named R, then you merge chromosomes in R, V will now have merged chromosomes. If you've already started adding mutations to V, then all the indexes used to store those mutations will be inaccurate. So when you do anything with V later, your R session will crash or have errors. The lesson here is that you shouldn't edit the reference genome after using it to create haplotypes.

• If a ref_genome object is used to create a haplotypes object, deleting the ref_genome object won't cause issues with the haplotypes object. However, the haplotypes class doesn't provide methods to edit chromosomes, so only remove the ref_genome object when you're done editing the reference genome.

Methods

Public methods

• haplotypes$new()

• haplotypes$print()

• haplotypes$ptr()

• haplotypes$n_chroms()

• haplotypes$n_haps()

• haplotypes$sizes()

• haplotypes$chrom_names()

• haplotypes$hap_names()

• haplotypes$chrom()

• haplotypes$gc_prop()

• haplotypes$nt_prop()

• haplotypes$set_names()

• haplotypes$add_haps()

• haplotypes$dup_haps()

• haplotypes$rm_haps()

• haplotypes$add_sub()

• haplotypes$add_ins()

• haplotypes$add_del()

Method new()

Do NOT use this; only use create_haplotypes to make new haplotypes.

Usage

haplotypes$new(genomes_ptr, reference_ptr)

Arguments

Method print()

Print a haplotypes object.

Usage

haplotypes$print()

Method ptr()

View pointer to underlying C++ object (this is not useful to end users).

Usage

haplotypes$ptr()

Returns

An externalptr object.

Method n_chroms()

View number of chromosomes.

Usage

haplotypes$n_chroms()

Returns

Integer number of chromosomes.

Method n_haps()

View number of haplotypes.

Usage

haplotypes$n_haps()

Returns

Integer number of haplotypes.

Method sizes()

View chromosome sizes for one haplotype.

Usage

haplotypes$sizes(hap_ind)

Arguments

Returns

Integer vector of chromosome sizes for focal haplotype.

Method chrom_names()

View chromosome names.

Usage

haplotypes$chrom_names()

Returns

Character vector of chromosome names.

Method hap_names()

View haplotype names.

Usage

haplotypes$hap_names()

Returns

Character vector of haplotype names.

Method chrom()

View one haplotype chromosome.

Usage

haplotypes$chrom(hap_ind, chrom_ind)

Arguments

Returns

A single string representing the chosen haplotype chromosome's DNA sequence.

Method gc_prop()

View GC proportion for part of one haplotype chromosome.

Usage

haplotypes$gc_prop(hap_ind, chrom_ind, start, end)

Arguments

Returns

A double in the range ⁠[0,1]⁠ representing the proportion of DNA sequence that is either G or C.

Method nt_prop()

View nucleotide content for part of one haplotype chromosome

Usage

haplotypes$nt_prop(nt, hap_ind, chrom_ind, start, end)

Arguments

Returns

A double in the range ⁠[0,1]⁠ representing the proportion of DNA sequence that is nt.

Method set_names()

Change haplotype names.

Usage

haplotypes$set_names(new_names)

Arguments

Returns

This R6 object, invisibly.

Method add_haps()

Add one or more blank, named haplotypes

Usage

haplotypes$add_haps(new_names)

Arguments

Returns

This R6 object, invisibly.

Method dup_haps()

Duplicate one or more haplotypes by name.

Usage

haplotypes$dup_haps(hap_names, new_names = NULL)

Arguments

Returns

This R6 object, invisibly.

Method rm_haps()

Remove one or more haplotypes by name.

Usage

haplotypes$rm_haps(hap_names)

Arguments

Returns

This R6 object, invisibly.

Method add_sub()

Manually add a substitution.

Usage

haplotypes$add_sub(hap_ind, chrom_ind, pos, nt)

Arguments

Returns

This R6 object, invisibly.

Method add_ins()

Manually add an insertion.

Usage

haplotypes$add_ins(hap_ind, chrom_ind, pos, nts)

Arguments

Returns

This R6 object, invisibly.

\item{`add_del(hap_ind, chrom_ind, pos, n_nts)`}{Manually add a deletion
    for a given haplotype (`hap_ind`), chromosome (`chrom_ind`), and position (`pos`).
    The designated number of nucleotides to delete (`n_nts`) will be deleted
    starting at `pos`, unless `pos` is near the chromosome end and doesn't have
    `n_nts` nucleotides to remove; it simply stops at the chromosome end in
    this case.}

Method add_del()

Manually add a deletion.

Usage

haplotypes$add_del(hap_ind, chrom_ind, pos, n_nts)

Arguments

Returns

This R6 object, invisibly.

See Also

create_haplotypes

Description

The following functions organize information that gets passed to create_haplotypes to generate haplotypes from a reference genome. Each function represents a method of generation and starts with "haps_". The first three are phylogenomic methods, and all functions but haps_vcf will use molecular evolution information when passed to create_haplotypes.

Details

haps_theta

Uses an estimate for theta, the population-scaled mutation rate, and a desired number of haplotypes.

haps_phylo

Uses phylogenetic tree(s) from phylo object(s) or NEWICK file(s), one tree per chromosome or one for all chromosomes.

haps_gtrees

Uses gene trees, either in the form of an object from the scrm or coala package or a file containing output in the style of the ms program.

haps_ssites

Uses matrices of segregating sites, either in the form of scrm or coala coalescent-simulator object(s), or a ms-style output file.

haps_vcf

Uses a haplotype call format (VCF) file that directly specifies haplotypes.

See Also

create_haplotypes

Description

This function organizes higher-level information for creating haplotypes from gene trees output from coalescent simulations. Note that all gene trees must be rooted and binary.

Usage

haps_gtrees(obj = NULL, fn = NULL)

write_gtrees(gtrees, out_prefix)

Arguments

Details

Using the obj argument is designed after the trees fields in the output from the scrm and coala packages. (These packages are not required to be installed when installing jackalope.) To get gene trees, make sure to add + sumstat_trees() to the coalmodel for coala, or make sure that "-T" is present in args for scrm. If using either of these packages, I encourage you to cite them. For citation information, see output from citation("scrm") or citation("coala").

If using an output file from a command-line program like ms/msms, add the -T option.

Value

A haps_gtrees_info object containing information used in create_haplotypes to create variant haplotypes. This class is just a wrapper around a list of NEWICK tree strings, one for each gene tree, which you can view (but not change) using the object's trees() method.

Functions

• write_gtrees(): Write gene trees to ms-style output file.

Description

This function organizes higher-level information for creating haplotypes from phylogenetic tree(s) output as phylo or multiPhylo objects (both from the ape package) or NEWICK files. Note that all phylogenetic trees must be rooted and binary. If using this function, I encourage you to cite ape. For citation information, see output from citation("ape").

Usage

haps_phylo(obj = NULL, fn = NULL)

Arguments

Details

See ?ape::write.tree for writing phylogenies to an output file.

Value

A haps_phylo_info object containing information used in create_haplotypes to create variant haplotypes. This class is just a wrapper around a list containing phylogenetic tree information for each reference chromosome, which you can view (but not change) using the object's phylo() method.

Description

This function organizes higher-level information for creating haplotypes from matrices of segregating sites output from coalescent simulations.

Usage

haps_ssites(obj = NULL, fn = NULL)

Arguments

Details

For what the seg_sites field should look like in a list, see output from the scrm or coala package. (These packages are not required to be installed when installing jackalope.) If using either of these packages, I encourage you to cite them. For citation information, see output from citation("scrm") or citation("coala").

Value

A haps_ssites_info object containing information used in create_haplotypes to create variant haplotypes. This class is just a wrapper around a list of matrices of segregating site info, which you can view (but not change) using the object's mats() method.

Description

This function organizes higher-level information for creating haplotypes from the population-scaled mutation rate and a desired number of haplotypes.

Usage

haps_theta(theta, n_haps)

Arguments

Value

A haps_theta_info object containing information used in create_haplotypes to create variant haplotypes. This class is just a wrapper around a list containing the phylogenetic tree and theta parameter, which you can view (but not change) using the object's phylo() and theta() methods, respectively.

Description

This function organizes higher-level information for creating haplotypes from Variant Call Format (VCF) files.

Usage

haps_vcf(fn, print_names = FALSE)

Arguments

Value

A haps_vcf_info object containing information used in create_haplotypes to create variant haplotypes. This class is just a wrapper around a list containing the arguments to this function, which you can view (but not change) using the object's fn() and print_names() methods.

Description

From either a reference genome or set of variant haplotypes, create Illumina reads from error profiles and write them to FASTQ output file(s). I encourage you to cite the reference below in addition to jackalope if you use this function.

Usage

illumina(obj,
         out_prefix,
         n_reads,
         read_length,
         paired,
         frag_mean = 400,
         frag_sd = 100,
         matepair = FALSE,
         seq_sys = NULL,
         profile1 = NULL,
         profile2 = NULL,
         ins_prob1 = 0.00009,
         del_prob1 = 0.00011,
         ins_prob2 = 0.00015,
         del_prob2 = 0.00023,
         frag_len_min = NULL,
         frag_len_max = NULL,
         haplotype_probs = NULL,
         barcodes = NULL,
         prob_dup = 0.02,
         sep_files = FALSE,
         compress = FALSE,
         comp_method = "bgzip",
         n_threads = 1L,
         read_pool_size = 1000L,
         show_progress = FALSE,
         overwrite = FALSE)

Arguments

Value

Nothing is returned.

Sequencing profiles

This section outlines how to use the seq_sys, profile1, and profile2 arguments. If all arguments are NULL (their defaults), a sequencing system is chosen based on the read length. If, however, one or more arguments has been provided, then how they're provided should depend on whether you want single- or paired-end reads.

For single-end reads

• profile2 should be NULL.

• Only seq_sys or profile1 should be provided, not both.

For paired-end reads

• If providing seq_sys, don't provide either profile1 or profile2.

• If providing profile1, you must also provide profile2 (they can be the same if you want) and you cannot provide seq_sys.

Sequencing systems

Sequencing system options are the following, where, for each system, "name" is the full name, "abbrev" is the abbreviated name, "max_len" indicates the maximum allowed read length, and "paired" indicates whether paired-end sequencing is allowed.

ID lines

The ID lines for FASTQ files are formatted as such:

⁠@<genome name>-<chromosome name>-<starting position>-<strand>[/<read#>]⁠

where the part in ⁠[]⁠ is only for paired-end Illumina reads, and where ⁠genome name⁠ is always REF for reference genomes (as opposed to haplotypes).

References

Huang, W., L. Li, J. R. Myers, and G. T. Marth. 2012. ART: a next-generation sequencing read simulator. Bioinformatics 28:593–594.

Examples

rg <- create_genome(10, 100e3, 100)
dir <- tempdir(TRUE)
illumina(rg, paste0(dir, "/illumina_reads"), n_reads = 100,
         read_length = 100, paired = FALSE)

Description

Construct necessary information for insertions and deletions (indels) that will be used in create_haplotypes.

Usage

indels(rate, max_length = 10, a = NULL, rel_rates = NULL)

Arguments

Details

All indels require the rate parameter, which specifies the overall indels rate among all lengths. The rate parameter is ultimately combined with a vector of relative rates among the different lengths of indels from 1 to the maximum possible length. There are three different ways to specify/generate relative-rate values.

1. Assume that rates are proportional to exp(-L) for indel length L from 1 to the maximum length (Albers et al. 2011). This method will be used if the following arguments are provided:

   • rate

   • max_length

2. Generate relative rates from a Lavalette distribution (Fletcher and Yang 2009), where the rate for length L is proportional to {L * max_length / (max_length - L + 1)}^(-a). This method will be used if the following arguments are provided:

   • rate

   • max_length

   • a

3. Directly specify values by providing a numeric vector of relative rates for each insertion/deletion length from 1 to the maximum length. This method will be used if the following arguments are provided:

   • rate

   • rel_rates

Value

An indel_info object, which is an R6 class that wraps the info needed for the create_haplotypes function. It does not allow the user to directly manipulate the info inside, as that should be done using this function. You can use the rates() method to view the indel rates by size.

References

Albers, C. A., G. Lunter, D. G. MacArthur, G. McVean, W. H. Ouwehand, and R. Durbin. 2011. Dindel: accurate indel calls from short-read data. Genome Research 21:961–973.

Fletcher, W., and Z. Yang. 2009. INDELible: a flexible simulator of biological sequence evolution. Molecular Biology and Evolution 26:1879–1888.

Examples

# relative rates are proportional to `exp(-L)` for indel
# length `L` from 1 to 5:
indel_rates1 <- indels(0.1, max_length = 5)

# relative rates are proportional to Lavalette distribution
# for length from 1 to 10:
indel_rates2 <- indels(0.2, max_length = 10, a = 1.1)

# relative rates are all the same for lengths from 1 to 100:
indel_rates3 <- indels(0.2, rel_rates = rep(1, 100))

Description

From either a reference genome or set of variant haplotypes, create PacBio reads and write them to FASTQ output file(s). I encourage you to cite the reference below in addition to jackalope if you use this function.

Usage

pacbio(obj,
       out_prefix,
       n_reads,
       chi2_params_s = c(0.01214, -5.12, 675, 48303.0732881,
                         1.4691051212330266),
       chi2_params_n = c(0.00189237136, 2.53944970, 5500),
       max_passes = 40,
       sqrt_params = c(0.5, 0.2247),
       norm_params = c(0, 0.2),
       prob_thresh = 0.2,
       ins_prob = 0.11,
       del_prob = 0.04,
       sub_prob = 0.01,
       min_read_length = 50,
       lognorm_read_length = c(0.200110276521, -10075.4363813,
                               17922.611306),
       custom_read_lengths = NULL,
       prob_dup = 0.0,
       haplotype_probs = NULL,
       sep_files = FALSE,
       compress = FALSE,
       comp_method = "bgzip",
       n_threads = 1L,
       read_pool_size = 100L,
       show_progress = FALSE,
       overwrite = FALSE)

Arguments

Value

Nothing is returned.

ID lines

The ID lines for FASTQ files are formatted as such:

⁠@<genome name>-<chromosome name>-<starting position>-<strand>⁠

where ⁠genome name⁠ is always REF for reference genomes (as opposed to haplotypes).

References

Stöcker, B. K., J. Köster, and S. Rahmann. 2016. SimLoRD: simulation of long read data. Bioinformatics 32:2704–2706.

Examples

rg <- create_genome(10, 100e3, 100)
dir <- tempdir(TRUE)
pacbio(rg, paste0(dir, "/pacbio_reads"), n_reads = 100)

Description

Accepts uncompressed and gzipped fasta files.

Usage

read_fasta(fasta_files, fai_files = NULL, cut_names = FALSE)

Arguments

Value

A ref_genome object.

Description

Interactive wrapper for a pointer to a C++ object that stores reference genome information.

Details

This class should NEVER be created using ref_genome$new. Only use read_fasta or create_genome. Because this class wraps a pointer to a C++ object, there are no fields to manipulate directly. All manipulations are done through this class's methods.

Methods

Public methods

• ref_genome$new()

• ref_genome$print()

• ref_genome$ptr()

• ref_genome$n_chroms()

• ref_genome$sizes()

• ref_genome$chrom_names()

• ref_genome$chrom()

• ref_genome$gc_prop()

• ref_genome$nt_prop()

• ref_genome$set_names()

• ref_genome$clean_names()

• ref_genome$add_chroms()

• ref_genome$rm_chroms()

• ref_genome$merge_chroms()

• ref_genome$filter_chroms()

• ref_genome$replace_Ns()

Method new()

Do NOT use this; only use read_fasta or create_genome to make a new ref_genome.

Usage

ref_genome$new(genome_ptr)

Arguments

Method print()

Print a ref_genome object.

Usage

ref_genome$print()

Method ptr()

View pointer to underlying C++ object (this is not useful to end users).

Usage

ref_genome$ptr()

Returns

An externalptr object.

Method n_chroms()

View number of chromosomes.

Usage

ref_genome$n_chroms()

Returns

Integer number of chromosomes.

Method sizes()

View chromosome sizes.

Usage

ref_genome$sizes()

Returns

Integer vector of chromosome sizes.

Method chrom_names()

View chromosome names.

Usage

ref_genome$chrom_names()

Returns

Character vector of chromosome names.

Method chrom()

View one reference chromosome.

Usage

ref_genome$chrom(chrom_ind)

Arguments

Returns

A single string representing the chosen chromosome's DNA sequence.

Method gc_prop()

View GC proportion for part of one reference chromosome.

Usage

ref_genome$gc_prop(chrom_ind, start, end)

Arguments

Returns

A double in the range ⁠[0,1]⁠ representing the proportion of DNA sequence that is either G or C.

Method nt_prop()

View nucleotide content for part of one reference chromosome

Usage

ref_genome$nt_prop(nt, chrom_ind, start, end)

Arguments

Returns

A double in the range ⁠[0,1]⁠ representing the proportion of DNA sequence that is nt.

Method set_names()

Change chromosome names.

Usage

ref_genome$set_names(new_names)

Arguments

Returns

This R6 object, invisibly.

Examples

ref <- create_genome(4, 10)
ref$set_names(c("a", "b", "c", "d"))

Method clean_names()

Clean chromosome names, converting ⁠" :;=%,\\|/\"\'"⁠ to "_".

Usage

ref_genome$clean_names()

Returns

This R6 object, invisibly.

Examples

ref <- create_genome(4, 10)
ref$set_names(c("a:", "b|", "c;", "d'"))
ref$clean_names()

Method add_chroms()

Add one or more chromosomes.

Usage

ref_genome$add_chroms(new_chroms, new_names = NULL)

Arguments

Returns

This R6 object, invisibly.

Examples

ref <- create_genome(4, 10)
ref$add_chroms("TCAGTCAG")

Method rm_chroms()

Remove one or more chromosomes by name

Usage

ref_genome$rm_chroms(chrom_names)

Arguments

Returns

This R6 object, invisibly.

Examples

ref <- create_genome(4, 10)
ref$set_names(c("a", "b", "c", "d"))
ref$rm_chroms("b")

Method merge_chroms()

Merge chromosomes into one.

Usage

ref_genome$merge_chroms(chrom_names)

Arguments

Returns

This R6 object, invisibly.

Examples

ref <- create_genome(4, 10)
ref$merge_chroms(ref$chrom_names()[1:2])
ref$merge_chroms(NULL)

Method filter_chroms()

Filter chromosomes by size or for a proportion of total bases.

Usage

ref_genome$filter_chroms(threshold, method)

Arguments

Returns

This R6 object, invisibly.

Examples

ref <- create_genome(4, 100, 50)
ref$filter_chroms(90, "size")
ref$filter_chroms(0.4, "prop")

Method replace_Ns()

Replace Ns in the reference genome.

Usage

ref_genome$replace_Ns(pi_tcag, n_threads = 1, show_progress = FALSE)

Arguments

Returns

This R6 object, invisibly.

See Also

read_fasta create_genome

Examples

## ------------------------------------------------
## Method `ref_genome$set_names`
## ------------------------------------------------

ref <- create_genome(4, 10)
ref$set_names(c("a", "b", "c", "d"))


## ------------------------------------------------
## Method `ref_genome$clean_names`
## ------------------------------------------------

ref <- create_genome(4, 10)
ref$set_names(c("a:", "b|", "c;", "d'"))
ref$clean_names()


## ------------------------------------------------
## Method `ref_genome$add_chroms`
## ------------------------------------------------

ref <- create_genome(4, 10)
ref$add_chroms("TCAGTCAG")

## ------------------------------------------------
## Method `ref_genome$rm_chroms`
## ------------------------------------------------

ref <- create_genome(4, 10)
ref$set_names(c("a", "b", "c", "d"))
ref$rm_chroms("b")


## ------------------------------------------------
## Method `ref_genome$merge_chroms`
## ------------------------------------------------

ref <- create_genome(4, 10)
ref$merge_chroms(ref$chrom_names()[1:2])
ref$merge_chroms(NULL)


## ------------------------------------------------
## Method `ref_genome$filter_chroms`
## ------------------------------------------------

ref <- create_genome(4, 100, 50)
ref$filter_chroms(90, "size")
ref$filter_chroms(0.4, "prop")

Description

For a more detailed explanation, see vignette("sub-models").

Usage

sub_JC69(lambda, mu = 1, gamma_shape = NULL, gamma_k = 5, invariant = 0)

sub_K80(alpha, beta, mu = 1, gamma_shape = NULL, gamma_k = 5, invariant = 0)

sub_F81(pi_tcag, mu = 1, gamma_shape = NULL, gamma_k = 5, invariant = 0)

sub_HKY85(
  pi_tcag,
  alpha,
  beta,
  mu = 1,
  gamma_shape = NULL,
  gamma_k = 5,
  invariant = 0
)

sub_F84(
  pi_tcag,
  beta,
  kappa,
  mu = 1,
  gamma_shape = NULL,
  gamma_k = 5,
  invariant = 0
)

sub_TN93(
  pi_tcag,
  alpha_1,
  alpha_2,
  beta,
  mu = 1,
  gamma_shape = NULL,
  gamma_k = 5,
  invariant = 0
)

sub_GTR(
  pi_tcag,
  abcdef,
  mu = 1,
  gamma_shape = NULL,
  gamma_k = 5,
  invariant = 0
)

sub_UNREST(Q, mu = 1, gamma_shape = NULL, gamma_k = 5, invariant = 0)

Arguments

Value

A sub_info object, which is an R6 class that wraps the info needed for the create_haplotypes function. It does not allow the user to directly manipulate the info inside, as that should be done using the sub_models functions. You can use the following methods from the class to view information:

Q()

View a list of substitution rate matrices, one for each Gamma category.

pi_tcag()

View the equilibrium nucleotide frequencies.

gammas()

View the discrete Gamma-class values.

invariant()

View the proportion of invariant sites.

model()

View the substitution model.

U()

View list of the U matrices (one matrix per Gamma category) used for calculating transition-probability matrices. This is empty for UNREST models.

Ui()

View list of the U^-1 matrices (one matrix per Gamma category) used for calculating transition-probability matrices. This is empty for UNREST models.

L()

View list of the lambda vectors (one vector per Gamma category) used for calculating transition-probability matrices. This is empty for UNREST models.

Functions

• sub_JC69(): JC69 model.

• sub_K80(): K80 model.

• sub_F81(): F81 model.

• sub_HKY85(): HKY85 model.

• sub_F84(): F84 model.

• sub_TN93(): TN93 model.

• sub_GTR(): GTR model.

• sub_UNREST(): UNREST model.

See Also

create_haplotypes

Examples

# Same substitution rate for all types:
obj_JC69 <- sub_JC69(lambda = 0.1)

# Transitions 2x more likely than transversions:
obj_K80 <- sub_K80(alpha = 0.2, beta = 0.1)

# Incorporating equilibrium frequencies:
obj_HKY85 <- sub_HKY85(pi_tcag = c(0.1, 0.2, 0.3, 0.4),
                       alpha = 0.2, beta = 0.1)

# 10-category Gamma distribution for among-site variability:
obj_K80 <- sub_K80(alpha = 0.2, beta = 0.1,
                   gamma_shape = 1, gamma_k = 10)

# Invariant sites:
obj_K80 <- sub_K80(alpha = 0.2, beta = 0.1,
                   invariant = 0.25)

Description

This file produces 1 FASTA file for a ref_genome object and one file for each haplotype in a haplotypes object.

Usage

write_fasta(
  obj,
  out_prefix,
  compress = FALSE,
  comp_method = "bgzip",
  text_width = 80,
  show_progress = FALSE,
  n_threads = 1,
  overwrite = FALSE
)

Arguments

Value

NULL

Description

Compression in this function always uses "bgzip" for compatibility with "tabix".

Usage

write_vcf(
  haps,
  out_prefix,
  compress = FALSE,
  sample_matrix = NULL,
  show_progress = FALSE,
  overwrite = FALSE
)

Arguments

Value

NULL