{"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\nSequences alignments are extremely important algorithms used in modern bioinformatics: https://en.wikipedia.org/wiki/Sequence_alignment\nSee e.g. very nice BioPython Kaggle tutorial notebook by ANDREY SHTRAUSS : https://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#5-|-PAIRWISE-SEQUENCE-ALIGNMENT (and other notebooks by the same author).\n\nHere we first compare the results of the alignments by different methods and params on the CAFA5 dataset. \n\nWe also compare \"new\" and \"standard\" (but depricated) BioPython syntax for the alignment functions. \n**NEW one is MUCH FASTER (!) up to 30-200 times!!!**\n\n\n**Outcomes:** we see that global and local alignment with default params agree not so bad with each other - \n\nBiopython implementations are slow. So we can process only not big lists.\n\nNotes:  simple alignment are quite faster than with substitutions ; \nglobal - faster than local\n\nOutcomes: we see not so bad agreement of different methods:\n\n    Let us consider topN by similarity proteins by different alignments methods and compare these lists:\n    \n    For 100 proteins (notebook version 12).\n    Local vs global with default params: \n    Top50 - intersection size is  48; intersection with top50 for Levenshtein  - 38, so Levenshtein quite agrees with the alighnments.\n    However when we consider intersection with top50 for local alignment with blossum62 (and different penalties),\n    intersection falls down to 25. \n    Thus - such parameters like substitution matrix and penalties - can change results quite strongly.\n    \n\n    For 10_000 proteins (notebook version 11 - 6 hour computation):\n    Local vs global with default params: \n    Top20 - intersection size is  13;  intersection with top20 for Levenshtein  - 0, so for large sets Levenshtein does NOT quite agrees with the alighnments. \n    Top100 - intersetction is 81\n    However when we consider intersection with results obtained with non-trivial substitution matrices (and penalties),\n    the itersections fall down: \n    Sizes of intersection for four methods+params (local/global and blossum62, penalties) \n    Intersection: 0 out of  200 top\n    Intersection: 5 out of  300 top\n    Intersection: 13 out of  500 top\n    Intersection: 58 out of  1000 top\n\nSo it seems thus we should be careful with choosing appropriate parameters for the alignment - substitution matrix, penalties - results might change depending on that. Though more tests are required to estimate how big are the changes due to these parameters. \n\n\nPS \n\nFor the competition purposes one can use alignments in a several ways:\nfind most similar proteins to the given one and try to transfer labels from them to the it. \nOr 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\nhttps://www.kaggle.com/code/alexandervc/cafa5-levenshtein-distance-features). \nOr one can use \"similarities\" to obtain groups for groupwise validation. \nThough the main cavet is that alignment are quite slow.\n\nPS\n\nThanks to : https://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#5-|-PAIRWISE-SEQUENCE-ALIGNMENT\n\n\"Diamond\" package works much faster and allows extremely fast alignment many x many :\nSee Liza Geraseva notebook: https://www.kaggle.com/code/geraseva/diamond\n\nAnother package \"scbio\" provides another Python packages for the alignment see examples:\nhttps://www.kaggle.com/alexandervc/cafa5-19-alignments-sсbio\nThey work much faster \n\n","metadata":{}},{"cell_type":"markdown","source":"## Key params:","metadata":{}},{"cell_type":"code","source":"# Alignments work quite slowly so we need to restrict number of samples to process: \nn_process_globalxx = 100 # simple global alignment - without substitution matrix, etc\nn_process_globalds = 100 # global alignment with substitution matrix , opening gap penalty -4, extension penalty -1\nn_process_localxx = 100 # simple local alignment - without substitution matrix, etc \nn_process_localds = 100 # local alignment with substitution matrix , opening gap penalty -4, extension penalty -1\n\nn_print_every_n_samples = 10\n\n# simple alignment are quite faster than with substitutions \n# global - faster than local","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:15:57.933513Z","iopub.execute_input":"2023-05-31T17:15:57.934586Z","iopub.status.idle":"2023-05-31T17:15:57.967493Z","shell.execute_reply.started":"2023-05-31T17:15:57.934529Z","shell.execute_reply":"2023-05-31T17:15:57.965986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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# 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-05-31T17:15:57.969971Z","iopub.execute_input":"2023-05-31T17:15:57.970412Z","iopub.status.idle":"2023-05-31T17:15:59.31543Z","shell.execute_reply.started":"2023-05-31T17:15:57.970378Z","shell.execute_reply":"2023-05-31T17:15:59.314763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-05-31T17:15:59.31643Z","iopub.execute_input":"2023-05-31T17:15:59.317739Z","iopub.status.idle":"2023-05-31T17:15:59.572342Z","shell.execute_reply.started":"2023-05-31T17:15:59.317669Z","shell.execute_reply":"2023-05-31T17:15:59.571696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Substitution matrices ","metadata":{}},{"cell_type":"code","source":"from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \nnames = substitution_matrices.load()\nprint(names )\nsm = substitution_matrices.load(\"BLOSUM62\") \nprint(type(sm), sm.shape )\nprint(np.array(sm)[:3,:5])\nprint( sm.alphabet )\nsm","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:15:59.573898Z","iopub.execute_input":"2023-05-31T17:15:59.574462Z","iopub.status.idle":"2023-05-31T17:15:59.590655Z","shell.execute_reply.started":"2023-05-31T17:15:59.574439Z","shell.execute_reply":"2023-05-31T17:15:59.589988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load Train data","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\nprint(\"Sequence example:\\n\\n\", next(iter(SeqIO.parse(fn, \"fasta\"))))\nsequences = SeqIO.parse(fn, \"fasta\")\nnum_sequences = sum(1 for seq in sequences)\nprint()\nprint(\"Number of sequences in super-test:\", num_sequences)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:15:59.592169Z","iopub.execute_input":"2023-05-31T17:15:59.592494Z","iopub.status.idle":"2023-05-31T17:16:01.957008Z","shell.execute_reply.started":"2023-05-31T17:15:59.592471Z","shell.execute_reply":"2023-05-31T17:16:01.955721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nsequences = SeqIO.parse(fn, \"fasta\")\nlist_seqs = [str(seq.seq) for seq in sequences]\nprint(len(list_seqs),  list_seqs[:3])\nsequences = SeqIO.parse(fn, \"fasta\")\nlist_ids = [str(seq.id) for seq in sequences]\nprint(len(list_ids),  list_ids[:3])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:01.958312Z","iopub.execute_input":"2023-05-31T17:16:01.958645Z","iopub.status.idle":"2023-05-31T17:16:04.546828Z","shell.execute_reply.started":"2023-05-31T17:16:01.958619Z","shell.execute_reply":"2023-05-31T17:16:04.546165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install python-Levenshtein\nfrom Levenshtein import distance\nedit_dist = distance(\"ah\", \"aho\")\nedit_dist","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:04.548017Z","iopub.execute_input":"2023-05-31T17:16:04.548262Z","iopub.status.idle":"2023-05-31T17:16:14.680319Z","shell.execute_reply.started":"2023-05-31T17:16:04.548239Z","shell.execute_reply":"2023-05-31T17:16:14.679357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Global Alignment","metadata":{}},{"cell_type":"markdown","source":"## Toy examples","metadata":{}},{"cell_type":"code","source":"%%time\nseq1 = list_seqs[0]\nseq2 = list_seqs[2]\n\n#lPRT = pairwise2.align.globalds(seq1,seq2 ,sm,-4,-1)\nlPRT = pairwise2.align.globalxx(seq1,seq2)#  ,sm,-4,-1)\nfor i in lPRT:\n    #print(i)\n    print('len1=%d,len2=%d'%(len(seq1), len(seq2) ) )\n    print('Levenshtein distance:', distance(seq1, seq2) )\n    print('score=%d,start=%d, end=%d'%(i[2],i[3],i[4]))\n    print(format_alignment(*i))\n    break\n\n#lPRT = pairwise2.align.globalds(seq1,seq2 ,sm,-4,-1)\nsm = substitution_matrices.load(\"BLOSUM62\")\nlPRT = pairwise2.align.globalds(seq1,seq2  ,sm,-4,-1)\nfor i in lPRT:\n    #print(i)\n    print('len1=%d,len2=%d'%(len(seq1), len(seq2) ) )\n    print('Levenshtein distance:', distance(seq1, seq2) )\n    print('score=%d,start=%d, end=%d'%(i[2],i[3],i[4]))\n    print(format_alignment(*i))\n    break\n    ","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:14.681565Z","iopub.execute_input":"2023-05-31T17:16:14.681811Z","iopub.status.idle":"2023-05-31T17:16:14.975693Z","shell.execute_reply.started":"2023-05-31T17:16:14.681788Z","shell.execute_reply":"2023-05-31T17:16:14.974826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"i","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:14.979451Z","iopub.execute_input":"2023-05-31T17:16:14.979726Z","iopub.status.idle":"2023-05-31T17:16:14.986302Z","shell.execute_reply.started":"2023-05-31T17:16:14.979705Z","shell.execute_reply":"2023-05-31T17:16:14.985138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Global alignment for many proteins ","metadata":{}},{"cell_type":"code","source":"dict_df = {}","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:14.988042Z","iopub.execute_input":"2023-05-31T17:16:14.988385Z","iopub.status.idle":"2023-05-31T17:16:14.997836Z","shell.execute_reply.started":"2023-05-31T17:16:14.988359Z","shell.execute_reply":"2023-05-31T17:16:14.996834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# globalxx - 3secs on 100 seqs:\n#\n# For globalds - quite slower:\n# 1min for 100\n# 100_000: 1000 minutes for one protein - 16 hours \n\nimport time\nt0 = time.time()\ncc = 0\nl = []\ndf_stat = pd.DataFrame(); IX_df_stat = 0\nfor seq1 in list_seqs:\n    # if  ('U' in seq1) or ('Y' in seq1): continue  \n    print(len(seq1), seq1)\n    print()\n    for ic2,seq2 in enumerate(list_seqs):\n        if seq1==seq2: continue \n        sm = substitution_matrices.load(\"BLOSUM62\")\n        # print(cc, len(seq1), len(seq2))\n        # if  ('U' in seq2) or ('Y' in seq2): continue\n        #if  ('U' in seq2): continue\n        #if  ('O' in seq2): continue\n        lPRT = pairwise2.align.globalxx(seq1,seq2)# ,sm,-4,-1)\n        s = lPRT[0][2]\n        _i = lPRT[0]\n        ss = s/(_i[4] - _i[3] )\n        ed =  distance(seq1, seq2)\n        if (cc %n_print_every_n_samples)==0:\n            print(cc,'similarity: %.2f'%(ss), len(seq1), len(seq2), 'score:',  s , 'Levenshtein distance:',ed,'start=%d, end=%d'%(_i[3],_i[4]),\n                 'time: %.1f'%(time.time() - t0))\n        df_stat.loc[IX_df_stat, 'Id'] = list_ids[ic2]    \n        df_stat.loc[IX_df_stat,'similarity'] = ss\n        df_stat.loc[IX_df_stat,'score'] = s\n        df_stat.loc[IX_df_stat,'Levenshtein'] = ed\n        df_stat.loc[IX_df_stat,'len alignment'] = _i[4] - _i[3]\n        df_stat.loc[IX_df_stat,'len'] = len(seq2)\n        IX_df_stat += 1\n        l.append(s)\n        cc += 1\n        if cc>=n_process_globalxx: # 142246:\n            break\n    break\n\nprint( '%.1f'%(time.time() - t0) )\n\n# display( pd.Series(l).describe() )\n# plt.hist(l, bins  = 20)\n# plt.show()\n\ndict_df['globalxx'] = df_stat.sort_values('similarity', ascending = False) \n\ndisplay(df_stat.sort_values('similarity', ascending = False).head(20))\ndisplay(df_stat.sort_values('score', ascending = False).head(20))\ndisplay(df_stat.sort_values('Levenshtein', ascending = True).head(20))\n\ndisplay(df_stat.describe() )\nplt.figure(figsize = (20,4))\nll = df_stat.columns[1:]\nfor i,col in enumerate(ll):\n    plt.subplot(1,len(ll),i+1 )\n    plt.hist(df_stat[col], bins = 100 )\n    plt.title(col,fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:14.999055Z","iopub.execute_input":"2023-05-31T17:16:14.999339Z","iopub.status.idle":"2023-05-31T17:16:18.506278Z","shell.execute_reply.started":"2023-05-31T17:16:14.999314Z","shell.execute_reply":"2023-05-31T17:16:18.505382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.sort_values('similarity', ascending = False).to_csv('df_stat_globalxx.csv')","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:18.507132Z","iopub.execute_input":"2023-05-31T17:16:18.507336Z","iopub.status.idle":"2023-05-31T17:16:18.520574Z","shell.execute_reply.started":"2023-05-31T17:16:18.507317Z","shell.execute_reply":"2023-05-31T17:16:18.519119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## - Use substitute matrix BLOSUM64, opening gap penalty -4, extension penalty -1","metadata":{}},{"cell_type":"code","source":"%%time\n# globalxx - 3secs on 100 seqs:\n#\n# For globalds - quite slower:\n# 1min for 100\n# 100_000: 1000 minutes for one protein - 16 hours \n\nimport time\nt0 = time.time()\ncc = 0\nl = []\ndf_stat = pd.DataFrame(); IX_df_stat = 0\nfor seq1 in list_seqs:\n    # if  ('U' in seq1) or ('Y' in seq1): continue  \n    print(len(seq1), seq1)\n    print()\n    for ic2,seq2 in enumerate(list_seqs):\n        if seq1==seq2: continue \n        sm = substitution_matrices.load(\"BLOSUM62\")\n        # print(cc, len(seq1), len(seq2))\n        # if  ('U' in seq2) or ('Y' in seq2): continue\n        if  ('U' in seq2): continue\n        if  ('O' in seq2): continue\n        #sm.alphabet\n        \n        # - Using a substitute matrix BLOSUM64, opening gap penalty -4, extension penalty -1\n        lPRT = pairwise2.align.globalds(seq1,seq2 ,sm,-4,-1)\n        s = lPRT[0][2]\n        _i = lPRT[0]\n        ss = s/(_i[4] - _i[3] )\n        ed =  distance(seq1, seq2)\n        if (cc %n_print_every_n_samples)==0:\n            print(cc,'similarity: %.2f'%(ss), len(seq1), len(seq2), 'score:',  s , 'Levenshtein distance:',ed,'start=%d, end=%d'%(_i[3],_i[4]),\n                 'time: %.1f'%(time.time() - t0))\n        df_stat.loc[IX_df_stat, 'Id'] = list_ids[ic2]    \n        df_stat.loc[IX_df_stat,'similarity'] = ss\n        df_stat.loc[IX_df_stat,'score'] = s\n        df_stat.loc[IX_df_stat,'Levenshtein'] = ed\n        df_stat.loc[IX_df_stat,'len alignment'] = _i[4] - _i[3]\n        df_stat.loc[IX_df_stat,'len'] = len(seq2)\n        IX_df_stat += 1\n        l.append(s)\n        cc += 1\n        if cc>=n_process_globalds: # 142246:\n            break\n    break\n\nprint( '%.1f'%(time.time() - t0) )\n\n# display( pd.Series(l).describe() )\n# plt.hist(l, bins  = 20)\n# plt.show()\n\ndict_df['globalds'] = df_stat.sort_values('similarity', ascending = False) \n\ndisplay(df_stat.sort_values('similarity', ascending = False).head(20))\ndisplay(df_stat.sort_values('Levenshtein', ascending = True).head(20))\ndisplay(df_stat.sort_values('score', ascending = False).head(20))\ndisplay(df_stat.describe() )\nplt.figure(figsize = (20,4))\nll = df_stat.columns[1:]\nfor i,col in enumerate(ll):\n    plt.subplot(1,len(ll),i+1 )\n    plt.hist(df_stat[col], bins = 100 )\n    plt.title(col,fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:16:18.522778Z","iopub.execute_input":"2023-05-31T17:16:18.5238Z","iopub.status.idle":"2023-05-31T17:17:01.11414Z","shell.execute_reply.started":"2023-05-31T17:16:18.52376Z","shell.execute_reply":"2023-05-31T17:17:01.11275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.sort_values('similarity', ascending = False).to_csv('df_stat_globalds.csv')","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:17:01.115585Z","iopub.execute_input":"2023-05-31T17:17:01.115863Z","iopub.status.idle":"2023-05-31T17:17:01.123483Z","shell.execute_reply.started":"2023-05-31T17:17:01.11584Z","shell.execute_reply":"2023-05-31T17:17:01.122214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Local alignment","metadata":{}},{"cell_type":"markdown","source":"## Toy Examples","metadata":{}},{"cell_type":"code","source":"''' Local Alignment of two DNA sequences ''' \n# Local Alignment Problem of DNA Sequences\n# Match score = 3, mismatch score = -2, constant gap penalty g=-3 x2\n\nseq1 = 'ATGGCAGATAGA'\nseq2 = 'ATAGAGAATAG'\n\nlDNA = pairwise2.align.localms(seq1,seq2, 3,-2,-3,-3)\nfor i in lDNA: \n    print(format_alignment(*i))\n    \nprint();    print();    print();    \nprint( ''' Local Alignment of Protein Sequences ''' )\n# - Local Alignment Problem of Protein Sequences\n# - Using a substitute matrix BLOSUM64, opening gap penalty -4, extension penalty -1\n\npseq1 = \"EVSAW\"\npseq2 = \"KEVLA\"\npPROT = pairwise2.align.localds(pseq1,pseq2,sm,-4,-1)\nfor i in pPROT:\n    print(format_alignment(*i))","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:17:01.124902Z","iopub.execute_input":"2023-05-31T17:17:01.125237Z","iopub.status.idle":"2023-05-31T17:17:01.144887Z","shell.execute_reply.started":"2023-05-31T17:17:01.12521Z","shell.execute_reply":"2023-05-31T17:17:01.14357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## local without substitution ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nt0 = time.time()\ncc = 0\nl = []\ndf_stat = pd.DataFrame(); IX_df_stat = 0\nfor seq1 in list_seqs:\n    # if  ('U' in seq1) or ('Y' in seq1): continue  \n    print(len(seq1), seq1)\n    print()\n    for ic2, seq2 in enumerate(list_seqs):\n        if seq1==seq2: continue \n        sm = substitution_matrices.load(\"BLOSUM62\")\n        # print(cc, len(seq1), len(seq2))\n        # if  ('U' in seq2) or ('Y' in seq2): continue\n        if  ('U' in seq2): continue\n        if  ('O' in seq2): continue\n        #sm.alphabet\n        \n        # - Using a substitute matrix BLOSUM64, opening gap penalty -4, extension penalty -1\n#         lPRT = pairwise2.align.globalds(seq1,seq2 ,sm,-4,-1)\n#         lPRT = pairwise2.align.localds(seq1,seq2,sm,-4,-1)\n        lPRT = pairwise2.align.localxx(seq1,seq2)# ,sm,-4,-1)\n        s = lPRT[0][2]\n        _i = lPRT[0]\n        ss = s/(_i[4] - _i[3] )\n        ed =  distance(seq1, seq2)\n        if (cc %n_print_every_n_samples)==0:\n            print(cc,'similarity: %.2f'%(ss), len(seq1), len(seq2), 'score:',  s , 'Levenshtein distance:',ed,'start=%d, end=%d'%(_i[3],_i[4]),\n                 'time: %.1f'%(time.time() - t0))\n        df_stat.loc[IX_df_stat, 'Id'] = list_ids[ic2]    \n        df_stat.loc[IX_df_stat,'similarity'] = ss\n        df_stat.loc[IX_df_stat,'score'] = s\n        df_stat.loc[IX_df_stat,'Levenshtein'] = ed\n        df_stat.loc[IX_df_stat,'len alignment'] = _i[4] - _i[3]\n        df_stat.loc[IX_df_stat,'len'] = len(seq2)\n        IX_df_stat += 1\n        l.append(s)\n        cc += 1\n        if cc>=n_process_localxx: # 142246:\n            break\n    break\n\nprint( '%.1f'%(time.time() - t0) )\n\n# display( pd.Series(l).describe() )\n# plt.hist(l, bins  = 20)\n# plt.show()\n\ndict_df['localxx'] = df_stat.sort_values('similarity', ascending = False) \n\ndisplay(df_stat.sort_values('similarity', ascending = False).head(20))\ndisplay(df_stat.sort_values('Levenshtein', ascending = True).head(20))\ndisplay(df_stat.sort_values('score', ascending = False).head(20))\n\ndisplay(df_stat.describe() )\nplt.figure(figsize = (20,4))\nll = df_stat.columns[1:]\nfor i,col in enumerate(ll):\n    plt.subplot(1,len(ll),i+1 )\n    plt.hist(df_stat[col], bins = 100 )\n    plt.title(col,fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:17:01.146019Z","iopub.execute_input":"2023-05-31T17:17:01.146262Z","iopub.status.idle":"2023-05-31T17:17:09.070562Z","shell.execute_reply.started":"2023-05-31T17:17:01.146241Z","shell.execute_reply":"2023-05-31T17:17:09.0693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.sort_values('similarity', ascending = False).to_csv('df_stat_localxx.csv')","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:17:09.072032Z","iopub.execute_input":"2023-05-31T17:17:09.072301Z","iopub.status.idle":"2023-05-31T17:17:09.082722Z","shell.execute_reply.started":"2023-05-31T17:17:09.072279Z","shell.execute_reply":"2023-05-31T17:17:09.081176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## local with substitution matrix ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nt0 = time.time()\ncc = 0\nl = []\ndf_stat = pd.DataFrame(); IX_df_stat = 0\nfor seq1 in list_seqs:\n    # if  ('U' in seq1) or ('Y' in seq1): continue  \n    print(len(seq1), seq1)\n    print()\n    for ic2, seq2 in enumerate(list_seqs):\n        if seq1==seq2: continue \n        sm = substitution_matrices.load(\"BLOSUM62\")\n        # print(cc, len(seq1), len(seq2))\n        # if  ('U' in seq2) or ('Y' in seq2): continue\n        if  ('U' in seq2): continue\n        if  ('O' in seq2): continue\n        #sm.alphabet\n        \n        # - Using a substitute matrix BLOSUM64, opening gap penalty -4, extension penalty -1\n#         lPRT = pairwise2.align.globalds(seq1,seq2 ,sm,-4,-1)\n        lPRT = pairwise2.align.localds(seq1,seq2,sm,-4,-1)\n        s = lPRT[0][2]\n        _i = lPRT[0]\n        ss = s/(_i[4] - _i[3] )\n        ed =  distance(seq1, seq2)\n        if (cc %n_print_every_n_samples)==0:\n            print(cc,'similarity: %.2f'%(ss), len(seq1), len(seq2), 'score:',  s , 'Levenshtein distance:',ed,'start=%d, end=%d'%(_i[3],_i[4]),\n                 'time: %.1f'%(time.time() - t0))\n        df_stat.loc[IX_df_stat, 'Id'] = list_ids[ic2]    \n        df_stat.loc[IX_df_stat,'similarity'] = ss\n        df_stat.loc[IX_df_stat,'score'] = s\n        df_stat.loc[IX_df_stat,'Levenshtein'] = ed\n        df_stat.loc[IX_df_stat,'len alignment'] = _i[4] - _i[3]\n        df_stat.loc[IX_df_stat,'len'] = len(seq2)\n        IX_df_stat += 1\n        l.append(s)\n        cc += 1\n        if cc>=n_process_localds: # 142246:\n            break\n    break\n\nprint( '%.1f'%(time.time() - t0) )\n\n# display( pd.Series(l).describe() )\n# plt.hist(l, bins  = 20)\n# plt.show()\n\ndict_df['localds'] = df_stat.sort_values('similarity', ascending = False) \n\ndisplay(df_stat.sort_values('similarity', ascending = False).head(20))\ndisplay(df_stat.sort_values('Levenshtein', ascending = True).head(20))\ndisplay(df_stat.sort_values('score', ascending = False).head(20))\n\n\ndisplay(df_stat.describe() )\nplt.figure(figsize = (20,4))\nll = df_stat.columns[1:]\nfor i,col in enumerate(ll):\n    plt.subplot(1,len(ll),i+1 )\n    plt.hist(df_stat[col], bins = 100 )\n    plt.title(col,fontsize = 20 )\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:17:09.084313Z","iopub.execute_input":"2023-05-31T17:17:09.084653Z","iopub.status.idle":"2023-05-31T17:18:09.7077Z","shell.execute_reply.started":"2023-05-31T17:17:09.084592Z","shell.execute_reply":"2023-05-31T17:18:09.706845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.sort_values('similarity', ascending = False).to_csv('df_stat_localds.csv')","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:18:09.708739Z","iopub.execute_input":"2023-05-31T17:18:09.708993Z","iopub.status.idle":"2023-05-31T17:18:09.71556Z","shell.execute_reply.started":"2023-05-31T17:18:09.708941Z","shell.execute_reply":"2023-05-31T17:18:09.714594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compare alignments ","metadata":{}},{"cell_type":"code","source":"N = 50\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:48:51.859666Z","iopub.execute_input":"2023-05-31T17:48:51.860136Z","iopub.status.idle":"2023-05-31T17:48:51.865271Z","shell.execute_reply.started":"2023-05-31T17:48:51.860102Z","shell.execute_reply":"2023-05-31T17:48:51.864199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Simple global and local seems to quite agree and partially agree with the Levenshtein \n","metadata":{}},{"cell_type":"code","source":"for ic,k in enumerate( ['globalxx', 'localxx'] ): # enumerate(dict_df):\n    df = dict_df[k]\n    print(k)\n    print(df['Id'].iloc[:N].tolist() )\n    if ic == 0:\n        s = set(df['Id'].iloc[:N].tolist() )\n    else:\n        s = s & set(df['Id'].iloc[:N].tolist() )\n    print('Intersection:', len(s), 'out of',N,'top')\n    if len(s) < 5: print(s)\n\n    \ndf = dict_df[list(dict_df.keys())[0]]\ndf = dict_df[k]\nprint('Levenshtein')\nprint(df.sort_values('Levenshtein')['Id'].iloc[:N].tolist() )\ns = s & set(df.sort_values('Levenshtein')['Id'].iloc[:N].tolist() )\nprint('Intersection to Levenshtein:', len(s), 'out of',N,'top')\nif len(s) < 5: print(s)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:48:53.22002Z","iopub.execute_input":"2023-05-31T17:48:53.22036Z","iopub.status.idle":"2023-05-31T17:48:53.230658Z","shell.execute_reply.started":"2023-05-31T17:48:53.220334Z","shell.execute_reply":"2023-05-31T17:48:53.229608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Compare methods with substitition matrix and penalties ","metadata":{}},{"cell_type":"code","source":"for ic,k in enumerate( ['globalxx', 'localxx'] ): # enumerate(dict_df):\n    df = dict_df[k]\n    print(k)\n    print(df['Id'].iloc[:20].tolist() )\n    if ic == 0:\n        s = set(df['Id'].iloc[:N].tolist() )\n    else:\n        s = s & set(df['Id'].iloc[:N].tolist() )\n    print('Intersection:', len(s),'out of', N,'top')\n    if len(s) < 5: print(s)\n\n    \ndf = dict_df[list(dict_df.keys())[0]]\ndf = dict_df[k]\nprint('Levenshtein')\nprint(df.sort_values('Levenshtein')['Id'].iloc[:N].tolist() )\ns = s & set(df.sort_values('Levenshtein')['Id'].iloc[:N].tolist() )\nprint('Intersection to Levenshtein:', len(s), 'out of', N,'top')\nif len(s) < 5: print(s)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:42:01.379655Z","iopub.execute_input":"2023-05-31T17:42:01.380047Z","iopub.status.idle":"2023-05-31T17:42:01.390931Z","shell.execute_reply.started":"2023-05-31T17:42:01.38002Z","shell.execute_reply":"2023-05-31T17:42:01.389751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## In general agreement is also not bad for about 50% for N=50","metadata":{}},{"cell_type":"code","source":"# N = 50\nfor ic,k in enumerate(dict_df):\n    df = dict_df[k]\n    print(k)\n    print(df['Id'].iloc[:20].tolist() )\n    if ic == 0:\n        s = set(df['Id'].iloc[:N].tolist() )\n    else:\n        s = s & set(df['Id'].iloc[:N].tolist() )\n    print('Intersection:', len(s), 'out of ', N, 'top')\n    if len(s) < 5: print(s)\n    \nprint()\nfor NN in [10,20,50,100,200,300,500,1000]:\n    for ic,k in enumerate(dict_df):\n        df = dict_df[k]\n        #print(k)\n        #print(df['Id'].iloc[:20].tolist() )\n        if ic == 0:\n            s = set(df['Id'].iloc[:NN].tolist() )\n        else:\n            s = s & set(df['Id'].iloc[:NN].tolist() )\n    print('Intersection:', len(s), 'out of ', NN, 'top')\n    if (len(s) < 5)and(len(s)>0): print(s)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:42:14.908763Z","iopub.execute_input":"2023-05-31T17:42:14.909122Z","iopub.status.idle":"2023-05-31T17:42:14.921195Z","shell.execute_reply.started":"2023-05-31T17:42:14.909096Z","shell.execute_reply":"2023-05-31T17:42:14.919942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# N = 80\ndf = dict_df[list(dict_df.keys())[0]]\ndf = dict_df[k]\nprint('Levenshtein')\nprint(df.sort_values('Levenshtein')['Id'].iloc[:N].tolist() )\ns = set(df.sort_values('Levenshtein')['Id'].iloc[:N].tolist() )\n\nfor ic,k in enumerate(dict_df):\n    df = dict_df[k]\n    print(k)\n    print(df['Id'].iloc[:N].tolist() )\n    s = s & set(df['Id'].iloc[:N].tolist() )\n    print('Intersection:', len(s),'óut of',N,'top')\n    if len(s) < 5: print(s)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:42:23.85429Z","iopub.execute_input":"2023-05-31T17:42:23.854644Z","iopub.status.idle":"2023-05-31T17:42:23.862932Z","shell.execute_reply.started":"2023-05-31T17:42:23.854616Z","shell.execute_reply":"2023-05-31T17:42:23.862351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# New syntax for BioPython alignments and speed tests \n\nThe syntax used above is deprecated and new syntax is introduced here we give some examples.\n\nNew functions work MUCH faster than the older ones !!! ","metadata":{}},{"cell_type":"markdown","source":"## New syntax","metadata":{}},{"cell_type":"code","source":"%%time \nseq1 = \"EVSAW\"\nseq2 = \"KEVLA\"\n\nprint(); \nprint('---------------- global alignment from BioPython ----------------------')\nprint();\n\nfrom Bio import Align\naligner = Align.PairwiseAligner()\n#alignments = aligner.align(\"TACCG\", \"ACG\")\nalignments = aligner.align(seq1, seq2 )\nprint(alignments[0])\n\nprint(); \nprint('---------------- local alignment from BioPython ----------------------')\nprint();\n\nfrom Bio import Align\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\nalignments = aligner.align(seq1, seq2 )\nprint(alignments[0])\n\n\n\nprint(); \nprint('---------------- local alignment from BioPython with BLOSUM62 ----------------------')\nprint();\n\nfrom Bio import Align\nfrom Bio.Align import substitution_matrices\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\naligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\nalignments = aligner.align(seq1, seq2 )\nprint(alignments[0])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:19:29.056355Z","iopub.execute_input":"2023-05-31T17:19:29.056708Z","iopub.status.idle":"2023-05-31T17:19:29.068247Z","shell.execute_reply.started":"2023-05-31T17:19:29.056681Z","shell.execute_reply":"2023-05-31T17:19:29.066535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Speed tests ","metadata":{}},{"cell_type":"code","source":"%%time\nseq1 = list_seqs[0]\nprint(len(seq1), seq1[:10])\n\nN = 1000\nprint('------------------- Global alignment --------------------- ')\nprint()\nprint('------------------- new syntax --------------------- ')\n\nfrom Bio import Align\naligner = Align.PairwiseAligner()\n#alignments = aligner.align(\"TACCG\", \"ACG\")\nt0 = time.time()\nl = [aligner.align(seq1, seq2 ) for seq2 in list_seqs[:N] ] \n#print(l[0][0])\nprint('%.1f seconds for %d alignments'%( (time.time() - t0), N ))\n\nprint('------------------- old syntax - slower --------------------- ')\nt0 = time.time()\nl = [pairwise2.align.globalxx(seq1,seq2) for seq2 in list_seqs[:N] ] \nprint('%.1f seconds for %d alignments'%( (time.time() - t0), N ))\n\n\nprint('------------------- Local alignment --------------------- ')\nprint()\nprint('------------------- new syntax --------------------- ')\n\nfrom Bio import Align\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\nt0 = time.time()\nl = [aligner.align(seq1, seq2 ) for seq2 in list_seqs[:N] ] \nprint('%.1f seconds for %d alignments'%( (time.time() - t0), N ))\n\nprint('------------------- old syntax - slower --------------------- ')\nt0 = time.time()\nl = [pairwise2.align.localxx(seq1,seq2) for seq2 in list_seqs[:N] ] \nprint('%.1f seconds for %d alignments'%( (time.time() - t0), N ))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:27:51.284263Z","iopub.execute_input":"2023-05-31T17:27:51.2846Z","iopub.status.idle":"2023-05-31T17:29:11.400743Z","shell.execute_reply.started":"2023-05-31T17:27:51.284575Z","shell.execute_reply":"2023-05-31T17:29:11.399874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint('------------------- Local alignment with substitution matrix BLOSUM62  --------------------- ')\nprint()\nprint('------------------- new syntax --------------------- ')\n\nN= 10\nfrom Bio import Align\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\naligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\nt0 = time.time()\nl = [aligner.align(seq1, seq2 ) for seq2 in list_seqs[:N] ] \nprint('%.3f seconds for %d alignments'%( (time.time() - t0), N ))\n\nprint('------------------- old syntax - slower --------------------- ')\nsm = substitution_matrices.load(\"BLOSUM62\")\nt0 = time.time()\nl = [pairwise2.align.localds(seq1,seq2,sm,-4,-1) for seq2 in list_seqs[:N] ] \nprint('%.3f seconds for %d alignments'%( (time.time() - t0), N ))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T17:32:43.609101Z","iopub.execute_input":"2023-05-31T17:32:43.609504Z","iopub.status.idle":"2023-05-31T17:32:51.473319Z","shell.execute_reply.started":"2023-05-31T17:32:43.609475Z","shell.execute_reply":"2023-05-31T17:32:51.471869Z"},"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":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{},"execution_count":null,"outputs":[]}]}