Teaching

As we are spending our time in quarantine, I thought we might try and consider reconstructing some phylogenetic trees that will tell us something about the evolutionary history of the coronavirus 2019-nCov. In particular, we will try and replicate (not completely, but the general idea) the panels b and c from Figure 2 in this paper:Identification_of_a_novel_coronavirus . The idea is to create and compare two phylogenetic trees created from the whole genome sequence of the coronavirus and from its spike protein (you might read more about spike proteins in this article, back from the previous SARS epidemic).

Our goal will be to write a script that creates these two plots using the following steps:

 

1. Downloading the sequences from Entrez using Bio.Entrez (23 pts)
2. Creating multiple alignments of both DNA (whole genome) and protein (spike protein) sequences (3 pts)
3. Creating phylogenetic trees for both protein and genome alignments using the UPGMA and neighbor joining algorithms (3pts)
4. Visualizing the trees (1 pt)

Now let us describe all these steps in more detail

1. We need to download whole genome sequences and spike protein sequences from the Genbank Entrez databases (“nucleotide” and “protein” databases, respectively), together with some other viruses. We will only look at 6 species:

• 2019-nCov virus (current coronavirus)
• SARS – another coronavirus back from 2002
• bat coronavirus
• MERS – middle eastern respiratory syndrome coronavirus
• influenza A –  flu virus, to see how far are coronaviruses from the flu
• hepatitis A virus – less related virus

You can use the Bio.Entrez interface. An example of how this works can be seen in a jupyter notebook here (please remember to use your e-mail and not mine in your code…). We will need to do this for all 6 species. You should provide a single script that does that for all species, using different query terms for each genome and protein. You will probably need some trial and error work here, but I want to see just the end result – a script that works.

2. Here we will want to run the ClustalW algorithm on the downloaded sequences. You will need to have the clustalw executable installed (it is the package clustalw in ubuntu). Your script should use the Bio.Align.Applications.Clustalw interface. in the result you should get two multiple sequence alignments (for genomes and for spike proteins).

3. Here you should use the two multiple alignments to create distance matrices and from these distance matrices to create phylogenetic trees (both with upgma and nj methods).

4. You should create images of the 4 trees (upgma or nj x protein or genome) using the draw function.

A full solution is a script (in a .py file sent as a real attachment -not some usos attachment, not an ipynb, not a link) together with the images for part 4 (in png format). All of this should be sent to me (bartek[at]mimuw.edu.pl) with the tag [WBO] in the e-mail subject by April 6th. The code should be reasonably easy to read (with comments and the code divided into sensible functions or classes. Not a huge spaghetti of code.)

 

Today we will hold the lecture remotely via google meet. I will send you the link via the chat server. If you cannot access the chat, please let me know by e-mail.

Polish version below…

Today we will discuss the simple methods for phylogenetic tree reconstruction based on a distance matrix. The slides are in two parts: My part regarding simple progressive heuristic methods: wyk4  and another part on ML methods by P. Górecki : gorecki-ml-tree

During the lab, we will try to build trees from sequences in practice. First, we should get ourselves acquainted with the Bio.Phylo module  tutorial  and wiki.

1. Let us first consider two lists of sequences containing paralogues and orthologues of phenylalanine hydroxylase gene as well as histone H2B paralogues and orthologues. Let us assume, that we have them already aligned (we will cover the multiple alignment process next week). See the files  Human_PAH_orthologues,  Human_PAH_paralogues and  Human_H2BFS_paralogues ,  Human_PAH_orthologues-91. Let  us use the read() function from the Bio.AlignIO module to read these files.

2.Calculate distance matrices for these groups of sequences using the  Bio.Phylo.TreeConstruction.DistanceCalculator  (using  BLOSUM62 substitution matrix, separately for each file).

3. Create the trees using Bio.Phylo.TreeConstruction.DistanceTreeConstructor both for UPGMA and nj (neighbor-joining)

4. look at the resulting trees using  draw_ascii() and draw() methods

5. Write out the resulting trees into newick and phyloxml formats.Compare resulting files.

Wersja polska:

Dziś rozmawiamy o prostych metodach kostrukcji drzew na podstawie macierzy odległości. Slajdy na dziś składają się z mojej prezentacji nt. metod heurystycznych:  wyk4 i slajdów P. Góreckiego nt. metod ML: gorecki-ml-tree

Na laboratorium będziemy konstruować drzewa w praktyce.

Warto zapoznać się z dokumentacją modułu Bio.Phylo w tutorialu jak i na stronie wiki.

1. Na początek ustalmy listy sekwencji. Rozważmy sekwencje paralogiczne i ortologiczne ludzkiej hydroksylazy fenyloalaniny.  Załóżmy, że mamy już sekwencje aminokwasowe uliniowione globalnie.  Ze względu na czasochłonność procesu uliniowienia, użyjemy plików zuliniowieniami  Human_PAH_orthologues i Human_PAH_paralogues oraz Human_H2BFS_paralogues i (jeśli ktoś liczy lokalnie – uwaga duży) plik Human_PAH_orthologues-91. Wczytaj te pliki przy pomocy metod z modułu Bio.AlignIO

2.Wylicz macierze odległości dla tych grup sekwencji przy pomocy klasy Bio.Phylo.TreeConstruction.DistanceCalculator  (dla macierzy BLOSUM62, osobno dla paralogów i osobno dla ortologów genów PAH – to  zajmie chwilę).

3. Stwórz drzewa filogenetyczne na podstawie macierzy przy pomocy klasy Bio.Phylo.TreeConstruction.DistanceTreeConstructor zarówno metodą UPGMA – hierarchiczną jak i nj (neighbor joining)

4. Wyświetl uzyskane drzewa przy pomocy metody draw_ascii() i draw()

5. Zapisz uzyskane drzewa do formatów newick i phyloxml. obejrzyj wyniki.

Dziś zajmujemy się uliniawianiem par sekwencji – wyk3  (english version below)

Continue reading “WBO 3 – Pairwise sequence alignment”

Today, we will discuss basic mathematical models of sequence evolution. (polska wersja poniżej)

The slides for the lecture can be found here: wyk2

Again, you can use our jupyter server. There are some tutorials on jupyter usage  here. Otherwise, you will need a python interpreter with biopython installed.

Assignments for today:

1. Take a look at the Bio.SubsMat.MatrixInfo module, where PAM and BLOSUM matrices are defined.

2. Implement a simple, discrete Markov model, that allows for simulation of mutations in a sequence over discrete time steps, according to substitution matrices defined using Bio.SubsMat. Try to implement Kimura and Felsenstein models.

3. Take the test_fasta file from last week and see how many steps of random mutations without selection does it take to change one of the sequences to another without selection pressure.

4. Implement a simple Markov process with selection, i.e. after making a random mutation that brings you closer (according to Hamming distance) to the target sequence, then accept it. Reject it otherwise. Compare run times of this implementation with the previous one (without selection).

 

Dziś będziemy mówić o modelach ewolucji sekwencji.

Slajdy do wykładu są tu: wyk2.

Podczas labów mogą Państwo korzystać z serwera Jupyter, który ma zainstalowany pakiet biopython. (kilka porad jak używać jupyter’a tu ). Zadania do wykonania na dziś:

1. Obejrzyj moduł Bio.SubsMat.MatrixInfo, gdzie zdefiniowane są modele PAM i BLOSUM

2. Zaimplementuj prosty łańcuch Markowa, który pozwoli na symulowanie mutacji w czasie dyskretnym, zgodnie z macierzami przejścia zdefiniowanymi wykorzystując moduł SubsMat (substitution table). Rozważ model analogiczny do modelu Kimury i Felsensteina

3. Rozważ sekwencje z pliku fasta z poprzednich zajęć. Zbadaj ile mutacji losowych potrzeba aby Łańcuch Markowa “przeszedł” pomiędzy dwiema sekwencjami bez selekcji.

4. Zaimplementuj prosty proces Markowa z selekcją – jeśli mutacja zbliża nas do “celu” w sensie odległości Hamminga, to ją akceptuj, w przeciwnym razie wracaj do punktu wyjścia.