2013-08-23 16 views
7

Podstawowy problem to:Porównanie wielu ciągów znaków arbitralnych ze znakami zorientowanymi

Szukam algorytmu do obliczenia maksymalnej oszczędnej odległości między zestawem ciągów. W przypadku odległości rozumiem coś podobnego do Damerau–Levenshtein distance, tj. Minimalną liczbę usunięć, wstawień, podstawień i transpozycji znaków lub bloków sąsiednich znaków. Ale zamiast zwykłych ciągów chcę zbadać ciągi znaków o zorientowanych znakach.

Zatem łańcuch może wyglądać tak:

  1. (A,1) (B,1) (C,1) (D,1)

i możliwych pochodnych mogą być:

  1. (A,1) (C,0) (B,0) (D,1)
  2. (A,1) (C,1) (B,1) (D,1)
  3. (A,1) (B,0) (C,0) (D,1)

Gdzie A,B,C,D są tożsamość i charakter 1 = forward i 0 = reverse.

Tutaj pochodna 1. miałaby odległość 2, ponieważ można wyciąć blok BC i ponownie wkleić go odwrócony (1 cięcie, 1 wklej). Pochodna 2. miałaby również 2, ponieważ można wyciąć C i ponownie wkleić go przed B (1 cięcie, 1 wklej), natomiast numer 3. wymagałoby 4 operacji (2 cięcia, 2 pasty) do przekształcenia. Analogiczne, delecję lub insercję bloków dałoby odległość 1.

Jeżeli można określić (X,0) i (X,1) w dwóch różnych zorientowanych znakami nie (X0, X1) dla wszystkich możliwych X Przykład 3. spowodowałoby odległości 2 ponieważ może wówczas wyciąć blok B1C1 i wstawić blok B0C0 w dwóch etapach.

przykładem prawdziwy świat:

genów w genomie bakterii można uznać za zorientowanych postać (A, 0), (B, 0) ... Dzięki określaniu odległości, sekwencji genomowego orientację Geny homologów w dwóch pokrewnych bakteriach można wykorzystać jako ewolucyjny ślad markerowy. Fakt, że genomy bakterii są cyklicznymi łańcuchami wprowadza dodatkowy warunek graniczny ABC jest równy BCA.

Prawdziwe genomy mają unikatowe geny bez odpowiednika u partnera, co powoduje powstanie znaku właściciela miejsca @. Te elementy zastępujące zmniejszają zawartość informacyjną porównania do niższej granicy, ponieważ np. (A, 1) (B, 1) @ (C, 1) można przekształcić do (A, 1) @@@ (B, 1) @ (C, 1) przez wstawienie bloku @@@. Jednak orientacja przywraca częściowo treść informacji, ponieważ może się okazać, że (A, 1) @@@ (B, 0) @ (C, 1) wskazuje minimalną odległość 3. Jeszcze lepiej byłoby algorytmem porównywania wielu powiązanych sekwencji (genomy) jednocześnie, ponieważ można wtedy znaleźć półprodukty w historii ewolucyjnej, co zwiększa rozdzielczość.

Zdaję sobie sprawę, istnieje kilka pytań już opublikowanych na temat porównania ciągów tekstowych. Ale nie można ich łatwo rozszerzyć, aby uwzględnić orientację. Ponadto istnieje wiele metod radzenia sobie z sekwencjami biologicznymi, w szczególności w przypadku analizy wielu sekwencji.Jednak są one ograniczone do sekwencji makrocząsteczek, które nie istnieją w alternatywnych orientacjach i zwykle wywołują określone wagi dla każdego konkretnego dopasowania postaci.

Jeśli istnieje już biblioteka Pythona, która pozwoliłaby na niezbędne dostosowanie w celu rozwiązania tego problemu, byłoby to fantastyczne. Ale bardzo przydatny mógłby być dowolny odpowiedni algorytm orientacji.

Odpowiedz

1

wierzę następujący kod może pomóc:

from itertools import permutations 
from random import randint 
from pprint import pprint 


def generate_genes(): 
    """ 
    Generates a boilerplate list of genes 
    @rtype : list 
    """ 
    tuple_list = [] 

    for i in range(16): 

     binary_var = bin(i)[2:] 

     if len(binary_var) != 4: 
      binary_var = "0" * (4 - len(binary_var)) + binary_var 

     tuple_list.append([('A', (1 if binary_var[0] == '1' else 0)), 
          ('B', (1 if binary_var[1] == '1' else 0)), 
          ('C', (1 if binary_var[2] == '1' else 0)), 
          ('D', (1 if binary_var[3] == '1' else 0))]) 

    return tuple_list 


def all_possible_genes(): 
    """ Generates all possible combinations of ABCD genes 
    @return: returns a list of combinations 
    @rtype: tuple 
    """ 
    gene_list = generate_genes() 
    all_possible_permutations = [] 
    for gene in gene_list: 
     all_possible_permutations.append([var for var in permutations(gene)]) 

    return all_possible_permutations 


def gene_stringify(gene_tuple): 
    """ 
    @type gene_tuple : tuple 
    @param gene_tuple: The gene tuple generated 
    """ 

    return "".join(str(var[0]) for var in gene_tuple if var[1]) 


def dameraulevenshtein(seq1, seq2): 
    """Calculate the Damerau-Levenshtein distance between sequences. 

    This distance is the number of additions, deletions, substitutions, 
    and transpositions needed to transform the first sequence into the 
    second. Although generally used with strings, any sequences of 
    comparable objects will work. 

    Transpositions are exchanges of *consecutive* characters; all other 
    operations are self-explanatory. 

    This implementation is O(N*M) time and O(M) space, for N and M the 
    lengths of the two sequences. 

    >>> dameraulevenshtein('ba', 'abc') 
    2 
    >>> dameraulevenshtein('fee', 'deed') 
    2 

    It works with arbitrary sequences too: 
    >>> dameraulevenshtein('abcd', ['b', 'a', 'c', 'd', 'e']) 
    2 
    """ 
    # codesnippet:D0DE4716-B6E6-4161-9219-2903BF8F547F 
    # Conceptually, this is based on a len(seq1) + 1 * len(seq2) + 1 matrix. 
    # However, only the current and two previous rows are needed at once, 
    # so we only store those. 
    oneago = None 
    thisrow = range(1, len(seq2) + 1) + [0] 
    for x in xrange(len(seq1)): 
     # Python lists wrap around for negative indices, so put the 
     # leftmost column at the *end* of the list. This matches with 
     # the zero-indexed strings and saves extra calculation. 
     twoago, oneago, thisrow = oneago, thisrow, [0] * len(seq2) + [x + 1] 
     for y in xrange(len(seq2)): 
      delcost = oneago[y] + 1 
      addcost = thisrow[y - 1] + 1 
      subcost = oneago[y - 1] + (seq1[x] != seq2[y]) 
      thisrow[y] = min(delcost, addcost, subcost) 
      # This block deals with transpositions 
      if (x > 0 and y > 0 and seq1[x] == seq2[y - 1] 
       and seq1[x - 1] == seq2[y] and seq1[x] != seq2[y]): 
       thisrow[y] = min(thisrow[y], twoago[y - 2] + 1) 
    return thisrow[len(seq2) - 1] 


if __name__ == '__main__': 
    genes = all_possible_genes() 
    list1 = genes[randint(0, 15)][randint(0, 23)] 
    list2 = genes[randint(0, 15)][randint(0, 23)] 

    print gene_stringify(list1) 
    pprint(list1) 
    print gene_stringify(list2) 
    pprint(list2) 
    print dameraulevenshtein(gene_stringify(list1), gene_stringify(list2)) 

kuponów

Michael Homer for the algorithm

Powiązane problemy