sim.seqs: Simulate DNA sequences according to DNA substitution models

Usage Arguments Value Note Examples

View source: R/sim.seqs.R

Usage

1
sim.seqs(num.seqs, num.haps, length.seqs, nucl.freqs, count.haps, subst.model = c("JC69", "K80", "F81", "HKY85"), mu.rate, transi.rate, transv.rate)

Arguments

num.seqs

Number of simulated DNA sequences

num.haps

Number of simulated unique species' haplotypes

length.seqs

Basepair length of DNA sequences

nucl.freqs

Nucleotide frequency distribution vector of A, C, G, and T respectively

count.haps

Haplotype frequency distribution vector

subst.model

Model of DNA substitution

mu.rate

Overall nucleotide mutation rate/site/generation

transi.rate

Nucleotide transition rate/site/generation

transv.rate

Nucleotide transversion rate/site/generation

Value

A FASTA file of DNA sequences

Note

num.seqs must be greater than or equal to num.haps.

Both num.seqs and num.haps must be greater than 1.

nucl.freqs must have a length of four and its elements must sum to 1.

count.haps must have a length of num.haps and its elements must sum to num.seqs.

subst.model must be one of "JC69" (Jukes Cantor corrected p-distance), "K80" (Kimura-2-Parameter (K2P), "F81" (Felenstein) or "HKY85"(Hasegawa-Kishino-Yano)

mu.rate must be specified for both "JC69" and "F81" models

transi.rate and transv.rate must be specified for both "K80" and "HKY85" models

All elements nucl.freqs must be equal to 0.25 when subst.model is either "JC69" or "K80"

All elements nucl.freqs must differ from 0.25 when subst.model is either "F81" or "HKY85"

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# Simulate DNA sequences from the 5'-COI DNA barcode region under a Jukes Cantor nucleotide substitution model

num.seqs <- 100 # number of DNA sequences
num.haps <- 15 # number of haplotypes
length.seqs <- 658 # length of DNA sequences
count.haps <- c(60, rep(10, 2), rep(5, 2), rep(1, 5)) # haplotype frequency distribution
nucl.freqs <- rep(0.25, 4) # nucleotide frequency distribution
subst.model <- "JC69" # desired nucleotide substitution model
mu.rate <- 1e-3 # mutation rate
transi.rate <- NULL # transition rate
transv.rate <- NULL # transversion rate

sim.seqs(num.seqs = num.seqs, num.haps = num.haps, length.seqs = length.seqs, nucl.freqs = nucl.freqs, count.haps = count.haps, subst.model = subst.model, transi.rate = transi.rate, transv.rate = transv.rate)

jphill01/HACSim.R documentation built on Jan. 7, 2021, 3:04 a.m.