{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ?\n\nExample to call blast web-server (ncbi) via BioPython for alignment. ( See  https://blast.ncbi.nlm.nih.gov/Blast.cgi )\n\nCAFA5 protein data is taken as an example. \n\nA retrieval for one sequence takes from 20 seconds to 50+ minutes. So can be quite long. \nThe median time is about 2 minutes, mean time 5 minutes. \nMoreover sometimes one can get \"<Iteration_message>[blastsrv4.REAL]: Error: CPU usage limit was exceeded, resulting in SIGXCPU (24).</Iteration_message>\" - thus getting no results at all. \n\nThe output is restricted to 50 proteins by default (can be changed by hitlist_size param). \n\nBasically we need just the two functions from BioPython: \n\n    result_handle = NCBIWWW.qblast('blastp','nr',seq1 ) # query the web server \n\nThen result can be saved to .xml file and that .xml file can be loaded and parsed with the help of \n\n    blast_records = NCBIXML.parse(result_handle2 )  # parse xml  \n\nParsing loop goes over three layers: queried proteins, list of returned alignments and adidditional \"hsp\" which is typically just one element (so - kind of - fake layer) , then we can access: various characteristics of the alignment: 'align_length', 'bits', 'expect',  etc...\n\n\nVersion 3 of the present notebook run 7.5 hours on 100 proteins: \nhttps://www.kaggle.com/code/alexandervc/cafa5-20-ncbiwww-blast-biopython?scriptVersionId=132171102\n\n\n\nPS\n\nThanks to:\n    \nhttps://www.kaggle.com/code/shtrausslearning/biological-sequence-operations\n    \nhttps://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#7-|-BLAST    \n\n\nPSPS\n\nFor the competition purposes one can use alignments in a several ways: find most similar proteins to the given one and try to transfer labels from them to  it. Or one can use \"similarities\" obained by alignment as features, for examples selecting several sequences and \"similarities\" to them - are features (similar to Levenshtein distance feature here https://www.kaggle.com/code/alexandervc/cafa5-levenshtein-distance-features). Or one can use \"similarities\" to obtain groups for groupwise validation. Though the main cavet  is that alignments are quite slow.\n\nOne can use alignments from BioPython: \nhttps://www.kaggle.com/code/alexandervc/cafa5-18-alignments-biopython-compare , some faster version from skbio: \nhttps://www.kaggle.com/code/alexandervc/cafa5-19-alignments-skbio , or very fast \"Diamond\": https://www.kaggle.com/code/geraseva/diamond\nHowever for any option there are certain concerns - either speed or quality.\n\n","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time()\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-04T10:50:58.626811Z","iopub.execute_input":"2023-06-04T10:50:58.627484Z","iopub.status.idle":"2023-06-04T10:51:00.449674Z","shell.execute_reply.started":"2023-06-04T10:50:58.627442Z","shell.execute_reply":"2023-06-04T10:51:00.448767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Imports ","metadata":{}},{"cell_type":"code","source":"from Bio.pairwise2 import format_alignment\n# from Bio.SubsMat import MatrixInfo \nfrom Bio import pairwise2\nfrom Bio import SeqIO, SearchIO\nfrom Bio.Seq import Seq\nfrom Bio.SeqRecord import SeqRecord\nfrom Bio.Blast import NCBIWWW\nfrom Bio.Blast import NCBIXML\n\nfrom Bio.Phylo.TreeConstruction import DistanceTreeConstructor\nfrom Bio.Phylo.TreeConstruction import DistanceCalculator\nfrom Bio.Phylo.PhyloXML import Phylogeny\nfrom Bio import Phylo\nfrom pprint import pprint","metadata":{"execution":{"iopub.status.busy":"2023-06-04T10:51:00.451572Z","iopub.execute_input":"2023-06-04T10:51:00.452475Z","iopub.status.idle":"2023-06-04T10:51:00.793707Z","shell.execute_reply.started":"2023-06-04T10:51:00.45244Z","shell.execute_reply":"2023-06-04T10:51:00.792343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Simple example - just one protein","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nseq1 = 'MNSVTVSHAPYTITYHDDWEPVMSQLVEFYNEVASWLLRDETSPIPDKFFIQLKQPLRNKRVCVCGIDPYPKDGTGVPFESPNFTKKSIKEIASSISRLTGVIDYKGYNLNIIDGVIPWNYYLSCKLGETKSHAIYWDKISKLLLQHITKHVSVLYCLGKTDFSNIRAKLESPVTTIVGYHPAARDRQFEKDRSFEIINVLLELDNKVPINWAQGFIY'\n# Just first seq example take from '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\n# ('P20536',\n#  'P20536 sp|P20536|UNG_VACCC Uracil-DNA glycosylase OS=Vaccinia virus (strain Copenhagen) OX=10249 GN=UNG PE=1 SV=1')\n\n\nprint(f'Length of Sequence: {len(seq1)}')\n\nprint()\nprint('Start calling blast web server ')\nprint('Pay attention on param \"blastP(!)\": \"p\"-for proteins')\nprint('\"nr\" - database for proteins, not \"nt\" - that is for DNA ')\nprint('  (can specify other databases like \"pdb\", etc - see ncbi website: https://blast.ncbi.nlm.nih.gov/Blast.cgi ) ')\nt0 = time.time()\nresult_handle = NCBIWWW.qblast('blastp','nr',seq1 )\nprint('call finished %.1f'%(time.time()-t0),'seconds passed' )\n\n\n# Save results to disk in xml format\nfn2 = \"/kaggle/working/test_blast.xml\"\nwith open(fn2, \"w\") as save_to:\n    save_to.write(result_handle.read()) # Save results to xml \n    result_handle.close()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-04T10:51:00.795593Z","iopub.execute_input":"2023-06-04T10:51:00.796082Z","iopub.status.idle":"2023-06-04T10:52:02.569085Z","shell.execute_reply.started":"2023-06-04T10:51:00.79604Z","shell.execute_reply":"2023-06-04T10:52:02.567957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Parse and show the results","metadata":{}},{"cell_type":"code","source":"print(result_handle)","metadata":{"execution":{"iopub.status.busy":"2023-06-04T10:52:02.572584Z","iopub.execute_input":"2023-06-04T10:52:02.57303Z","iopub.status.idle":"2023-06-04T10:52:02.580353Z","shell.execute_reply.started":"2023-06-04T10:52:02.572997Z","shell.execute_reply":"2023-06-04T10:52:02.578818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result_handle2 = open(fn2,'r')\nblast_records = NCBIXML.parse(result_handle2 )  # parse xml   \n\ne_thresh = 1.00; # 1 - display ALL matches; set zero for only perfect matches\nn_results_to_show = 3 # show only  several results\n\nprint('show only first ',n_results_to_show, 'alignments')\nfor blast_record in blast_records: # go through all blast records\n    jj=-1\n    for alignment in blast_record.alignments:\n        jj+=1;\n        if(jj < n_results_to_show):\n            print(f'\\nAlignment #{jj}')\n            print(f'sequence: {alignment.title}')\n            print(f'length: {alignment.length}')\n            for hsp in alignment.hsps:\n                if(hsp.expect <= e_thresh):\n                    print('***** Alignment *****')\n                    if 1: # Show all results\n                        for t in ['align_length', 'bits', 'expect', 'frame', 'gaps', 'identities', 'match', 'num_alignments', 'positives', 'query', 'query_end', 'query_start', 'sbjct', 'sbjct_end', 'sbjct_start', 'score', 'strand']:\n                            if hasattr(hsp, t):\n                                print(t, getattr(hsp,t))\n                    else: # Show only some fields\n                        print(f'sequence: {alignment.title}')\n                        print(f'length: {alignment.length}')\n                        print(f'e value: {hsp.expect}')\n","metadata":{"execution":{"iopub.status.busy":"2023-06-04T10:52:02.581715Z","iopub.execute_input":"2023-06-04T10:52:02.582078Z","iopub.status.idle":"2023-06-04T10:52:02.612548Z","shell.execute_reply.started":"2023-06-04T10:52:02.582049Z","shell.execute_reply":"2023-06-04T10:52:02.611453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loop over several proteins and save results ","metadata":{}},{"cell_type":"code","source":"%%time\nfrom numbers import Number\nimport time\n\nn_proteins_for_query = 3 # cannot be big since one query takes 5-7 minutes\n\nhitlist_size = 100 # Number of hits to return. Default 50\n\ne_thresh = 1.00; # 1 - display ALL matches; set zero for only perfect matches\nn_results_to_show = 2 # show only  several results\n\nfile_fasta = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'#  '/kaggle/input/biopython-genbank/NC_005816.fna'\n\nsequences = SeqIO.parse(file_fasta, \"fasta\")\n# record = SeqIO.read(open(file_fasta),format='fasta')\n# record = next((sequences))\n\ndf = pd.DataFrame(); IX = -1\ncc = 0\nfor record in sequences:\n    cc += 1\n    if cc > n_proteins_for_query: break\n\n    print()\n    print(record)\n    print(f'Length of Sequence: {len(record)}')\n    t0 = time.time()\n    result_handle = NCBIWWW.qblast('blastp','nr',record.format('fasta'), hitlist_size = hitlist_size)\n    #result_handle = NCBIWWW.qblast('blastp','nr', sequence = record.seq) # record.format('fasta'))\n\n    time_query = np.round( time.time() - t0, 1)\n    print(time_query , 'secs passed')\n\n\n    fn2 = \"/kaggle/working/protein_\"+str(cc)+ \"_\"+str(record.id) + \".xml\"\n    with open(fn2, \"w\") as save_to:\n        save_to.write(result_handle.read())\n        result_handle.close()\n\n    result_handle2 = open(fn2,'r')\n    blast_records = NCBIXML.parse(result_handle2)    \n\n    for blast_record in blast_records: # go through all blast records\n        jj=-1\n        for alignment in blast_record.alignments:\n            jj+=1;\n            if 1:\n                for i_hsp,hsp in enumerate(alignment.hsps):\n                    IX+=1 \n                    df.loc[IX,'query_id'] = record.id\n                    df.loc[IX,'query_len'] = len(record.seq)\n                    df.loc[IX,'time'] = time_query\n                    df.loc[IX,'n_alignment'] = jj\n                    df.loc[IX,'n_hsp'] = i_hsp\n                    for t in ['align_length', 'bits', 'expect', 'frame', 'gaps', 'identities', 'num_alignments', 'positives', \n                               'query_end', 'query_start', 'sbjct_end', 'sbjct_start', 'score', 'strand']: # , 'query', 'sbjct', 'match', ]:\n                        if hasattr(hsp, t):\n                            if isinstance( getattr(hsp,t) , Number ):\n                                df.loc[IX,t] = getattr(hsp,t)\n                            else:\n                                df.loc[IX,t] = str( getattr(hsp,t) )\n                    \n                    t = 'title'\n                    if hasattr(alignment, t):\n                        if isinstance( getattr(alignment,t) , Number ):\n                            df.loc[IX,t] = getattr(alignment,t)\n                        else:\n                            df.loc[IX,t] = str( getattr(alignment,t) )\n                                \n            if(jj < n_results_to_show):\n                print(f'\\nAlignment #{jj}')\n                for i_hsp,hsp in enumerate(alignment.hsps):\n                    if(hsp.expect <= e_thresh):\n                        print('***** Alignment *****')\n                        if 1: # Show all results\n                            print(f'sequence: {alignment.title}')\n                            print(f'length: {alignment.length}')\n                            for t in ['title', 'align_length', 'bits', 'expect', 'frame', 'gaps', 'identities', 'match', 'num_alignments', 'positives', 'query', 'query_end', 'query_start', 'sbjct', 'sbjct_end', 'sbjct_start', 'score', 'strand']:\n                                if hasattr(hsp, t):\n                                    print(t, getattr(hsp,t))\n                                    \n                        else: # Show only some fields\n                            print(f'sequence: {alignment.title}')\n                            print(f'length: {alignment.length}')\n                            print(f'e value: {hsp.expect}')\n                            \n                        \nprint(df.shape)\ndf.to_csv('df_alignment_results.csv')\ndf.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-06-04T10:52:02.614361Z","iopub.execute_input":"2023-06-04T10:52:02.615221Z","iopub.status.idle":"2023-06-04T11:01:04.042161Z","shell.execute_reply.started":"2023-06-04T10:52:02.61518Z","shell.execute_reply":"2023-06-04T11:01:04.041181Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Statistics ","metadata":{}},{"cell_type":"code","source":"print(df.shape, df.columns)","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.043707Z","iopub.execute_input":"2023-06-04T11:01:04.044156Z","iopub.status.idle":"2023-06-04T11:01:04.050568Z","shell.execute_reply.started":"2023-06-04T11:01:04.044115Z","shell.execute_reply":"2023-06-04T11:01:04.049395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.051798Z","iopub.execute_input":"2023-06-04T11:01:04.052115Z","iopub.status.idle":"2023-06-04T11:01:04.130263Z","shell.execute_reply.started":"2023-06-04T11:01:04.052088Z","shell.execute_reply":"2023-06-04T11:01:04.129129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.groupby(df.columns[0]).mean()","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.13178Z","iopub.execute_input":"2023-06-04T11:01:04.132141Z","iopub.status.idle":"2023-06-04T11:01:04.16615Z","shell.execute_reply.started":"2023-06-04T11:01:04.132111Z","shell.execute_reply":"2023-06-04T11:01:04.165087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.groupby(df.columns[0]).min()","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.169279Z","iopub.execute_input":"2023-06-04T11:01:04.170012Z","iopub.status.idle":"2023-06-04T11:01:04.2068Z","shell.execute_reply.started":"2023-06-04T11:01:04.169979Z","shell.execute_reply":"2023-06-04T11:01:04.205821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.groupby(df.columns[0]).max()","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.208433Z","iopub.execute_input":"2023-06-04T11:01:04.208784Z","iopub.status.idle":"2023-06-04T11:01:04.244493Z","shell.execute_reply.started":"2023-06-04T11:01:04.208756Z","shell.execute_reply":"2023-06-04T11:01:04.243239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df.head(20) )","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.246077Z","iopub.execute_input":"2023-06-04T11:01:04.246696Z","iopub.status.idle":"2023-06-04T11:01:04.30789Z","shell.execute_reply.started":"2023-06-04T11:01:04.246657Z","shell.execute_reply":"2023-06-04T11:01:04.306958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df.tail(20) )","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.309173Z","iopub.execute_input":"2023-06-04T11:01:04.309466Z","iopub.status.idle":"2023-06-04T11:01:04.372478Z","shell.execute_reply.started":"2023-06-04T11:01:04.309441Z","shell.execute_reply":"2023-06-04T11:01:04.371342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2023-06-04T11:01:04.373945Z","iopub.execute_input":"2023-06-04T11:01:04.374254Z","iopub.status.idle":"2023-06-04T11:01:04.379587Z","shell.execute_reply.started":"2023-06-04T11:01:04.374228Z","shell.execute_reply":"2023-06-04T11:01:04.378599Z"},"trusted":true},"execution_count":null,"outputs":[]}]}