2

私は次のようなコードを持っています:

import HTSeq
reference = open('genome.fa','r')
sequences = dict( (s.name, s) for s in HTSeq.FastaReader(reference))
out = open('homopolymers_in_ref','w')

def find_all(a_str,sub):
    start = 0
    while True:
        start = a_str.find(sub, start)
        if start == -1: return
        yield start
        start += len(sub)
homa = 'AAAAAAAAAA'
homc = 'CCCCCCCCCC'
homg = 'GGGGGGGGGG'
homt = 'TTTTTTTTTT'
for key,line in sequences.items():
    seq = str(line)
    a= list(find_all(seq,homa))
    c = list(find_all(seq,homc))
    g = list(find_all(seq,homg))
    t = list(find_all(seq,homt))
    for i in a:
##        print i,key,'A'
        out.write(str(i)+'\t'+str(key)+'\t'+'A'+'\n')
    for i in c:
        out.write(str(i)+'\t'+str(key)+'\t'+'C'+'\n')
##        print i,key,'C'
    for i in g:
        out.write(str(i)+'\t'+str(key)+'\t'+'G'+'\n')
    for i in t:
        out.write(str(i)+'\t'+str(key)+'\t'+'T'+'\n')
out.close()

HTSeq を使用してリファレンスを開きました。機能 - 長さ 10 の単純なホモポリマーを探し、開始位置、染色体、タイプ (A、C、T、G) を出力します。

THE シーケンスは常に次のようになります: ACCGCTACGATCGATCGAAAAAAAAAAAAAAAAAACGATCGAC N が含まれることがあります。

したがって、探しているホモポリマーは次のとおりです: AAAAAAAAAA (または C、G、T のみで構成されるその他)

基本的に、あなたからのヘルプは find_all 関数のみに関するものです。ここで変更したいのは、各ホモポリマーの長さを見つけることです。現在、ホモポリマーの長さが 15 である場合、私のスクリプトではそれを判断できません。つまり、現在のように少なくとも 10 bp を見つけ、次の塩基がホモポリマーのようにならなくなるまで +1 を加えて長さを計算します。

Pythonで正規表現を使用する方法について何か提案はありますか?

4

1 に答える 1

6

正規表現でこれを行いたい場合は、次のようなものを試すことができます。

>>> import re
>>> seq = 'ACCGCTACGATCGATCGAAAAAAAAAAAAAAAAAACGATCGAC'
>>>
>>> [(m.group(), m.start())
...     for m in re.finditer(r'([ACGT])\1{9,}', seq)
...         if len(m.group()) >= 10]
[('AAAAAAAAAAAAAAAAAA', 17)]

(sequence, start_index)これにより、タプルのリストが生成されます。

于 2013-09-09T16:05:44.503 に答える