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

Find a Motif in DNA: Exploring Sequence Similarity

onouyek
January 21, 2022

Find a Motif in DNA: Exploring Sequence Similarity

onouyek

January 21, 2022
Tweet

More Decks by onouyek

Other Decks in Programming

Transcript

  1. 【第8回】ゼロから始めるゲノム解析
    (Python編)
    Find a Motif in DNA: Exploring Sequence Similarity
    @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/subs/
    DNA配列sとtからsの部分配列になりうる tの位置をすべて求めよ

    View full-size slide

  6. 前提知識編

    View full-size slide

  7. 部分配列
    文字列sの中に含まれる文字列tと同一の配列
    s:GATATATGCATATACTT
    t:ATAT
    の例ではsの2,4,10番目にtの開始位置が存在する
    G A T A T A T G C A T A T A C T T
    2 4 10

    View full-size slide

  8. in演算子
    >>> seq = 'GATATATGCATATACTT'
    >>> subseq = 'ATAT'
    seqとsubseqを定義する
    >>> subseq in seq
    1
    seqにsubseqが存在するかはinで確認できるが位置情報は分からない

    View full-size slide

  9. str.find()
    第2引数で開始位置を指定できる
    >>> seq.find(subseq)
    1
    部分配列を調べる場合、str.find()でマッチする最小の位置が求められる
    >>> seq.find(subseq, 2)
    3
    存在しなければ-1が返される
    >>> seq.find(subseq, 10)
    -1

    View full-size slide

  10. str.index()
    第2引数で開始位置を指定できる
    >>> if subseq in seq:
    >>> seq.index(subseq)
    1
    str.find()と同様にマッチする最小の位置が求められる
    >>> if subseq in seq[2:]:
    >>> seq.index(subseq, 2)
    3
    存在しなければValueErrorが返される
    >>> seq.index(subseq, 10)
    ValueError: substring not found

    View full-size slide

  11. reモジュール
    >>> import re
    >>> re.findall(subseq, seq)
    ['ATAT', 'ATAT']
    re.findall()でseqに存在するsubseqを見つけることができる
    >>> list(re.finditer(subseq, seq))
    [,
    ]
    re.finditer()でseqに存在するsubseqの位置も分かる
    1つめのsubseqと重なる開始位置4のsubseqは見つけられない

    View full-size slide

  12. Solution 1: str.find()
    ①subseqが見つからなくなったら終了
    ②lastを開始位置の+1の位置に更新
    ③*でリストを展開



    View full-size slide

  13. Solution 2: str.index()
    >>> ' '.join(found)
    Traceback (most recent call last):
    File "", line 1, in
    ' '.join(found)
    TypeError: sequence item 0: expected str instance, int found
    ①subseqが見つからなくなったら終了
    ②lastを開始位置の+1の位置に更新
    ③str.join()を使う場合はstringに変換



    View full-size slide

  14. Solution 3: 高階関数
    ① lambda n: 0 <= n
    ② lambda n: seq.find(subseq, n)
    ③ lambda n: n + 1



    View full-size slide

  15. Solution 4: k-mers
    seqからすべてのkmerを作成してsubseqにマッチするものだけに絞り込む
    >>> kmers = [seq[i:i + k] for i in range(len(seq) - k + 1)]
    >>> kmers
    ['GATA', 'ATAT', 'TATA', 'ATAT', 'TATG', 'ATGC', 'TGCA', 'GCAT', 'CATA',
    'ATAT', 'TATA', 'ATAC', 'TACT', 'ACTT']

    View full-size slide

  16. Solution 5: re.finditer()
    先読みアサーション(?=()を使うことでオーバーラップする配列の位置も分かる
    >>> list(re.finditer('(?=(ATAT))', 'GATATATGCATATACTT'))
    [,
    ,
    ]

    View full-size slide

  17. ベンチマーキング
    hyperfine -m 1000 -L prg ./solution1_str_find.py,./solution2_str_index.py,\
    ./solution3_functional.py,./solution4_kmers_functional.py,\
    ./solution4_kmers_imperative.py,./solution5_re.py \
    '{prg} GATATATGCATATACTT ATAT' --prepare 'rm -rf __pycache__'

    View full-size slide

  18. 追加課題
    プログラムを拡張して、IUPACコードを使った部分配列を扱えるようにしてください。(例
    えばRはAまたはGを表すので、ARCはAACとAGCにマッチする。)

    View full-size slide

  19. 本日学んだこと
    ● str.find()
    ● str.index()
    ● k-mer
    ● 正規表現モジュール

    View full-size slide