Creating a version of SequenceAlignment that accepts patterns

Question

SequenceAlignment seems like a function that is not (yet?) fully integrated into Mathematica. I want a function that accepts general patterns instead.

I naively fed a StringPattern instead of a literal string to SequenceAlignment and met silence, though MatchQ is positive:

{StringMatchQ["A", Alternatives["A", "B"]], SequenceAlignment["A", Alternatives["A", "B"]]}
(* Out= {True, SequenceAlignment["A", "A" | "B"]} *)

One approach would be to write out all alternatives and repeatedly feed them to SequenceAlignment. The question is strongly related to this one: How to generally match, unify and merge patterns?.

Background: I want to answer questions like "What DNA sequences give rise to a given protein sequence and where can you find them?". Then one naturally stumbles across Alternatives.

For example, a transmembrane protein that is supposed to selectively filter K ions contains an amino acid sequence "TVGYG". I want to find all its DNA parents, ie DNA code of all proteins that may contain such a filter. Due to the redundancy of the genetic code, that subsequence can be coded by a diversity of DNA sequences, giving rise to a StringExpression with alternatives:

In[80]:= StringReplace["TVGYG", backlist]    
Out[80]= "ACA" | "ACC" | "ACG" | "ACT" ~~ "GTA" | "GTC" | "GTG" | "GTT" ~~ 
 "GGA" | "GGC" | "GGG" | "GGT" ~~ "TAC" | "TAT" ~~ 
 "GGA" | "GGC" | "GGG" | "GGT"

The first residue, "T" for Threonine, may thus be coded by either "ACA" or "ACC" or "ACG" or "ACT", etc.

Answer 1

A brute force approach, generating a literal list of all alternatives be explained as follows:

I use a substitution list for the genetic code:

In[3]:= CodonRules

Out[3]= {"TCA" -> "S", "TCC" -> "S", "TCG" -> "S", "TCT" -> "S", "TTC" -> "F", 
 "TTT" -> "F", "TTA" -> "L", "TTG" -> "L", "TAC" -> "Y", "TAT" -> "Y", 
 "TAA" -> "_", "TAG" -> "_", "TGC" -> "C", "TGT" -> "C", "TGA" -> "_", 
 "TGG" -> "W", "CTA" -> "L", "CTC" -> "L", "CTG" -> "L", "CTT" -> "L", 
 "CCA" -> "P", "CCC" -> "P", "CCG" -> "P", "CCT" -> "P", "CAC" -> "H", 
 "CAT" -> "H", "CAA" -> "Q", "CAG" -> "Q", "CGA" -> "R", "CGC" -> "R", 
 "CGG" -> "R", "CGT" -> "R", "ATA" -> "I", "ATT" -> "I", "ATC" -> "I", 
 "ATG" -> "M", "ACA" -> "T", "ACC" -> "T", "ACG" -> "T", "ACT" -> "T", 
 "AAC" -> "N", "AAT" -> "N", "AAA" -> "K", "AAG" -> "K", "AGC" -> "S", 
 "AGT" -> "S", "AGA" -> "R", "AGG" -> "R", "GTA" -> "V", "GTC" -> "V", 
 "GTG" -> "V", "GTT" -> "V", "GCA" -> "A", "GCC" -> "A", "GCG" -> "A", 
 "GCT" -> "A", "GAC" -> "D", "GAT" -> "D", "GAA" -> "E", "GAG" -> "E", 
 "GGA" -> "G", "GGC" -> "G", "GGG" -> "G", "GGT" -> "G"}

to produce a reverse translation list, (demonstrating the redundancy of the code):

backlist = 
 Rule @@@ ({#[[1, 2]], Alternatives @@ #[[All, 1]]} & /@ 
    GatherBy[CodonRules, Last])
backlist::usage = 
  "backlist gives all DNA codes for each amino acid residue. It is the \
re(con?)verse of 'CodonRules'";

Out[4]= {"S" -> "TCA" | "TCC" | "TCG" | "TCT" | "AGC" | "AGT", "F" -> "TTC" | "TTT", 
 "L" -> "TTA" | "TTG" | "CTA" | "CTC" | "CTG" | "CTT", "Y" -> "TAC" | "TAT", 
 "_" -> "TAA" | "TAG" | "TGA", "C" -> "TGC" | "TGT", 
 "W" -> Alternatives["TGG"], "P" -> "CCA" | "CCC" | "CCG" | "CCT", 
 "H" -> "CAC" | "CAT", "Q" -> "CAA" | "CAG", 
 "R" -> "CGA" | "CGC" | "CGG" | "CGT" | "AGA" | "AGG", 
 "I" -> "ATA" | "ATT" | "ATC", "M" -> Alternatives["ATG"], 
 "T" -> "ACA" | "ACC" | "ACG" | "ACT", "N" -> "AAC" | "AAT", 
 "K" -> "AAA" | "AAG", "V" -> "GTA" | "GTC" | "GTG" | "GTT", 
 "A" -> "GCA" | "GCC" | "GCG" | "GCT", "D" -> "GAC" | "GAT", 
 "E" -> "GAA" | "GAG", "G" -> "GGA" | "GGC" | "GGG" | "GGT"}

Now all DNA quintuplets that code for the 5 amino acid residues in "TVGYG" may be generated:

In[6]:= (allparents = Tuples[List @@@ StringReplace["TVGYG", backlist]]) // Length

Out[6]= 512

Then finally those 512 can be compared with Mathematica 's curated data and elicit a list of exact matches:

alles = GenomeLookup[#] & /@ allparents (* Took 150 seconds! *)

It turns out that 216 different parents of the 512 are actually present in the human genome, distributed over in total 608 places.

The drawback of this brute force approach is of course that combinatorial explosion may stop the procedure at the Tuples stage.