kuroの覚え書き

96の個人的覚え書き

塩基配列(文字列)のどこにリストに挙げた配列が位置しているかを示すプログラム

モチーフの抽出ができたらそれが一体どこにあるのか列記したい。

import re

text = """GACTTTAGATGGCTTCTTCCTTTATAACCAATTGATATTGCATACTCTGATGAGATTTAT 
AATTAAAGAAGCAGAAACAAAAACAAGTAAAACAGAAACAATCAACACAGAGAAACCACC 
CCGAGAATATCTCCATTGGTTGGTGACTGATATCCCTGCTACAACTGGAACAACCTTTGG"""

list1 = text.split()

motif = """AAACA AAACAA AACAA AGAAA AGAAAC GAAAC"""

list2= motif.split()

for i, w in enumelate(list1):
    for k in list2:
        m = re.finditer(k,w)
        for match in m:
            print(i, match)

実行すると

$ python3 motif1.py
1 <_sre.SRE_Match object; span=(14, 19), match='AAACA'>
1 <_sre.SRE_Match object; span=(20, 25), match='AAACA'>
1 <_sre.SRE_Match object; span=(29, 34), match='AAACA'>
1 <_sre.SRE_Match object; span=(35, 40), match='AAACA'>
1 <_sre.SRE_Match object; span=(14, 20), match='AAACAA'>
1 <_sre.SRE_Match object; span=(20, 26), match='AAACAA'>
1 <_sre.SRE_Match object; span=(35, 41), match='AAACAA'>
1 <_sre.SRE_Match object; span=(15, 20), match='AACAA'>
1 <_sre.SRE_Match object; span=(21, 26), match='AACAA'>
1 <_sre.SRE_Match object; span=(36, 41), match='AACAA'>
1 <_sre.SRE_Match object; span=(12, 17), match='AGAAA'>
1 <_sre.SRE_Match object; span=(33, 38), match='AGAAA'>
1 <_sre.SRE_Match object; span=(50, 55), match='AGAAA'>
1 <_sre.SRE_Match object; span=(12, 18), match='AGAAAC'>
1 <_sre.SRE_Match object; span=(33, 39), match='AGAAAC'>
1 <_sre.SRE_Match object; span=(50, 56), match='AGAAAC'>
1 <_sre.SRE_Match object; span=(13, 18), match='GAAAC'>
1 <_sre.SRE_Match object; span=(34, 39), match='GAAAC'>
1 <_sre.SRE_Match object; span=(51, 56), match='GAAAC'>
2 <_sre.SRE_Match object; span=(48, 53), match='AACAA'>

こんな感じ。
span=()のところに存在位置が表示されている
このあとの処理がやりやすいようにもうちょっと整形しておきたいが。

import re

text = """GACTTTAGATGGCTTCTTCCTTTATAACCAATTGATATTGCATACTCTGATGAGATTTAT
AATTAAAGAAGCAGAAACAAAAACAAGTAAAACAGAAACAATCAACACAGAGAAACCACC
CCGAGAATATCTCCATTGGTTGGTGACTGATATCCCTGCTACAACTGGAACAACCTTTGG"""

list1 = text.split()

motif = """AAACA AAACAA AACAA AGAAA AGAAAC GAAAC"""

list2= motif.split()

for i, w in enumerate(list1):
    for k in list2:
        m = re.finditer(k,w)
        for match in m:
            print(i, match.group(), match.start())

これだ。

$ python3 motif2.py
1 AAACA 14
1 AAACA 20
1 AAACA 29
1 AAACA 35
1 AAACAA 14
1 AAACAA 20
1 AAACAA 35
1 AACAA 15
1 AACAA 21
1 AACAA 36
1 AGAAA 12
1 AGAAA 33
1 AGAAA 50
1 AGAAAC 12
1 AGAAAC 33
1 AGAAAC 50
1 GAAAC 13
1 GAAAC 34
1 GAAAC 51
2 AACAA 48

なお、match.span()とすると(start, end)のタプルが返る。end()だとendだけ。

検索したい配列とモチーフをそれぞれファイルで用意しておいて(各1行に1配列)スクリプトに食わせる

import sys
import argparse
import re

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set seq file", type=str)
        parser.add_argument("motif", help="please set motif list", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'motif'):
        motif = args.motif

    with open(file) as list1:
        for i, w in enumerate(list1):
            with open(motif) as list2:
                for k in list2:
                    m = re.finditer(k.strip(),w)
                    for match in m:
                        print(i, match.group(), match.start())

if __name__ == '__main__':
    main()

第一引数に配列ファイル、第二引数にモチーフファイルをつけて実行。出力はそのままだと標準出力。


2/11追記
reのfinditerでは重複したモチーフは飛ばされてしまうことが判明したので修正

import sys
import argparse
import regex

def get_args():
    parser = argparse.ArgumentParser()

    if sys.stdin.isatty():
        parser.add_argument("file", help="please set seq file", type=str)
        parser.add_argument("motif", help="please set motif list", type=str)

    args = parser.parse_args()

    return(args)

def main():
    args = get_args()

    if hasattr(args, 'file'):
        file = args.file
    if hasattr(args, 'motif'):
        motif = args.motif

    with open(file) as list1:
        for i, w in enumerate(list1):
            with open(motif) as list2:
                for k in list2:
                    m = regex.finditer(k.strip(),w,overlapped=True)
                    for match in m:
                        print(i, match.group(), match.start(), match.end())

if __name__ == '__main__':
    main()

これでOKのはず。