In this lecture, we'll talk about Biopython. We cannot do our Python course on genomics without at least mentioning Biopython. Biopython contains tons of freely available tools for bioinformatics. And the purpose of this lecture is not to teach you about how to use all of these tools. We will need an entire course only for this. But to introduce you to Biopython. The Biophython Project was founded in 1999. Since then it has growing to a large collection of modules and scripts for bioinformatics, which you can download easily from biopython.org. The goal of Biopython is to make it as easy as possible to use Python for bioinformatics. It provides parsers for lots of file formats such as FASTA, Genbank, SwissProt and many others. It also provides access to bioinformatics online services and databases, it also gives the interfaces for lots of common and not so common bioinformatics programs, various clustering methods, as well as other tools for performing common bioinformatics task. And lots of documentation on how to use these tools. But first you might want to check if Biopython is installed on your computer, if you work or study in an institution that does bioinformatics chances are that Biophyton is already on your system. To check if Biopython is installed just type import bio, the Python prompter and if Python doesn't return an error then you are lucky, Biopython is on your system. You might also want to check the version of the Biopython. So by typing print (Bio.__version__), you can find out what version of Biopython is installed on your computer. On my system it's 1.65, which is the latest at this moment. But, if yours is not this fashion, if it's an older one and you want the latest fashion, or if the input command fails, then you might want to install Biopython again. To install Biopython, make sure your follow the instructions from the download link at biopython.work. You can install Biopython for both versions of Python two and Python three. You can install it for many platforms. As I mentioned before, I will not give you the nuts and bolts of Bioython, but I just want to show you a simple example of how can use Biopythons to solve an interesting problem for bioinformatics, we are now at a point where genome sequencing is becoming a reliable tool for clinical diagnostics. When you are sick, you go to your doctor and your doctor will often request a sample from you to find out what makes you sick. Maybe is that bacteria in your urine that gives you that nasty urinary tract infection, or maybe your meningitis is caused by that virus in your cerebral spinal fluid. Whatever it is, next generation sequencing is a fast and cheap way to find out the pathogens in a culture from a patient's sample. Speaking about samples, I finally got my sequence, the sample, from a very sick patient, and that's when I have to find out what makes him so sick. So, I mapped most of the sequences is to the human genome, but some of them didn't map anywhere. So I have no idea what species they are coming from. One of those sequences is shown on this slide. Steven, do you think we can use Biopython to find out what species this sequence is coming from? >> Sure Ella, we can use Biopython to find out where the sequence comes from. We're not exactly going to use Python itself. We're going to use a program called BLAST. BLAST is one of the most successful and most widely used programs ever developed in bioinformatics. It takes a short sequence or maybe a slightly longer sequence of DNA or protein and aligns it to a database, which can be very large, of other sequences, and tells us how the two are related. So what we want to do here is we want to take our query sequence that we've taken from a patient and we want to align it to every bacterium, every virus, every sequence by scientists. We can do that with BLAST and Biopython provides a very nice handle to the BLAST program. So the server is a BLAST server that's run at NCBI. The National Center for Biotechnology Information. The main issue is everyday you can go to the server on the web, and you can paste your query in. We don't want to do that. We want to call it directly over the internet. So, BioPython gives us a handle to let us to do that. So, first we're going to use the BioPython module called BLAST. So, the method called BLAST on the BioPython module. And we're going to import a set of methods from NCBI. So, their web methods. Those are called NCBIWWW so it will say from BIO.BLAST import NCBIWWW. That gives us a large collection of NCBI build programs, that we can call over the web. First we open a file, and we read the sequence into a variable. We'll call that fasta_string. So we say fasta_string = open("myseq.fa"), and we call the read method. .read. So, on that input string. So now we've got a handle on our fasta sequence, and our DNA sequence, and we want to just call blast, and for that, we'll call NCBIWWW.qblast. So that's our blast query method that's attached to this part of this NCBIWWW package. And it's going to return just a handle that we're going to be able to look at later as a way of getting out the results. But this call right here is actually calling BLAST. So we say NCBIWWW.qblast and we give it the arguments blastn. There's actually several BLAST programs. Blastn searches nucleotides against nucleotides. You have to give it a database to search against. NT is the database. That's the non redundant nucleotide database of pretty much everything that's stored at NCBI. So it's an enormous database. But again, this is at NCBI. We don't need to have it locally. Now, we give it our fastest string. So our short little query sequence. So this little call here calls the blast program. So if you want to find out more and other ways to use BLAST, you can always type help and give it the argument NCBIWWW.qblast. It will tell you other things about qblast. Here's a bit of the output you would get. And you see that when you do that the qblast function has many, many parameters. All sorts of things, dozens of different parameters. That BLAST can take and you can provide all those to this qblast function and modify how you're running BLAST. So we were running BLAST with all default parameters just to keep things simple, but you could modify many of those. The default output format by the way is XML, which is a web type format that lets you display things nicely in a web browser. So let's look at this BLAST Record that we got, so we got this handle which is all we have so far. And by the way when you ran BLASTt, I hope you waited a minute because sometimes it may take a minute or even a few minutes for your results to comeback as you're launching a pretty sophisticated search against a large database over the web. So we're going to now import another package from NCBI called XML which that says format things a little more nicely so we're going to use NCBI XML. So, we say from Bio.Blast import NCBIXML. So, again, this is another part of Biopython and we're going to take that handle that we just got back and we're going to put it into a blast record that is going to be a little more nicely formatted so we can look at it, so we're going to say our BLAST record, that's our variable, equals NCBIXML.read. This is a special read method that handles XML. And we're going to call that on result_handle. Which is the results we just got back from our BLAST query. So this is a. And here I've just shown you a little bit about the structure of this BLAST record that we're getting back. And you can find more about this if you look more in to the Biopython. So the BLAST results themselves are going to give you lots of information, actually. You can get an illumanous output from one call to BLAST. So first you get some descriptions of what BLAST found. That is what species it hit against. And you also get a list of descriptions. For each thing each species or each sequence that a particular query matched you'll get the alignments of those sequences to the references. And then for each of these individual matches you're going to have a score. And a score is a bit score. And there's also an e value which is just shown here as E, is supposed to be a representation of the likelihood that you would get this same result by chance. So, if that value is very, very small, that means the search was not accidental, but it's a genuine match. And then you get some number representing how many alignments were found against that particular query. And on the alignments themselves, you get a number of things. Each particular alignment you have something that's called the HSP record, which stands for High Scoring Pairs. So for all the alignments of your query against the target, which might be a genome, it'll show you the score of each of those separate alignments, and it'll show you the actual matches and mismatches all capture this in this elaborate record. So, you don't want to Parse that yourself, you want to use these tools to Parse that so you can put things out in a nice human readable way. So with the Parsing BLAST Output, first let's look and see what we've got, we have our BLAST record we have captured. So one of the things in that BLAST record are the alignments. So I can just print those out to the screen by typing len(blast_record.alignments) and that will show me how long that list is. That's got 50 items in it, so I got 50 matches back. And by the way, BLAST by default limits your output to 50 matches. That's why there's 50 here. There probably were more than 50. It will give me the 50 best ones. So now I want to print out just the matches that gave me a significant hit. So I am going to set a threshold for that evalue. That's that expectation value which is a probability. And I want to make it pretty small so we'll say e_value_thresh equals 0.01 meaning I don't want to see any matches. Which don't have an e value threshold. Which is one percent or smaller. So they are unlikely to be chance matches. So what I'd like to do is loop through all the alignments of last return, and just print out the ones that have a good evalue. So here's how we'll do that. We'll get the alignments out of our BLAST record. Which is blast_record.alignments. That gets all of them. We want to loop through them, so we'll say for alignment in blast_record.alignments, okay? So now we've got the variable alignment is now bound as we go through the loop to one of our alignments. And we then want to look at the high scoring pairs at HSPs, the subrecord of alignment. So we're going to look at each of the HSPs, so that sort of serve each of the subalignments in alignment. So we do that by looping again with four HSP end alignment.hsp's. And then each of those are going to have these evalues attached. We write if hsp.expect is less than evalue threshold, so if the evalue that we got back from BLAST is smaller than the threshold we set, then we are going to print it out. So here's how we print it out. We just print some text say the word alignment and then we print the sequence. The sequence, that's called the title in the BLAST record. So, we print alignment.title. The length is alignment.length and the e value is hsp.expect. So, we're going to print out each of those with some plain text so that we can recognize in their outputs and print out the sequence, the length, the Evalue and then we print out the query sequence itself, that was the sequence that we provided to BLAST as our query. We want to print out how it matches and you'll see in the next slide how that looks. And we'll print out the target, the subject, which is hsp.subject. So this is going to print out, for all of our significant matches, it's going to print them out in a nice format. So, let's see what we got back. So, here we are for these alignments, each alignment the text ****alignments prints out first. And then you see that we have sequence which has a fancy [INAUDIBLE] number and some other information and and Evalue which is a very small number you see the first one is less than 2 times 10 to the minus 29. So very, very small number. So highly significant match. Then you see the two sequences, one above the other, with little vertical bars in between showing where they match. So the vertical bars means there's a match. If there's no vertical bar there, that would be a mismatch. But as you see, this sequence is a perfect match. And you also see on this sequence name on the first line there you can see that the name of this is Ebola virus there. The Ebola virus, isolated Ebola virus. So that's what we have. So our patient was in fact very sick and what the patient had now with the sequence we see right away that it's matching Ebola virus and all of the matches, if you scroll down through them, all the matches are to different strains of the Ebola virus. So, Ella, we found out what the patient had. If you want more help with Biopython, this is just one example. There are, literally, thousands of programs now in Biopython, all pre written for you, there's no reason to write any of them yourself, is much more useful for you to learn about what programs are available and then you can simply call them from within your Python program. So there's a Biopython tutorial available at this link I'm putting here on the screen. And then there's a FAQ I also recommend you look at.
