【第8回】ゼロから始めるゲノム解析 （Python編） Find a Motif in DNA: Exploring Sequence Similarity @onouyek

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

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

環境構築 - 必要パッケージ群のインストール # 公開されているレポジトリからファイル群を取得 $ 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

本日のお題 https://rosalind.info/problems/subs/ DNA配列sとtからsの部分配列になりうる tの位置をすべて求めよ

前提知識編

部分配列 文字列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

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

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

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

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は見つけられない

解法編

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

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を開始位置の＋１の位置に更新 ③str.join()を使う場合はstringに変換 ① ② ③

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

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']

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

ベンチマーキング 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__'

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

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