benchmark/bm_so_fasta.rb

   1 # The Computer Language Shootout
   2 # http://shootout.alioth.debian.org/
   3 # Contributed by Sokolov Yura
   4
   5 $last = 42.0
   6 def gen_random (max,im=139968,ia=3877,ic=29573)
   7     (max * ($last = ($last * ia + ic) % im)) / im
   8 end
   9
  10 alu =
  11    "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG"+
  12    "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA"+
  13    "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT"+
  14    "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA"+
  15    "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG"+
  16    "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC"+
  17    "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
  18
  19 iub = [
  20     ["a", 0.27],
  21     ["c", 0.12],
  22     ["g", 0.12],
  23     ["t", 0.27],
  24
  25     ["B", 0.02],
  26     ["D", 0.02],
  27     ["H", 0.02],
  28     ["K", 0.02],
  29     ["M", 0.02],
  30     ["N", 0.02],
  31     ["R", 0.02],
  32     ["S", 0.02],
  33     ["V", 0.02],
  34     ["W", 0.02],
  35     ["Y", 0.02],
  36 ]
  37 homosapiens = [
  38     ["a", 0.3029549426680],
  39     ["c", 0.1979883004921],
  40     ["g", 0.1975473066391],
  41     ["t", 0.3015094502008],
  42 ]
  43
  44 def make_repeat_fasta(id, desc, src, n)
  45     puts ">#{id} #{desc}"
  46     v = nil
  47     width = 60
  48     l = src.length
  49     s = src * ((n / l) + 1)
  50     s.slice!(n, l)
  51     puts(s.scan(/.{1,#{width}}/).join("\n"))
  52 end
  53
  54 def make_random_fasta(id, desc, table, n)
  55     puts ">#{id} #{desc}"
  56     rand, v = nil,nil
  57     width = 60
  58     chunk = 1 * width
  59     prob = 0.0
  60     table.each{|v| v[1]= (prob += v[1])}
  61     for i in 1..(n/width)
  62         puts((1..width).collect{
  63             rand = gen_random(1.0)
  64             table.find{|v| v[1]>rand}[0]
  65         }.join)
  66     end
  67     if n%width != 0
  68         puts((1..(n%width)).collect{
  69             rand = gen_random(1.0)
  70             table.find{|v| v[1]>rand}[0]
  71         }.join)
  72     end
  73 end
  74
  75
  76 n = (ARGV[0] or 250_000).to_i
  77
  78 make_repeat_fasta('ONE', 'Homo sapiens alu', alu, n*2)
  79 make_random_fasta('TWO', 'IUB ambiguity codes', iub, n*3)
  80 make_random_fasta('THREE', 'Homo sapiens frequency', homosapiens, n*5)
  81