🔖

バイオインフォマティクスコンテスト2021-Problem3 Diagnosis

2021/08/02に公開

algorithm

bioinformatics

コンテスト参加記

idea

バイオインフォマティクスコンテスト2021参加記の問題編です。
この問題は部分点のみで満点回答までたどり着けていないためideaにしています。

問題概要

URL : https://stepik.org/lesson/541855/

与えられる情報

以下の情報が与えられる。

各ノードが表現型異常を表している木構造（表現型ツリー）
- より低レベルのノードのほうがより具体的な表現型異常になっている
- ノードの特異性の情報はIC(v)で表す
いくつかの病気の情報
- それぞれの病気はいくつかの表現型異常を持っている
患者の情報
- それぞれの患者は同様にいくつかの表現型異常を持っている

課題

患者の最も可能性の高い疾患を特定する。

具体的には、表現型 Qp を持つ各患者 p に対して、以下の式を最大化する表現型 Dm を持つ病気 m を見つける。

\sum_{\begin{subarray}{l}q\in\\Qp\\\end{subarray}}(\max_{\begin{subarray}{l}d\in\\Dm\\\end{subarray}} IC(LCA(q,d)))

LCA(q, d) はノード q, d の Lowest common ancestor を表す

入力

行数	入力内容	説明
1	整数 N	表現型ツリーのノード数を表す
2	N-1 要素の配列(スペース区切りの整数)	配列のインデックスは表現型ツリーの root 以外のノードを表し、要素の値はそのノードの親を表す
3	N 要素の配列(スペース区切りの整数)	配列のインデックスは表現型ツリーのノードを表し、要素の値は表現型ツリーの特異性の値を表す( IC(v) )
4	整数 M	病気の数を表す
続くM行	いくつかの要素の配列スペース区切りの整数	各行が病気を表し、1つ目の値がその病気に含まれる表現型異常の数を表す。続く値はその病気に含まれる表現型異常の内容（表現型ツリーのノードのインデックス）を表す
続く行	整数 nq	患者の数を表す
続く nq 行	いくつかの要素の配列スペース区切りの整数	各行が患者を表し、1つ目の値がその患者の表現型異常の数を表す。続く値は患者の表現型異常の内容（表現型ツリーのノードのインデックス）を表す

出力

nq 行を出力する
各行は1つの整数を出力する
各行は患者を表し、出力内容は患者の最も可能性が高い疾患の番号を表す

サンプル

インプット

10                            // ノードの数
1 1 3 3 4 4 5 5 5             // 2つ目以降のノードの親子関係
5 7 8 13 18 14 15 21 20 29    // 各ノードのICの値
2                             // 病気の数
2 4 2                         // 病気1 の表現型情報
1 10                          // 病気2 の表現型情報
4                             // 患者の数
3 5 9 8                       // 患者1 の表現型情報
1 6                           // 患者1 の表現型情報
2 7 10                        // 患者3 の表現型情報
1 10                          // 患者4 の表現型情報

アウトプット

サンプル解説

10個の表現型異常があり、それぞれの親子関係は下のようになる（カッコ内は特異性ICの値）

         1(5)
      ／      ＼
   2(7)         3(8)
              ／     ＼
         4(13)         5(18)
        ／  ＼        ／   |   ＼
   6(14)    7(15)  8(21)  9(20)  10(29)

病気の表現型はそれぞれ
病気1 : 4, 2
病気2 : 10
なので下のようになる

         1(5)
      ／       ＼
    2①           3(8)
              ／      ＼
          4①            5(18)
        ／  ＼        ／   |   ＼
   6(14)    7(15)  8(21)  9(20)  10②

患者①は 5,9,8 の表現型異常を持っているため、ツリーに次のようにマップできます。

         1(5)
      ／      ＼
    2(7)         3(8)
              ／     ＼
         4(13)          5①
        ／  ＼        ／  |  ＼
   6(14)    7(15)  8①    9①    10(29)

患者1とそれぞれの病気の値は以下のようになります。
病気①：
MAX(IC(LCA(5, 2)), IC(LCA(5, 4)) + MAX(IC(LCA(8, 2)), IC(LCA(8, 4))) + MAX(IC(LCA(9, 2)), IC(LCA(9, 4)))
= 24

病気②：
MAX(IC(LCA(5, 10))) + MAX(IC(LCA(8, 10))) + MAX(IC(LCA(9, 10)))
= 54

-> 値が大きい病気②が患者①の病気と予測できる。

これらをそれぞれの患者と病気で実施していくと

患者2:
病気① :
MAX(IC(LCA(6, 2))) + MAX(IC(LCA(6, 4)))
= 18
病気② :
MAX(IC(LCA(6, 10)))
= 8

-> 病気①と予測

患者3:
病気① :
MAX(IC(LCA(7, 2)), IC(LCA(7, 4))) + MAX(IC(LCA(10, 2)), IC(LCA(10, 4)))
= 21
病気② :
MAX(IC(LCA(7, 10))) + MAX(IC(LCA(10, 10)))
= 37

-> 病気②と予測

患者4:
病気① :
MAX(IC(LCA(10, 2))) + MAX(IC(LCA(10, 4)))
= 13
病気② ： MAX(IC(LCA(10, 10)))
= 29

-> 病気②と予測

考えたこと

1 : 愚直にシミュレーション

まずは簡単な問題を回答できるようにするため、愚直にシミュレーションを行いました。
具体的には、以下のようなイメージで、全患者と全病気で、それぞれの表現型異常を突き合せてLCAを計算し値を計算していくというものです。
また、LCAの計算についてもqの親を列挙し、その後dの親に含まれるかを確認していく全探索で実装していきました。

def findParent(child):
    parent = set()
    # child で指定したノードの親すべて探し出して返す。
    ...
    return parent

def findLCA(node,target):
    # node の親をたどっていき、target に存在するかを確認していく
    return lca

def LCA(q,d):
    # q の親を列挙する
    qParent = findParent(q)
    # d の親を探していき、q の親を探し当てるまで続ける
    return findLCA(d,qParent)

def IC(node):
    return IC[node]

def solve(patient,disease):
    for q in patient:
        for d in disease:
            IC(LCA(q,d))

for patient in patients:
    score = 0
    maxIdx = 0
    for i,disease in enumerate(diseases):
        dScore = solve(patient,disease)
        if score < dScore:
            score = dScore
            maxIdx = i
    print(maxIdx)

しかしこの方法では最悪で O(N * M * qn^2) (患者数×病気数×木の深さの2乗)の計算量がかかってしまいます。
この問題の小課題1でもN = 19830, M = 5000, qn = 5000 となっているため、おそらくまともに動作しません。
そこで次の手を考えていきました。

2 : 一度計算したものは保存しておく

まず考えたのはLCAの値をキャッシュしておく方法です。
また LCA(a,b) = LCA(b,a) が成り立つため、それぞれを計算せずに一度計算したらどちらも残しておくことで無駄な計算を省いていくことで高速化していくというアプローチをとりました。
具体的には、function LCA(q,d) 内でfindParent(q) を呼び出す前に、まずはキャッシュ内に(q,d) または(d,q) が存在することを確認し、存在する場合はそれが持つ値を返す、キャッシュに存在しない場合は探索後にキャッシュに新たに値を保存する単純な方法で実装しました。
このアプローチで（おそらく）高速化できるであろうと思っていたのですが、データ量が多いからなのか実装で無駄があったからなのか、PCのメモリ使用量がいっぱいになりディスクIOもパツパツになってしまっている状態で、結局のところこれではあまりいい結果が得られませんでした。

memo = set()
...
def LCA(q,d):
    # q の親を列挙する
    if (q,d) in memo or (q,d) in memo:
        return memo[(q,d)]
    qParent = findParent(q)
    # d の親を探していき、q の親を探し当てるまで続ける
    lca = findLCA(d,qParent)
    memo[(q,d)] = lca
    memo[(d,q)] = lca
    return lca

3 : LCA を高速化する（未実装）

また愚直にLCAを計算していくと木の構造が下記のようになっている場合に非常に低速になってしまいます。

               1
             ／  ＼
           2        3
         ／           ＼
        4               5
     ／                   ＼
    6                       7
   ...

そこでLCAをどのように求めれば高速にできるかを検討しました。
しかしこれはキャッシュした際のバグを対処していて結局手が付けられずでした。

感想

完答できなかったものの問題の内容も面白く、色々試行錯誤しながら楽しめました。
この問題を解いているときに、例えばプログラミング言語のみでどうこうするのではなく、データベースなど別の何かをりようしたらうまく行かないだろうかと考えていました。特にこの問題は木構造の探索が重要になるのでグラフデータベースなどに与えられたデータを取り込んであげればもしかしたら計算も早く終わり、可視化してみたりもできるのでは?と考えていました。