Intro to jackalope

Based on a reference

The examples here produce FASTQ files from the known reference assembly that could test strategies for how to assemble a similar genome using HTS data.

The first strategy is to use only PacBio sequencing. The PacBio Sequel system produces up to 500,000 reads per Single Molecule, Real-Time (SMRT) cell, so you could run the following for two cells (with the file pacbio_R1.fq as output):

pacbio(ref, out_prefix = "pacbio", n_reads = 2 * 500e3)

An alternative, hybrid strategy uses 1 SMRT cell of PacBio sequencing and 1 lane (\(\sim 500\) million reads) of \(2 \times 100\)bp Illumina sequencing on the HiSeq 2500 system (the default Illumina system in jackalope):

pacbio(ref, out_prefix = "pacbio", n_reads = 500e3)
illumina(ref, out_prefix = "illumina", n_reads = 500e6, paired = TRUE,
         read_length = 100)

The last strategy combines 1 lane of \(2 \times 100\)bp Illumina HiSeq 2500 sequencing with 1 flow cell of \(2 \times 250\)bp mate-pair sequencing on an Illumina MiSeq v3. The mate-pair sequencing uses longer fragments (defaults are mean of 400 and standard deviation of 100) to better cover highly repetitive regions.

illumina(ref, out_prefix = "ill_pe", n_reads = 500e6, paired = TRUE,
         read_length = 100)
illumina(ref, out_prefix = "ill_mp", seq_sys = "MSv3",
         read_length = 250, n_reads = 50e6, matepair = TRUE, 
         frag_mean = 3000, frag_sd = 500)

These data could then be used to compare genome assembly performance between the strategies above, or between programs within a given strategy. Extensions of these tests include adjusting sequencing depth, sequencing platform, and error rates.

Based on a diploid individual

For diploid species, scientists won’t be sampling a haploid reference, and heterozygosity can be a real problem for the assembly process. So below I’ll show you how to simulate a diploid individual and create reads based on that. I’ll just show the hybrid assembly strategy, but the others simply differ in the sequencing step as shown above.

First, we’ll simulate two haplotypes based on the \(\theta\) (population-scaled mutation rate) and the other molecular evolution information for D. melanogaster we got from the evo_rates object earlier. I’m also renaming the haplotypes.

haps <- create_haplotypes(ref, haps_theta(theta = theta, n_haps = 2), 
                          sub, ins, del)
haps$set_names(c("A", "B"))

This results in the following haplotypes object:

#>                               << haplotypes object >>
#> # Haplotypes: 2
#> # Mutations: 21
#> 
#>                           << Reference genome info: >>
#> < Set of 4 chromosomes >
#> # Total size: 961 bp
#>   name                             chromosome                             length
#> 2          TCTCACGGAAAACGGGAAGTTTTTGGCC...CAAGGCTTCGACGTAGAGCCGGTGGATT       232
#> 3          GCTCCGGCAGCATCGGCTTCCGTGATCA...GGACTGCGGGGGTTCTAACGGGCGTCGG       245
#> 4          GAATGCTCGCATGCGAACACTGAACGTT...GAATGCCAGTAAGCTCGAGCGTATAAGG       247
#> X          CCCCATTTCAGCAACGAAGTCGATCGCA...GTTCAATCACGCGAAAAAAAACACGTCA       237

We generate 1 SMRT cell of PacBio sequencing and 1 lane of \(2 \times 100\)bp Illumina sequencing on the HiSeq 2500 system as before, except using the new haps object instead of ref:

pacbio(haps, out_prefix = "pacbio", n_reads = 500e3)
illumina(haps, out_prefix = "illumina", n_reads = 500e6, paired = TRUE,
         read_length = 100)

If we want to save the information for these haplotypes, we can output FASTA files (one per haplotype) or a VCF file:

write_fasta(haps, "haps")
write_vcf(haps, "haps", sample_matrix = cbind(1, 2))

This would result in the following files being generated: haps__A.fa, haps__B.fa, and haps.vcf. The sample_matrix argument to write_vcf states that the first two (i.e., all) haplotypes are from the first (and only) sample. These data could test different assembly strategies, similar to the above section, except that, for diploid organisms, it includes heterozygosity. Increasing the \(\theta\) parameter will result in more heterozygosity, which is an obvious extension of these tests.

From one phylogenetic tree

This section shows how jackalope can generate haplotypes from a phylogeny, then simulate sequencing data from those haplotypes to test phylogeny reconstruction methods. First, I simulated a random species tree of 10 species, then scaled it to have a maximum tree depth of \(4 N_0\) generations:

tree <- rcoal(10)
tree$edge.length <- 4 * N0 * tree$edge.length / max(node.depth.edgelength(tree))

Function haps_phylo organizes and checks the tree object, and including it with the mutation-type information allowed me to create haplotypes based on this phylogeny:

haps <- create_haplotypes(ref, haps_phylo(tree), sub, ins, del)

This results in the following haplotypes object:

#>                               << haplotypes object >>
#> # Haplotypes: 10
#> # Mutations: 547
#> 
#>                           << Reference genome info: >>
#> < Set of 4 chromosomes >
#> # Total size: 961 bp
#>   name                             chromosome                             length
#> 2          TCTCACGGAAAACGGGAAGTTTTTGGCC...CAAGGCTTCGACGTAGAGCCGGTGGATT       232
#> 3          GCTCCGGCAGCATCGGCTTCCGTGATCA...GGACTGCGGGGGTTCTAACGGGCGTCGG       245
#> 4          GAATGCTCGCATGCGAACACTGAACGTT...GAATGCCAGTAAGCTCGAGCGTATAAGG       247
#> X          CCCCATTTCAGCAACGAAGTCGATCGCA...GTTCAATCACGCGAAAAAAAACACGTCA       237

Now I can generate data for 1 flow cell of \(2 \times 250\)bp sequencing on an Illumina MiSeq v3, where haplotype_barcodes is a character string that specifies the barcodes for each haplotype. I also wrote the true phylogenetic tree to a NEWICK file.

haplotype_barcodes <- c("CTAGCTTG", "TCGATCCA", "ATACCAAG", "GCGTTGGA",
                        "CTTCACGG", "TCCTGTAA", "CCTCGGTA", "TTCTAACG", 
                        "CGCTCGTG", "TATCTACA")
illumina(haps, out_prefix = "phylo_tree", seq_sys = "MSv3",
         paired = TRUE, read_length = 250, n_reads = 50e6,
         barcodes = haplotype_barcodes)
ape::write.tree(tree, "true.tree")

The true phylogenetic tree would then be compared to the final tree output from the program(s) the user chooses to test.

From gene trees

Similar to the section above, the ultimate goal here is to test phylogeny reconstruction methods. The difference in this section is that instead of using a single, straightforward phylogeny, I use multiple gene trees per chromosome. In the species used in these simulations, species 1 diverged from 2 and 3 at \(t = 1.0\), where \(t\) indicates time into the past and is in units of \(4 N_0\) generations. Species 2 and 3 diverged at \(t = 0.5\). I assume a recombination rate of \(1 / (4 N_0)\) recombination events per chromosome per generation. There are 4 diploid individuals sampled per species. I used the following code to call scrm to simulate the gene trees and create the object gtrees. And like the section using segregating sites, we have to wrap lapply inside list(trees = ...) to replicate the data structure from scrm for jackalope to use later.

# Run scrm for one chromosome size:
one_chrom <- function(.size) {
    sims <- scrm(
        paste("24 1",
              # Output gene trees:
              "-T",
              # Recombination:
              "-r 1", .size,
              # 3 species with no ongoing migration:
              "-I 3 8 8 8 0",
              # Species 2 derived from 1 at time 1.0:
              "-ej 1.0 2 1",  
              # Species 3 derived from 2 at time 0.5:
              "-ej 0.5 3 2"
        ))
    trees <- sims$trees[[1]]
    # scrm outputs branch lengths in units of 4*N0 generations, but we want just
    # generations:
    adjust_tree <- function(.p) {
        # Read to phylo object and adjust branch lengths:
        .tr <- read.tree(text = .p)
        .tr$edge.length <- .tr$edge.length * 4 * N0
        # "prefix" from `.p` showing how large the region this gene tree refers to is
        prefix <- paste0(strsplit(.p, "\\]")[[1]][1], "]")
        # Put back together into NEWICK text
        return(paste0(prefix, write.tree(.tr)))
    }
    trees <- sapply(trees, adjust_tree)
    names(trees) <- NULL
    return(trees)
}
# For all chromosomes:
gtrees <- list(trees = lapply(ref$sizes(), one_chrom))

We can write the true gene trees using write_gtrees to an ms-style file:

write_gtrees(haps_gtrees(gtrees), "gtrees")

The create_haplotypes function uses these gene trees to create variant haplotypes. As for the other haplotype-creation methods, function haps_gtrees checks and organizes information from the gtrees object.

haps <- create_haplotypes(ref, haps_gtrees(gtrees),
                          sub, ins, del)

This results in the following haplotypes object:

#>                               << haplotypes object >>
#> # Haplotypes: 24
#> # Mutations: 2,782
#> 
#>                           << Reference genome info: >>
#> < Set of 4 chromosomes >
#> # Total size: 961 bp
#>   name                             chromosome                             length
#> 2          TCTCACGGAAAACGGGAAGTTTTTGGCC...CAAGGCTTCGACGTAGAGCCGGTGGATT       232
#> 3          GCTCCGGCAGCATCGGCTTCCGTGATCA...GGACTGCGGGGGTTCTAACGGGCGTCGG       245
#> 4          GAATGCTCGCATGCGAACACTGAACGTT...GAATGCCAGTAAGCTCGAGCGTATAAGG       247
#> X          CCCCATTTCAGCAACGAAGTCGATCGCA...GTTCAATCACGCGAAAAAAAACACGTCA       237

To store mutation information by diploid sample, the write_vcf function writes the haplotypes object to a VCF file. It assigns every other haplotype to a new diploid sample using a matrix for the sample_matrix argument:

write_vcf(haps, out_prefix = "hap_gtrees",
          sample_matrix = matrix(1:haps$n_haps(), ncol = 2, byrow = TRUE))

Next I generated data for 1 flow cell of \(2 \times 250\)bp sequencing on an Illumina MiSeq v3, with barcodes in the object haplotype_barcodes.

# 2 of each barcode bc it's diploid
haplotype_barcodes <- rep(c("TCGCCTTA", "CTAGTACG", "TTCTGCCT", "GCTCAGGA", "AGGAGTCC", 
                            "CATGCCTA", "GTAGAGAG", "CCTCTCTG", "AGCGTAGC", "CAGCCTCG", 
                            "TGCCTCTT", "TCCTCTAC"), each = 2)
illumina(haps, out_prefix = "phylo_gtrees", seq_sys = "MSv3",
         read_length = 250, n_reads = 50e6, paired = TRUE, 
         barcodes = haplotype_barcodes)

Topologies of the gene trees would then be compared to final phylogenies output from software the user is interested in testing. Varying recombination rates or adding gene flow after separation of species would be natural extensions of these simulations.