 #### 联系方式

• QQ：99515681
• 邮箱：99515681@qq.com
• 工作时间：8:00-21:00
• 微信：codinghelp #### 您当前位置：首页 >> Database作业Database作业

###### 日期：2022-06-05 08:38

See Canvas for due dates

In the first part of this assignment, we use a Hidden Markov Model to model secondary

structure in protein sequences and implement a couple of algorithms we saw in lectures.

In the second part, we simulate sequences down a tree according to the Jukes-Cantor

model then use distance methods to try to reconstruct the tree.

Write your code in Python and present your code embedded in a report in a Jupyter

Notebook. Make sure you test your code thoroughly and write clear, commented code

that others can understand.

Submit two files to Canvas: the .ipynb and .html both showing code and results by 10pm

on the due date.

There are 30 marks in total for this assessment.

1. [14 marks total] Suppose we wish to estimate basic secondary structure in protein

(amino acid) sequences. The model we consider is a simplistic rendition of the

model discussed in S C. Schmidler et al. (2004) Bayesian Segmentation of Protein

Secondary Structure, doi:10.1089/10665270050081496

We assume that at each point of the sequence, the residue is associated with one

of three secondary structures: α-helix, β-strand and loops which we label H, S

and T , respectively. To simplify the problem, we classify the amino acids as either

hydrophobic, hydrophilic or neutral (B, I or N , respectively) so a sequence can be

represented by this 3-letter alphabet.

In a α-helix, the residues are 15% neutral, 20% hydrophobic and 65% hydrophilic.

In a β-strand, they are 30%, 60%, 10% and in a loop they are 70%, 15%, 15%.

Assume that all secondary structures have geometrically distributed length with

α-helices having mean 15 residues, β-strands having a mean of 8 residues and loops

a mean of 6 residues. A β-strand is followed by an α-helix 40% of the time and a

loop 60% of the time. An α-helix is followed by a β-strand 30% of the time and a

loop 70% of the time and a loop is equally likely to be followed by a strand or a

helix. At the start of a sequence, any structure is equally likely.

When writing code below, work in natural logarithms throughout to make your

calculations robust to numerical error.

(a) [3 marks] Sketch a diagram of the HMM (a hand-drawn and scanned picture

is fine). In your diagram, show only state nodes and transitions. Show the

emission probabilities using a separate table.

Note that the transition probabilities of states to themselves (e.g., aHH) are

not given. Derive them by noticing that you are given the expected lengths

of α-strands, β-helices and loops, and that if a quantity L is geometrically

distributed with parameter p then the expected value of L is E[L] = 1/p.

Make sure you use the correct parametrisation of the geometric distribution

1

(noting that you can’t have a secondary structure of length 0) and remember

that

l akl = 1 for any state k.

(b) [3 marks] Write a method to simulate state and symbol sequences of arbitrary

length from the HMM. Your method should take sequence length, and model

parameters (a and e) as arguments. Simulate and print out a state and symbol

sequence of length 200.

(c) [3 mark] Write a method to calculate the natural logarithm of the joint prob-

ability P (x, pi). Your method should take x, pi, and model parameters as

arguments.

Use your method to calculate P (x, pi) for pi and x given below and for the

sequences you simulated in Q1b.

pi = S,S,H,H,H,T,T,S,S,S,H,T,T,H,H,H,S,S,S,S,S,S

x = B,I,B,B,N,I,N,B,N,I,N,B,I,N,B,I,I,N,B,B,N,N

(d) [5 marks] Implement the forward algorithm for HMMs to calculate the natural

logarithm of the probability P (x). Your method should take x as an argument.

Note that we don’t model the end state here.

Use your method to calculate log(P (x)) for pi and x given in Q1c and for the

sequences you simulated in Q1b.

How does P (x) compare to P (x, pi) for the examples you calculated? Does

2

2. [16 marks total] In this question you will write a method that simulates random

trees, simulates sequences using a mutation process on these trees, calculate a

distance matrix from the simulated sequences and then, using existing code, re-

construct the tree from this distance matrix.

(a) [5 marks] Write a method that simulates trees according to the Yule model

(described below) with takes as input the number of leaves, n, and the branch-

ing parameter, λ. Use the provided Python classes.

The Yule model is a branching process that suggests a method of constructing

trees with n leaves. From each leaf, start a lineage going back in time. Each

lineage coalesces with others at rate λ. When there k lineages, the total rate

of coalescence in the tree is kλ. Thus, we can generate a Yule tree with n

leaves as follows:

Set k = n, t = 0.

Make n leaf nodes with time t and labeled from 1 to n. This is the set of

available nodes.

While k > 1, iterate:

Generate a time tk ～ Exp (kλ). Set t = t+ tk.

Make a new node, m, with height t and choose two nodes, i and j,

uniformly at random from the set of available nodes. Make i and j

the child nodes of m.

Add m to the set of available nodes and remove i and j from this set.

Set k = k-1.

Simulate 1000 trees with λ = 0.5 and n = 10 and check that the mean height

of the trees (that is, the time of the root node) agrees with the theoretical

mean of 3.86.

Use the provided plot tree method to include a picture of a simulated tree

with 10 leaves and λ = 0.5 in your report. To embed the plot in your report,

include in the first cell of your notebook the command %matplotlib inline

(b) [5 marks] The Jukes-Cantor model of DNA sequence evolution is simple:

each site mutates at rate μ and when a mutation occurs, a new base is chosen

uniformly at random from the four possible bases, {A,C,G, T}. If we ignore

mutations from base X to base X, the mutation rate is 34μ. All sites mutate

independently of each other. A sequence that has evolved over time according

to the Jukes-Cantor model has each base equally likely to occur at each site.

The method mutate is provided to simulate the mutation process.

Write a method to simulate sequences down a simulated tree according to the

Jukes-Cantor model.

Your method should take a tree with n leaves, sequence length L, and a

mutation rate μ. It should return either a matrix of sequences corresponding

to nodes in the tree or the tree with sequences stored at the nodes.

3

Your method should generate a uniform random sequence of length L at the

root node and recursively mutate it down the branches of the tree, using the

node heights to calculate branch length.

In your report, include a simulated tree with n = 10 and λ = 0.5 and a set

of sequences of length L = 20 and mutation parameter μ = 0.5 simulated on

that tree.

(c) [3 marks] Write a method to calculate the Jukes-Cantor distance matrix, d,

from a set of sequences, where dij is the distance between the ith and the

jth sequences. Recall that the Jukes-Cantor distance for sequences x and y

is defined by

dxy = ?3

4

log (1? 4fxy

3

)

where fxy is the fraction of differing sites between x and y. Since we will be

dealing with short sequences, use the following definition of fxy so that the

distances are well-defined:

fxy = min

(

Dxy

L

, 0.75? 1

L

)

where Dxy is the number of differing sites between x and y and L is the length

of x.

Include a simulated set of sequences of length L = 20 from the tree leaves and

corresponding distance matrix in your report for a tree with n = 10, λ = 0.5

and mutation parameter μ = 0.5.

(d) [3 marks] Now simulate a tree with n = 10 and λ = 0.5 and on that tree,

simulate three sets of sequences with lengths L = 20, L = 50 and L = 200,

respectively, with fixed μ = 0.1. For each simulated set of sequences, calculate

the distance matrix and print it out.

Then reconstruct the tree using the provided compute upgma tree method.

Use the plot tree method to include a plot of the original tree and a plot of

the reconstructed tree for each distance matrix.

Comment on the quality of the reconstructions and the effect that increasing

the sequence length has on the accuracy of the reconstruction.

【上一篇】：到头了
【下一篇】：没有了