While the discrete-loci, discrete-time spatial structured Wright Fisher model has its own output format, at times it may be beneficial to work with “traditional” genomic data formats. In this tutorial, we will use the model to simulate a infections - with multiple parasites, or with a COI > 1, in ten hosts. We will then will simulate “reads” from a subset of respective host “samples”. Reads will produce allele depths for our variants (i.e. in VCF terminology, the AD and DP fields). For the purpose of this tutorial, we will force all sites to be biallelic.
‘polySimiIBD’ allows the user to interact with individual parasites. As a result, this requires some thought when it comes to combining data within-hosts. Below, we simulate data (as before) for a population of parasites among ten hosts. We then create the ARG for those parasites that are within Host 1 and Host 2. Finally, we collapse the parasite haplotypes from host 1 and 2 into “sequencing read counts”, respectively.
set.seed(1)
# define some parameters parameters
pos <- seq(1, 1e3, 1e2)
N <- 10
m <- 0.5
rho <- 1e-3
mean_coi <- 2
# run forwards simulation
swfsim <- polySimIBD::sim_swf(pos = pos,
N = N,
m = m,
rho = rho,
mean_coi = mean_coi,
migr_mat = 1, # nonspatial
tlim = 10)
#> completed in 0.000099 seconds
# get ARG
ARG <- polySimIBD::get_arg(swfsim)
#> completed in 0.000054 seconds
# layer on mutations
parasite.hapmat <- polySimIBD::layer_mutations_on_ARG(arg = ARG, mutationrate = 0.1)We now have created a matrix of haplotypes that has dimensions Loci x
Parasites. We know which parasites belong to which hosts (i.e. which
columns) by our host1parasites and
host2parasites vectors. Using this information, we are
going to simulate “reads” for the parasites within our host and collapse
them to form a (non-referent) Within Sample Allele Frequency Matrix.
WSAF.list <- polySimIBD::sim_biallelic(COIs = swfsim$coi[1:2], # how many haplotypes in each sample
haplotypematrix = parasite.hapmat[, 1:sum(swfsim$coi[1:2])],
coverage = 100,
alpha = 1,
overdispersion = 0.1,
epsilon = 0.05)The return form the sim_billalelic function is a list of
phased and unphased items. Phased items include the strain proportions,
the biallelic haplotype matrix, the non-referent within sample allele
counts, and the per-haplotype coverage. We then collapse all of this
information to form convoluted, or unphased, non-referent within sample
allele frequencies for each of the hosts that were specified above.
These non-referent within sample allele frequencies can now be used for
downstream analyses or converted to a more standard file format.
Below, we show how this data could be used to create a variant call
file (VCF; for more details, see VCF
specifciations) with the vcfR package. It is important
to note, that we continue in a biallelic SNP framework after using the
polySimIBD::sim_biallelic function.
# VCF files (and vcfR Objects) contain three pieces:
# 1. Meta (Header information)
# 2. Fixed (Loci and Filtering Information)
# 3. Genotype Matrix
# Here we will assume that we are using the VCFv4.3 specifications
# Also assuming user called the VCF with a ploidy of 2 specification (allowed for homozyg-ref, heterozyg, homozyg-alt calls)
# get meta
meta <- paste("##fileformat=VCFv4.3", "##Simulated with polySimIBD", "##ploidy=2", collapse = "\n")
# get Fix: must be in this order and only these
fix <- data.frame(CHROM = "CHROM1",
POS = swfsim$pos,
ID = NA,
REF = "A", # arbitrary choice
ALT = "T", # arbitrary choice
QUAL = NA,
FILTER = "PASS",
INFO = NA)
# get GT
# here we will need to make a choice on the "threshold" for "calling" a homozyg-ref, heterozyg, or homozyg-alt allele
# selecting 10% as 2x the simulated error rate
gtmatsim <- matrix("0/1", nrow = nrow(WSAF.list$NRWSAcounts), ncol = ncol(WSAF.list$NRWSAcounts))
gtmatsim[WSAF.list$NRWSAF > 0.9] <- "1/1"
gtmatsim[WSAF.list$NRWSAF < 0.1] <- "0/0"
# now that we have genotype calls, we can combine these into depths
gt <- matrix(NA, nrow = nrow(gtmatsim), ncol = ncol(gtmatsim))
# quick for loop to protect against vector setting
for (i in 1:nrow(gtmatsim)) {
for (j in 1:ncol(gtmatsim)) {
AD <- paste(WSAF.list$WS.coverage[i,j] - WSAF.list$NRWSAcounts[i,j], WSAF.list$NRWSAcounts[i,j], sep = ",")
gt[i, j] <- paste(gtmatsim[i,j], AD, WSAF.list$WS.coverage[i,j], sep = ":")
}
}
gt <- cbind(FORMAT = "GT:AD:DP", gtmatsim)
# write out new vcfRobj
require(vcfR) # need this for class
newvcfR <- new("vcfR", meta = meta, fix = as.matrix(fix), gt = gt)Reminder, we leave it to the user to decide how to wrangle and combine this data per their project specifications!
Sepearately, polySimIBD provides the functionality to
convert a from the into a Newick tree