Constant number of mutation instead of rate

[shaharMazia]

Hi all,
I am trying to simulate a simple simulation. I want to simulate a constant number (say N) of SNPs, and to generate a tree sequence for each SNP separately (I assume the SNPs are very far away on the sequence, so there is no linkage at all).

My attempt is to use sim_ancestry with sequence_length=N, recombination_rate=0.5
This is running relatively slow, because of the recombination, I don't think there is a way to speed it up, but if there are suggestions I will be glad to hear.

Later I want to simulate a single SNP (a binary mutation) for each point on the sequence, with union distribution on the corresponding tree "branches". So I do something like this:

ts = msprime.sim_ancestry(....., sequence_length=N, recombination_rate=0.5)
mts = msprime.sim_mutations(ts, model=msprime.BinaryMutationModel())

But, I am not sure what rate should I specify. The doc says "The rate of mutation per unit of sequence length per unit time, as either a single number (for a uniform rate) or as a RateMap"
I have in mind I might need to compute the sum on all branches (in time unit) , and rate should be sequence_length / total time
So I will have a single mutation in each tree? Does it make sense?

Btw, what if I have more than 1 mutation in the same nucleotide? Can it create reverse mutations? I don't want my simulation to have more than 1 mutation, or reverse mutation in the same nucleotide.

Wonderful lib guys, well done.
Thanks a lot,
Shahar

[jeromekelleher]

Hey @shaharMazia 👋

If you're assuming unlinked loci (i.e with free recombination) then you're almost certainly better off just simulating those loci separately without recombination. This will be much faster.

Simulating a fixed number of mutations is very subtle, and we don't recommend doing it usually. See #1806 and #68 for previous discussions. Why do you need exactly one mutation per site?

If you want to be certain of no recurrent mutations, you can specify "infinite sites" using the discrete_genome=False option -- see here for discussion.

[shaharMazia]

Hi, thanks for the fast reply and the good references!
I will make my homework and will get back in here in a day or two to mark as an answer, or to ask a bit more to clarify it all :)

[shaharMazia]

Ok I understood most of it. Using discrete_genome=False seems like what I need to avoid reverse mutations.
Regards simulating loci separately, I do it using num_replicates in sim_ancestry, and that should not be a problem.
Then I can simulate mutations on top of any of them. I can let down the assumption of a single mutation in any unlinked area, and with enough replicates that will be fine.
What eventually I couldn't understand, is how to combine all the sequences, so I can write a single VCF file. My current code (random seeds are for easy debugging, should I remove in order to have variance in the mutations?):

        ts_iterator = msprime.sim_ancestry(num_replicates=100, random_seed=1)
        mts = []
        for ts in ts_iterator:
            mts.append(msprime.sim_mutations(ts, model=msprime.BinaryMutationModel(), random_seed=1, discrete_genome=False))
        

Not sure how to continue. Previously I would use mts.write_vcf()
Thanks a lot!
Shahar

[jeromekelleher]

Great, glad it's working out for you @shaharMazia!

Merging the outputs into a single VCF is a bit tricky, and probably the simplest thing to do would be to combine the VCFs afterwards, adding some value to all the output positions so you can think of them as being from different chunks along the genome. I don't have any specific suggestions on how to do that, though.

[shaharMazia]

Ok, specifically for my use I can analyse the genotype matrices, and I don't need an actual VCF.
Just in case someone will use it somewhen in the future, I made a numpy array of all my mutations by doing:

        mts = np.empty(0)
        for ts in ts_iterator:
            mt = msprime.sim_mutations(ts, .......)
            mt_matrix = np.array([e.genotypes for e in mt.variants()])
            mts = np.concatenate((mts, mt_matrix), axis=0) if mts.size else mt_matrix

Thanks again!