Feature request: SNP data

[Schraiber]

Simulate a fixed number of SNPs by simulating trees and placing one mutation on them. Trees need to be rejection sampled, and a reasonable thing to do is to fix theta<<1 and accept trees with probability 1-exp(-theta/2T_total) ~ theta/2T_total. I think the right acceptance rate in the limit as theta goes down to 0 is something like T_total/E(T_total), but that seems weird because it will frequently be greater than 1?

[jeromekelleher]

Thanks @Schraiber, I'll have a think about this. I think this can be fast, as we can use the (super sneaky!) linear time simulation algorithm from Hudson 1990 to simulate the trees. Hopefully we don't have to reject too much for reasonable values of theta.

[jeromekelleher]

I've had a think about this @Schraiber, and this is the most efficient approach I can think of without coding stuff up specifically:

def snp_example(sample_size, num_snps):
    theta = 1 # TODO calculate correct value of theta s.t the E[#mutations] ~ 10
    replicates = msprime.simulate(
        sample_size=sample_size, Ne=0.5, mutation_rate=theta,
        num_replicates=100 * num_snps)
    j = 0
    variants = np.empty((num_snps, sample_size))
    for ts in replicates:
        for variant in ts.variants():
            variants[j] = variant.genotypes
            j += 1
            break
        if j == num_snps:
            break
    assert j == num_snps
    return variants

Is this correct, or are we biasing things by doing it this way?

It takes about 2 minutes to generate 100K unlinked SNPs for 1000 samples. I've chosen theta such that we seem to generate ~10 mutations, so a rejection is unlikely (generating a few extra mutations is much less costly than throwing away a replicate). Looking at the profile, we're not going to do a huge amount better than this without having a specialised single-locus main-loop. I'm guessing this would probably halve the time, but there might be more in it. The question then is, is it worth it?

[Schraiber]

So I don't think this algorithm produces the right thing. The basic reasoning is that you have to take all the SNPs from a tree, or else you are undercounting the big trees. Unfortunately, taking all the SNPs defeats the idea of getting a bunch of unlinked data.

To demonstrate that taking only the first SNP biases the result, I'll show that it gets the wrong frequency spectrum. Under boring conditions, the probability that a SNP occurs in i out of n chromosomes is proportional to 1/i. So,

In [5]: theory = 1/arange(1,11.)
In [6]: theory /= sum(theory)
In [7]: theory
Out[7]: 
array([ 0.34141715,  0.17070858,  0.11380572,  0.08535429,  0.06828343,
        0.05690286,  0.04877388,  0.04267714,  0.03793524,  0.03414172])

is what we're looking for in a sample of size 10.

Using a very slight modification of the code you wrote:

In [8]: test = snp.snp_example(10, 10000, theta = 1)
In [9]: histogram(sum(test,axis=1),bins=range(1,11),normed=True)
Out[9]: 
(array([ 0.3774,  0.1732,  0.1182,  0.0888,  0.0685,  0.0551,  0.044 ,
         0.0381,  0.0367]), array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10]))

which is obviously not right. In particular, you can see that you have way more low frequency alleles, as expected if you are under representing big trees, on which you could have more of a chance of high frequency mutations.

If we reduce theta to something small, such that you get <1 mutation per tree on average, you get closer:

In [10]: test = snp.snp_example(10, 10000, theta = .01)
In [11]: histogram(sum(test,axis=1),bins=range(1,11),normed=True)
Out[11]: 
(array([ 0.3584,  0.1738,  0.1198,  0.0882,  0.0703,  0.057 ,  0.0485,
         0.0442,  0.0398]), array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10]))

But even this is biased toward singletons.

"Proper" SNP data is something like the limit of this approach as theta goes down to zero, which is obviously going to result in throwing away lots and lots of trees. As argued in the initial post, in the limit you want to simulate a tree, and accept it with probability proportional to T_total on the tree, and then just throw one mutation on it. One suggestion, from Ben Peter, is to just set some arbitrary large-but-not-too-large T_max, and accept trees with probability T_total/T_max. This algorithm would be approximate because T_total has support on (0, infinity) but it is probably good enough for many practical cases.

[jeromekelleher]

Ah OK, that's making sense now. Thanks for taking the time to explain it. This is subtle!

[BenjaminPeter]

I think there are two ways to do this: One is a rejection sampler as mentioned by @Schraiber. However, T_max can be picked adaptively using a simple update scheme, using an algorithm like this (pseudocode, I am unfortunately not familiar with msprime). The update scheme works because if a simulation was accepted with probability T/T_max_old already, we can change T_max_old to T_max_new by rejecting again with probability T_max_old/T_max_new.

from random import random                         

def simulate_snp(n_individuals, n_snps, **kwargs):                    
    trees = []                                                        
    t_max = 0.                                                        

    while len(trees) < n_snps:                                        
        tree = simulate_single_tree(n_individuals, **kwargs)          
        if tree.tree_length > t_max:                                  
            new_t_max = tree.tree_length                              
            trees = [t for t in trees if random() <= t_max/new_t_max] 
            t_max = new_t_max             
            trees.append(tree)
        else:                                                         
            if random() <= tree.tree_length / t_max:                  
                trees.append(tree)                                    

    snps = [tree.get_single_snp() for tree in trees]
    return snps

An alternative that is harder to implement is to create a table of the allele frequeny spectrum from the trees, and then draw from that.

[jeromekelleher]

Thanks for this @BenjaminPeter. I've coded your algorithm up using msprime and numpy here:

def simulate_snp(sample_size, num_snps):
    replicates = msprime.simulate(
        sample_size=sample_size, Ne=0.5, mutation_rate=10,
        num_replicates=100 * num_snps)
    t_max = 0
    variants = np.empty((num_snps, sample_size), dtype="u1")
    total_branch_length = np.empty(num_snps)
    j = 0
    num_adaptive_updates = 0
    num_rejected_trees = 0
    for ts in replicates:
        tree = next(ts.trees())
        tbl = tree.get_total_branch_length()
        if tbl > t_max:
            new_t_max = tbl
            new_variants = np.empty((num_snps, sample_size), dtype="u1")
            new_total_branch_length = np.empty(num_snps)
            keep = np.where(np.random.random(j) < t_max / new_t_max)[0]
            j = keep.shape[0]
            new_variants[:j] = variants[keep]
            new_total_branch_length[:j] = total_branch_length[keep]
            variants = new_variants
            total_branch_length = new_total_branch_length
            t_max = new_t_max
            num_adaptive_updates += 1
        else:
            if np.random.random() < tbl / t_max:
                total_branch_length[j] = tbl
                for variant in ts.variants():
                    variants[j] = variant.genotypes
                    break
                else:
                    raise Exception("Must have at least one mutation")
                j += 1
                if j == num_snps:
                    break
            else:
                num_rejected_trees += 1
    assert j == num_snps
    print("num adaptive updates: ", num_adaptive_updates)
    print("num rejected trees", num_rejected_trees)
    return variants

I haven't checked this against the expectations, so it's very likely wrong somewhere in the details!

It takes about 3 minutes to simulate 100K SNPs for a sample of 100 individuals, which seems slower than it should be (most of the time is being spent in Python, I'm not sure what it's doing though). In this example there were 17 adaptive updates (which are very fast) and 220K rejected trees.

[petrelharp]

This just came up again in the mailing list. I'm not totally clear why people want to do this (please explain?), but, my proposals are:

1. if you want a fixed number of linked snps, set the mutation rate high enough you'll be certain to get at least as many as you need, then randomly subsample these.
2. if you want unlinked snps, run a small-ish simulation and pick a random snp, if there are any, until you have enough.

I think people often want to do this to match their dataset, and so the FAQ should strongly ask people to consider how their lici were ascertained.

[jeromekelleher]

I'm going to close this one as a "wontfix" - we are unlikely to ever actually implement something like this. Happy to reopen if someone has a compelling use-case though.