人工的なたんぱく質データを生成する

はじめに †

開発中の相同性検索ツール（配列が類似しているたんぱく質を検索するツール）の性能を評価するために，シロイヌナズナのたんぱく質データを人工的に生成しました．

たんぱく質の数は30,000，長さはすべて1,000としています．

生起確率のみに基づく生成プログラム †

FASTAフォーマットで保存されているAt_GB_all_prot.gzを読み込んで，各アミノ酸の生起確率に基づいて人工的に生成した配列をFASTAフォーマットで出力し，At_GB_all_prot_art.gzに保存します．

Ruby 1.9用です．

require 'zlib'

h = Hash.new  # ハッシュ
"ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c|
  h[c] = 0  # カウントを0に初期化
end

sum = 0  # 合計
Zlib::GzipReader.open('At_GB_all_prot.gz') do |gz| # 入力ファイル
  gz.each do |l|
    next if (l[0] == '>')  # ヘッダ行はスキップ
    l.chomp!.each_char do |c|  # データ行から一文字ずつ取り出す
      h[c] += 1
      sum += 1
    end
  end
end

h.each do |key,value|
  h[key] /= sum.to_f  # 確率に変換
end

orig_name = "At_GB_all_prot_art"  # 出力ファイルの名前
Zlib::GzipWriter.open(orig_name+".gz", Zlib::BEST_COMPRESSION) do |gz|  # gzip圧縮
  gz.mtime = Time.now
  gz.orig_name = orig_name
  30000.times do |i|  # たんぱく質の数
    gz.printf(">id|%08d\n", i);
    1000.times do |j|  # たんぱく質の長さ
      r = rand
      sum = 0
      h.each do |key, value|
        sum += value
        if r <= sum then
          gz.print key  # アミノ酸を出力
          gz.print "\n" if j % 80 == 79  # 80アミノ酸ごとに改行
          break
        end
      end
    end
    gz.print "\n"
  end
end

開始アミノ酸の生起確率と条件付き生起確率に基づく生成プログラム †

FASTAフォーマットで保存されているAt_GB_all_prog.gzを読み込んで，開始アミノ酸の生起確率と直前のアミノ酸を条件とする条件付き生起確率に基づいて人工的に生成した配列をFASTAフォーマットで出力し，At_GB_all_prot_art.gzに保存します．

これも，Ruby 1.9用です．

require 'zlib'

# 初期化
h = Hash.new  # ハッシュ
"ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c1|
  h[c1] = 0  # カウントを0に初期化
  "ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c2|
    h[c1+c2] = 0  # カウントを0に初期化
  end
end

# カウント
sum = 0  # 合計
first = false;
Zlib::GzipReader.open('At_GB_all_prot.gz') do |gz|
  gz.each do |l|
    c = l[0]
    if (c == '>')  # ヘッダ行はスキップ
      first = true
      next
    end
    if first then  # ヘッダ行の次の行
      h[c] += 1  # 開始アミノ酸のカウント
      sum += 1
      first = false
    end
    for i in 0..(l.chomp!.length-2) do
      token = l[i]+l[i+1]  # バイグラムのトークン
      h[token] += 1
    end
  end
end

# 確率計算
"ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c|
  h[c] /= sum.to_f  # 開始アミノ酸の生起確率に変換
end
"ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c1|
  sum = 0  # 合計
  "ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c2|
    sum += h[c1+c2]  # 合計
  end
  "ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c2|
    h[c1+c2] /= sum.to_f  # 条件付き確率に変換
  end
end

# データ生成
orig_name = "At_GB_all_prot_art"  # 出力ファイルの名前
Zlib::GzipWriter.open(orig_name+".gz", Zlib::BEST_COMPRESSION) do |gz|  # gzip圧縮
  gz.mtime = Time.now
  gz.orig_name = orig_name
  30000.times do |i|  # たんぱく質の数
    gz.printf(">id|%08d\n", i);
    r = rand
    sum = 0
    c1 = "c"
    "ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c|
      sum += h[c]
      if r <= sum then
        gz.print c  # 開始アミノ酸を出力
        c1 = c
        break
      end
    end
    999.times do |j|  # たんぱく質の長さ（開始アミノ酸を除く）
      r = rand
      sum = 0
      "ABCDEFGHIKLMNPQRSTUVWXYZ".each_char do |c2|
        sum += h[c1+c2]
        if r <= sum then
          gz.print c2  # アミノ酸を出力
          c1 = c2
          gz.print "\n" if j % 80 == 78  # 80アミノ酸ごとに改行
          break
        end
      end
    end
    gz.print "\n"
  end
end

参考情報 †

• アミノ酸の生起確率を調べる | とうごろうぃき
• アミノ酸の条件付き生起確率を調べる | とうごろうぃき