2022/01/28 (金) 【第9回】ゼロから始めるゲノム解析(Python編) 資料
【第9回】ゼロから始めるゲノム解析(Python編)Overlap Graphs: Sequence Assembly Using Shared K-mers@kimoton
View Slide
本勉強会の概要・目的書籍名対象者/目的Mastering Python for BioinformaticsPython・バイオインフォ知識ほぼゼロの人を対象に、正しい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_0498AAATAAA>Rosalind_2391AAATTTT>Rosalind_2323TTTTCCC>Rosalind_0442AAATCCC>Rosalind_5013GGGTGGGRosalind_0498AAATAAAAAATTTTRosalind_2391Rosalind_2391AAATTTTTTTTCCCRosalind_2323Rosalind_0498AAATAAAAAATTTTRosalind_2391k=3 のオーバーラップグラフ
ログファイルの作成今回のプログラムでは--debug引数を付けることでランタイムメッセージが出力される仕様としている。このような工夫によりプログラムの異常時にデバッグが容易となる。$ ./grph.py tests/inputs/1.fa --debugRosalind_2391 Rosalind_2323Rosalind_0498 Rosalind_2391Rosalind_0498 Rosalind_0442実行時に--debug引数を指定$ cat .logDEBUG:root:STARTSdefaultdict(,{'AAA': ['Rosalind_0498', 'Rosalind_2391','Rosalind_0442'],'GGG': ['Rosalind_5013'],'TTT': ['Rosalind_2323']})DEBUG:root:ENDSdefaultdict(,{'AAA': ['Rosalind_0498'],'CCC': ['Rosalind_2323', 'Rosalind_0442'],'GGG': ['Rosalind_5013'],'TTT': ['Rosalind_2391']})出力ログ
ランタイムメッセージの制御デフォルトでは、print()関数は標準出力に引数で与えた文字列を出力するが、その出力先はfile引数で変更することができる。$ cat log.py#!/usr/bin/env python3import sysprint('This is STDOUT.')print('This is also STDOUT.', file=sys.stdout)print('This is STDERR.', file=sys.stderr)$ ./log.py 1>out 2>errbashでは標準出力はファイルハンドル1、標準エラー出力は2を用いて制御できる
loggingモジュールの使用標準出力、標準エラー出力の2種類以上の出力のコントロールを行いたい場合、loggingモジュールを使用するimport loggingdef 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モジュールの挙動にグローバルに影響を与える設定122 ファイル名とファイルモードを指定3 debug引数が付いた場合はDEBUGレベルまで出力、そうでない場合はCRITICALレベルまで出力344 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 + 1return [] if n < 1 else [seq[i:i + k] for i in range(n)]1 seq中のk文字の長さの文字列数を算出2 nが負の場合は空リストを返し、そうでなければk-merを返す12
>>> 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')]
解法編
解法1
解法2グラフを描画して探索