Full credit will only be given to correct answers with a clear explanation of how they are obtained. Use additional paper as necessary.

1. Look up the helpfile of the S-Plus command "sample".
2. Generate a random nucleotide sequence of length 1000 using the independent model (i.e. the tetrahedral die model) with probabilities P(A) = 0.3, P(C) = 0.2, P(G) = 0.25, P(T) = 0.25. Name your sequence "myseq".
   Solution: myseq<-sample(c("A", "C", "G", "T"),size=1000,replace=T,prob=c(0.2,0.3,0.25,0.25))
3. Let N_A (respectively N_C, N_G, and N_T) be the random variable representing the number of A's (respectively C's, G's, and T's) in such a random sequence as generated in problem #2.
   
   (a). What kind of probability distribution does N_A (respectively N_C, N_G, and N_T) have? What are the parameters in that distribution?
   
   (b). What is the expected value, of N_A (respectively N_C, N_G, and N_T)?
   
   

   x	Probability Distribution	Parameters	E(x)
   NA
   NC
   NG
   NT

4. Figure out the observed counts of the bases in myseq you obtained in problem #2. Fill in the table below.
   
   (Hint: In Exercise3A, you have created an S-Plus object called "mysequence". Find out what happens when you do mysequence[mysequence=="A"]. Do the same for "C", "G", and "T". This might give you some idea.)
   
   

   Base	Expected Count (from #3)	Observed Count
   A
   C
   G
   T

   Show the S-Plus command(s) that you used.
   

   cta<-length(myseq[myseq=="A"]) ctc<-length(myseq[myseq=="C"]) ctg<-length(myseq[myseq=="G"]) ctt<-length(myseq[myseq=="T"]) cta ctc ctg ctt
5. Generate a random nucleotide sequence of length 1000 according to the Markov model with the transition probability matrix P given in Problem #1 in Exercise 2A. Assume that the sequence has already attained its stationary distribution at the beginning. Repeat problem 4 for this sequence.
   
   Hint: You'll need to use "for"-loops and some "if-else" statements. The "dimnames" function for matrices may also be useful. Check out some of the introductory S-Plus documents.)
   
   

   Base	Expected Count	Observed Count
   A
   C
   G
   T

   
   
   The expected values were obtained as follows:
   Let I^(j)_A = 1 if the j^th base is A; 0 if the j^th base is not A.
   Then,
   N_A = I⁽¹⁾_A + I⁽²⁾_A + ... + I⁽¹⁰⁰⁰⁾_A
   Since the expected value of a sum is equal to the sum of the expected values, we have
   E[N_A] = sum(j=1 to 1000)[E(I^(j)_A)] = 0.2691466 + 0.2691466 + ... + 0.2691466 = 269.1466
   The expected counts for C, G, and T are obtained similarly.

   Show the S-Plus commands that you used.
   

   dimnames(P)<-list(c("A", "C", "G", "T"),c("A", "C", "G", "T")) DNA<-c("A", "C", "G", "T") markovseq<-rep(" ", 1000) markovseq[1]<-sample(DNA,size=1,replace=T,prob=c(0.2691466,0.2319475,0.1816193,0.3172867)) for (i in 2:1000) {markovseq[i]<-sample(DNA,size=1,replace=T,prob=P[markovseq[i-1], ]) if (i==1000) { cta<-(length(markovseq[markovseq=="A"])) ctc<-(length(markovseq[markovseq=="C"])) ctg<-(length(markovseq[markovseq=="G"])) ctt<-(length(markovseq[markovseq=="T"])) print (cta) print (ctc) print (ctg) print (ctt) } }