4

我有一个带~700键的默认字典。密钥格式为A_B_STRING。我需要做的是用'_'分割键,并且比较每个键的'STRING'之间的距离,如果A和B是相同的。如果距离为< = 2,我想将这些键的列表分组到单个默认值键:值组中。将有多个键应该匹配并分组。我也想保留新的组合defaultdict,所有没有成为一个组的key:value对。Python - 按键汉明距离对defaultdict值进行分组

输入文件采用FASTA格式,其中标题为关键字,值为序列(因为多个序列基于来自原始fasta文件的blast报告具有相同的头部,所以使用defaultdict)。

这是我到目前为止有:

!/usr/bin/env python 

import sys 
from collections import defaultdict 
import itertools 

inp = sys.argv[1]              # input fasta file; format '>header'\n'sequence' 

with open(inp, 'r') as f: 
     h = [] 
     s = [] 
     for line in f: 
       if line.startswith(">"): 
         h.append(line.strip().split('>')[1])   # append headers to list 
       else: 
         s.append(line.strip())       # append sequences to list 

seqs = dict(zip(h,s))             # create dictionary of headers:sequence 

print 'Total Sequences: ' + str(len(seqs))        # Numb. total sequences in input file 

groups = defaultdict(list) 

for i in seqs: 
     groups['_'.join(i.split('_')[1:])].append(seqs[i])      # Create defaultdict with sequences in lists with identical headers 

def hamming(str1, str2): 
     """ Simple hamming distance calculator """ 
     if len(str1) == len(str2): 
       diffs = 0 
       for ch1, ch2 in zip(str1,str2): 
         if ch1 != ch2: 
           diffs += 1 
       return diff 

keys = [x for x in groups] 

combos = list(itertools.combinations(keys,2))       # Create tupled list with all comparison combinations 

combined = defaultdict(list)           # Defaultdict in which to place groups 

for i in combos:              # Combo = (A1_B1_STRING2, A2_B2_STRING2) 
     a1 = i[0].split('_')[0] 
     a2 = i[1].split('_')[0] 

     b1 = i[0].split('_')[1]           # Get A's, B's, C's 
     b2 = i[1].split('_')[1] 

     c1 = i[0].split('_')[2] 
     c2 = i[1].split('_')[2] 

     if a1 == a2 and b1 == b2:          # If A1 is equal to A2 and B1 is equal to B2 
       d = hamming(c1, c2)          # Get distance of STRING1 vs STRING2 
       if d <= 2:            # If distance is less than or equal to 2 
         combined[i[0]].append(groups[i[0]] + groups[i[1]])  # Add to defaultdict by combo 1 key 

print len(combined) 
for c in sorted(combined): 
     print c, '\t', len(combined[c]) 

的问题是,因为预计代码是行不通的。在组合defaultdict中打印键时;我清楚地看到有很多可以合并的东西。但是,组合defaultdict的长度大约是原始大小的一半。

编辑

替代没有itertools.combinations:

for a in keys: 
     tocombine = [] 
     tocombine.append(a) 
     tocheck = [x for x in keys if x != a] 
     for b in tocheck: 
       i = (a,b)            # Combo = (A1_B1_STRING2, A2_B2_STRING2) 
       a1 = i[0].split('_')[0] 
       a2 = i[1].split('_')[0] 

       b1 = i[0].split('_')[1]           # Get A's, B's, C's 
       b2 = i[1].split('_')[1] 

       c1 = i[0].split('_')[2] 
       c2 = i[1].split('_')[2] 

       if a1 == a2 and b1 == b2:          # If A1 is equal to A2 and B1 is equal to B2 
         if len(c1) == len(c2):           # If length of STRING1 is equal to STRING2 
           d = hamming(c1, c2)          # Get distance of STRING1 vs STRING2 
           if d <= 2: 
             tocombine.append(b) 
     for n in range(len(tocombine[1:])): 
       keys.remove(tocombine[n]) 
       combined[tocombine[0]].append(groups[tocombine[n]]) 

final = defaultdict(list) 
for i in combined: 
     final[i] = list(itertools.chain.from_iterable(combined[i])) 

然而,这些方法,我仍然失踪,少数不符合任何人。

+0

你错过了你的汉明机制的文档串一“,这是造成格式错误,你可以把在那里我会编辑就为你,但堆栈溢出至少需要6个字符编辑:?/ –

+0

改变了它。在我遇到? –

回答

1

我想我看到一个问题与您的代码考虑此方案:

0: A_B_DATA1 
1: A_B_DATA2  
2: A_B_DATA3 

All the valid comparisons are: 
0 -> 1 * Combines under key 'A_B_DATA1' 
0 -> 2 * Combines under key 'A_B_DATA1' 
1 -> 2 * Combines under key 'A_B_DATA2' **opps 

我可以想象你想所有这三个下1个键组合。但是考虑这样的场景:

0: A_B_DATA111 
1: A_B_DATA122  
2: A_B_DATA223 

All the valid comparisons are: 
0 -> 1 * Combines under key 'A_B_DATA111' 
0 -> 2 * Combines under key 'A_B_DATA111' 
1 -> 2 * Combines under key 'A_B_DATA122' 

现在,它变得有点棘手,因为行0是距离2,从第1行和行1是距离2第2行,但你可能不希望他们都在一起,行0距第2排距离3!

这里是一个可行的解决方案假设这的一个例子是你想要的输出看起来像:

def unpack_key(key): 
    data = key.split('_') 
    return '_'.join(data[:2]), '_'.join(data[2:]) 

combined = defaultdict(list) 
for key1 in groups: 
    combined[key1] = [] 
    key1_ab, key1_string = unpack_key(key1) 
    for key2 in groups: 
     if key1 != key2: 
      key2_ab, key2_string = unpack_key(key2) 
      if key1_ab == key2_ab and len(key1_string) == len(key2_string): 
       if hamming(key1_string, key2_string) <= 2: 
        combined[key1].append(key2) 

在我们的第二个例子中,这将导致下面的单词表,如果这不是你正在寻找的答案,你能准确地输出这个例子的最终字典应该是什么吗?

A_B_DATA111: ['A_B_DATA122'] 
A_B_DATA122: ['A_B_DATA111', 'A_B_DATA223'] 
A_B_DATA223: ['A_B_DATA122'] 

请记住这是一个O(n^2)算法,这意味着当您的密钥集越来越大时它不可伸缩。

+0

这个问题有什么想法我明白你的意思与关键问题是什么。我想我已经想出了一个解决方案,不使用itertoosl组合,看看编辑 –

+0

你能给我你会期待什么输出成为上述情况 –

+0

输出应该是默认字典,其中键包含通过上述标准值:所有值原始报头必须具有相同的A和B,并且每个字符串应该由<= 2个字符是不同的。 –