Simple Mardi Gras King Cake

.... as they say in Louisiana, laissez le bon temps roulet!  Let the good times roll!

As Mardi Gras season is in full swing, I set out to make a king cake from a kit my sister brought me from New Orleans.  It seemed pretty authentic, so I was very excited to try it out.

Ingredients:
The kit contained the 'king cake mix', yeast, praline sugar, glaze, three packets of colored sprinkles, and of course the baby!  I added 1/3 cup of warm water, 3 eggs and 2 egg yolks, and 1 1/2 sticks of butter.

Recipe for king cake (following a kit):
1) Preheat oven to 350. Mix king cake mix, yeast, eggs, warm water and room temperature butter.  Approx 3-4 mins

2) Cover dough with plastic wrap and set in warm place for 30 minutes to rise.  Next, knead dough on a lightly floured surface.  This dough was very sticky and not easy to work with, so I kept adding more flour, probably more than I should have. This made the dough a bit drier than I would have liked.

3) Brush 2 tablespoons of melted butter over the surface and sprinkle the praline sugar topping all over. Next, start to roll up the dough like a jelly roll.  Pinch ends closed.

 4) Place dough on a cookie sheet and bring the two ends together to form a circle or oval. Let rise for 30 minutes until the dough has doubled in size.

5) Bake for 20-30 minutes until lightly browned. Remove from oven and let cool completely.  At this point put the baby somewhere under the cake.

6) To decorate, mix glaze with a bit of water until the consistency is of a thick icing.  Add decorative sprinkle or sugars.

Happy Mardi Gras! 

Hardy-Weinberg Genotype Frequencies

An important principle in population genetics is called the Hardy-Weinberg principle (or H-W equilibrium or H-W law) which "describes the equilibrium state of a single locus in a randomly mating diploid population that is free of other evolutionary forces such as mutation, migration, and genetic drift" [Population Genetics: A Concise Guide by John Gillespie, 2nd edition]. This means the allele frequencies and genotype frequencies of a single locus are in equilibrium and stay constant generation to generation (in an ideal world) until outside forces act on them such as mutation, migration, genetic drift (see post on simulating genetic drift) or natural selection (in the real world).

Let's try to put some notation to all this.  If we have two alleles $A$ and $a$ with genotype frequencies
\begin{eqnarray*}
P(AA) & = & u \\
P(Aa) & = & v \\
P(aa) & = & w
\end{eqnarray*}
then without any more assumptions, the allele frequencies are easily estimated:
\[ P(A) = u + 0.5v = p \hspace{2in} P(a) = w + 0.5v = q \]
If we have the allele frequencies, we need a few more assumptions to calculate the genotype frequencies.  This is where the Hardy-Weinberg principle comes into play.

H-W Assumptions:

1. diploid organism, infinite population, discrete generations
2. random mating
3. no outside forces at play (e.g. no selection, no migration, no mutation)
4. equal (or unequal)* initial genotype frequencies in two sexes

When a population is in H-W equilibrium, the alleles that make up a genotype are thought of as randomly sampling alleles from the population.  Thus, we can estimate genotype frequencies from allele frequencies as followed:
\begin{eqnarray*}
P(AA) & = & P(A) P(A) = p^2 \\
P(Aa) & = & 2P(A) P(a) = 2*p*q\\
P(aa) & = & P(a) P(a) = q^2
\end{eqnarray*}
*In dioecious species (individual is either male or female), H-W equilibrium can be reached in just two generations assuming unequal genotype frequencies in two sexes (will be achieved in one generation assuming equal genotype frequencies).

Rcode:
Assume equal genotype frequencies in two sexes. Show HW equilibrium in one generation.

> # Start with two alleles A, a with respective allele frequencies p.0, q.0

> p.0 = 0.2
> q.0 = 1 - p.0
> 
> # After one round of mating
> # Calculate genotype frequencies (Assumes HWE)
> AA = p.0*p.0
> AA
[1] 0.04
> 
> Aa = 2*p.0*q.0
> Aa
[1] 0.32
> 
> aa = q.0*q.0
> aa
[1] 0.64
> 
> 
> # Calulate new allele frequencies (Does not assume HWE)
> p.1 = AA + 0.5*Aa 
> p.1
[1] 0.2
> 
> q.1 = 1 - p.1
> q.1
[1] 0.8
> 
> # After two rounds of mating 
> # Calculate genotype frequencies (Assumes HWE)
> AA = p.1*p.1
> AA
[1] 0.04
> 
> Aa = 2*p.1*q.1
> Aa
[1] 0.32
> 
> aa = q.1*q.1
> aa
[1] 0.64

Therefore, we see the genotype frequencies after one round of mating were equal to the genotype frequencies after two rounds of mating (i.e. H-W equilibrium is attained).

Now this is all nice, but the H-W assumptions are almost never met in the real world.  For example, there may be small population sizes, deviations from random mating (e.g. assortative mating, inbreeding) and there more than likely outside forces at play such as mutation, migration and selection.  To test for deviations from H-W equilibrium, we can use the $\chi^2$ Goodness-of-Fit test or Exact tests (better for small sample sizes) in standard software tools such as R or PLINK.

Example: If we have two alleles $A$ and $a$, we can compare the observed genotype counts with values expected under H-W equilibrium:

\begin{eqnarray*}
Genotype &  Observed &  Expected \\
AA & n_{AA} & np^2 \\
Aa & n_{Aa} & 2np(1-p) \\
aa & n_{aa} & n(1-p)^2
\end{eqnarray*}
where $p = P(A) = (n_{Aa} + 2n_{AA}) / 2n$

Markov Chains: Stationary Distributions

In my first post on Markov chains, I introduced a particular type of stochastic process with the basic first-order Markov property (in discrete time).  We call this a Markov chain. Examples explaining what are 1-step and 2-step transition probabilities and how to compute a marginal or unconditional probability of being in a particular state at time $n$ were all provided. This post will be a bit more technical, but ultimately I hope to end with a concrete example (with R code) to help make it more clear.

In this post, I will introduce stationary distributions by answering the following questions:

1. What is a stationary distribution? 
2. What properties are needed to guarantee the unique existence of the stationary distribution? 
3. How to compute the stationary distribution? 

Q1: A stationary distribution $\pi$ is a probability distribution such that when the Markov chain reaches the stationary distribution, then it remains in that probability distribution forever. This means whenever you are interested in asking questions like "In the long-run, what it the proportion of time the Markov chain will be in a particular state?" (such as the rain or no rain example), you're interested in computing the stationary distribution (if it's possible).

Let's start with a bit of notation.  Define $\pi_j$ as the limiting probability that the process will be in state $j$ at time $n$ or
\[ \pi_j = \lim_{n \rightarrow \infty} P_{ij}^n \]
Consider:
\begin{eqnarray*}
\lim_{n \rightarrow \infty} P[X_{n+1} = j] & = & lim_{n \rightarrow \infty}  \sum_{i = 0}^{\infty} P[X_{n+1} = j | X_n = i] P[X_n = i] \\
& = & \sum_{i = 0}^{\infty}  lim_{n \rightarrow \infty}  P_{ij} P[X_n = i] \hspace{.2in} (Fubini's \hspace{.05in}Theorem)
\end{eqnarray*}
which defines the stationary distribution:
\[ \pi_j = \sum_{i=0}^{\infty} P_{ij} \pi_i \]

Q2: This question is more difficult to answer.
Short answer: The Markov chain needs to be a (1) finite (2) aperiodic and (3) irreducible chain.
Longer answer: Assume $P$ is the transition probability matrix for a Markov chain with the following properties:

1. finite = finite number of states
2. aperiodic = state $i$ has period $d$ if $P_{ii}^n = 0$ whenever $n$ is not divisible by $d$ and $d$ is the largest integer with this property
3. irreducible = two states that communicate are in the same class. classes are disjoint or identical. If there exists only 1 class, then the Markov chain is irreducible

Then, there exists a unique probability distribution satisfying
\[ \pi_j = \sum_{i = 0}^{\infty} \pi_i P_{ij} \hspace{1in} j \geq 0, \hspace{.5in} \sum_j \pi_j = 1 \]
In vector notation $\pi^{'} = \pi^{'} P$ and $\pi^{'} 1 = 1$ where $\pi^{'} = (\pi_0, \ldots, \pi_s)$.

This theorem also gives us the fact that if $P$ is finite, aperiodic and irreducible, then not only does the unique stationary distribution exist, but also $\pi_j$ = the limiting probability that the process will be in state $j$ at time $n$ equals the long-run proportion of time that the process will be in state $j$ [Sheldon Ross, 2007]. Or

\[ P^n \rightarrow \left[ \begin{array}{ccc}

\pi_0 & \ldots & \pi_s \\
\pi_0 & \ldots & \pi_s \\
\vdots & & \\
\pi_0 & \ldots & \pi_s
\end{array} \right]  \]

Q3: Back to our rain example. The transition probability matrix $P$ was given by

\[ P = \left[ \begin{array}{cc}

\alpha & 1-\alpha  \\
\beta & 1- \beta \\
\end{array} \right]  \]
where $\alpha$ = 0.7 and $\beta$= 0.4.  As we consider the 2-step transition probability matrix ($P_{ij}^2$), we can consider the 4-step ($P_{ij}^4$)or $n$-step transition probability matrix ($P_{ij}^n$) as $n \rightarrow \infty$.

Rcode:

trans.mat <- matrix(c(0.7, 0.3, 0.4, 0.6), 2,2, byrow = TRUE)
> trans.mat %*% trans.mat %*% trans.mat
      [,1]  [,2]
[1,] 0.583 0.417
[2,] 0.556 0.444

> trans.mat %*% trans.mat %*% trans.mat %*% trans.mat
       [,1]   [,2]
[1,] 0.5749 0.4251
[2,] 0.5668 0.4332
> trans.mat %*% trans.mat %*% trans.mat %*% trans.mat %*% trans.mat
        [,1]    [,2]
[1,] 0.57247 0.42753
[2,] 0.57004 0.42996
> trans.mat %*% trans.mat %*% trans.mat %*% trans.mat %*% trans.mat %*% trans.mat
         [,1]     [,2]
[1,] 0.571741 0.428259
[2,] 0.571012 0.428988

We see each row in starting to converge to what's called the 'stationary distribution' as $n \rightarrow \infty$.  Thus, let's find the stationary distribution for Markov chain.  Another way of asking this is in the long run, what is the proportion of time the process is in the rain or no rain state?

Because our Markov chain is finite, aperiodic, and irreducible, we can use the set of equations $\pi_j = \sum_{i = 0}^{\infty} \pi_i P_{ij}$ and $ \sum_j \pi_j = 1$.

\begin{eqnarray*}
\pi_0 & = & \alpha \pi_0 + \beta \pi_1 \\
\pi_1 & = & (1-\alpha) \pi_0 + (1-\beta) \pi_1 \\
\pi_0 + \pi_1 & = & 1 \\
\end{eqnarray*}
Solving for $\pi_0$ and $\pi_1$ yields the solutions
\begin{eqnarray*}
\pi_0 & = &  \frac{\beta}{ 1 + \beta - \alpha} \\
\pi_1 & = & \frac{1-\alpha}{1 + \beta - \alpha}
\end{eqnarray*}
Plugging in $\alpha = 0.7$ and $\beta = 0.4$, we see the long-run proportion of time it will rain is $\pi_0 = 0.571$ and the long-run proportion of time it won't rain is $\pi_1 = 0.428$. This matches the limiting probabilities we calculated.

Hopefully this was helpful in understanding what a stationary distribution is and how to compute one.  Stay tuned for second-order Markov chains and a great example using the RHmm package in R.

Easy Introduction to Markov Chains in R

Markov chains are an important concept in probability and many other areas of research.  In this post, I provide the basic Markov property and then a few examples including R code to give an example of how they work.  Most of the background can be found in one of my favorite books on this topic which is called Introduction to Probability Models by Sheldon Ross.

We define a stochastic process $\{X_n, n = 0,1,2, \ldots \}$ that takes on a finite or countable number of possible values.  Let the possible values be nonnegative integers (i.e. $X_n \in \mathbb{Z}_+ $).  If $X_n = i$ then the process is said to be in state $i$ at time $n$.

Now, we define the 'Markov property' in discrete time:
\[ P_{ij} = P[X_{n+1} = j | X_n = i_n, X_{n-1} = i_{n-1}, \ldots, X_0 = i_0] = P[X_{n+1} = j | X_n = i_n] \hspace{.2in} \forall i, j \in \mathbb{Z}_+ \]
Such a stochastic process is known as Markov chain.  We call $P_{ij}$ a 1-step transition probability because we moved from time $n$ to time $n+1$.  It is a first-order Markov chain because the probability of being in state $j$ at the $(n+1)$ time point only depends on the state at the the $n$th time point.  Since proabilities are nonnegative and the process must make a transition into some state, we have these properties:

1. $P_{ij} \geq 0, \hspace{.2in} \forall i,j \geq 0$
2. $\sum_{j \geq 0} P_{ij} = 1, \hspace{.2in} i = 0, 1, 2, \ldots$

The $n$-step transition probability $P_{ij}^n$ is defined as
\[ P_{ij}^n = P[X_{n+k} = j | X_k = i] \hspace{.2in} \forall n \geq 0, \hspace{.2in} i,j \geq 0 \]
An important set of equations called the Chapman-Kolmogorov equations allow us to compute these $n$-step transition probabilities.  It states
\[ P_{ij}^{n+m} = \sum_{k=0}^{\infty} P_{ik}^n P_{kj}^m \hspace{.5in} \forall n, m \geq 0, \hspace{.2in} \forall i, j \geq 0\]
Let's try an example.  Consider two states: 0 = rain and 1 = no rain.  Define two probabilities $\alpha$ = the probability it will rain tomorrow given it rained today and $\beta$ = the probability it will rain tomorrow given it didn't rain today.  We are interested in asking what is the probability it will rain the day after tomorrow given it rained today?

First, we must set up the transition probability matrix. Translating $\alpha$ and $\beta$ to probabilities leads to $\alpha = P_{00} = P[X_{n+1} = 0 | X_n = 0]$ and $\beta = P_{01} = P[X_{n+1} = 1 | X_n = 0]$.  Because the rows must sum to 1, this leads to the transition probability matrix:
\[ P = \left[ \begin{array}{cc}
\alpha & 1-\alpha  \\
\beta & 1- \beta \\
\end{array} \right]  \]
If we assign numbers to $\alpha$ = 0.7 and $\beta$= 0.4, then we can compute the probability it will rain the day after tomorrow (aka 2 steps) given it rained today.
\[ P^2 = \left[ \begin{array}{cc}
\alpha & 1-\alpha  \\
\beta & 1- \beta \\
\end{array} \right]  \left[ \begin{array}{cc}
\alpha & 1-\alpha  \\
\beta & 1- \beta \\
\end{array} \right]  = \left[ \begin{array}{cc}
.61 & .39 \\
.52 & .48 \\
\end{array} \right]  \]
R code:
> trans.mat <- matrix(c(0.7, 0.3, 0.4, 0.6), 2,2, byrow = TRUE)
> trans.mat %*% trans.mat
     [,1] [,2]
[1,] 0.61 0.39
[2,] 0.52 0.48

These probabilities are conditional probabilities because they depend on what happened the previous day.  We can also ask: 'What is the unconditional probability of rain at time $n$? '   This can be easily found by computing the unconditional or marginal distribution of the state at time $n$:
\[  P[X_n = j] = \sum_{i = 0}^{\infty} P[X_n = j| X_0 = i] P[X_0 = i] = \sum_{i = 0}^{\infty} P_{ij}^n \alpha_i \]
where $\alpha_i = P[X_0 = i]$, $\forall i \geq 0$ and note $\sum_{i = 0}^{\infty} \alpha_i = 1$ because we must start in some state.

Back to the example. If we let $\alpha_0 = .4$ and $\alpha_1 = 0.6$ in our rain example, then the unconditional or marginal probability it will rain in two days is
\[ P[X_2 = 0] = 0.4 P_{00}^2 + 0.6 P_{10}^2 = 0.4 (0.61) + 0.6 (0.52) = 0.566 \]
Rcode:
> twostep.trans <- trans.mat %*% trans.mat 
> init.dist = matrix(c(0.4, 0.6), 1, 2)
> init.dist %*% twostep.trans
      [,1]  [,2]
[1,] 0.556 0.444

There are many more fun questions to ask, but for now I hope this was a basic introduction to discrete time Markov chains. Stay tuned for a blogpost about stationary distributions and second-order Markov chains!

Spanish Wine Tasting Party and Tapas

There are so many great wines out there that I have yet to explore.  I recently read Spain has one of the largest grape varieties in the world (wikipedia says over 600).  In an effort to learn about wines from another country, I decided to have a Spanish wine tasting party!  I read through many websites and it seemed the best way to expose your palette to wines from another country is to explore the range of available wines starting from the bubbly to the white wines and finally comparing multiple reds (one younger and one older).

Many spanish laws related to the classification of how long a wine must be aged have been put in place. The most common three are as follows:
Crianza:
     red - aged for at least 2 years (> 6 months in oak)
     white & rosé - aged for at least 1 years (> 6 months in oak)
Reserva:
     red - aged for at least 3 years (> 12 months in oak)
     white & rosé - aged for at least 2 years (> 6 months in oak)
Gran Reserva:
     red - aged for at least 5 years (> 18 months in oak)
     white & rosé - aged for at least 4 years (> 6 months in oak)

When I went to the store, I didn't have particular bottles of wine in mind. Instead, I asked for wines in the $10-15 range starting with a cava (Spain's version of champagne), one albariño, two riojas (one crianza, one reserva) and finally a sherry.

Wines
Here are the wines we picked and in order of tasting (from left to right):

1) Segura Viudas Aria Brut, Penedes region, Spain (blend of macabre, parellada, and xarel-lo grapes)
2) 2010 Vionta Albariño, Riax Baixas, Spain (albariño grape)
3) 2008 Cune Rioja Crianza, Rioja, Spain (tempranillo grape)
4) 2006 Cune Rioja Reserva, Rioja, Spain (tempranillo grape) - This one was the only over $15, but it was highly suggested and close enough to the budget.
5) Lustau Solera Reserva Pedro Ximenez San Emilio, Andalusia region, Spain (pedro ximenez grape)

Tapas
Since we were tasting spanish wines, we needed spanish tapas of course! :)  Here is the line up of the three tapas.

Manchego cheese (aged 1 year), chorizo, olives, crusty bread
This was by far the easiest to make. Just chop everything up and place on plate.  This was perfect to snack on at the beginning of the tasting with the cava.

Shrimp in garlic sauce (Gambas al Ajillo)

This recipe is very simple too.  I didn't follow a particular recipe because I used my mom's recipe for shrimp scampi, but just omitted the pasta.   The albariño was nice and light (good pairing with the shrimp). 

Ingredients: 

- 1 lb of shrimp

- 1/3 cup of fresh squeezed lemon juice

- 1/3 cup of white wine

- 1/3 cup of canola oil

- 2 garlic cloves

- handful of parsley

- salt and pepper

Recipe: 

1) Finely chop garlic and parsley. Mix garlic, parsley, lemon juice, white wine, canola oil and seasoning together. Set aside marinade (can be done early).

2) Clean and de-veine shrimp.  20 minutes before you are ready to cook the shrimp, let shrimp marinate in lemon garlic marinade in refrigerator.  Preheat oven to 350 degrees. 

3) After 20 minutes marinating, place shrimp and marinade in a baking dish in the oven for 5 minutes. After 5 minutes, turn on broil in oven and cook for additional 4-5 minutes.  Keep a close eye on the shrimp.   This recipe can also be used for shrimp scampi too! 

Patas Bravas

This was slightly more complicated than the first two tapas, but well worth the effort!  Bravas refers to the spicy or "fierce" tomato sauce with a bit of tabasco.  Moving along in the wine tasting, I thought the potatoes were hearty enough to pair with the riojas (red wines). 

Ingredients: 

- 2 lbs of baking potatoes (russet) 

- canola oil, salt, pepper

- For the sauces you can pick whatever you'd like, but I made a spicy tomato sauce and simple mayonnaise

Recipe: 

1) Peel potatoes and cut into 1/2'' chunks.  Place in large pot with cold water. Bring potatoes to a boil (~7-8 minutes) then immediately drain and cool potatoes under cold water to stop the cooking process. This is just to parboil the potatoes. 

2) After cooling the potatoes, pour oil in a large pan and place a single layer of potatoes in it to cook on medium to high heat. If the pan is not hot enough the potato will just absorb the oil instead of turning a nice golden brown.  Once potatoes seem brown on one side, flip them over and let them cook on another side until you have achieved the golden brown color.  After, remove potatoes to a paper towel to drain any excess oil and lightly coat with salt.  This may be need to be done in batches depending on the size of your pan.

3) Sauces. I made two sauces: a quick mayonnaise (or aioli) and a spin off of these two spicy tomato sauces (1) and (2). The tomato sauce began by sauteing 1/2 diced onion until golden and pureeing 1 can of diced tomatoes in a food processor.  Add tomatoes to onion mixture. Add apple cider vinegar, tabasco sauce, spanish paprika, sugar, salt, pepper. Bring to a boil and let simmer for 15 minutes.

We finished the tasting with the sherry (very sweet as noted on bottle).  Overall, I really enjoyed learning about the different grape varieties and their respective regions in Spain.  Hopefully it will come in handy if I ever make it to Spain someday!

Savory Roasted Pumpkin Seeds

With fall in full swing and Halloween around the corner, we always have lots of squash and/or pumpkins around. Whether you are carving a jack-o-lantern or making a pumpkin pie, there are leftover winter pumpkin seeds which can be turned into delicious, savory roasted snacks!

Ingredients:
- pumpkin seeds (or any winter squash seeds such as my personal favorite: acorn squash)
- olive oil
- salt, pepper
- my seasoning blend: garlic powder, paprika, chipotle powder, worchestire sauce

Recipe:
1) Rinse seeds under water and remove any squash pulp. Pat dry.

2) Add olive oil.

3) Add seasoning blend. Mix well. Spread seasoned seeds on a baking sheet.

4) Bake at 350 degrees for 15 minutes. Halfway through, move seeds around so they don't stick. Let cool and enjoy!

And in case you were curious, here is our jack-o-lantern this year: pumpkin 'pi' :)

UCSC Genome Browser: A few useful tips

Most of my research revolves around analyzing next-generation sequencing data such as whole-exome sequencing data.  As a statistician, I always appreciate finding useful bioinformatic tricks/tips from various tools that make me more efficient.  Here are three examples of using the UCSC Genome Browser that I've found helpful.

Tip #1:  How do you find a list of chromosome positions given a list of dbSNP identifiers? (Taken from the Guide to the UCSC Genome Browser FAQ by Nature)
Use the 'Variation and Repeats' group in Table Browser and the SNPs track of choice.  Just specify the genome (e.g. Human) and assembly (e.g. hg19).  For 'Region', if you want the chromosomal positions for a specific regions, click position and specify the region OR click genome and upload a list of dbSNP identifiers. Finally, choose your output format (e.g. GTF, BED) and click 'get output'.

Tip #2:  How do you lift over a set of genomic coordinates from hg18 to hg19? 
Use the Batch Coordinate Conversion (liftOver) tool in Utilities. Selected Original and New assemblies. Upload the original genomic coordinates (in a BED format) and submit.  

Tip #3: How can you extract data from the UCSC browser and use it in R? 
For this, we need to install the package rtracklayer.  Here is an example on how to extract recombination rates: 

library(rtracklayer)
my.session <- browserSession()
genome(my.session) <- "hg19"
recomb.rates <- getTable(ucscTableQuery(my.session, "recombRate"))