In [27]:

    

import os, os.path
import random
os.chdir('/home/will/DeepSequencing/testfix/')


    

In [28]:

    

with open('RefSequence_single.fasta') as handle:
    _ = handle.next()
    seq = handle.next().strip()[:600]


    

In [29]:

    

haplotypes = []
for _ in range(100):
    tmp = range(len(seq))
    random.shuffle(tmp)
    nseq = list(seq)
    for n in tmp[:50]:
        nseq[n] = 'A'
    haplotypes.append(''.join(nseq))


    

In [30]:

    

tfracs = [0.4, 0.2, 0.1, 0.05]
nleft = len(haplotypes) - len(tfracs)
fracs = tfracs + [(1-sum(tfracs))/(nleft)]*nleft
cfracs = [sum(fracs[:n]) for n in range(1,len(fracs))]
sum(fracs)
print cfracs


    

    




[0.40000000000000002, 0.60000000000000009, 0.70000000000000007, 0.75000000000000011, 0.75260416666666674, 0.75520833333333337, 0.7578125, 0.76041666666666663, 0.76302083333333326, 0.76562499999999989, 0.76822916666666652, 0.77083333333333315, 0.77343749999999978, 0.77604166666666641, 0.77864583333333304, 0.78124999999999967, 0.7838541666666663, 0.78645833333333293, 0.78906249999999956, 0.79166666666666619, 0.79427083333333282, 0.79687499999999944, 0.79947916666666607, 0.8020833333333327, 0.80468749999999933, 0.80729166666666596, 0.80989583333333259, 0.81249999999999922, 0.81510416666666585, 0.81770833333333248, 0.82031249999999911, 0.82291666666666574, 0.82552083333333237, 0.828124999999999, 0.83072916666666563, 0.83333333333333226, 0.83593749999999889, 0.83854166666666552, 0.84114583333333215, 0.84374999999999878, 0.84635416666666541, 0.84895833333333204, 0.85156249999999867, 0.8541666666666653, 0.85677083333333193, 0.85937499999999856, 0.86197916666666519, 0.86458333333333182, 0.86718749999999845, 0.86979166666666508, 0.87239583333333171, 0.87499999999999833, 0.87760416666666496, 0.88020833333333159, 0.88281249999999822, 0.88541666666666485, 0.88802083333333148, 0.89062499999999811, 0.89322916666666474, 0.89583333333333137, 0.898437499999998, 0.90104166666666463, 0.90364583333333126, 0.90624999999999789, 0.90885416666666452, 0.91145833333333115, 0.91406249999999778, 0.91666666666666441, 0.91927083333333104, 0.92187499999999767, 0.9244791666666643, 0.92708333333333093, 0.92968749999999756, 0.93229166666666419, 0.93489583333333082, 0.93749999999999745, 0.94010416666666408, 0.94270833333333071, 0.94531249999999734, 0.94791666666666397, 0.95052083333333059, 0.95312499999999722, 0.95572916666666385, 0.95833333333333048, 0.96093749999999711, 0.96354166666666374, 0.96614583333333037, 0.968749999999997, 0.97135416666666363, 0.97395833333333026, 0.97656249999999689, 0.97916666666666352, 0.98177083333333015, 0.98437499999999678, 0.98697916666666341, 0.98958333333333004, 0.99218749999999667, 0.9947916666666633, 0.99739583333332993]

In [31]:

    

seq_len = 80
nreads = 500000
pos_sites = range(600-80)
with open('sim_reads.fa', 'w') as handle:
    for _ in range(nreads):
        val = random.uniform(0,1)
        for hap, cval in enumerate(cfracs):
            if val <= cval:
                break
        pos = random.randint(0, 600-80)
        read = haplotypes[hap][pos:pos+seq_len]
        rname = 'hap-%i-start-%i' % (hap, pos)
        handle.write('>%s\n%s\n' % (rname, read))
    
with open('known_haps.fasta', 'w') as handle:
    for hap, (seq, freq) in enumerate(zip(haplotypes, fracs)):
        handle.write('>hap-%i-%f\n%s\n' % (hap, freq, seq))


    

In [ ]: