2 - Sequences_and_BLAST - v9.ipynb

{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Working with Sequence files using Biopython\n","\n","The code boxes below will increase in complexity as we go on. Comments in the code begin with #s. Read these if you want help understanding the code.\n","\n","In the first code box below, the first line \"turns on\" the SeqIO function of Biopython, a package (set of tools) built for biology and biochemistry! You can learn more at https://biopython.org/\n","\n","The second line uses SeqIO (think: sequence input and output) to read a fasta file and stores the information as a list of records.\n","\n","<font color=blue><b>STEP 1:</b></font> The command below won't run correctly (you can try it!) unless you enter in the file to read. <b> Between the quotes, where it says <<\\<your file here>>>, change it to : files/uniprot-dtxr.fasta.</b>\n","    \n","Then shift+enter to run the code. Be sure to wait for a number to appear in the brackets to the left of the code box. When the asterisk appears, that code is actively running!\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["from Bio import SeqIO # imports the SeqIO function from Biopython\n","\n","records = list(SeqIO.parse(\"<<<your file here>>>\", \"fasta\"))     # reads the fasta file into a list of records \n","print(\"Finished storing the FASTA file in the list called \\\"records\\\".\")"]},{"cell_type":"markdown","metadata":{},"source":["While this <em>did</em> read the fasta file into a list of records, it isn't obvious how this helps us. Let's see how we can access different information about the sequence records stored in the list, records.\n","\n","<font color=blue><b>STEP 2:</b></font> The first line below will tell us how many sequences were in the fasta file, and then the information stored in the first item. Note that python numbering starts with the number 0. Go ahead and run the code box below."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["print(\"There are %i sequences in your file.\\n\" % len(records))   # prints the number of sequences, that is, the length of the list, named records\n","\n","print(records[0]) #prints the information in the first record"]},{"cell_type":"markdown","metadata":{},"source":["You should be able to see that first record includes several things: an ID, Name, Description, and a Sequence.\n","\n","We can access each of those items specifically\n","\n","<font color=blue><b>STEP 3:</b></font>  Run the code below to print only the sequence of the first record."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["print(records[0].seq)"]},{"cell_type":"markdown","metadata":{},"source":["The next code box shows us how to list the first 10 ids. Then it lists the first record id and its sequence.\n","\n","<font color=blue><b>STEP 4:</b></font> Run the code box below."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["print(\"The first 10 sequence record ids are:\\n\")\n","for i in range(10):                                              # this creates a variable i and counts to 10\n","    print(records[i].id)                                         # prints the id for record i\n","    \n","print(\"\\nThe record: %s has a sequence of: %s\\n\" % (records[0].id, records[0].seq))  # prints the record id and its sequence!"]},{"cell_type":"markdown","metadata":{},"source":["***\n","Great! The code above finds the first record (recall in Python we start counting at zero), so records[0].id gets the identification of the first record. \n","\n","<font color=blue><b>STEP 5:</b></font> Edit the last print statement in the code above to give the id and sequence of the 100th record.\n","\n","We can also look at the last record id. You could put in the number of records (less one), but -1 is easier. The -1 starts counting from the end of the records list and you don't need to know how many records you have.\n","\n","<font color=blue><b>STEP 6:</b></font> Edit the last print statement in the code above to give the id and sequence of last record.\n","\n","***"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":[]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 7:</b></font> Before you move on, check your knowledge by using the empty code box above to enter some code to answer these questions.\n","\n","    1. What is the sequence id for the 3rd record in the file?\n","    2. What is the sequence id for the 250th record in the file?\n","    3. What is the sequence id for the penultimate record in the file?\n","   \n","<font color=blue><b>STEP 8:</b></font> It might be quite clear that simply locating records by their number isn't always the best way to interact with the data. Perhaps we would like to search by name, ID, or description. Let's examine the comments below to see how we can do this. You will need to edit the line that sets the variable search_term. In the code box below, change the text <b>\\<<<search term here\\>>></b> to one of the search terms below. Then run the code box. If the output list is too long, comment the print statement inside the loop.\n","\n","    Search term ideas: DTXR, IDER, MNTR, dtxr, ider, mntr\n","\n","    Where to search ideas: id, description, name"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import re  # imports the re (Regular Expressions) functions\n","\n","search_term = \"<<<search term here>>>\"   # in between the quotes we can add a search term.\n","                                         \n","\n","counter = 0   # a variable to keep track of how many times we find the term\n","\n","for item in records:                       # this creates a loop to iterate over all the records               \n","    if re.search(search_term, item.description, re.IGNORECASE):    # Here we use re.search to find the search term in the item description, returns true if yes\n","        print(item.description)                     # If the above is true we print the item description...\n","        counter = counter + 1              # and increment our counter\n","    else: continue                         # if the search term wasn't found we do nothing, just continue on\n","\n","        \n","# the next two lines print a summary for us - either no hits or how many hits.\n","        \n","if not counter: print(\"We didn't find any results using the search term: %s\" % search_term)\n","else: print(\"\\nWe found %i records using the search term: %s\" % (counter, search_term))"]},{"cell_type":"markdown","metadata":{},"source":["***\n","\n","<font color=blue><b>STEP 9:</b></font> Edit the code above to search the ID rather than the description. Make sure to change it in both places!\n","\n","***\n","\n","Using Biopython's SeqIO function allows us to store lots of data in a way that is rapidly accessible.\n","\n","As you have seen above, the IDs are a little long and redundant with the name. The code below simplifies the record ID and writes a new, simpler file.\n","\n","<font color=blue><b>STEP 10:</b></font> Read through the code below and run it."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["original_file = \"files/uniprot-dtxr.fasta\"         # original file path\n","simple_file = \"files/uniprot-dtxr_simple.fasta\"   # new file path\n","\n","with open(original_file) as original, open(simple_file, 'w') as simple:  #opens the file to read and one to write\n","    records = SeqIO.parse(original, 'fasta')   # here we use the Biopython SeqIO again to parse the file into records\n","\n","    for item in records:                       # iterate through each item in the list called records\n","        tmp_name = str.split(item.id, \"|\")[1]      # sets the variable tmp_name to the second item generated by splitting at the | character \n","        item.id = tmp_name                         # now this sets the id of that item to the name\n","        item.description = \"\"\n","        SeqIO.write(item, simple, 'fasta')  # writes out the item information to the new, corrected fasta file"]},{"cell_type":"markdown","metadata":{},"source":["There is no output given for the above code. Let's consider how we changed this file by looking at the simplfied file.\n","\n","<font color=blue><b>STEP 11:</b></font> You can find the new file in the file browser in the left of your screen and double click to open it."]},{"cell_type":"markdown","metadata":{},"source":["***\n","## Big Data Strategies - Filtering and Reducing Redundancy in Datasets\n","\n","Now that we have a simpler fasta file, we still have a lot of sequences to deal with. It is often best to use fast and easy methods to pare down a dataset before using more computationally intensive ones. \n","\n","We will start by filtering out small fragments and large proteins, then move on to removing redundancy.\n","***"]},{"cell_type":"markdown","metadata":{},"source":["## Filtering by length\n","\n","Let's visualize the distribution of sequence lengths in our dataset. The code below will generate a histogram of the lengths of sequences. Charlie Weiss will give a thorough demonstration of this tool on Thursday!\n","\n","<font color=blue><b>STEP 12:</b></font> Run the code below to generate a histogram of sequence lengths.\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import matplotlib.pyplot as plt\n","\n","records = list(SeqIO.parse(\"files/uniprot-dtxr_simple.fasta\", \"fasta\"))\n","\n","lengths = []\n","\n","for i in records:\n","    #print(i.seq)\n","    lengths.append(len(i.seq))\n","\n","lengths.sort()\n","#print(lengths)\n","\n","plt.hist(lengths, bins = 100)\n","plt.xlabel('seq length')\n","plt.ylabel('count')\n","plt.show()"]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 13:</b></font> Answer the question:\n","\n","    1. Using the histogram, what are the most common lengths of proteins in this dataset?\n","\n","***\n","\n","We will further reduce the complexity of this dataset by removing sequences that are much larger and much smaller than the most common lengths shown. Those variables can be determined using a histogram and easily changed for different datasets. Using the histogram we can identify appropriate values for the small_len and large_len variables.\n","\n","<font color=blue><b>STEP 14:</b></font>   Run the code below to remove short and long sequences and generate a new fasta file called \"files/uniprot-dtxr_simple_trim.fasta\".\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["from Bio import SeqIO\n","\n","original_file = \"files/uniprot-dtxr_simple.fasta\"\n","trimmed_file = \"files/uniprot-dtxr_simple_trim.fasta\"\n","\n","small_len = 115\n","large_len = 252\n","\n","with open(original_file) as original, open(trimmed_file, 'w') as trimmed:\n","    records = SeqIO.parse(original_file, 'fasta')\n","\n","    for record in records:\n","        if len(record.seq) > small_len and len(record.seq) < large_len: \n","            #print(len(record.seq))\n","            SeqIO.write(record, trimmed, 'fasta')\n","            \n","records2 = list(SeqIO.parse(\"files/uniprot-dtxr_simple_trim.fasta\", \"fasta\"))\n","print(\"We now have %i sequences in our fasta file.\" % len(records2))"]},{"cell_type":"markdown","metadata":{},"source":["The above procedure removed a few thousand sequences. This is a good start, but still a lot of sequences to visualize.\n","***\n","## Reducing sequence redundancy using CD-HIT\n","\n","We will use the program CD-HIT to remove sequence within a given sequence similarity to another sequence. For example, the flag \"-c 0.4\" given below will keep only sequences with less than 40% sequence identity. The -n flag selects the word size used in processing. Less sequence identity requires the use of smaller words in comparing the sequences. From the CD-HIT manual:\n","\n","    Choice'of'word'size:\n","        -n 5 for thresholds 0.7 ~ 1.0\n","        -n 4 for thresholds 0.6 ~ 0.7\n","        -n 3 for thresholds 0.5 ~ 0.6\n","        -n 2 for thresholds 0.4 ~ 0.5\n","\n","CD-HIT was installed when you launched the Binder. This is one of the benefits of using Binder.org\n","\n","For this exercise we will use a cutoff of 40%. <b>In practice you might select higher values, but this will increase runtimes for the following steps and the visualization in the next Jupyter Notebook.</b>\n","\n","<font color=blue><b>STEP 15:</b></font> Let's run the code below. This runs CD-HIT using the simplified and trimmed fasta as input (-i) and creates a new fasta output (-o) named uniprot-dtxr_40.fasta that will be much smaller. The bang (!) at the beginning tells the notebook this is not a Python command, but rather a program to be run from the linux environment. <b>Note running CD-HIT takes a few minutes..</b>\n","We do this in two steps as this produces a better selection because of the way CD-HIT is implemented. It is a practical, not theoretical, issue. CD-HIT manual: http://www.bioinformatics.org/cd-hit/cd-hit-user-guide.pdf"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["!cd-hit -i files/uniprot-dtxr_simple_trim.fasta -o files/uniprot-dtxr_70.fasta -c 0.7 -n 5\n","!cd-hit -i files/uniprot-dtxr_70.fasta -o files/uniprot-dtxr_40.fasta -c 0.4 -n 2\n"]},{"cell_type":"markdown","metadata":{},"source":["CD-HIT has reduced over 30,000 sequences to less than 1500 at 40% sequence identity! This is one strategy for dealing with very large datasets - they can be reduced. However, reduction needs to be done in a logical and reproducible manner. CD-HIT does just that. Again, depending on how closely related these proteins are, we might need to try different percent identity cutoffs. We will stay with 40%."]},{"cell_type":"markdown","metadata":{},"source":["***\n","## Adding in sequences of known function\n","\n","We are almost done processing our FASTA dataset. In order to ensure that we can identify sequences of known function, we have created a small set of FASTA sequences of known function - obtained from the Protein Databank (https://rcsb.org). We will add these knowns to the FASTA file using the cat command (it stands for concatenate). \n","\n","<font color=blue><b>STEP 16:</b></font> Run the code below which combines the two files into a new file called \"files/final_40.fasta."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["!cat \"files/uniprot-dtxr_40.fasta\" \"files/dtxr_pdbs.fasta\" > \"files/final_40.fasta\""]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 17:</b></font> Run the code below to determine the final number of sequences in your file! "]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["records = list(SeqIO.parse(\"files/final_40.fasta\", \"fasta\"))    # use SeqIO to process the file\n","\n","print(\"There are %i sequences in your file.\\n\" % len(records))  # print the number of sequences"]},{"cell_type":"markdown","metadata":{},"source":["***\n","\n","## Making a BLAST-able database and performing an all-by-all BLAST search\n","\n","Now that we have our dataset we are going to use some of the NCBI BLAST tools to 1) create a searchable database using our processed set of sequences and 2) use protein BLAST to complete an all-by-all search of that database using each of those sequences.\n","\n","We will perform an analysis of the evolutionary links between these sequences by determining which sequences are most closely related to which others.\n","\n","<font color=blue><b>STEP 18:</b></font> Let's see what the command makeblastdb does by using the -h (help) flag after the program name. Go ahead and run the code below. "]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["!makeblastdb -help"]},{"cell_type":"markdown","metadata":{},"source":["The output above tells us the application of the program and how to indicate input and outout files.\n","\n","<font color=blue><b>STEP 19:</b></font> In the code box below we have added input and output files and told the program it is a protein database. The input is the files/final_40.fasta file. We will name the output files files/finalpro_40. Go ahead and run this code."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["!makeblastdb -in files/final_40.fasta -dbtype prot -out files/finalpro_40"]},{"cell_type":"markdown","metadata":{},"source":["That last step was really quite fast given the smaller number of sequences. You can try this with more sequences, it will be slower, but you might need to use more sequences to find connections among your proteins. We should be good with this number for the exercise today.\n","\n","<font color=blue><b>STEP 20:</b></font> Next we want to use blastp (Basic Local Alignment Search Tool - Protein) to take each of our protein sequences and search the database for connections. Let's start by finding out the usage of blastp by running the command below."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["!blastp -help"]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 21:</b></font> <b>This is the final step in generating new data to visualize in the next Jupyter Notebook.</b> We will run the code below to use blastp to search the files/finalpro_40 database using each of the sequences in the files/final_40.fasta file. <b>To be very clear - this will run almost 2000 BLAST searches!!!</b>\n","\n","The outfmt controls the output formatting and the -evalue is the expectation value. Briefly evalues tell us how often we expect a search to return a given result by chance. An evalue of 1 would mean that the blast result would always be found by chance (and not by an evolutionary link). Smaller evalues help us to ensure that the results are not found simply by chance. In the code below we are using a value of 10e-40. Less stringent values of 10e-22 could also be used; however, this results in more data to visualize and longer processing times. "]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["!blastp -db files/finalpro_40 -query files/final_40.fasta -outfmt \"6 qseqid sseqid evalue\" -out files/BLASTe40_out -num_threads 4 -evalue 10e-40"]},{"cell_type":"markdown","metadata":{},"source":["Once the above code has completed running you should see a file in the file browswer at the left called \"BLASTe40_out\". This is simply a text file, but contains the information we want to visualize in the next Jupyter Notebook. <b>Let's download \"BLASTe40_out\" and \"final_40.fasta\" to our personal computers so we can use them later.<b> \n","    \n","<font color=blue><b>STEP 22:</b></font> Download the file by clicking on the file \"BLASTe40_out\" and then go up to the Jupyter file menu and Download. Repeat for \"final_40.fasta\" Please note where you are downloading these files on your computer.\n","\n","<b>We will take a 10 min break here.</b>\n","\n","When we return, we will most likely need to either 1) restart the kernel or 2) relaunch the binder. Then we will open Jupyter Notebook \"3 - Visualizing_BLAST_data\" and upload our BLASTe40 file if needed."]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.10"}},"nbformat":4,"nbformat_minor":4}