Month: March 2020

The lecture is again online. I will send the link before. Labs are again after the lecture. English version below:

Dziś rozmawialiśmy o algorytmach przybliżonego wyszukiwania sekwencji w bazach danych. Slajdy są tu: wyk6-blast.

Na ćwiczeniach spróbujemy zająć się wykorzystaniem programu blast: 0. Wczytaj plik w formacie FastQ microbial_reads.fastq (jest już na naszym serwerze jupyter w katalogu WBO, nie trzeba go tam ładować), przy pomocy SeqIO.parse(). Są to odczyty z mikrobiomu jelitowego myszy.

Continue reading “WBO 6 – Searching for distant relatives- BLAST”

English version below:

Dzisiaj mówiliśmy o metodach uliniowienia wielu sekwencji. Slajdy są tu: wyk5

Na zajęciach naszym celem jest zapoznanie się z metodami uliniowienia wielu sekwencji w pakiecie Biopython oraz metodami uliniowienia progresywnego. Warto przeczytać odpowiedni rozdział tutorialu biopythona

0. Tutaj warto, aby każdy z Państwa korzystających z serwera jupyter, założył sobie na serwerze Jupyter własny folder, żeby pliki nam się nie mieszały. Proponuję numer indeksu, albo inny identyfikator jednoznaczny.

1.Weźmy na początek naszą rodzinę sekwencji białek histonowych. Jest dostępna na naszym serwerze jupyter w pliku histones.fa Na początek potrzebujemy wczytać je z pliku (SeqIO.parse) przetłumaczyć je na sekwencje białkowe (metoda translate()) i zapisać do pliku fasta (SeqIO.write) wraz z odpowiednimi identyfikatorami i usuniętymi kodonami stopu (*, która powstaje w wyniku translate)

2. Wykorzystajmy program clustalw (Bio.Align.Applications.ClustalwCommandline) do wykonania mulituliniowienia tych sekwencji na naszym serwerze plik wykonywalny jest zainstalowany w katalogu /usr/bin. Po wykonaniu zadania wczytaj uliniowienie z pliku .aln (AlignIO.read format “clustal”) oraz drzewo filogenetyczne z pliku .dnd (Phylo.read format=”newick”) (osoby pracujące na mac’ach mogą zainstalować sobie clustal Omega)

3. Wykorzystaj teraz program muscle (Bio.Align.Applications.MuscleCommandline) do wykonania multiuliniowienia tych samych sekwencji. Wczytaj je z pliku w formacie fasta

4. Stwórz drzewo filogenetyczne na podstawie uzyskanego uliniowienia metodą neighbor joining jak na poprzednich zajęciach. Porównaj to drzewo do drzewa otrzymanego z programu clustalw

5. Powtórz ćwiczenia 2-4 dla większego zbioru sekwencji (np. Human_PAH_homologs). Czy tutaj widoczne są większe różnice między metodami?

6(*). Praca domowa za 2 punkty. Napisz prostą implementację uliniowienia progresywnego, która wykorzystuje globalne uliniowienie z modułu Bio.pairwise2. To będzie wymagać implementacji funkcji call-back dla substytucji profili i niewielkich zmian w samym module pairwise2

English version

0. If you are using our Jupyter server, please make sure that you work in a separate folder that identifies you by name.

1.Firstly, let us take the histones.fa file. We need to parse it (using Bio.SeqIO.parse method), then translate each of them into protein alphabet (using the  .translate() method) and save to another file using SeqIO.write( remember about removing the * symbols that represent translations of a stop codon)

2. Let us then use the clustalw (Bio.Align.Applications.ClustalwCommandline) to create a multiple alignment of these sequences (the clustalw executable should be in  /usr/bin on the server). After you are able to run it, you need to read the  .aln file with the alignment (AlignIO.read(file, format=”clustal”)) and the guide tree from the  .dnd file (Phylo.read format=”newick”) (students using apple macs can install clustal Omega)

3. Now, for comparison, you can use the  muscle  program (Bio.Align.Applications.MuscleCommandline) to create a multiple alignment of the same sequences. It should give you the alignment in a  fasta format.

4. You can use neighbor joining algorithm we discussed last week to create the tree based on the muscle alignment and compare it visually to the guide tree from clustalw. Do you see differences?

5. If you still have time, you can repeat assignments 2-4 for a larger set of sequences (e.g. Human_PAH_homologs). Are the differences between trees more visible here?

6(*). Homework assignment(2pts). Please implement a simple progressive alignment method that uses Bio.pairwise2 module. It is not as simple as it seems, as you need to implement profile substitution matrix using a callback function and possibly some changes in bio.pairwise2 module.

 

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.