Confronto di sequenze multiple di stringhe arbitrarie con caratteri orientati

Question

Il problema principale è:

Sto cercando un algoritmo per calcolare una distanza massima parsimoniosa tra un insieme di stringhe. Con la distanza intendo qualcosa di simile allo Damerau–Levenshtein distance cioè il numero minimo di cancellazioni, inserzioni, sostituzione e trasposizione di caratteri o blocchi di caratteri adiacenti. Ma invece di stringhe regolari voglio indagare su stringhe con caratteri orientati.

Così, una stringa potrebbe sembrare:

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

e possibili derivati ​​potrebbe essere:

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)

Dove A,B,C,D sono le identità dei personaggi e 1 = forward e 0 = reverse.

Qui, la derivata 1. avrebbe la distanza 2, dal momento che è possibile ritagliare il blocco BC e incollarlo nuovamente invertito (1 taglio, 1 incolla). Derivativo 2. avrebbe anche 2, poiché è possibile ritagliare C e incollarlo di fronte a B (1 taglio, 1 incolla) mentre il numero 3. richiederebbe 4 operazioni (2 tagli, 2 paste) da trasformare. Analogamente, delezione o inserzione di blocchi produrrebbe una distanza 1.

Se definirebbe (X,0) e (X,1) come due differenti caratteri non orientato (X0, X1) per tutte le possibili X, esempio 3. comporterebbe una distanza di 2 in quanto si potrebbe allora ritaglia il blocco B1C1 e inserisci il blocco B0C0 in due passaggi.

Un esempio pratico:

I geni in un genoma batterico può essere considerato il carattere orientato (A, 0), (B, 0) ... Per determinare la distanza sequenza, l'orientamento della genomica i geni omologhi in due batteri correlati potrebbero essere utilizzati come traccia di un marker evolutivo. Il fatto che i genomi batterici siano stringhe circolari introduce la condizione di bordo aggiuntiva ABC uguale a BCA.

I genomi reali hanno geni univoci senza equivalenti in un partner che danno origine a un carattere di segnaposto @. Quei titolari di posto riducono il contenuto di informazione del confronto ad un limite inferiore, poiché ad es. (A, 1) (B, 1) @ (C, 1) può essere trasformato in (A, 1) @@@ (B, 1) @ (C, 1) inserendo il blocco @@@. Tuttavia, l'orientamento ripristina parzialmente il contenuto delle informazioni poiché è possibile trovare (A, 1) @@@ (B, 0) @ (C, 1) che indica una distanza minima di 3. Ancora meglio sarebbe un algoritmo per confrontare più sequenze correlate (genomi) simultaneamente, dal momento che è possibile trovare intermedi nella storia evolutiva, che aumenta la risoluzione.

Mi rendo conto, ci sono diverse domande già pubblicate sul confronto delle stringhe di testo. Ma non riescono ad essere facilmente espandibili per includere l'orientamento. Inoltre, esiste una vasta gamma di metodi per trattare le sequenze biologiche, in particolare per l'analisi di sequenze multiple.Tuttavia, quelli sono limitati a sequenze di macromolecole che non esistono in orientamenti alternati e di solito invocano pesi specifici per una particolare corrispondenza di caratteri.

Se esiste già una libreria Python che consentirebbe la personalizzazione necessaria per risolvere il problema, sarebbe fantastico. Ma qualsiasi algoritmo adatto all'orientamento sarebbe molto utile.

Answer 1

Credo che il seguente codice può aiutarti:

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)) 

Credits

Michael Homer for the algorithm