Increase the maximum SegSites number of runMacs

[jamonterotena]

Hello, I encounter this error when trying to create a 100 founder population with 34k SNPs. The maximum seg sites produced with 'NULL' is 18k, which is too low for my aim. Is there any way to skip this error and force any number of loci to be simulated?

Error in runMacs(nInd = numberFounders, nChr = 1, segSites = numberSNPs,  : 
  MaCS did not return enough segSites, use segSites=NULL to return all sites generated by MaCS

[gregorgorjanc]

@jamonterotena how many segregating sites you get depends on the used demographic/species model. Which one are you using now? You can tweak that choice or use runMacs2() instead.

@gaynorr this is more of a discussion point, so consider moving it to disucssions.

[jamonterotena]

Thanks for the fast reply. I'm using species="GENERIC" now

I tried runMacs2(), but it still depends on the chromosome base pair size. I wanted to simulate 34k SNPs in a 39Mbp chromosome. With runMacs2(), I would need to simulate a 200Mbp chromosome in order to get that number of SNPs 😆

The reason why I need such high density is because I would like to simulate WGS data. I suppose this is based on SNP chips. Is there any way to overturn that?

[gaynorr]

@jamonterotena the number of base pairs is highly dependent on effective population size. To simplify, we often model agricultural species with a small current effective population size and a large historic effective population size. This combination means that there will be lots of loci (due to large historic population size) with strong LD (due to small current population size). You can certainly manipulate these parameters to get the number of sites you want, but there's no telling if the LD is appropriate without carefully selecting the parameters.

Here is some example code that will give you more than enough sites. However, I caution that the LD pattern might not be desirable.

library(AlphaSimR)

# Simulating a single chromosome for 100 individuals
# Assuming current effective population size is 100
# Modeling an extreme bottleneck of 100,000 Ne 100 generations ago
founderPop = runMacs2(nInd=100, 
                      nChr=1, 
                      Ne=100, 
                      histNe=10000, 
                      histGen=100,
                      bp=39e6)

founderPop@nLoci # I got 75K