Find a Motif in DNA: Exploring Sequence Similarity

【第8回】ゼロから始めるゲノム解析（Python編）

Transcript

1. 【第8回】ゼロから始めるゲノム解析 （Python編） Find a Motif in DNA: Exploring Sequence Similarity @onouyek
2. 本勉強会の概要・目的 書籍名 対象者/目的 Mastering Python for Bioinformatics Python・バイオインフォマティクス知識ほぼ ゼロの人を対象に、正しい Pythonのコー ディング手法について学ぶ 頻度 毎週〜隔週開催予定 登壇者 募集中！
3. Rosalindとは • 問題解決を通じてバイオインフォマティク ス、プログラミング、およびアルゴリズムを 学習するためのプラットフォーム • 大学やハッカソン、就職の面接にも 600回 以上の採用実績あり 参考：https://qiita.com/_kimoton/items/d534d0fa9b83dd7dc412 概要
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
5. 本日のお題 https://rosalind.info/problems/subs/ DNA配列sとtからsの部分配列になりうる tの位置をすべて求めよ
6. 前提知識編
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
8. in演算子 >>> seq = 'GATATATGCATATACTT' >>> subseq = 'ATAT' seqとsubseqを定義する >>> subseq in seq 1 seqにsubseqが存在するかはinで確認できるが位置情報は分からない
9. str.find() 第2引数で開始位置を指定できる >>> seq.find(subseq) 1 部分配列を調べる場合、str.find()でマッチする最小の位置が求められる >>> seq.find(subseq, 2) 3 存在しなければ-1が返される >>> seq.find(subseq, 10) -1
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
11. reモジュール >>> import re >>> re.findall(subseq, seq) ['ATAT', 'ATAT'] re.findall()でseqに存在するsubseqを見つけることができる >>> list(re.finditer(subseq, seq)) [<re.Match object; span=(1, 5), match='ATAT'>, <re.Match object; span=(9, 13), match='ATAT'>] re.finditer()でseqに存在するsubseqの位置も分かる 1つめのsubseqと重なる開始位置4のsubseqは見つけられない
12. 解法編
13. Solution 1: str.find() ①subseqが見つからなくなったら終了 ②lastを開始位置の＋１の位置に更新 ③*でリストを展開 ① ② ③
14. Solution 2: str.index() >>> ' '.join(found) Traceback (most recent call last): File "<input>", line 1, in <module> ' '.join(found) TypeError: sequence item 0: expected str instance, int found ①subseqが見つからなくなったら終了 ②lastを開始位置の＋１の位置に更新 ③str.join()を使う場合はstringに変換 ① ② ③
15. Solution 3: 高階関数 ① lambda n: 0 <= n ② lambda n: seq.find(subseq, n) ③ lambda n: n + 1 ① ② ③
16. 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']
17. Solution 5: re.finditer() 先読みアサーション(?=(<pattern>)を使うことでオーバーラップする配列の位置も分かる >>> list(re.finditer('(?=(ATAT))', 'GATATATGCATATACTT')) [<re.Match object; span=(1, 1), match=''>, <re.Match object; span=(3, 3), match=''>, <re.Match object; span=(9, 9), match=''>]
18. ベンチマーキング 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__'
19. 追加課題 プログラムを拡張して、IUPACコードを使った部分配列を扱えるようにしてください。（例 えばRはAまたはGを表すので、ARCはAACとAGCにマッチする。）
20. 本日学んだこと • str.find() • str.index() • k-mer • 正規表現モジュール