Module 7 Exercises
Overview
In todays workshop we will be looking to:
- Explore various tools for predicting protein secondary structure, function, and localization from primary sequence.
We will also see how to:
- Use R to perform repetative structure prediction of proteins
Distance-based trees
There are many ways to build a phylogenetic tree, and each one has its pros and cons. Here, we will cover two of these. UPGMA and NJ Trees can be built from any sequences, whether DNA, RNA, or protein. What you are tracking is the history of mutations. Today we will build trees using an available data set for mammalian DNA.
Download the data file GGH.dna (https://bit.ly/3n0ql1n). Upload this data to your RStudio server account or store this data in your RStudio/R working directory on your personal computer. Go ahead and click on this file. You should see this open a new text tab. You will notice that this file contains aligned gamma-glutamyl hydrolase genes from 40 mammalian species. (If you are unfamiliar with the Latin names of mammals, they are listed here: http://evolution.gs.washington.edu/gs570/2016/data/species)
The file starts with two numbers, which are the dimensions of the data matrix. The first number is the number of sequences (lines in this case) and the second number is the number of aligned base pairs.
To perform our exercise today, we will need to load a number of libraries using the following commands:
library(ape)
library(phangorn)
library(seqinr)This loads the programs that we will be using and sets your working directory to the relevant folder. You may get an error message regarding seqinr and ape but you can ignore it. However, always make note of such conflicts in your own work. Which function, from which library, you use can have large (often unintentional) effects on your analysis.
Now execute:
mammals<-read.dna("GGH.dna",format="sequential")
mammals_phyDat<-phyDat(mammals,type="DNA",levels=NULL)The first line of code reads the data file into the program and gives it a handy name (mammals). The second line of code converts the data file into a format that the program can manipulate. In other words, you have renamed your data file twice, first calling it mammals and then calling it mammals_phyDat.
Next type:
dna_dist<-dist.ml(mammals_phyDat)This converts the original data matrix to a matrix of distances. To see what this looks like, type:
dna_dist Dasypus Echinops Pongo Ictidomys Loxodonta Canis Macaca Myotis Felis
Echinops 0.292734977
Pongo 0.202036851 0.242208060
Ictidomys 0.202149558 0.263548411 0.138162071
Loxodonta 0.215647557 0.218454552 0.197325452 0.211277630
Canis 0.211008250 0.296077362 0.181132821 0.194594452 0.223778310
Macaca 0.207406803 0.244066797 0.032021939 0.143236086 0.204196778 0.186473298
Myotis 0.215647557 0.301572723 0.206213919 0.213304262 0.220346567 0.165917747 0.206213919
Felis 0.212008981 0.281664453 0.191186094 0.187130807 0.228885559 0.120256988 0.188485821 0.169134214
Macropus 0.419932722 0.453185956 0.368837072 0.415875515 0.420353141 0.418281907 0.368837072 0.438553411 0.436256333
Callithrix 0.226710654 0.262651822 0.073439912 0.162381089 0.221060294 0.178847830 0.087180754 0.223690833 0.218165946
Erinaceus 0.251809411 0.359679810 0.221378793 0.260640871 0.280574587 0.224666470 0.216474143 0.257765811 0.230613525
Sorex 0.333407422 0.368640004 0.249168106 0.283237814 0.299767731 0.270843305 0.267186991 0.277647726 0.266591015
Choloepus 0.149477029 0.287969177 0.200656474 0.200761054 0.219559031 0.226876372 0.209915851 0.243495807 0.220372588
Nomascus 0.241055104 0.272937079 0.038786551 0.166639617 0.213362474 0.198120341 0.053403335 0.220180678 0.212511900
Ornithorhy 0.421638950 0.435200874 0.405268478 0.421662818 0.398871008 0.411424464 0.397138318 0.442985923 0.434415977
Vicugna 0.197134441 0.282224679 0.168984723 0.170339775 0.199277374 0.147096010 0.174931235 0.166029054 0.147674567
Monodelphi 0.434221089 0.455741150 0.409138573 0.431316060 0.435111135 0.435147372 0.410966728 0.452393500 0.442753265
Sarcophilu 0.439186648 0.465738550 0.403795432 0.427582077 0.437028046 0.438943664 0.396601132 0.460133667 0.435147372
Cavia 0.271115497 0.337071172 0.244971614 0.242430634 0.290423016 0.265530904 0.247874375 0.298615778 0.290851733
Equus 0.211008250 0.292098572 0.164034316 0.191873361 0.213923640 0.154986383 0.162735052 0.181777369 0.167226648
Rattus 0.310797600 0.339180878 0.281639745 0.264048248 0.326918846 0.304885407 0.286231596 0.315622237 0.310832438
Ochotona 0.199741145 0.281553781 0.179568939 0.165696750 0.233043121 0.185501732 0.189483162 0.222119174 0.225428994
Gorilla 0.197361664 0.238885006 0.015666366 0.147350977 0.197860943 0.175080008 0.043308650 0.204367069 0.191348420
Procavia 0.220349941 0.254397775 0.198684723 0.219273726 0.146071694 0.220320839 0.210627651 0.234113452 0.216828736
Dipodomys 0.268532897 0.339180878 0.225196750 0.193262907 0.271586000 0.238084895 0.222362547 0.262549746 0.250668310
Macropus Callithrix Erinaceus Sorex Choloepus Nomascus Ornithorhy Vicugna Monodelphi
Echinops
Pongo
Ictidomys
Loxodonta
Canis
Macaca
Myotis
Felis
Macropus
Callithrix 0.401529692
Erinaceus 0.449219974 0.253200679
Sorex 0.414170522 0.250015485 0.309926366
Choloepus 0.393233479 0.195993787 0.304098831 0.279588900
Nomascus 0.376971642 0.106431711 0.243831353 0.245615701 0.234360890
Ornithorhy 0.389557025 0.428594827 0.410685176 0.454476738 0.378008466 0.410509271
Vicugna 0.420507428 0.209047292 0.212473106 0.201999700 0.215761554 0.190175009 0.417631396
Monodelphi 0.169836211 0.414159738 0.463483228 0.489049362 0.394184204 0.442933510 0.384913477 0.430239115
Sarcophilu 0.147117966 0.418940630 0.449781982 0.486381143 0.399220436 0.430495007 0.367411837 0.419859610 0.167688424
Cavia 0.443133578 0.251942761 0.310074958 0.318117930 0.294651540 0.265177965 0.458181812 0.260196454 0.467992285
Equus 0.413850580 0.184089552 0.236288174 0.269012920 0.223074083 0.184812944 0.436566483 0.164555576 0.427623797
Rattus 0.483404611 0.285303396 0.349567297 0.391715300 0.323985342 0.298088496 0.478320244 0.297455741 0.543038658
Ochotona 0.422875068 0.223375833 0.242361867 0.218932806 0.213218563 0.210580442 0.432436549 0.193485840 0.435253381
Gorilla 0.381424780 0.076471750 0.221996548 0.231121234 0.209251611 0.049197962 0.401192382 0.167505432 0.413077411
Procavia 0.410958273 0.191729771 0.315824865 0.322038969 0.203855211 0.196318791 0.439455596 0.179268368 0.477404741
Dipodomys 0.459460039 0.232240625 0.304098831 0.313325234 0.263325118 0.243069327 0.458244696 0.234884261 0.498599118
Sarcophilu Cavia Equus Rattus Ochotona Gorilla Procavia Dipodomys Otolemur
Echinops
Pongo
Ictidomys
Loxodonta
Canis
Macaca
Myotis
Felis
Macropus
Callithrix
Erinaceus
Sorex
Choloepus
Nomascus
Ornithorhy
Vicugna
Monodelphi
Sarcophilu
Cavia 0.438943664
Equus 0.425801901 0.261077650
Rattus 0.498218001 0.340337764 0.303313079
Ochotona 0.452086695 0.296578521 0.230505006 0.305267966
Gorilla 0.411230890 0.240815488 0.163973836 0.273082442 0.185501732
Procavia 0.439731761 0.294613610 0.213039228 0.350505384 0.209220315 0.186759698
Dipodomys 0.470064475 0.285790989 0.255621009 0.300940932 0.236022390 0.225078444 0.259518082
Microcebus Mus Pteropus Tursiops Tarsius Mustela Tupaia Ailuropoda Pan
Echinops
Pongo
Ictidomys
Loxodonta
Canis
Macaca
Myotis
Felis
Macropus
Callithrix
Erinaceus
Sorex
Choloepus
Nomascus
Ornithorhy
Vicugna
Monodelphi
Sarcophilu
Cavia
Equus
Rattus
Ochotona
Gorilla
Procavia
Dipodomys
Bos Homo Oryctolagu
Echinops
Pongo
Ictidomys
Loxodonta
Canis
Macaca
Myotis
Felis
Macropus
Callithrix
Erinaceus
Sorex
Choloepus
Nomascus
Ornithorhy
Vicugna
Monodelphi
Sarcophilu
Cavia
Equus
Rattus
Ochotona
Gorilla
Procavia
Dipodomys
[ reached getOption("max.print") -- omitted 14 rows ]Which writes the matrix to the screen. The matrix is large so it wraps around and is a little hard to read. Scroll to the top of the output so you can see the beginning of the matrix. The first column of numbers is labeled Dasypus, the second is Echinops, the third is Pongo, etc. The species names are in the same order down the left side of the matrix. This shows you that the distance between Dasypus (armadillo) and Echinops (tenrec) is 0.292605364, that between Dasypus (armadillo) and Pongo (orangutan) is 0.202628039, etc.
Question
What sequences are most distant from Dasypus? What do those animals have in common?
Now use the distance matrix to produce a tree. We will compare two trees, one produced by the Unweighted Pair Group Method with Averaging (UPGMA), and one produced by the Neighbor Joining (NJ) method.
mammals_UPGMA<-upgma(dna_dist)To see this tree:
plot(mammals_UPGMA)
To save the tree, do it two ways.
First, just save as a pdf, using the Export menu in your Plots menu.
Second, write it to a file:
write.tree(mammals_UPGMA,file="UPGMA_tree")You can download this file to your desktop via the Files panel (same window that you can visualize your plots). Files -> More -> Export. You can use this tree file to view your tree in tree viewing software such as FigTree (http://tree.bio.ed.ac.uk/software/figtree/).
Now produce a neighbor-joining tree:
mammals_NJ<-NJ(dna_dist)
plot(mammals_NJ)
Save the plot as a pdf and
write.tree(mammals_NJ,file="NJ_tree")Now, let’s compare the two trees.
Question
What is different? What is the same?
Sometimes, it is useful to redefine the root of a tree (or reroot a tree). Let’s take a look at how we can do this in R.
First, plot your tree as we’ve done before. Except, now we will use an additional command from the ape package called nodelabels(). This will label all tree nodes with their numerical node ID in the tree. We can then use this node ID to define the root of the tree at this node.
plot(mammals_NJ)
nodelabels()
Now, let’s reroot the tree around node 77 (the clade with marsupials). It is often best to export these labeled trees to a PDF to see them more clearly.
marsupial_rooted_NJ <- root(mammals_NJ, node = 77)
plot(marsupial_rooted_NJ)
See if you can do the same for your UPGMA tree!
Maximum likelihood trees, maximum parsimony trees, and bootstrapping
We just built phylogenetic trees using the clustering algorithms UPGMA and neighbor joining. Now, we will use the same data set, GGH.dna, to build a tree using maximum likelihood and maximum parsimony. We will also do a bootstrap analysis.
Maximum likelihood
To get started, let’s do a preliminary neighbor-joining tree. This is exactly as we did above. We need to make this tree as a starting point for a maximum likelihood optimization of the tree.
dna_dist<-dist.ml(mammals_phyDat)
mammals_NJ<-NJ(dna_dist)Now we can calculate the likelihood of the NJ tree.
fit<-pml(mammals_NJ, mammals_phyDat)
print(fit)
loglikelihood: -15999.05
unconstrained loglikelihood: -6087.627
Rate matrix:
a c g t
a 0 1 1 1
c 1 0 1 1
g 1 1 0 1
t 1 1 1 0
Base frequencies:
0.25 0.25 0.25 0.25 The output will give you a loglikelihood. Write down this number somewhere.
Now see if you can improve the likelihood of the NJ tree by adjusting the tree topology and branch lengths. This, in effect, is the process of maximum likelihood optimization. Let’s do this first using a Jukes-Cantor model of evolution.
fitJC<-optim.pml(fit,model="JC",rearrangement="stochastic")optimize edge weights: -15999.05 --> -15909.49
optimize edge weights: -15909.49 --> -15909.49
optimize topology: -15909.49 --> -15881.9
optimize topology: -15881.9 --> -15876.25
optimize topology: -15876.25 --> -15866.6
8
optimize edge weights: -15866.6 --> -15866.6
optimize topology: -15866.6 --> -15854.15
optimize topology: -15854.15 --> -15848.86
optimize topology: -15848.86 --> -15841.21
5
optimize edge weights: -15841.21 --> -15841.2
optimize topology: -15841.2 --> -15837.22
optimize topology: -15837.22 --> -15833.67
optimize topology: -15833.67 --> -15830.68
5
optimize edge weights: -15830.68 --> -15830.68
optimize topology: -15830.68 --> -15821.98
optimize topology: -15821.98 --> -15818.52
optimize topology: -15818.52 --> -15815.55
3
optimize edge weights: -15815.55 --> -15815.55
optimize topology: -15815.55 --> -15815.55
0
[1] "Ratchet iteration 1 , best pscore so far: -15815.5473627609"
[1] "Ratchet iteration 2 , best pscore so far: -15815.5472449143"
[1] "Ratchet iteration 3 , best pscore so far: -15815.546991391"
[1] "Ratchet iteration 4 , best pscore so far: -15815.5468588054"
[1] "Ratchet iteration 5 , best pscore so far: -15815.5467412791"
[1] "Ratchet iteration 6 , best pscore so far: -15815.5465696573"
[1] "Ratchet iteration 7 , best pscore so far: -15815.5463980959"
[1] "Ratchet iteration 8 , best pscore so far: -15815.54622666"
[1] "Ratchet iteration 9 , best pscore so far: -15815.54622666"
[1] "Ratchet iteration 10 , best pscore so far: -15815.5460670799"
[1] "Ratchet iteration 11 , best pscore so far: -15815.5457601723"
[1] "Ratchet iteration 12 , best pscore so far: -15815.5455363778"
[1] "Ratchet iteration 13 , best pscore so far: -15815.545311109"
[1] "Ratchet iteration 14 , best pscore so far: -15815.5451938252"
[1] "Ratchet iteration 15 , best pscore so far: -15815.5450498017"
[1] "Ratchet iteration 16 , best pscore so far: -15815.5447466052"
[1] "Ratchet iteration 17 , best pscore so far: -15815.5446028576"
[1] "Ratchet iteration 18 , best pscore so far: -15815.5446028576"
[1] "Ratchet iteration 19 , best pscore so far: -15815.5446028576"
[1] "Ratchet iteration 20 , best pscore so far: -15815.544441575"
[1] "Ratchet iteration 21 , best pscore so far: -15815.5442711599"Is this analysis faster or slower than the NJ analysis you did earlier?
Get the log likelihood by typing:
logLik(fitJC)'log Lik.' -15815.54 (df=77)Question
What is this log likelihood and is it higher or lower than the log likelihood of the NJ analysis?
Look at the tree topology by typing:
plot(fitJC)
How does this tree compare to the neighbor joining tree you got earlier?
Save the tree as a pdf if you want. Write the tree to a file by typing:
write.tree(fitJC$tree,file="fitJC")Now change the model of evolution to the General Time Reversible model (GTR). Type:
fitGTR<-optim.pml(fit,model="GTR",rearrangement="stochastic")optimize edge weights: -15999.05 --> -15909.49
optimize base frequencies: -15909.49 --> -15899.83
optimize rate matrix: -15899.83 --> -15473.94
optimize edge weights: -15473.94 --> -15471
optimize topology: -15471 --> -15439.66
optimize topology: -15439.66 --> -15432.77
optimize topology: -15432.77 --> -15419.57
9
optimize base frequencies: -15419.57 --> -15409.19
optimize rate matrix: -15409.19 --> -15405.07
optimize edge weights: -15405.07 --> -15405.04
optimize topology: -15405.04 --> -15390.4
optimize topology: -15390.4 --> -15381.54
optimize topology: -15381.54 --> -15380.04
5
optimize base frequencies: -15380.04 --> -15378.7
optimize rate matrix: -15378.7 --> -15378.06
optimize edge weights: -15378.06 --> -15378.04
optimize topology: -15378.04 --> -15377.69
optimize topology: -15377.69 --> -15376.78
optimize topology: -15376.78 --> -15365.14
3
optimize base frequencies: -15365.14 --> -15364.75
optimize rate matrix: -15364.75 --> -15364.57
optimize edge weights: -15364.57 --> -15364.57
optimize topology: -15364.57 --> -15359.14
optimize topology: -15359.14 --> -15355.32
optimize topology: -15355.32 --> -15353.9
3
optimize base frequencies: -15353.9 --> -15353.85
optimize rate matrix: -15353.85 --> -15353.81
optimize edge weights: -15353.81 --> -15353.81
optimize topology: -15353.81 --> -15353.81
0
[1] "Ratchet iteration 1 , best pscore so far: -15353.8080846898"
[1] "Ratchet iteration 2 , best pscore so far: -15353.8071873916"
[1] "Ratchet iteration 3 , best pscore so far: -15344.747840337"
[1] "Ratchet iteration 4 , best pscore so far: -15344.7469727731"
[1] "Ratchet iteration 5 , best pscore so far: -15344.7469727731"
[1] "Ratchet iteration 6 , best pscore so far: -15344.7469727731"
[1] "Ratchet iteration 7 , best pscore so far: -15344.7433412151"
[1] "Ratchet iteration 8 , best pscore so far: -15344.7433412151"
[1] "Ratchet iteration 9 , best pscore so far: -15344.7433412151"
[1] "Ratchet iteration 10 , best pscore so far: -15344.7433412151"
[1] "Ratchet iteration 11 , best pscore so far: -15344.7433412151"
[1] "Ratchet iteration 12 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 13 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 14 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 15 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 16 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 17 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 18 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 19 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 20 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 21 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 22 , best pscore so far: -15344.7433405812"
[1] "Ratchet iteration 23 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 24 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 25 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 26 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 27 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 28 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 29 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 30 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 31 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 32 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 33 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 34 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 35 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 36 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 37 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 38 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 39 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 40 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 41 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 42 , best pscore so far: -15344.7433405709"
[1] "Ratchet iteration 43 , best pscore so far: -15344.7433405709"
optimize base frequencies: -15344.74 --> -15344.71
optimize rate matrix: -15344.71 --> -15344.69
optimize edge weights: -15344.69 --> -15344.69
optimize topology: -15344.69 --> -15344.69
0
optimize base frequencies: -15344.69 --> -15344.69
optimize rate matrix: -15344.69 --> -15344.69
optimize edge weights: -15344.69 --> -15344.69
optimize base frequencies: -15344.69 --> -15344.69
optimize rate matrix: -15344.69 --> -15344.69
optimize edge weights: -15344.69 --> -15344.69 When the analysis finishes, get the loglikelihood of the tree by typing:
logLik(fitGTR)'log Lik.' -15344.69 (df=85)Is this log likelihood higher or lower than the log likelihood of the NJ analysis? How does it compare to the log likelihood of the JC analysis?
Look at the tree topology by typing:
plot(fitGTR)
Is this tree different from the one using the JC model?
Save the tree as a pdf if you want. Write the tree to a file by typing:
write.tree(fitGTR$tree,file="fitGTR")Maximum parsimony
Another method for constructing trees is maximum parsimony. In this case, the MP approach looks to minimize the number of changes in character states between branches. This contrasts to ML that looks to find the tree with the highest probability of each character state in their respective positions in the tree. Let’s see how we can implement MP in R using phangorn and ape.
Like before, it is best to start with an existing model to optimize. We can use the NJ tree we made above (mammals_NJ). Now, let’s determine the parsimony score for this tree.
parsimony(mammals_NJ, mammals_phyDat)[1] 3165Similar to ML, we can perform an optimization of the tree except optimizing for parsimony.
fitNJ_MP <- optim.parsimony(mammals_NJ, mammals_phyDat, perturbation="stochastic")Final p-score 3119 after 5 nni operations Let’s see what we’ve got!
plot(fitNJ_MP)
Boostrapping
Finally, let’s assess the support for a tree using the bootstrap. Here, we’ll examine the ML tree we produced, above. Run the bootstrap analysis using the command below. We are setting the number of iterations to 10. Typically, you would set this value much higher to produce more accurate branch support estimates. If you are using the RStudio server, be aware this process will take some time to execute when you choose a much larger bootstrap value (bs).
bs<-bootstrap.pml(fitGTR, bs=10, optNni=TRUE, multicore=TRUE, control=pml.control(trace=0))Then summarize all the bootstrap trees with:
plotBS((fitGTR$tree),bs,p=50,type="p")
Save this tree as a pdf.
Which branches are poorly supported? Which ones are strongly supported? How do the poorly supported branches relate to branches that vary among the NJ and ML trees?
Write the set of bootstrap trees to a file using write.tree:
write.tree(bs,file="bootstrap.GTR.tree")