$30 off During Our Annual Pro Sale. View Details »

Overlap Graphs: Sequence Assembly Using Shared K-mers

nkimoto
January 29, 2022

Overlap Graphs: Sequence Assembly Using Shared K-mers

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

nkimoto

January 29, 2022
Tweet

More Decks by nkimoto

Other Decks in Programming

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. 解法1

    View Slide

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

    View Slide