3 - Visualizing_BLAST_data_v9.ipynb

{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Visualizing SSNs using an internactive network\n","\n","Networks have nodes and edges. A node is a connection point and edges are the connections between them. In this exercise each protein is a node and each edge indicates an evalue (expectation value) smaller than 10e-40.\n","\n","## Examining the BLAST output\n","\n","<font color=blue><b>STEP 1:</b></font> Let's first look at the output of our BLASTP search by double clicking the file in the file browser. <b>N.B. if you had to relaunch the binder, you will need to upload this file from your personal computer!</b> If you can't find it, you should be able to download the file here: https://drive.google.com/drive/folders/1Qjx3u6T-LIoUV4d2xtyagPN2pXPqHpx_?usp=sharing. "]},{"cell_type":"markdown","metadata":{},"source":["***\n","You should see three columns, with each column separated by a tab, they aren't labeled, but we dicated them using the outfmt command in our blast search. To refresh your memory here is the command we ran:\n","\n","~~~python\n","!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\n","~~~\n","\n","    - qseqid is the query sequence id (we will call this the source)\n","    - sseqid is the subject sequence id (we will call this the target)\n","    - evalue is the expectation value\n","\n","<font color=blue><b>STEP 2:</b></font> Answer the following questions:\n","\n","    1. What is the range of expectations value when the source and target are the same (find a few in the BLASTe40_out)?\n","    2. Given that the expectation value refers to the chance that a BLAST hit is found by chance, does your answer to question 1 make sense?\n","    3. In the first ten results in the file, what is the expectation value of the closest non-identical match (give the source, target, and evalue)?\n","    \n","***\n","\n","## Creating a Dataframe from the BLAST output\n","\n","A little lingo here. An API is an Application Programming Interface. APIs are pieces of software that allow applications to talk to each other. A dataframe is a popular API that resembles a spreadsheet. The BLAST output is a tab separated file, and we will use pandas - a \"powerful Python data analysis toolkit\" to read our file and convert it to a dataframe.\n","\n","<font color=blue><b>STEP 3:</b></font> Edit the code to replace the <b>\\<<<your file here\\>>></b> with the BLAST output file path (remember it is in the files directory). Then run the code below to convert the BLASTe40_out file into a dataframe. Since our BLAST output did not contain any headers (column labels), we can add them in.\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import pandas as pd # imports the pandas functions\n","\n","headerList = ['source','target','evalue']\n","\n","blast_data_con = pd.read_csv('<<<your file here>>>', sep='\\t', header=None)  # reads the BLAST output and looks for 'tab' to separate the values\n","\n","blast_data_con.columns = ['source', 'target', 'evalue']  # assigns names to the columns\n","\n","blast_data_con    # show what is in the dataframe."]},{"cell_type":"markdown","metadata":{},"source":["Note that the complete dataframe is not shown, but that it contains 3000-4000 connections, which we call edges in a network.\n","***\n","## Creating dataframes of edges and nodes\n","\n","## Removing duplicates and self-references from edges\n","\n","The code below removes duplicates (e.g. if <font color=\"blue\">a</font> finds <font color=\"blue\">b</font> and <font color=\"blue\">b</font> finds <font color=\"blue\">a</font>, we only need to keep one of them) and self-references (e.g. remove all instances of <font color=\"blue\">a</font> finds <font color=\"blue\">a</font>).\n","\n","The code is a bit complicated and uses another function called numpy. Briefly, the code uses pandas and the numpy.sort function to create another dataframe with only the duplicates. We then \"subtract\" the dataframe containing duplicates from the original dataframe.\n","\n","We call the new dataframe \"edges\".\n","\n","<font color=blue><b>STEP 4:</b></font> Run the code below to create a dataframe of unique edges. We have named this dataframe \"edges\"."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["import numpy as np\n","\n","df = blast_data_con       # we assign the variable df to our blast_data_con so we can retain the original dataframe.\n","\n","#remove duplicates\n","m=pd.DataFrame(np.sort(df[['source','target']])).duplicated()\n","df = df[~m]\n","\n","#removes self-reference\n","df = df[df.source != df.target]\n","\n","edges = df    # this is a unique set of edges\n","edges         # show us the edges dataframe"]},{"cell_type":"markdown","metadata":{},"source":["## Creating a unique list of nodes\n","\n","We will use the numpy.unique function to read through the sources and targets in the dataframe and make a list (called uniq_list) of nodes.\n","\n","Then we will use pandas to convert this list into a simple dataframe of nodes.\n","\n","<font color=blue><b>STEP 5:</b></font> Run the code below to create dataframe of unique nodes."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["uniq_list = np.unique(df[['source', 'target']].values)   # find the unique values and put them in a list\n","\n","nodes = pd.DataFrame(uniq_list, columns = ['id'])  # make a node dataframe with the column header id\n","\n","nodes     # show us the nodes dataframe"]},{"cell_type":"markdown","metadata":{},"source":["***\n","\n","## Visualizing nodes and edges\n","\n","You now have a set of nodes and edges that you can visualize. We will import ipycytoscape, A Cytoscape widget for displaying interactive networks. You can find more about Cytoscape (a really cool stand alone software package here: https://cytoscape.org/ and ipycytoscape here: https://github.com/cytoscape/ipycytoscape.\n"," \n","This code is borrowed and edited from https://github.com/joseberlines, who has done some neat work with ipycytoscape.\n","\n","<font color=blue><b>STEP 6:</b></font> Run the code box below to visualize the graph (this might take a few minutes and might take a few seconds to appear even after the asterisk disappears).\n","\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["# There isn't a real need to edit any of this since it is just making the network graph.\n","\n","import json                     # json stands for JavaScript Object Notation and ipycytoscape reads the data in this format\n","import ipycytoscape             # the widget to visualize interactive networks.\n","from ipywidgets import Output\n","\n","def transform_into_ipycytoscape(nodes_df,edges_df):\n","    \n","    nodes_dict = nodes_df.to_dict('records')    # converts the nodes to a dictionary\n","    edges_dict = edges_df.to_dict('records')    # convertst the edges to a dictionary \n","\n","    # building nodes\n","\n","    data_keys = ['id']  #this is a list of keys in stations (nodes)\n","    position_keys = ['position_x','position_y']\n","    rest_keys = ['score','idInt','name','score','group','removed','selected','selectable','locked','grabbed'\n","                 'grabbable']\n","    \n","    \n","    nodes_graph_list=[] #an empty list for making the json-like? file\n","    for node in nodes_dict: #iterating over each node\n","        dict_node = {}\n","        data_sub_dict = {'data':{el:node[el] for el in data_keys}}\n","        rest_sub_dict = {el:node[el] for el in node.keys() if el in rest_keys}\n","        posi_sub_dict = {}\n","        if 'position_x' in node.keys() and 'position_y' in node.keys():\n","            posi_sub_dict = {'position':{el:node[el] for el in node.keys() if el in position_keys}}\n","        \n","        dict_node = {**data_sub_dict,**rest_sub_dict,**posi_sub_dict}\n","        nodes_graph_list.append(dict_node)\n","    \n","    # building edges\n","    \n","    data_keys  = ['source','target','evalue'] #this is a list of keys in edges\n","    data_keys2 = ['label','classes'] \n","    rest_keys  = ['score','weight','group','networkId','networkGroupId','intn','rIntnId','group','removed','selected','selectable','locked','grabbed','grabbable','classes']\n","    position_keys = ['position_x','position_y']\n","    \n","    edges_graph_list = []\n","    for edge in edges_dict:\n","        dict_edge = {}\n","        data_sub_dict = {el:edge[el] for el in data_keys}\n","        data_sub_dict2 = {el:edge[el] for el in edge.keys() if el in data_keys2}\n","        rest_sub_dict = {el:edge[el] for el in edge.keys() if el in rest_keys}\n","        \n","        dict_edge = {'data':{**data_sub_dict,**data_sub_dict2},**rest_sub_dict}\n","        edges_graph_list.append(dict_edge)\n","    \n","    #print(edges_graph_list)\n","    \n","    total_graph_dict = {'nodes': nodes_graph_list, 'edges':edges_graph_list}\n","    \n","    #print(total_graph_dict)\n","    \n","    # building the style\n","    all_node_style = ['background-color','background-opacity',\n","                     'font-family','font-size','label','width',\n","                     'shape','height','width','text-valign','text-halign']\n","    all_edge_style = ['background-color','background-opacity',\n","                     'font-family','font-size','label','width','line-color', \n","                     ]\n","    \n","    total_style_dict = {}\n","    style_elements=[]\n","    for node in nodes_dict:\n","        node_dict = {'selector': f'node[id = \\\"{node[\"id\"]}\\\"]'}\n","        style_dict ={\"style\": { el:node[el] for el in node.keys() if el in all_node_style}}\n","        node_dict.update(style_dict)\n","        style_elements.append(node_dict)\n","    \n","    for edge in edges_dict:\n","        edge_dict = {'selector': f'edge[id = \\\"{edge[\"source\"]}\\\"]'}\n","        style_dict ={\"style\": { el:edge[el] for el in edge.keys() if el in all_edge_style}}\n","        edge_dict.update(style_dict)\n","        style_elements.append(edge_dict)\n","    \n","    # the graph\n","    data_graph = json.dumps(total_graph_dict)\n","    json_to_python = json.loads(data_graph)\n","    result_cyto = ipycytoscape.CytoscapeWidget()\n","    result_cyto.graph.add_graph_from_json(json_to_python)    \n","    result_cyto.set_style(style_elements)\n","    result_cyto.set_layout(name='grid')   #concentric, cola, grid, cose https://blog.js.cytoscape.org/2020/05/11/layouts/\n","\n","    out = Output()\n","\n","    \"\"\"\n","    def log_clicks(node):\n","        with out:\n","            print(f'clicked: {(node)}')\n","\n","    #def log_mouseovers(node):\n","        #with out:\n","            #print(f'mouseover: {pformat(node)}')\n","\n","    result_cyto.on('node', 'click', log_clicks)\n","    #result_cyto.on('node', 'mouseover', log_mouseovers)\n","    \n","    \"\"\"  \n","    return result_cyto, out\n","\n","network, out =transform_into_ipycytoscape(nodes,edges)\n","display(network)\n","display(out)"]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 7:</b></font> Try zooming in and out of the network. You can also grab and move nodes.\n","\n","If you are underimpressed - you should be. We did make an interactive network of nodes and edges. <b>But</b>...the problem is that this network doesn't contain enough information to give us insight into the connections among DtxR-like sequences.\n","\n","<font color=blue><b>STEP 8:</b></font> The code below adds some information to the dataframe for our nodes. It includes color and labels and sets a default size for the nodes. Run the code below to update the dataframe and then run the code to visualize the network again. In case you would like to see the full range of named colors, check out this site: https://matplotlib.org/stable/gallery/color/named_colors.html.\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["# Let's add some new columns to our dataframe\n","nodes['label'] = nodes['id']        # creates a label using the id\n","nodes['background-color']='cyan'    # our default color is cyan, but could be anything.\n","nodes['width']='24'                \n","nodes['height']='24'\n","nodes['font-size']='20'\n","nodes['text-valign']='center'\n","nodes['text-halign']='center'\n","nodes['count'] = '1'\n","\n","nodes"]},{"cell_type":"markdown","metadata":{},"source":["***\n","The network is still very difficult to use. Knowledge of the number of connections might help us to find nodes that make many connections and nodes that make only a few connections easily.\n","\n","<font color=blue><b>STEP 9:</b></font> The code below creates a list using our node IDs. It goes through the list and counts the number of edges to targets. Then it changes the variable count in the dataframe for that item. Lastly, the height and width are determined using the size * 10. I just made that up, you could try any multiple or even an exponent (e.g. \\**2).\n"]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["#Let's change the size of the nodes based on the number of connections. \n","col_one_list = nodes['id'].tolist()     # make a list from the dataframe\n","\n","for item in col_one_list: \n","    size = len(edges[edges['source']==item]) + len(edges[edges['target']==item])\n","    nodes.loc[nodes['id'] == item, 'count']=size\n","    size = size*2         # 2 is an arbitrary scale factor, you can adjust it also depending on the size and resolution of your screen; try 4, 7 or 10\n","    nodes.loc[nodes['id'] == item, 'width']=size\n","    nodes.loc[nodes['id'] == item, 'height']=size\n","\n","nodes['font-size']='20'   # let's also adjust the size of the font here.\n","    \n","nodes"]},{"cell_type":"markdown","metadata":{},"source":["This is finally starting to give us some information. \n","\n","<font color=blue><b>STEP 10:</b></font> Let's try another last amendment to our dataframe by giving a color to our knowns. Note that these are the IDs in the \"dtxr_pdbs.fasta\" file, and I just picked a different color for each. Run the code below and then rerun the network visualization."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["\n","#Here we can assign colors to nodes that connect to one of our knowns!\n","\n","nodes.loc[nodes['id'] == '1U8R_IDER','background-color']  = 'red'\n","nodes.loc[nodes['id'] == '1C0W_DTXR','background-color']  = 'orange'\n","nodes.loc[nodes['id'] == '6O5C_MTSR','background-color']  = 'yellow'\n","nodes.loc[nodes['id'] == '3HRT_SCAR','background-color']  = 'green'\n","nodes.loc[nodes['id'] == '5CVI_SLOR','background-color']  = 'blue'\n","nodes.loc[nodes['id'] == '3R60_MNTR','background-color']  = 'magenta'\n","\n","nodes"]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 11:</b></font> Answer the question:\n","\n","    1. Was adding color to your knowns as helpful as you thought it might be? Why or why not?\n","    \n","***\n","\n","Let's see if adding more color can help us to generate clusters.\n","\n","<font color=blue><b>STEP 12:</b></font> In this step, if a node is identified by a known, we will color it the same as the known. Note that the last section recolors the known nodes just in case they were identified by another known directly. Run the code below and then rerun the network visualization."]},{"cell_type":"code","execution_count":null,"metadata":{"trusted":true},"outputs":[],"source":["records = edges.to_records(index=False)\n","result = list(records)\n","\n","for item in result:\n","    #print(item)\n","    if item[1] == '1U8R_IDER':\n","        nodes.loc[nodes['id'] == item[0],'background-color']  = 'red'\n","    if item[1] == '5CVI_SLOR':\n","        nodes.loc[nodes['id'] == item[0],'background-color']  = 'blue'\n","    if item[1] == '3HRT_SCAR':\n","        nodes.loc[nodes['id'] == item[0],'background-color']  = 'green'\n","    if item[1] == '1C0W_DTXR':\n","        nodes.loc[nodes['id'] == item[0],'background-color']  = 'orange'\n","    if item[1] == '3R60_MNTR':\n","        nodes.loc[nodes['id'] == item[0],'background-color']  = 'magenta'\n","    if item[1] == '6O5C_MTSR':\n","        nodes.loc[nodes['id'] == item[0],'background-color']  = 'yellow'\n","\n","nodes.loc[nodes['id'] == '1U8R_IDER','background-color']  = 'red'\n","nodes.loc[nodes['id'] == '1C0W_DTXR','background-color']  = 'orange'\n","nodes.loc[nodes['id'] == '6O5C_MTSR','background-color']  = 'yellow'\n","nodes.loc[nodes['id'] == '3HRT_SCAR','background-color']  = 'green'\n","nodes.loc[nodes['id'] == '5CVI_SLOR','background-color']  = 'blue'\n","nodes.loc[nodes['id'] == '3R60_MNTR','background-color']  = 'magenta'\n","        \n","nodes"]},{"cell_type":"markdown","metadata":{},"source":["***\n","\n","<font color=blue><b>STEP 13:</b></font> Interact with the graph to create clusters. Briefly, cluster by first moving the larger known clusters to the outside of the grid graph. Then bring the similarly colored groups together. Finally, arrange the remaining cyan nodes nearer to their connections as dictated by the edges.\n","\n","<font color=blue><b>STEP 14:</b></font> Use your clustered graph to answer the following questions:\n","\n","    1. Which DtxR-like proteins are most closely related to each other?\n","    2. What are the MNTRs connected to? What does this mean? What step or steps might you need to change to identify sequence connections to MNTRs?\n","    3. What sequence (or sequences) connect the 5CVI_SLOR protein to the 1C0W_DTXR protein?"]},{"cell_type":"markdown","metadata":{},"source":["<font color=blue><b>STEP 15:</b></font> <b>Challenge Question:</b> Having identified the sequence links between 5CVI_SLOR protein and 1C0W_DTXR protein (Step 14 question 3), create a fasta file that contains these linking sequences and the knowns (in dtxr_pdbs.fasta). Then create a multiple alignment from this fasta. Using the mutliple alignment output, assign functions to the linking sequences - choices are: DtxR/IdeR-like, SloR/ScaR-like, MtsR-like or unknown function.   "]},{"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}