Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Finding the Longest Shared Subsequence: Finding K-mers, Writing Functions, and Using Binary Search

onouyek
February 04, 2022

Finding the Longest Shared Subsequence: Finding K-mers, Writing Functions, and Using Binary Search

onouyek

February 04, 2022
Tweet

More Decks by onouyek

Other Decks in Programming

Transcript

  1. 【第10回】ゼロから始めるゲノム解析
    (Python編)
    Finding the Longest Shared Subsequence: Finding K-mers, Writing Functions, and
    Using Binary Search
    @onouyek

    View full-size slide

  2. 本勉強会の概要・目的
    書籍名
    対象者/目的
    Mastering Python for Bioinformatics
    Python・バイオインフォマティクス知識ほぼ
    ゼロの人を対象に、正しい Pythonのコー
    ディング手法について学ぶ
    頻度 毎週〜隔週開催予定
    登壇者 募集中!

    View full-size slide

  3. Rosalindとは
    ● 問題解決を通じてバイオインフォマティク
    ス、プログラミング、およびアルゴリズムを
    学習するためのプラットフォーム
    ● 大学やハッカソン、就職の面接にも
    600回
    以上の採用実績あり
    参考:https://qiita.com/_kimoton/items/d534d0fa9b83dd7dc412
    概要

    View full-size slide

  4. 環境構築 - 必要パッケージ群のインストール
    # 公開されているレポジトリからファイル群を取得
    $ git clone https://github.com/kyclark/biofx_python
    $ cd biofx_python
    # requirements.txtに記載のパッケージをインストール
    $ pip3 install -r requirements.txt
    # pylintの設定ファイルをホームディレクトリに移動
    $ cp pylintrc ~/.pylintrc
    # mypyの設定ファイルをホームディレクトリに移動
    $ cp mypy.ini ~/.mypy.ini

    View full-size slide

  5. 本日のお題 https://rosalind.info/problems/lcsm/
    k個のDNA配列(FASTA形式)から最も長い共通部分文字列を求めよ

    View full-size slide

  6. 前提知識編

    View full-size slide

  7. 最長共通部分文字列
    $ cat tests/inputs/1.fa
    >Rosalind_1
    GATTACA
    >Rosalind_2
    TAGACCA
    >Rosalind_3
    ATACA
    以下の例ではCA、TA、ACが最長共通部分文字列
    >Rosalind_1
    GATTACA
    >Rosalind_2
    TAGACCA
    >Rosalind_3
    ATACA
    >Rosalind_1
    GATTACA
    >Rosalind_2
    TAGACCA
    >Rosalind_3
    ATACA
    >Rosalind_1
    GATTACA
    >Rosalind_2
    TAGACCA
    >Rosalind_3
    ATACA

    View full-size slide

  8. FASTAファイルの読み込み
    >>> from Bio import SeqIO
    >>> fh = open('./tests/inputs/1.fa')
    >>> recs = SeqIO.parse(fh, 'fasta')
    >>> type(recs)

    Bio.SeqIO.parse()でイテレータを作る
    >>> fh = open('./tests/inputs/1.fa')
    >>> seqs = [str(rec.seq) for rec in SeqIO.parse(fh, 'fasta')]
    >>> seqs
    ['GATTACA', 'TAGACCA', 'ATACA']
    リスト内包表記で配列を取り出す
    >>> seqs = list(map(lambda rec: str(rec.seq), SeqIO.parse(fh, 'fasta')))
    map()を使って取り出す

    View full-size slide

  9. K-mersの取得
    第9回で作ったfind_kmers()
    def find_kmers(seq: str, k: int) -> List[str]:
    """ Find k-mers in string """
    n = len(seq) - k + 1
    return [] if n < 1 else [seq[i:i + k] for i in range(n)]
    def test_find_kmers() -> None:
    """ Test find_kmers """
    assert find_kmers('', 1) == []
    assert find_kmers('ACTG', 1) == ['A', 'C', 'T', 'G']
    assert find_kmers('ACTG', 2) == ['AC', 'CT', 'TG']
    assert find_kmers('ACTG', 3) == ['ACT', 'CTG']
    assert find_kmers('ACTG', 4) == ['ACTG']
    assert find_kmers('ACTG', 5) == []

    View full-size slide

  10. range()によるカウントダウン
    range()の第3引数でstepを-1にすればカウントダウンできる
    >>> shortest = min(map(len, seqs))
    >>> shortest
    5
    最短の配列の長さを求める
    >>> list(range(shortest, 0, -1))
    [5, 4, 3, 2, 1]
    カウントアップしてからreversed()しても同じ結果は得られる
    >>> list(reversed(range(1, shortest + 1)))
    [5, 4, 3, 2, 1]

    View full-size slide

  11. collections.Counter()によるK-mersのカウント
    >>> from collections import Counter
    >>> counts = Counter()
    第1回で扱ったCounter()でK-mersをカウントできる
    >>> kmers = [set(find_kmers(seq, shortest)) for seq in seqs]
    >>> for group in kmers:
    … counts.update(group)

    >>> pprint(counts)
    Counter({'ATTAC': 1,
    'GATTA': 1,
    'TTACA': 1,
    'AGACC': 1,
    'TAGAC': 1,
    'GACCA': 1,
    'ATACA': 1})
    Counter.update()でcountsを更新できる

    View full-size slide

  12. itertools.chain()によるリスト内リストの結合
    >>> from itertools import chain
    >>> list(chain.from_iterable(kmers))
    ['ATTAC', 'GATTA', 'TTACA', 'AGACC', 'TAGAC', 'GACCA', 'ATACA']
    chain.from_iterable()でkmersを1つのリストにできる
    >>> counts = Counter(chain.from_iterable(kmers))
    >>> pprint(counts)
    Counter({'ATTAC': 1,
    'GATTA': 1,
    'TTACA': 1,
    'AGACC': 1,
    'TAGAC': 1,
    'GACCA': 1,
    'ATACA': 1})
    Counter()と組み合わせて1行でcountsが作れる

    View full-size slide

  13. Solution 1: 線形探索
    ①入力配列の中で最も短い配列の長さを求める
    ②最短の配列長からkをカウントダウンしていく
    ③common_kmers()で共通配列があったら終了する



    View full-size slide

  14. (参考)二分探索
    第4回で扱った再帰を使った方法
    ①highがlowより小さくなったときを終了条件にする
    ②中間地点の値と探している値が同じ時にmidを返す
    ③中間地点の値が探している値より大きければhighをmid-1
    にしてbinary_search()する
    ④中間地点の値が探している値より小さければlowをmid+1に
    してbinary_search()する




    View full-size slide

  15. Solution 2: 二分探索
    今回の問題用に修正した二分探索
    ①引数で受け取ったf()とlowとhighから共通k-mersのリスト
    hiとloを作る
    ②hiとloがどちらも共通k-mersが存在したらhighの位置を
    返す
    ③共通k-mersがloにあってhiになければhighをmidにして
    binary_search()する
    ④hiもloも空の場合は-1を返す





    View full-size slide

  16. Solution 2: 二分探索
    ①lowを1にhighを最短の配列長で二分探索して開始位置を決
    める
    ②開始位置から共通k-mersが見つからなくなるまでカウントアッ
    プして候補k-mersをcandidatesに追加していく
    ③max()のkeyオプションで最大長のk-merを出力する



    View full-size slide

  17. ベンチマーキング
    hyperfine -L prg ./solution1_kmers_functional.py,./solution2_binary_search.py \
    '{prg} tests/inputs/2.fa'

    View full-size slide

  18. 本日学んだこと
    ● K-mersによる共通配列の検索
    ● itertools.chain()によるリスト内リストの結合
    ● 二分探索による検索の効率化
    ● min()、max()のkeyオプション

    View full-size slide