Matching limited elements using Postgres IN clause

So you have two tables, one a sequence table holds unique sequence ids.  The other table hold sequence mutations and has a 1 to many relationship with the sequence table.  So for a given sequence id in the mutation table, there are many records, one for each mutation.

Now you have two sets of mutations you are interested in, a primary set A and a secondary set B.  You want to know how many sequences contain 1 and only 1 mutation from set A and 1 and only 1 from set B.


select count(*) from (
SELECT distinct a.seq_id
FROM sequence a
JOIN seq_mutations b ON b.seq_id = a.seq_id
JOIN seq_mutations c ON c.seq_id = a.seq_id
WHERE b.mutation IN ('mutA','mutB','mutC')
AND c.mutation IN ('mut1', 'mut2', 'mut3', 'mut4', 'mut5')
GROUP BY a.seq_id
HAVING 
  count(distinct b.mutation) = 1 
  AND count(distinct c.mutation) = 1) a

This query gives you exact, fine control over the conditions you want to test for.  If you then want to test for 1 and only 1 mutation in set A and 2 and only 2 mutations in set B, just change

count(distinct c.mutation) = 2

in the HAVING clause.  This also lets you use > and <, so you can ask for 3 or more mutations by using

count(distinct c.mutation) >= 3

You can also vary the count for b.mutation to adjust the number of mutations allowed from set A, so any combination can now be extracted.

Thanks to Jake for the code help.

Generate all possible proteins from ambiguous DNA

This had me stumped for awhile, but this works pretty well.  Does NOT handle stop codons or gap characters like '-'.  Requires BioPython

import itertools
from Bio.Seq import Seq
from Bio.Data import CodonTable
from Bio.Data import IUPACData</pre>

# Takes Bio.Seq.Seq object as input
# Returns list of all possible proteins
# Assumes sequence is in frame +1
def generateProtFromAmbiguousDNA(s):
   std_nt = CodonTable.unambiguous_dna_by_name["Standard"]
   nonstd = IUPACData.ambiguous_dna_values
   aa_trans = []
   for i in range(0,len(s),3):
      codon = s.tostring()[i:i+3]
      aa = CodonTable.list_possible_proteins(codon,std_nt.forward_table,nonstd) 
      aa_trans.append(aa)
   proteins = list(itertools.product(*aa_trans))
   possible_proteins = []
   for x in proteins:
      possible_proteins.append("".join(x))
   return possible_proteins

def main():
   a = Seq('ATGGCARTTGTAHAC')
   print "DNA: ",a.tostring()
   print "Proteins:"
   foo = generateProtFromAmbiguousDNA(a)
   for s in foo: print s

if __name__ == '__main__':
   main()

Creating a quick codon table

I didn't think this up, the code comes from Peter Collingridge here.  But it is rather elegant.

bases = ['t', 'c', 'a', 'g']
codons = [a+b+c for a in bases for b in bases for c in bases]
amino_acids = "F F L L S S S S Y Y stop stop C C stop W L L L L P P P P H H Q Q R R R R I I I M T T T T N N K K S S R R V V V V A A A A D D E E G G G G".split(' ')
codon_table = dict(zip(codons, amino_acids))