<jupyter_start><jupyter_text>1. Parameters<jupyter_code># Defaults

## Random seed
random_seed <- 25524

## Directories
simulation_dir <- "simulations/unset"
reference_file <- "simulations/reference/reference.fa.gz"
initial_tree_file <- "input/salmonella.tre"

## Simulation parameters
sub_lambda <- 1e-2
sub_pi_tcag <- c(0.1, 0.2, 0.3, 0.4)
sub_alpha <- 0.2
sub_beta <- sub_alpha/2
sub_mu <- 1
sub_invariant <- 0.3
ins_rate <- 1e-4
ins_max_length <- 60
ins_a <- 1.6
del_rate <- 1e-4
del_max_length <- 60
del_a <- 1.6

## Read simulation information
read_coverage <- 30
read_length <- 250

## Other
ncores <- 48
# Parameters
read_coverage = 20
mincov = 10
simulation_dir = "simulations/cov-20"
iterations = 3
sub_alpha = 0.2

output_dir <- file.path(simulation_dir, "simulated_data")
output_vcf_prefix <- file.path(output_dir, "haplotypes")
reads_data_initial_prefix <- file.path(output_dir, "reads_initial", "data")
set.seed(random_seed)
print(output_dir)
print(output_vcf_prefix)<jupyter_output>[1] "simulations/cov-20/simulated_data"<jupyter_text>2. Generate simulated dataThis simulates *Salmonella* data using a reference genome and a tree.<jupyter_code>library(jackalope)

# Make sure we've complied with openmp
jackalope:::using_openmp()
reference <- read_fasta(reference_file)
reference_len <- sum(reference$sizes())
reference
library(ape)

tree <- read.tree(initial_tree_file)
tree <- root(tree, "reference", resolve.root=TRUE)
tree
sub <- sub_HKY85(pi_tcag = sub_pi_tcag, mu = sub_mu,
                 alpha = sub_alpha, beta = sub_beta, gamma_shape=1, gamma_k = 5,
                 invariant = sub_invariant)
ins <- indels(rate = ins_rate, max_length = ins_max_length, a = ins_a)
del <- indels(rate = del_rate, max_length = del_max_length, a = del_a)

ref_haplotypes <- create_haplotypes(reference, haps_phylo(tree), sub=sub, ins=ins, del=del)
ref_haplotypes<jupyter_output><empty_output><jupyter_text>3. Write simulated data<jupyter_code>write_vcf(ref_haplotypes, out_prefix=output_vcf_prefix, compress=TRUE)
assemblies_prefix = file.path(output_dir, "assemblies", "data")

write_fasta(ref_haplotypes, out_prefix=assemblies_prefix,
            compress=TRUE, n_threads=ncores, overwrite=TRUE)
n_samples <- length(tree$tip)
n_reads <- round((reference_len * read_coverage * n_samples) / read_length)
print(sprintf("Number of reads for coverage %sX and read length %s over %s samples with respect to reference with length %s: %s", 
              read_coverage, read_length, n_samples, reference_len, n_reads))

illumina(ref_haplotypes, out_prefix = reads_data_initial_prefix, sep_files=TRUE, n_reads = n_reads,
         frag_mean = read_length * 2 + 50, frag_sd = 100,
         compress=TRUE, comp_method="bgzip", n_threads=ncores,
         paired=TRUE, read_length = read_length)
# Remove the simulated reads for the reference genome since I don't want these in the tree
ref1 <- paste(toString(reads_data_initial_prefix), "_reference_R1.fq.gz", sep="")
ref2 <- paste(toString(reads_data_initial_prefix), "_reference_R2.fq.gz", sep="")
if (file.exists(ref1)) {
    file.remove(ref1)
    print(sprintf("Removing: %s", ref1))
}
if (file.exists(ref2)) {
    file.remove(ref2)
    print(sprintf("Removing: %s", ref2))
}
# Remove the new reference assembly genome since I don't need it
ref1 <- paste(toString(assemblies_prefix), "__reference.fa.gz", sep="")
if (file.exists(ref1)) {
    file.remove(ref1)
    print(sprintf("Removing: %s", ref1))
}<jupyter_output>[1] "Removing: simulations/cov-20/simulated_data/assemblies/data__reference.fa.gz"
{'fraction_non_alphanumeric': 0.12535211267605634, 'fraction_numerical': 0.023098591549295774, 'mean_word_length': 5.318505338078292, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '263395', 'license': 'Apache-2.0', 'n_tokens_mistral': 1376, 'n_tokens_neox': 1258, 'n_words': 271, 'path': 'evaluations/simulation/1-simulate-reads.simulation-cov-20.ipynb', 'repo_name': 'apetkau/genomics-data-index-evaluation', 'size': '3550'}