【第9回】ゼロから始めるゲノム解析 （Python編） Overlap Graphs: Sequence Assembly Using Shared K-mers @kimoton 

本勉強会の概要・目的 書籍名 対象者/目的 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 

本日のお題 FASTA形式で与えられたDNA配列セットからO₃のオーバーラップグラフを求めよ https://rosalind.info/problems/prot/ 

● オーバーラップグラフの作成にk-merを使用する方法 ● ファイルベースのランタイムログの取り方 ● collections.defaultdictの使い方 ● 積集合の求め方 ● itertools.productを使ったリストのデカルト積導出 ● iteration_utilities.starfilter() メソッドの使い方 ● グラフ構造を可視化するためのGraphvizの使い方 本日学ぶこと 

前提知識編 

オーバーラップグラフとは https://ja.wikipedia.org/wiki/%E3%82%B3%E3%83%89%E3%83%B3 >Rosalind_0498 AAATAAA >Rosalind_2391 AAATTTT >Rosalind_2323 TTTTCCC >Rosalind_0442 AAATCCC >Rosalind_5013 GGGTGGG Rosalind_0498 AAATAAA AAATTTT Rosalind_2391 Rosalind_2391 AAATTTT TTTTCCC Rosalind_2323 Rosalind_0498 AAATAAA AAATTTT Rosalind_2391 k=3 のオーバーラップグラフ 

ログファイルの作成 今回のプログラムでは--debug引数を付けることでランタイムメッセージが出力される仕様としている。 このような工夫によりプログラムの異常時にデバッグが容易となる。 $ ./grph.py tests/inputs/1.fa --debug Rosalind_2391 Rosalind_2323 Rosalind_0498 Rosalind_2391 Rosalind_0498 Rosalind_0442 実行時に--debug引数を指定 $ cat .log DEBUG:root:STARTS defaultdict(, {'AAA': ['Rosalind_0498', 'Rosalind_2391', 'Rosalind_0442'], 'GGG': ['Rosalind_5013'], 'TTT': ['Rosalind_2323']}) DEBUG:root:ENDS defaultdict(, {'AAA': ['Rosalind_0498'], 'CCC': ['Rosalind_2323', 'Rosalind_0442'], 'GGG': ['Rosalind_5013'], 'TTT': ['Rosalind_2391']}) 出力ログ 

ランタイムメッセージの制御 デフォルトでは、print()関数は標準出力に引数で与えた文字列を出力するが、その出力先は file引数で変更すること ができる。 $ cat log.py #!/usr/bin/env python3 import sys print('This is STDOUT.') print('This is also STDOUT.', file=sys.stdout) print('This is STDERR.', file=sys.stderr) $ ./log.py 1>out 2>err bashでは標準出力はファイルハンドル1、標準エラー出力は2を用いて制御できる 

loggingモジュールの使用 標準出力、標準エラー出力の2種類以上の出力のコントロールを行いたい場合、 loggingモジュールを使用する import logging def main() -> None: args = get_args() logging.basicConfig( filename='.log', filemode='w', level=logging.DEBUG if args.debug else logging.CRITICAL) logging.debug('input file = "%s"', args.file.name) 1 loggingモジュールの挙動にグローバルに影響を与える設定 1 2 2 ファイル名とファイルモードを指定 3 debug引数が付いた場合はDEBUGレベルまで出力、そうでない場合はCRITICALレベルまで出力 3 4 4 logging.debugによる出力はログレベルの設定がdebug以下の時に出力 

loggingモジュールの使用 - ログレベルについて loggingモジュールでは、設定を変えることでどこまでのログレベルを出力するかを制御することができる Level Numeric value 出力コマンド CRITICAL 50 logging.critical() ERROR 40 logging.error() WARNING 30 logging.warning() INFO 20 logging.info() DEBUG 10 logging.debug() NOTSET 0 logging.notset() 

K-merの取得 第7回で扱ったk-merを取得する関数を再利用する 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)] 1 seq中のk文字の長さの文字列数を算出 2 nが負の場合は空リストを返し、そうでなければk-merを返す 1 2 

>>> from pprint import pprint >>> from Bio import SeqIO >>> from collections import defaultdict >>> from grph import find_kmers >>> k = 3 >>> start, end = defaultdict(list), defaultdict(list) >>> for rec in SeqIO.parse('tests/inputs/1.fa', 'fasta'): ... if kmers := find_kmers(str(rec.seq), k): ... start[kmers[0]].append(rec.id) ... end[kmers[-1]].append(rec.id) ... >>> pprint(start) {'AAA': ['Rosalind_0498', 'Rosalind_2391', 'Rosalind_0442'], 'GGG': ['Rosalind_5013'], 'TTT': ['Rosalind_2323']} >>> pprint(end) {'AAA': ['Rosalind_0498'], 'CCC': ['Rosalind_2323', 'Rosalind_0442'], 'GGG': ['Rosalind_5013'], 

intersection関数による積集合の抽出 intersection関数を使うとset同士の積集合を抽出することができる >>> set(start) {'TTT', 'GGG', 'AAA'} >>> set(end) {'TTT', 'CCC', 'AAA', 'GGG'} >>> set(start).intersection(set(end)) {'TTT', 'GGG', 'AAA'} 

product関数によるデカルト積の取得 product関数を使うとリスト内のデカルト積（総当たり組合わせ）を作成することができる >>> from itertools import product >>> kmer = 'AAA' >>> pairs = list(product(end[kmer], start[kmer])) >>> pprint(pairs) [('Rosalind_0498', 'Rosalind_0498'), ('Rosalind_0498', 'Rosalind_2391'), ('Rosalind_0498', 'Rosalind_0442')] 

starfilterモジュール starfilterを使うとfilter関数に複数の引数を取ることができる >>> from iteration_utilities import starfilter >>> list(starfilter(lambda a, b: a != b, pairs)) [('Rosalind_0498', 'Rosalind_2391'), ('Rosalind_0498', 'Rosalind_0442')] 

解法編 

解法１ 

解法２グラフを描画して探索