如何扩展模糊的DNA序列

Question

假设你有一个像这样的DNA序列：

AATCRVTAA

其中R和V是DNA核苷酸的不确定值，其中R代表A或G，而V代表A、C或G。

是否有一种Biopython方法可以生成由上述不确定序列表示的所有不同序列组合？

例如，在此处输出将为：

AATCAATAA
AATCACTAA
AATCAGTAA
AATCGATAA
AATCGCTAA
AATCGGTAA

- jrjc

Answer 1

也许有一种更短更快的方法，因为毫无疑问，这个函数将在非常大的数据上使用：

from Bio import Seq
from itertools import product

def extend_ambiguous_dna(seq):
   """return list of all possible sequences given an ambiguous DNA input"""
   d = Seq.IUPAC.IUPACData.ambiguous_dna_values
   return [ list(map("".join, product(*map(d.get, seq)))) ]

使用map可以让您的循环在C中执行，而不是在Python中执行。这比使用普通循环甚至列表推导式要快得多。

现场测试

使用简单的字典d代替由ambiguous_na_values返回的字典

from itertools import product
import time

d = { "N": ["A", "G", "T", "C"], "R": ["C", "A", "T", "G"] }
seq = "RNRN"

# using list comprehensions
lst_start = time.time()
[ "".join(i) for i in product(*[ d[j] for j in seq ]) ]
lst_end = time.time()

# using map
map_start = time.time()
[ list(map("".join, product(*map(d.get, seq)))) ]
map_end = time.time()

lst_delay = (lst_end - lst_start) * 1000
map_delay = (map_end - map_start) * 1000

print("List delay: {} ms".format(round(lst_delay, 2)))
print("Map delay: {} ms".format(round(map_delay, 2)))

输出：

# len(seq) = 2:
List delay: 0.02 ms
Map delay: 0.01 ms

# len(seq) = 3:
List delay: 0.04 ms
Map delay: 0.02 ms

# len(seq) = 4
List delay: 0.08 ms
Map delay: 0.06 ms

# len(seq) = 5
List delay: 0.43 ms
Map delay: 0.17 ms

# len(seq) = 10
List delay: 126.68 ms
Map delay: 77.15 ms

# len(seq) = 12
List delay: 1887.53 ms
Map delay: 1320.49 ms

显然，map 更优秀，但只有 2 或 3 倍的差距。可以肯定地说，它可以进一步优化。

- Jivan

---

由于Biopython的更改，可能与此有关，因此需要更改此代码才能正常工作。将第一行更改为import Bio.Data.IUPACData as bdi，然后将d = Seq.IUPAC.IUPACData.ambiguous_dna_values更改为d = bdi.ambiguous_dna_values。 - Wayne

Answer 2

我最终编写了自己的函数:

from Bio import Seq
from itertools import product

def extend_ambiguous_dna(seq):
   """return list of all possible sequences given an ambiguous DNA input"""
   d = Seq.IUPAC.IUPACData.ambiguous_dna_values
   r = []
   for i in product(*[d[j] for j in seq]):
      r.append("".join(i))
   return r 

In [1]: extend_ambiguous_dna("AV")
Out[1]: ['AA', 'AC', 'AG']

它允许您为给定尺寸生成每种图案。

In [2]: extend_ambiguous_dna("NN")

Out[2]: ['GG', 'GA', 'GT', 'GC',
         'AG', 'AA', 'AT', 'AC',
         'TG', 'TA', 'TT', 'TC',
         'CG', 'CA', 'CT', 'CC']

希望这能为其他人节省时间！

- jrjc

Answer 3

我不确定是否有一种Biopython的方法来实现这个，但是这里有一个使用itertools的方法：

s = "AATCRVTAA"
ambig = {"R": ["A", "G"], "V":["A", "C", "G"]}
groups = itertools.groupby(s, lambda char:char not in ambig)
splits = []
for b,group in groups:
    if b:
        splits.extend([[g] for g in group])
    else:
        for nuc in group:
            splits.append(ambig[nuc])
answer = [''.join(p) for p in itertools.product(*splits)]

输出：

In [189]: answer
Out[189]: ['AATCAATAA', 'AATCACTAA', 'AATCAGTAA', 'AATCGATAA', 'AATCGCTAA', 'AATCGGTAA']

- inspectorG4dget

Answer 4

另一个itertools解决方案：

from itertools import product
import re

lu = {'R':'AG', 'V':'ACG'}

def get_seqs(seq):
    seqs = []
    nrepl = seq.count('R') + seq.count('V')
    sp_seq = [a for a in re.split(r'(R|V)', seq) if a]
    pr_terms = [lu[a] for a in sp_seq if a in 'RV']

    for cmb in product(*pr_terms):
        seqs.append(''.join(sp_seq).replace('R', '%s').replace('V', '%s') % cmb)
    return seqs

seq = 'AATCRVTAA'

print 'seq: ', seq
print '\n'.join(get_seqs(seq))

seq1 = 'RAATCRVTAAR'
print 'seq: ', seq1
print '\n'.join(get_seqs(seq1))

输出：

seq:  AATCRVTAA
AATCAATAA
AATCACTAA
AATCAGTAA
AATCGATAA
AATCGCTAA
AATCGGTAA
seq:  RAATCRVTAAR
AAATCAATAAA
AAATCAATAAG
AAATCACTAAA
AAATCACTAAG
AAATCAGTAAA
AAATCAGTAAG
AAATCGATAAA
AAATCGATAAG
AAATCGCTAAA
AAATCGCTAAG
AAATCGGTAAA
AAATCGGTAAG
GAATCAATAAA
GAATCAATAAG
GAATCACTAAA
GAATCACTAAG
GAATCAGTAAA
GAATCAGTAAG
GAATCGATAAA
GAATCGATAAG
GAATCGCTAAA
GAATCGCTAAG
GAATCGGTAAA
GAATCGGTAAG

- vikramls

---

在特殊情况下，如果存在两个或更多相同的相邻模糊代码（如“RRATCGGTAAA”），则会产生错误的输出。 - Zingo