Overlap Graphs: Sequence Assembly Using Shared K-mers

2022/01/28 (金) 【第9回】ゼロから始めるゲノム解析（Python編） 資料

Transcript

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

   View Slide

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

   View Slide

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

   View 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 Slide

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

   View Slide

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

   View Slide

7. 前提知識編
   

   View Slide

8. オーバーラップグラフとは
   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 のオーバーラップグラフ
   

   View Slide

9. ログファイルの作成
   今回のプログラムでは--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']})
   出力ログ
   

   View Slide

10. ランタイムメッセージの制御
    デフォルトでは、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を用いて制御できる
    

    View Slide

11. 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以下の時に出力
    

    View Slide

12. 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()
    

    View Slide

13. 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
    

    View Slide

14. >>> 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'],
    

    View Slide

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

    View Slide

16. 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')]
    

    View Slide

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

    View Slide

18. 解法編
    

    View Slide

19. 解法１
    

    View Slide

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

    View Slide