<jupyter_start><jupyter_text>PRACTICAL PHYLOGENETICS NOTEBOOKDr <NAME> <EMAIL> Goals of these experiments This Jupyter notebook will take you through two case studies using phylogenetic analysis to understand biological questions. This will involve aligning sequences and building maximum likelihood phylogenetic trees, followed by annotation an interpretation.We hope that this will give you - experience in analysing DNA sequence data- Understanding of the steps involved in phylogenetics- Knowledge about the compleities of the specific case studies we are usingYou will write up one of the case study analyses you perform today for your assessment. Introduction to Jupyter computational notebooks **FIRSTLY, DO NOT PANIC. EVERYTHING YOU NEED TO KNOW ABOUT COMPUTERS AND CODE WILL BE TAUGHT HERE. YOU WILL BE ABLE TO DO THIS EVEN IF YOU HAVE LITTLE EXPERIENCE WITH COMPUTERS**This class of students has mixed prior experience however, so if you have not done the bioinformatics practicals in Genetic Analysis last semester then please make yourself known and we will give you a 5 minute catch-up to make your life easier.If you are familiar with Jupyter notebooks then you can skip this section and move to "A NEW SPECIES oF APE?" below.This document that you are reading now is a **Jupyter Notebook**. It is a web browser based text editor that is also able to execute scripts ie  code. Today we are using the programming language `python`, probably the most used language in bioinformatics, but we could also run `R`, `bash` or many other languages. Scripts are found in the grey cells (see below) and have something like `In [ ]:` or `[1]` to their left in the margin..  To execute a script, click the cell below and then press SHIFT+ENTER, or instead the triangular "Run" button in the tool bar above. Try running this code below now<jupyter_code>print('Hey there, good job in running the python print command!')<jupyter_output><empty_output><jupyter_text>This code cell (containing "print('Hey...") should have executed when you pressed SHIFT-ENTER and it's output was printed below the cell ("Hey there, good job..."). All Jupyter commands run in a similar fashion Can you identify which parts of this notebook are code, which parts output, and which parts documentation like this sentence? Discuss with us if you are in doubt.1. Try editing the code below and re-running. Replace "Good job" with "Even better job"2. Instead of the `Run` button at the top you can click in the cell and press Shift-Enter to run the code. Most people find this faster, edit the cell below then give it a try:<jupyter_code>print('Hey there, good job!')<jupyter_output><empty_output><jupyter_text>ACTION:Now edit the cell above to have two print statements. On a new line type `print('Your new phrase')` and then run it. It might be easier to copy/paste and just change the pasted phrase. If it doesn't run well, you have a typo. Yes, its always a typo.**Congratulations, you have now run, copy/pasted and edited cells. Those are all the skills you will need today**This iterative edit-and-run approach is how much of modern biological data is explored and analysed. This mix of code and explanation you are seeing in this Jupyter notebook is called "literate programming" This notebook will take you through the anaysis of the two case studies found in the practical handbook. To make this notebook concise, background information is excluded from it and only available in the practical handbook, and you will need to work with both documents. For each case study you will need to run several cells just as you did above. The programs will then align and clean the DNA sequences, build a tree and annotate it. **In most cases you will only need to run the cell just as you did above. In a few cases you will be able to tweak the script just a bit following clear instructions**. Good luck! STUDY1: A NEW SPECIES OF APE? ![orangutan males](images/Bornean,_Sumatran_&_Tapanuli_orangs.jpg) _Figure 1:_ Male Bornean, Sumatran and Tapanuli orangutans, three suggested species [wikipedia](https://en.wikipedia.org/wiki/Orangutan).  The first aim of today is to investigate what phylogenetics can tell us about different species of great ape. It is, of course, complex. You might like to think how you would conceptually go about trying to get information using a phylogenetic approach. **Table 1: Latin names and common names of species in this practical.** As always, Googling is encouraged. | Name             | Common name           | Name  | Common Name || ----------------|:----------------------| --------------|:---------- || Macaca macaca | Macaque (outgroup)|  Homo sapiens sapiens | Modern humans| Hylobates lar      | Gibbon (outgroup)     |  Homo sapiens neanderthalis  | Neanderthals (extinct)|| Gorilla gorilla | Western Gorilla      | Homo sapiens denisovan | Denisovans (extinct)| Gorilla beringei | Eastern/mountain Gorilla | Pongo abelii | Sumatran orangutan     || Pan troglodytes | Chimp      |   Pongo pygmaeus | Bornean orangutan   || Pan paniscus | Bonobo      |     Pongo tapanuliensis | Tapanuli orangutan|| How much data do you have?Your working directory has some DNA sequence files in fasta format. There are a number of ways to determine the number of sequences in a file, here is a quick one-liner. Edit the cell to replace `name.fas` with the correct file `data/ape.fas`. Shift-Enter to run the cell as usual<jupyter_code>!echo "Number of sequences: "; grep -c ">" data/ape.fas<jupyter_output><empty_output><jupyter_text>It should have displayed the number of sequences in the `ape.fas` file in the `data` directory  Some python code for processing sequences, just FYI Below we use some python code to allow more complex analyses. In the next example we are going to find the number and total length of sequences using a useful code package called BioPython [1].Remember: The code below has explanations of what each section does (explanations begin with the  or 'comment' symbol) as some people are interested in seeing bioinformatics code in action. I want you to see scripting like this so you can realise it's not magic **but you do not have to know python or understand this code. Just run the cell as usual.**<jupyter_code># --------------------------------------------
# Python code to report on number of sequences 
# in a file by using BioPython
# --------------------------------------------

# import BioPython code so we can use it
from Bio.SeqIO.FastaIO import SimpleFastaParser

# set counts to zero before starting
count = 0
total_len = 0

# open the data file and give it a handle (nickname)
with open("data/ape.fas") as in_handle:
    
# for each title line add 1 to count of records, 
# and add length of sequence to a count called total_len
     for title, seq in SimpleFastaParser(in_handle):
         count += 1
         total_len += len(seq)
            
# print the results in a readable format
print("The file contained %i records with total sequence length of %i nucleotides" % (count, total_len))<jupyter_output><empty_output><jupyter_text>QUESTIONS:- Can you see which part of the above code specifies the fasta file `ape.fas`?- How could you run this on a different file in the data directory called `testseqs.fasta`? You don't need any python knowledge to answer these. The idea here is that in much of bioinformatics you can modify someone else's code to point at your data file and everything will work.  ACTION:Try it, just change the name above and re-run the cell, or ask for help if you can't quite see it. Remember that the file is within the `data` directory. If you've done it correctly (watch for typos) then the number and length of sequence reported will change. Aligning the sequencesIn order to carry out a valid analysis you have to align the DNA sequences. If you're not quite sure why, look at the images below and discuss with a demonstrator.  ![Aligned DNA sequence](./images/aligned.png "A DNA sequence alignment")_A DNA sequence alignmnet. Each character (column) can be directly compared across the different species_ ![Un-aligned DNA sequence](./images/unaligned.png "An incomplete  DNA sequence alignment")_A set of DNA sequences not completely aligned. Each character (column) cannot be directly compared across the different species as some are 'shifted' so even though they are very similar, they look enormously different when just comparing down each column (character)_ To align the sequences we will use a program called MAFFT [2]. What piece of information will we have to add to the code? Yes, the name of the input DNA sequence file to be aligned.ACTIONS:- Change the name of the file in the following code to be `ape.fas`- run the cell<jupyter_code># ---------------------------
# Align sequences using MAFFT
# ---------------------------

!mafft --auto --quiet data/ape.fas > ape.afa<jupyter_output><empty_output><jupyter_text>In the command above can you see where we name the output file? You may remember that the greater than symbol `>` is used to specify where the output of a command goes, in this case to `ape.afa`.  Did it work? Can you find the `ape.afa` file? The ".afa." extension stands for 'aligned fasta' QC the alignmentTrimal [3] quality controls the alignment, removing badly aligned regions and alignment artefacts.<jupyter_code># ------------------------------------------
# Quality control the alignment using trimal
# ------------------------------------------

!trimal -in ape.afa -out ape_trimmed.afa -gappyout -keepheader<jupyter_output><empty_output><jupyter_text>Tree reconstructionThis section will reconstruct a maximum likelihood phylogenetic tree using the sequence alignment you have produced. We will use the program FastTree [4].<jupyter_code># -------------------------
# Build tree using FastTree
# -------------------------

!FastTree -gtr -nt ape_trimmed.afa > ape.nwk<jupyter_output><empty_output><jupyter_text>The treebuilding program "FastTree" gives lots of output, you can just ignore all those details. When it is finished a file called `ape.nwk` will appear. This is a Newick tree file (.nwk) containing the tree as bracket notation text. Double click it if you want to have a look. Tree Annotation and ViewingThe tree alone (example below) is in bracket notation format (called Newick) and not very meaningful to examine.```((A,B),(C,D));```Instead we are going to display it as a graphic, and then annotate it to be easier to interpret. To do this we are going to use a tree graphics program called ToyTree [5]. QUESTION:What treefile (.nwk) has just been written by the build tree cell above?ACTION:Take your newick treefile name and enter it into the cell below to replace "tree.nwk"<jupyter_code># -----------------------------------
# Drawing the phylogeny using ToyTree
# -----------------------------------
# import the code so we can use it here
import toytree       # a tree plotting library
import toyplot       # a general plotting library
# import numpy as np   # a numerical library, give it the shorthand 'np'

# read the newick format tree file, give it the name 'newick'
newick = "ape.nwk" # change this to point at your .nwk treefile
tre = toytree.tree(newick, tree_format=1)

tre.draw();<jupyter_output><empty_output><jupyter_text>If you see a graphic image of a phylogenetic tree, congratulations! If not please ask for a little help, its probably a quick fix for a demonstrator.NOW ROOT THE TREEYour tree will probably look very odd because it isn't yet rooted correctly. Use the next cell to root it by entering "Macaca" (Macacque) instead of "outgroup"<jupyter_code># ----------------
# Root and re-draw
# ----------------
# root and draw the tree
rtre = tre.root(wildcard="Macaca") # specify the outgroup taxon
rtre.draw(height=600, tip_labels_align=True); # draw the tree<jupyter_output><empty_output><jupyter_text>You should now have a tree that reveals a lot about the relationships betwen these species. Remember that there is a table with the species name earlier in this document. It will be easier to interpret though when you put it into a report if you annotate and colour it by taxon. NOW ANNOTATE THE TREE:Although you now have 'the answer' it is not so easy to study this tree. You will need to compare the divergences between the two species of orangutan and compare those to the divergences between the two species of chimpanzee. In this simple tree its not too hard, but in general phylogeneticists label and colour to maintain focus on the correct comaprisons. You are now going to use the script below to colour in the tips by their species identity. Run the cell and examine the tree<jupyter_code># -----------------------------
# Colouring the tree label text
# -----------------------------

# set list of colours depending on the taxon label text
# numbers like "#5384a3" are colour 'hex' codes, cyan in
# this case (google "hex colour codes" for other options)

colorlist = ["blue" if "Pan_paniscus" in tip
             else "darkblue" if "Pan_troglodytes" in tip 
             else "red" if "Pongo_abelii" in tip 
             else "brown" if "Pongo_pygmaeus" in tip
             else "#5384a3" for tip in rtre.get_tip_labels()]

# draw the tree
canvas = rtre.draw(
    width=600,  # set dimensions of the figure
    height=600,
    scalebar=True,  # scale bar of divergence levels
    tip_labels_align=True,
    tip_labels=True,
    tip_labels_colors=colorlist,
    node_labels=None,
    node_sizes=[0 if i else 8 for i in rtre.get_node_values(None, 1, 0)],
    node_markers="s", # use "o" for circles instead of squares
    node_colors=toytree.colors[0],
)<jupyter_output><empty_output><jupyter_text>You now have all the skills to edit this script and change colours. Pick some ones you like and rerun. IMPORTANT, SAVE YOUR FILEMake sure that you save and take away a copy of you tree file image in a format suitable to insert into your final report. Run the cell below and then find the file in your working directory and save it somewhere accesible.<jupyter_code># --------------------------------
# Save the tree as a graphics file
# --------------------------------

# import code to draw graphics files
import toyplot.pdf
import toyplot.svg
import toyplot.html

# draw graphics files
# toyplot.svg.render(canvas[0], "ape.svg")
toyplot.pdf.render(canvas[0], "ape.pdf")
toyplot.html.render(canvas[0], "ape.html")<jupyter_output><empty_output><jupyter_text>SPECIES DIFFERENTIATION Looking at the tree, it would seem the two *Pan* clades are as distant from each other as the two *Pongo* clades The scale along the bottom is genetic distances from 0 to 1, so 0.06 would be 6%. Find the common ancestor node of each genus. What is the distance between the common ancestor of both Pan species and the tree tips? What is that value for Pongo? Do you think it is very different?When this analysis was first done orangutans were all classified as the same species, whereas there were two Pan species. Do you think that classification accurately reflects evolutionary history? Does the classification in this way have any negatives for conservation planning? Discuss in your group.Soon after these analyses on single mitochondrial genes, exactly as you have done here, orangutans were divided into two species.  Pause You have just loaded a data file, aligned it, quality controled the alignment, constructed a maximum likelihood phyogenetic tree, and created an annotated figure of the phylogeny. Well done! Oh yes, some of you may also have modified a bioinformatics python script to characterise a DNA sequence file.  This was quite a lot of work to do the first time, trying to understand how to pass a specific data file through the analytical satges to create a phylogeny. Fortunately, as you learned above, doing it again on a different data file just requires a simple change, ie specifying a different file. Below you can quite rapidly analyse a "big ape" data set, containing a lot more sequences, by rerunning the same commands with different data. It should not take long. In the second case study (below) we are going to swap from apes to HIV, but again it should be rapid because the commands will be very similar. A big data analysis of great apesOne very useful aspect of using code to carry out analyses is that once you have written it, and it works, its very little effort to re-run it again on any number or any size of other data sets.Here I have collected from GenBank whole mitochondrial genomes (about 16,000 nucleotides per genome) from a lot of great apes including humas, neanderthals, and species of gorillas in addition to the species you have just analysed. The file is large but we can just run the same code again. If you want to find out how much data you have, you could insert a cell and paste in the code to quantify sequences (from the ape example) and run it for the big_ape dataset. This is optional. For efficiency reasons I've compressed the code below a little, but its the same as you have just run.This big analysis gives you the opportunity to decide whether the similarity of divergence between groups that you have just observed is true more widely. When you have produced their big trees you should discuss what the divergence levels might mean with your demonstrators. If you want to find out how much data you are analysing you can copy and paste the "how much data do you have?" cell from above, and run it here (on the new fasta file). But that is optional. Run this cell to align, trim the alignment, and then build a tree. It might take a few minutes to complete. When there is a number not an asterisk in the left margin then it is complete.<jupyter_code># Align
!mafft --auto --quiet data/big_ape.fas > big_ape.afa
print("\nThe sequence alignment has finished")
# Trim
!trimal -in big_ape.afa -out big_ape_trimmed.afa -gappyout -keepheader
print("The alignment trimming has finished")
# Tree build
print("The phylogenetic tree construction has started\n")
!FastTree -gtr -nt big_ape_trimmed.afa > big_ape.nwk
print("\nThe phylogenetic analysis has finished")<jupyter_output><empty_output><jupyter_text>Expect this to take a couple of minutes. Remember if there is an asterisk in the top left "`In [*]:`" then it is still working, when it is a number it is finished. If this completed without errors then you can just run the cell below and see the output tree. You may want to adjust colours and re-run a few times. If you had an error, see if you can spot what went wrong, but seek assistance if not. **Or maybe the cell below. Merge them somehow**<jupyter_code># import the tree plotting software
import toytree
import toyplot
import toyplot.svg
import toyplot.html
import toyplot.pdf

# read the tree file
newick = "big_ape.nwk"
bigtre = toytree.tree(newick, tree_format=1)

# specify the rooted tree with outgroup as Macaque, Hylobates lar
bigrtre = bigtre.root(wildcard="Hylobates")

# colour the tree, change these options until you are happy with the design
colorlist = ["#d6557c" if "Pan" in tip # pink
             else "blue" if "Gorilla_b" in tip #
             else "#4169E1" if "Gorilla_g" in tip #royalblue
             else "#008000" if "Homo_sapiens_sapiens" in tip #green
             else "#32CD32" if "neanderthal" in tip #limegreen
             else "#006400" if "denisovan" in tip #darkgreen
             else "red" if "Pongo_abelii" in tip #darkgreen
             else "orange" if "Pongo_pygmaeus" in tip #darkgreen
             else "brown" if "Pongo_tapanuliensis" in tip #darkgreen
             else "#5384a3" for tip in rtre.get_tip_labels()] # cyan


# create the tree using these colours and some other standard options
canvas = bigrtre.draw(
    width=700,
    height=1200,
    scalebar=True,
    node_labels=None,
    node_sizes=[0 if i else 8 for i in bigrtre.get_node_values(None, 1, 0)],
    node_markers="s",
    node_colors=toytree.colors[0],
    tip_labels_align=True,
    tip_labels_colors=colorlist
);
# --------------------------------
# Save the tree as a graphics file
# --------------------------------
#you can specify the name if you wish

toyplot.pdf.render(canvas[0], "bigape.pdf")
toyplot.html.render(canvas[0], "bigape.html")<jupyter_output><empty_output><jupyter_text>Check the files list to make sure that the tree has saved and looks how you intended. You will need this tree for your report. You can alter and recreate the tree under the same, or a different name. It's a bit more complicated than that... There is no level of percentage sequence divergence that determines what is a species. So what can this tree really tell us? Well, amount of sequence divergence from closest relatives gives us a good initial picture of the system, and can give a framework on which to carry out more sophisticated analyses. Evolutionary biologists have sequenced the entire nuclear genomes of these species (you have analysed the much smaller mitochondrial genomes). This allows them to look at estinates of coalescent date, evidence for hybridization and many other biological aspects.  Well DoneYou are now finished with case study 1, the apes. Case study two, the origins of HIV, will be much faster now you have experience.Please feel free to take a short break here. STUDY2: WHAT ARE THE ORIGINS OF HIV? ![HIV](images/HIV.png) **QUESTION: Has HIV (human immunodeficiency virus) coevolved with humans or does it have a recent zoonotic origin?**How can you test this? You now have all the skills required. We are going to repeat some of the work described by <NAME> (2011) in their paper "Origins of HIV and the AIDS Pandemic". Their figure 4 is very informative, and Zimmer and Emlen redraw it in their (3rd edition) Figure 8.12. Today you are going to reanalyse the HIV and SIV sequence data from great apes to produce a simlar figure and answer the question set above regarding zoonotic transfer.  Sequence data I have prepared aligned fasta files for you containing SIV and HIV sequence data from the *env* gene (Google it). They are called `SIVHIVspecies_ENV.afa`It would be useful for you to know how many sequences were in each file, you calculated this earlier today for different files. Paste in the code to the cell below and describe the sequences. Tree reconstruction<jupyter_code>#HIV
print("\nThe analysis has begun, this will take a few minutes, please be patient")

# # Alignment has already been done, so this is commented out
# !mafft --auto --quiet data/SIVHIVspecies_ENV.fasta > SIVHIVspecies_ENV.afa
# print("\nThe sequence alignment has finished")
# Trim
!trimal -in data/SIVHIVspecies_ENV.afa -out SIVHIVspecies_ENV_trimmed.afa -gappyout -keepheader
print("The alignment trimming has finished")
# Tree build
print("The phylogenetic tree construction has started\n")
!FastTree -gtr -nt SIVHIVspecies_ENV_trimmed.afa > SIVHIVspecies_ENV.nwk
print("\nThe phylogenetic analysis has finished")
# import the code so we can use it later, just run this cell
import toytree       # a tree plotting library
import toyplot       # a general plotting library
import toyplot.pdf
import toyplot.svg
newick = "SIVHIVspecies_ENV.nwk"
hivtre = toytree.tree(newick, tree_format=1)
hivrtre = hivtre.root(wildcard="schwein")

# change these options until you are happy with the design
colorlist = ["green" if "P.t.schweinfurthii" in tip 
             else "blue" if "P.t." in tip 
             else "blue" if "G.g" in tip
             else "blue" if "H.sapiens" in tip 
             else "blue" for tip in hivrtre.get_tip_labels()]

canvas = hivrtre.draw(
    node_labels=None,
    width=600,
    height=600,
    node_sizes=[0 if i else 8 for i in hivrtre.get_node_values(None, 1, 0)],
    node_markers="s",
    node_colors=toytree.colors[0],
    tip_labels_align=True,
    tip_labels_colors=colorlist);<jupyter_output><empty_output><jupyter_text>You can use modify the code to colour it in more appropriately if you wish.The name contains a lot of metadata. The letter after H.sapiens indicates the HIV1 group to which that sequence belongs. Then the country from which it was isolated. Then the year. Then some identifiers and accession number. So `H.sapiens.O.CM.96.LA51YBF35.KU168294` is HIV1 group O, from Cameroon, isolated in 1996, with sample code LA51YBF35 and genbank accession number KU168294. It may help you to talk about HIV1 group O, or group P, or group N in your explanations and in your annotations of the tree. Make sure you then download the annotated tree and save it somewhere safe for your report.<jupyter_code># --------------------------------
# Save the tree as a graphics file
# --------------------------------
# remember to save these files and take them away

toyplot.pdf.render(canvas[0], "hiv.pdf")
toyplot.html.render(canvas[0], "hiv.html")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08395898926293695, 'fraction_numerical': 0.006216194397352062, 'mean_word_length': 3.765339488363147, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 72, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '160323', 'license': 'Unlicense', 'n_tokens_mistral': 6781, 'n_tokens_neox': 6406, 'n_words': 3520, 'path': 'planB/backup.ipynb', 'repo_name': 'davelunt/phylo_practical', 'size': '24797'}