Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
2  
As for me, the phrase a transmembrane protein that is supposed to selectively filter K ions contains an amino acid sequence "TVGYG" in this site is worth a +1. –  belisarius May 27 '12 at 22:45
    
The syntax information on the function's doc page clearly lists its arguments as strings, not patterns. If your question really only is "Is it naive to expect more functionality or am I too demanding on Mathematica?" I'm afraid I have to vote to close it as either "not constructive" or "not a real question". –  Sjoerd C. de Vries May 27 '12 at 22:52
1  
If you could rephrase your question to something like "how can I mimic the dreamed behavior of SequenceAlignment["A", Alternatives["A", "B"]]?" then I think you should be OK. Otherwise, your question doesn't adhere to the rules as delineated in the FAQ. –  Sjoerd C. de Vries May 27 '12 at 23:07
1  
"BLAST" is specifically listed as a permissible value for the SimilarityRules option for SequenceAlignment. Perhaps it might be what you are looking for? –  WReach May 27 '12 at 23:54
1  
@WReach The SimilarityRules option specifies how acceptable a candidate candidate character is for matching another. In mma "BLAST" seems to me just the name of a specific square table, it did not help aligning a pattern. you might browse <en.wikipedia.org/wiki/BLAST>; –  wilbert van meerwijk May 28 '12 at 4:49
show 4 more comments

1 Answer

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.

share|improve this answer
    
@Mr.Wizard thanks for your editing, I must learn that too! –  wilbert van meerwijk May 28 '12 at 11:03
    
allparents may also be obtained as follows: Tuples[StringExpression @@ StringReplaceList[Characters@"TVGYG", Reverse /@ cr]] although I suppose that is not too much use to you (where cr is the list of codon rules) –  TomD May 28 '12 at 15:51
    
Or: StringJoin @@@ Tuples[StringReplaceList[Characters@"TVGYG", Reverse /@ cr]] –  TomD May 28 '12 at 16:39
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.