「最短ヌクレオチド連鎖問題」が解けない

高校生のとき、友人にこんな問題を出されました。適当に再構成してみます。

タンパク質を構成するアミノ酸は 20種類あるよね。DNA または RNA はそのアミノ酸の設計図で、いわゆる「コドン」を指定している。コドンというのはヌクレオチドの塩基3個から成るもので、ヌクレオチドの塩基は 4種類ある。RNA の場合なら U(ウラシル)、C(シトシン)、A(アデニン)、G(グアニン) だ。これら 4種類の 3個づつの重複順列でコドンは出来ているから、全部で 4×4×4 = 64 種類あることになる。これらコドンによってアミノ酸が指定されるのだよね。

さて、ここからは遊びだ。64種類のコドンをすべて一度づつ使って、それらを適当な順番ですべて連結したとする。このとき、実際にはあり得ないのだが、重複するものは圧縮できるものとする。例えば UUU と UUA というコドンをこの順につなげたとき、UUUUUA の UU の部分が重複しているから、圧縮できて UUUA になる。この次に例えば UAG が来れば、また 2つ圧縮できて UUUAG となる。重複は一度つなげるごとに最大でこのような 2つで、1つの場合もまたまったく重複しない場合もある。

このとき、適当な順番ですべての種類のコドンを連結して圧縮したなら、そのヌクレオチド連鎖の長さの最小値はいくらになるだろう?

 
たぶん、高校の生物の先生が冗談で出した問題だったと思う。そのとき自分はすでに PC(富士通 FM-7!)をもっていたので、機械語でそれを解くプログラムを考えてみたのだが、そのときはうまくいかなかった。学校ではかしこい生徒がかなり短い例を作ったらしいが、その答えは我々にはわかっていない。そしてそれが正解かどうかもわからない。

あれからもう何十年も過ぎたけれど、この問題はいまでも時々思い出す。そしていま素人プログラミングで遊んでいるので、つらつら考えてみるのだが、どうもわからない。これを仮に「最短ヌクレオチド連鎖問題」ということにして、この問題が自分には解けないのだ。


で、じつはそれを解く Ruby プログラムを考えてみたのだが、それはたぶんうまくない。誰かが正解を求めるコードを書いてくれないだろうか、そんな気もあってここに記した。コドンの種類は 64通りということはわかっているから、答えは 66 より小さくはならないことはわかっている。さて、どうであろうか。


なお、時間がどれだけかかってもよいというなら簡単である。用意されたコドンは完全に対称的なので、始まりは何でもよい。なので 64 - 1 = 63 の階乗とおりの順列を考え、それを圧縮して比較し、最短を出すことができる。しかし、63 の階乗というのは

1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000

なのである。これを仮に Ruby で解くとすると、自分の PC では 100万回のループでほぼ 1秒かかるというのが大まかな目安だから、おおよそ

62868097266756724509010233015027211363038232186303546045988064422017052054

年もかかってしまうのだ(細かい数値に意味はないけれど)。さて、それをどうするかということである。


おそらく、いちばん有望な解き方はこうである。まず、与えられた幾らかのコドンを連結圧縮して、それが最短になっているとする。このとき、まだ使われていないコドンをそこにひとつ追加して、そこで最短になるような連結の仕方(アルゴリズム)を与えることができるか。もしそれが可能ならば、再帰的な処理で答えを確実に求めることができる。問題は、そのようなアルゴリズムが存在するか、だ。

追記

いまのところ自分の見つけた最短連鎖の長さは 80 である。その順は例えば

UUU UUC UCU CUU UUA UAU AUU UUG UGU GUU UCC CCU CUC UCA CAU AUC
UCG CGU GUC UAC ACU CUA UAA AAU AUA UAG AGU GUA UGC GCU CUG UGA
GAU AUG UGG GGU GUG CCC CCA CAC ACC CCG CGC GCC CAA AAC ACA CAG
AGC GCA CGA GAC ACG CGG GGC GCG AAA AAG AGA GAA AGG GGA GAG GGG

である。これはたぶん最短でない。考え方としては、できるだけ 2個の圧縮が続くように、でなければ 1個の圧縮、となるように適当に並べただけである。最短の保証はない。

再追記(8/25)

Haskell で brute force するコードを書いてみた。これで求められる筈だけれど、もちろん永遠に終わらない。
shortest_chain_of_nucleotide.hs

import Data.List

main :: IO ()
main = do
    print $ head $ f $ map joinCodons $ permutations codons
        where nucleotide = "UCAG"
              codons = [[a, b, c] | a <- nucleotide, b <- nucleotide, c <- nucleotide]
              f a = [(length x, x) | x <- a, (minimum $ map length a) == length x]

joinCodons :: (Eq a) => [[a]] -> [a]
joinCodons (x:y:ys) = if null ys
                      then head st
                      else joinCodons st
                          where st = [compress x y] ++ ys

compress :: (Eq a) => [a] -> [a] -> [a] 
compress x y = if double x y
               then x ++ drop 2 y
               else if single x y
                    then x ++ drop 1 y
                    else x ++ y

double :: (Eq a) => [a] -> [a] -> Bool
double x y = drop (length x - 2) x == take 2 y

single :: (Eq a) => [a] -> [a] -> Bool
single x y = last x == head y

 

追記

これは本質的に「最短ハミルトン路」問題。

上の長さ 80 の例から単純な乱択アルゴリズムで長さ 72 まで短縮した例。

i = 2677, length = 72
UUU UUC UCU CUU UUA UAU AUU UUG UGU GUC UCC CCU CUC UCA CAU AUC
UCG CGU GUA UAC ACU CUA UAA AAU AUA UAG AGU GUG UGC GCU CUG UGA
GAU AUG UGG GGU GUU CCC CCA CAC ACC CCG CGC GCA CAA AAC ACA CAG
AGC GCG CGA GAC ACG CGG GGC GCC AAA AAG AGA GAA AGG GGA GAG GGG

Ruby コード。

def dist(a, b)
  case
  when a[1, 2] == b[0, 2] then 1
  when a[-1] == b[0] then 2
  else 3
  end
end

L = 64

table = %w(UUU UUC UCU CUU UUA UAU AUU UUG UGU GUU UCC CCU CUC UCA CAU AUC
           UCG CGU GUC UAC ACU CUA UAA AAU AUA UAG AGU GUA UGC GCU CUG UGA
           GAU AUG UGG GGU GUG CCC CCA CAC ACC CCG CGC GCC CAA AAC ACA CAG
           AGC GCA CGA GAC ACG CGG GGC GCG AAA AAG AGA GAA AGG GGA GAG GGG)
edge = table.each_cons(2).map {|a, b| dist(a, b)} + [0]

1.step do |i|
  a, b = rand(1...L), rand(1...L)
  next if a >= b
  prev_d = edge[a - 1] + edge[a] + edge[b - 1] + edge[b]
  table[a], table[b] = table[b], table[a]
  a_d1, a_d2 = dist(table[a - 1], table[a]), dist(table[a], table[a + 1])
  if b == L - 1
    b_d1, b_d2 = dist(table[b - 1], table[b]), 0
  else
    b_d1, b_d2 = dist(table[b - 1], table[b]), dist(table[b], table[b + 1])
  end
  if prev_d > a_d1 + a_d2 + b_d1 + b_d2
    edge[a - 1], edge[a] = a_d1, a_d2
    edge[b - 1], edge[b] = b_d1, b_d2
    puts "i = #{i}, length = #{3 + edge.sum}"
    p table
  else
    table[a], table[b] = table[b], table[a]
  end
end

 
焼きなまし法」でも挑戦してみたけれど、全然ダメ。一回での距離の動きの幅が小さすぎる。
「最短ヌクレオチド連鎖問題」焼きなまし法 · GitHub
(2019/2/14)