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）。

序列总是看起来像： 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 回答