Note: _ denotes subscript, so x_i is x_i.

1. Let x_1,...,x_n, x_{n+1}, x_{n+2} be independent, with Pr(x_r = i) = f_i for i = 1...k. Let n_i be the number of x_1...x_n that are equal to i, and assume a Dirichlet(alpha_1,...,alpha_k) prior for (f_1,...,f_k).

a) Show that Pr(x_{n+1} = i | x_1,...,x_n) = (alpha_i + n_i) / (n+sum_i alpha_i)

b) Find a similar expression for Pr(x_{n+1} = i and x_{n+2} = j | x_1,...,x_n), in both the case i not equal to j, and the case i =j.

2. Let g_1,...,g_n be the (observed) genotypes of n individuals, let h_1,...,h_n be the (unobserved) haplotype pairs for the same individuals, and let f=(f_1,...f_N) be the population frequencies of the N possible haplotypes. Assume a Dirichlet prior for f.

The following algorithm describes a Gibbs sampling strategy for sampling from the conditional distribution Pr(h_1,...,h_n | g_1,...,g_n): Start with some initial guess for h_1,...,h_n, and then iterate the following steps:

a) Choose i uniformly at random from i = 1...n.

b) Sample h_i from Pr(h_i | h_{-i}, g_1, ..., g_n) where h_{-i} denotes the set of current haplotype guesses, excluding individual i. Using results from 1 above, describe in detail how the step b) above would be performed.

3. Try applying Clark's algorithm (by hand) to estimate the haplotypes of each individual given the following genotype data on 5 SNPs in 6 individuals (each column is a SNP locus, each pair of rows is an individual). Does the algorithm give a unique solution?

00000
11111
00011
00010
00000
00010
01110
10010
11110
11110
00001
00011

4.

• a) Use the PHASE program (on abacus.biostat.washington.edu) to estimate the haplotype frequencies and haplotypes for each individual in problem 3 above. Compare your estimates with the estimates you get from Clark's algorithm.
• b) use PHASE on the data for a single individual at 2 loci that you considered in Q3 in homework 3. What do you find?

Instructions for using PHASE:

• 1. log-in to abacus
• 2. type "source ~statgen/.statgen.cshrc" (without the "").
• 3. copy the input files for the two datasets from my home directory by typing "cp ~ms1024/hw3*.phase.inp ." (note that the final period follows a space).
• 4. run PHASE by typing, for example "PHASE hw3a.phase.inp hw3a.phase.out"
• 5. examine the output in the file "hw3a.phase.out_pairs" (which lists the most probable haplotype pairs for each individual, along with their estimated probabilities), and the file "hw3a.phase.out_freqs" which gives estimates of the population frequencies of each haplotype.
• 6. Here are the full instructions for the PHASE software, if you need them, or are interested in understanding the input format etc. (If you have problems, the first thing to do is to email me!)