Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

fasta

  • Loading branch information...
commit 65fa625fac3965ab2e9f24f74c94678b89071875 1 parent ee1d821
Robert Bradshaw robertwb authored

Showing 2 changed files with 158 additions and 0 deletions. Show diff stats Hide diff stats

  1. +79 0 fasta.cython.pyx
  2. +79 0 fasta.python.py
79 fasta.cython.pyx
... ... @@ -0,0 +1,79 @@
  1 +# The Computer Language Benchmarks Game
  2 +# http://shootout.alioth.debian.org/
  3 +#
  4 +# modified by Ian Osgood
  5 +# modified again by Heinrich Acker
  6 +
  7 +import sys, bisect
  8 +
  9 +alu = (
  10 + 'GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG'
  11 + 'GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA'
  12 + 'CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT'
  13 + 'ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA'
  14 + 'GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG'
  15 + 'AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC'
  16 + 'AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA')
  17 +
  18 +iub = zip('acgtBDHKMNRSVWY', [0.27, 0.12, 0.12, 0.27] + [0.02]*11)
  19 +
  20 +homosapiens = [
  21 + ('a', 0.3029549426680),
  22 + ('c', 0.1979883004921),
  23 + ('g', 0.1975473066391),
  24 + ('t', 0.3015094502008),
  25 +]
  26 +
  27 +
  28 +def genRandom(lim, ia = 3877, ic = 29573, im = 139968):
  29 + seed = 42
  30 + imf = float(im)
  31 + while 1:
  32 + seed = (seed * ia + ic) % im
  33 + yield lim * seed / imf
  34 +
  35 +Random = genRandom(1.)
  36 +
  37 +def makeCumulative(table):
  38 + P = []
  39 + C = []
  40 + prob = 0.
  41 + for char, p in table:
  42 + prob += p
  43 + P += [prob]
  44 + C += [char]
  45 + return (P, C)
  46 +
  47 +def repeatFasta(src, n):
  48 + width = 60
  49 + r = len(src)
  50 + s = src + src + src[:n % r]
  51 + for j in xrange(n // width):
  52 + i = j*width % r
  53 + print s[i:i+width]
  54 + if n % width:
  55 + print s[-(n % width):]
  56 +
  57 +def randomFasta(table, n):
  58 + width = 60
  59 + r = xrange(width)
  60 + gR = Random.next
  61 + bb = bisect.bisect
  62 + jn = ''.join
  63 + probs, chars = makeCumulative(table)
  64 + for j in xrange(n // width):
  65 + print jn([chars[bb(probs, gR())] for i in r])
  66 + if n % width:
  67 + print jn([chars[bb(probs, gR())] for i in xrange(n % width)])
  68 +
  69 +
  70 +n = int(sys.argv[1])
  71 +
  72 +print '>ONE Homo sapiens alu'
  73 +repeatFasta(alu, n*2)
  74 +
  75 +print '>TWO IUB ambiguity codes'
  76 +randomFasta(iub, n*3)
  77 +
  78 +print '>THREE Homo sapiens frequency'
  79 +randomFasta(homosapiens, n*5)
79 fasta.python.py
... ... @@ -0,0 +1,79 @@
  1 +# The Computer Language Benchmarks Game
  2 +# http://shootout.alioth.debian.org/
  3 +#
  4 +# modified by Ian Osgood
  5 +# modified again by Heinrich Acker
  6 +
  7 +import sys, bisect
  8 +
  9 +alu = (
  10 + 'GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG'
  11 + 'GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA'
  12 + 'CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT'
  13 + 'ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA'
  14 + 'GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG'
  15 + 'AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC'
  16 + 'AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA')
  17 +
  18 +iub = zip('acgtBDHKMNRSVWY', [0.27, 0.12, 0.12, 0.27] + [0.02]*11)
  19 +
  20 +homosapiens = [
  21 + ('a', 0.3029549426680),
  22 + ('c', 0.1979883004921),
  23 + ('g', 0.1975473066391),
  24 + ('t', 0.3015094502008),
  25 +]
  26 +
  27 +
  28 +def genRandom(lim, ia = 3877, ic = 29573, im = 139968):
  29 + seed = 42
  30 + imf = float(im)
  31 + while 1:
  32 + seed = (seed * ia + ic) % im
  33 + yield lim * seed / imf
  34 +
  35 +Random = genRandom(1.)
  36 +
  37 +def makeCumulative(table):
  38 + P = []
  39 + C = []
  40 + prob = 0.
  41 + for char, p in table:
  42 + prob += p
  43 + P += [prob]
  44 + C += [char]
  45 + return (P, C)
  46 +
  47 +def repeatFasta(src, n):
  48 + width = 60
  49 + r = len(src)
  50 + s = src + src + src[:n % r]
  51 + for j in xrange(n // width):
  52 + i = j*width % r
  53 + print s[i:i+width]
  54 + if n % width:
  55 + print s[-(n % width):]
  56 +
  57 +def randomFasta(table, n):
  58 + width = 60
  59 + r = xrange(width)
  60 + gR = Random.next
  61 + bb = bisect.bisect
  62 + jn = ''.join
  63 + probs, chars = makeCumulative(table)
  64 + for j in xrange(n // width):
  65 + print jn([chars[bb(probs, gR())] for i in r])
  66 + if n % width:
  67 + print jn([chars[bb(probs, gR())] for i in xrange(n % width)])
  68 +
  69 +
  70 +n = int(sys.argv[1])
  71 +
  72 +print '>ONE Homo sapiens alu'
  73 +repeatFasta(alu, n*2)
  74 +
  75 +print '>TWO IUB ambiguity codes'
  76 +randomFasta(iub, n*3)
  77 +
  78 +print '>THREE Homo sapiens frequency'
  79 +randomFasta(homosapiens, n*5)

0 comments on commit 65fa625

Please sign in to comment.
Something went wrong with that request. Please try again.