Python Programming for Biology: Bioinformatics and Beyond

Figure 21.4.  The range of all possible outcomes for two dice. The outcomes which all

contribute to each numeric total are grouped together. The probability of getting a total of

seven is then calculated from the sum of the outcomes for that event.

When using dice rolls or DNA bases as examples it is clear that individual outcomes are

mutually exclusive, in that a die provides only one number and that only one of the DNA

nucleotides  is  possible  in  a  given  position.  However,  it  is  also  possible  to  define  events

that  are  not  mutually  exclusive.  For  example,  given  two  DNA  positions  we  could

investigate the event that there is at least one A nucleotide present (7 of 16 possibilities)

and  the  separate  event  that  the  two  nucleotides  are  different  (12  of  16  possibilities).

Clearly  these  are  not  exclusive  events  because  some  outcomes  are  present  in  both  sets.

The set of common outcomes, which apply to both one event and another, is referred to as

the  intersection.  Furthermore  knowing  what  the  intersection  is  allows  us  to  calculate

things like the probability that one event but not the other occurs.

We  can  use  Python  to  calculate  the  probabilities  for  the  illustrated  events  and  their

intersection  in  the  case  where  the  probabilities  of  the  outcomes  are  not  equal.  Here  the

probabilities  of  the  two-letter  outcome  are  obtained  by  multiplying  the  individual

probabilities for a nucleotide recorded in letterProbs, i.e. the probability of having x at the

first position times the probability of having y at the second. The outcomes are then tested

and added to the appropriate event, which we store as Python sets. We can then use the set

operation  ‘&’,  which  generates  a  new  set  with  the  common  elements.  Finally,  the

probabilities of the events are simply calculated from the summation of the probabilities

within each.

probs = {}

letters = ['G','C','A','T']

event1 = set()

event2 = set()

for x in letters:

for y in letters:

outcome = (x,y)

probs[outcome] = letterProbs[x] * letterProbs[y]

if 'A' in outcome:

event1.add(outcome)

if x != y:

event2.add(outcome)

---

intersection = event1 & event2

pEvent1 = sum([probs[xy] for xy in event1]) # 0.43156

pEvent2 = sum([probs[xy] for xy in event2]) # 0.74994

pEvent1and2 = sum([probs[xy] for xy in intersection]) # 0.37102