{"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 \n\nHere we explore sequence aligner \"Parasail\":  https://github.com/jeffdaily/parasail-python\nDaily, Jeff. (2016). Parasail: SIMD C library for global, semi-global, and local pairwise sequence alignments. BMC Bioinformatics, 17(1), 1-11. doi:10.1186/s12859-016-0930-z  http://dx.doi.org/10.1186/s12859-016-0930-z\n\nOutcomes: \n\n    It is really very fast, but strangely enough the quality is not good enough. Despite the claim that package implements the classical  Smith-Waterman (local) algorithm . It gives different and worse results than BioPython alignment Smith-Waterman  , however BioPython is much much slower. (T5 embeddings give similar results (in many, but not all cases) to Biopython slow-hard alignment ). \n    \n    \nTechnical notes:\n\n    (!) Be very careful with \"_8\" mode like: parasail.sw_stats_striped_8 , parasiail.sw_scan_8 (\"_8\" at the end) - it wraps the scores to 255 or 0 (depending on mode). Thus scores are incorrect. And different for \"_stats\" and \"_scan\" one sends high scores to 255, another to 0.  So do not use \"_8\", use \"_16\".\n  \n    \"_stats\" is desribed to return only statistics, comparing \"_scan\" with detailed alignments. So might be thought to work faster. It is NOT the case.     That is strange. So use \"_scan\", not \"_stats\"\n\n\nSome discussion and comparison can be found at: https://www.biostars.org/p/9505903/\n\n\nPS \n\nLists of aligners can be found: https://github.com/danielecook/Awesome-Bioinformatics#pairwise , \nhttps://en.wikipedia.org/wiki/List_of_sequence_alignment_software\n\nSee  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). (But pay attention that it uses slower and depricated alignment ways of BioPython see for newer ones:  https://www.kaggle.com/code/alexandervc/cafa5-18-alignments-biopython-compare )\n\nSome use-examples of other aligners can be found: \nweb-server NCBI blast: https://www.kaggle.com/code/alexandervc/cafa5-20-ncbiwww-blast-biopython\nskbio: https://www.kaggle.com/code/alexandervc/cafa5-19-alignments-skbio\nBioPython: https://www.kaggle.com/code/alexandervc/cafa5-18-alignments-biopython-compare\nDiamond: https://www.kaggle.com/code/geraseva/diamond\nLevenshtein distance can be thought as \"simplified\" global alignment score, there are many notebooks on it, e.g.: https://www.kaggle.com/code/alexandervc/cafa5-levenshtein-distance-features , \nor GPU way from top-grandmaster Chris Deotte:  https://www.kaggle.com/code/cdeotte/train-data-contains-mutations-like-test-data https://docs.rapids.ai/api/cudf/stable/api_docs/api/cudf.core.column.string.stringmethods.edit_distance_matrix/#cudf.core.column.string.StringMethods.edit_distance_matrix\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 the 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.","metadata":{}},{"cell_type":"markdown","source":"# Install import","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\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-07T13:02:43.747894Z","iopub.execute_input":"2023-06-07T13:02:43.748474Z","iopub.status.idle":"2023-06-07T13:02:43.784603Z","shell.execute_reply.started":"2023-06-07T13:02:43.748437Z","shell.execute_reply":"2023-06-07T13:02:43.783243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \n!pip install parasail\nimport parasail\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:02:43.786987Z","iopub.execute_input":"2023-06-07T13:02:43.788218Z","iopub.status.idle":"2023-06-07T13:02:57.232645Z","shell.execute_reply.started":"2023-06-07T13:02:43.788178Z","shell.execute_reply":"2023-06-07T13:02:57.231014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Toy examples ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T19:50:20.760009Z","iopub.execute_input":"2023-06-06T19:50:20.760961Z","iopub.status.idle":"2023-06-06T19:50:20.804981Z","shell.execute_reply.started":"2023-06-06T19:50:20.760916Z","shell.execute_reply":"2023-06-06T19:50:20.803812Z"}}},{"cell_type":"code","source":"%%time \nimport parasail\nresult = similarity_name\nprint(result.score)\n\nresult = parasail.sw_stats_striped_8(\"asdf\", \"asdf\", 11, 1, parasail.pam100)\nprint(result.score)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:02:57.234953Z","iopub.execute_input":"2023-06-07T13:02:57.23551Z","iopub.status.idle":"2023-06-07T13:02:57.246486Z","shell.execute_reply.started":"2023-06-07T13:02:57.235454Z","shell.execute_reply":"2023-06-07T13:02:57.244832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Tecnical problem: Score truncation to 255,0 problem for \"_8\" mode. Solution: use \"_16\" - mode , do not use \"_8\"-mode","metadata":{}},{"cell_type":"markdown","source":"## In \"_stats\" mode truncation goes to 255","metadata":{}},{"cell_type":"code","source":"%%time \nimport parasail\nprint(\"_16 mode works okay: \")\nresult = parasail.sw_scan_16(\"asdf\"*100, \"asdf\"*100, 11, 1, parasail.blosum62)\nprint(result.score)\nresult = parasail.sw_scan_16(\"eklmn\"*100, \"klmnxyz\"*100, 11, 1, parasail.pam100)\nprint(result.score)\nprint()\n\nprint('WARNING(!) the score truncation for \"_8\" mode - different results are truncated to the same 255 ')\nresult = parasail.sw_stats_striped_8(\"asdf\"*100, \"asdf\"*100, 11, 1, parasail.pam100)\nprint('WARNING(!) the score truncation for \"_8\" mode: ')\nprint(result.score)\nresult = parasail.sw_stats_striped_8(\"eklmn\"*100, \"klmnxyz\"*100, 11, 1, parasail.pam100)\nprint('WARNING(!) the score truncation for \"_8\" mode: ')\nprint(result.score)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:02:57.251918Z","iopub.execute_input":"2023-06-07T13:02:57.252499Z","iopub.status.idle":"2023-06-07T13:02:57.266344Z","shell.execute_reply.started":"2023-06-07T13:02:57.25245Z","shell.execute_reply":"2023-06-07T13:02:57.264762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## In \"_scan\" mode truncation goes to 0","metadata":{}},{"cell_type":"code","source":"print('WARNING(!) the score truncation for \"_8\" mode - different results are truncated to the same 255 ')\nresult = parasail.sw_scan_8(\"asdf\"*100, \"asdf\"*100, 11, 1, parasail.pam100)\nprint('WARNING(!) the score truncation for \"_8\" mode: ')\nprint(result.score)\nresult = parasail.sw_scan_8(\"eklmn\"*100, \"klmnxyz\"*100, 11, 1, parasail.pam100)\nprint('WARNING(!) the score truncation for \"_8\" mode: ')\nprint(result.score)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:02:57.268361Z","iopub.execute_input":"2023-06-07T13:02:57.268828Z","iopub.status.idle":"2023-06-07T13:02:57.28136Z","shell.execute_reply.started":"2023-06-07T13:02:57.268795Z","shell.execute_reply":"2023-06-07T13:02:57.280243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Speed tests on random data: 13.6 second for 100_000 alignments of both length 536  , 3.93 secs for 100_000 length 200","metadata":{}},{"cell_type":"code","source":"%%time\nprot_name = 'SRC human P12931'\nsel_seq = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\nprint(len(sel_seq))\nN = 100_000\nlist_seqs = []\nfor i in range(N):\n    characters = list(sel_seq)\n    np.random.shuffle(characters)\n    s = ''.join(characters ) \n    list_seqs.append( s )\n    \nprint(len(list_seqs), list_seqs[0][:10])","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:02:57.282753Z","iopub.execute_input":"2023-06-07T13:02:57.283253Z","iopub.status.idle":"2023-06-07T13:03:01.863104Z","shell.execute_reply.started":"2023-06-07T13:02:57.283207Z","shell.execute_reply":"2023-06-07T13:03:01.861771Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nl = [  parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs ]","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:01.865023Z","iopub.execute_input":"2023-06-07T13:03:01.867043Z","iopub.status.idle":"2023-06-07T13:03:15.399382Z","shell.execute_reply.started":"2023-06-07T13:03:01.867003Z","shell.execute_reply":"2023-06-07T13:03:15.398405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure( figsize = (20,6))\nplt.hist(l, bins = 100)\nplt.title('Score distribution for random strings alignment of length '+str(len(sel_seq)) , fontsize = 20 )\nplt.show()\ndisplay( pd.Series(l).describe() )","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:15.401079Z","iopub.execute_input":"2023-06-07T13:03:15.401739Z","iopub.status.idle":"2023-06-07T13:03:16.542938Z","shell.execute_reply.started":"2023-06-07T13:03:15.401706Z","shell.execute_reply":"2023-06-07T13:03:16.541725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Second test - other length","metadata":{}},{"cell_type":"code","source":"%%time\n#prot_name = 'SRC human P12931'\nsel_seq = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\nsel_seq = sel_seq[:200] # Truncate\nprint(len(sel_seq))\nN = 100_000\nlist_seqs = []\nfor i in range(N):\n    characters = list(sel_seq)\n    np.random.shuffle(characters)\n    s = ''.join(characters ) \n    list_seqs.append( s )\n    \nprint(len(list_seqs), list_seqs[0][:10])\n\nimport time\nt0 = time.time()\nl = [  parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs ]\nprint('%.4f seconds on %d alignments both lenghth %d'%(time.time()-t0, N, len(sel_seq) ))\n\nplt.figure( figsize = (20,6))\nplt.hist(l, bins = 100)\nplt.title('Score distribution for random strings alignment of length '+str(len(sel_seq)) , fontsize = 20 )\nplt.show()\ndisplay( pd.Series(l).describe() )","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:16.544512Z","iopub.execute_input":"2023-06-07T13:03:16.545112Z","iopub.status.idle":"2023-06-07T13:03:23.459685Z","shell.execute_reply.started":"2023-06-07T13:03:16.54505Z","shell.execute_reply":"2023-06-07T13:03:23.45845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load CAFA5 protein data \n\nAbout 140_000 in train and the same in test. But they intersect in about 70k proteins.\nTrain - 3600+ organisms, test 91 organisms. \n\nFor EDA see: \nhttps://www.kaggle.com/code/alexandervc/cafa5-towards-eda","metadata":{}},{"cell_type":"code","source":"%%time\n\n# --------------------------------------------------------------------------------------------------------------\n# ############################################# Import BioPython ###############################################\n# --------------------------------------------------------------------------------------------------------------\n\nfrom 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\n\n# --------------------------------------------------------------------------------------------------------------\n# ############################################# Load CAFA5 train and test data ##################################\n# --------------------------------------------------------------------------------------------------------------\n\nfile_fasta = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'#  '/kaggle/input/biopython-genbank/NC_005816.fna'\nsequences = SeqIO.parse(file_fasta, \"fasta\")\nlist_ids = [seq.id for seq in sequences]\nsequences = SeqIO.parse(file_fasta, \"fasta\")\nlist_lens = [len(seq.seq) for seq in sequences]\nsequences = SeqIO.parse(file_fasta, \"fasta\")\nlist_seqs = [str(seq.seq) for seq in sequences]\nsequences = SeqIO.parse(file_fasta, \"fasta\")\ndict_ids2description = {seq.id:seq.description for  seq in sequences }\nprint('Train:')\nprint('Mean len: %.1f, median %.1f'%(np.mean(list_lens), np.median(list_lens)) )\nprint('Seqs examples:', list_seqs[:2])\nprint()\n\n\nfile_fasta2 = '/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta'\nsequences = SeqIO.parse(file_fasta2, \"fasta\")\nlist_ids2 = [seq.id for seq in sequences]\nsequences = SeqIO.parse(file_fasta2, \"fasta\")\nlist_lens2 = [len(seq.seq) for seq in sequences]\nsequences = SeqIO.parse(file_fasta2, \"fasta\")\nlist_seqs2 = [str(seq.seq) for seq in sequences]\nprint('Test data:')\nprint('Mean len: %.1f, median %.1f'%(np.mean(list_lens2), np.median(list_lens2)) )\n\n# Our sequences contain small number of examples with  \"O\" and \"U\" symbols \n# they are NOT present in blossum62 alphabet\n# so we will create corrected sequence lists \n# \tO\t4\tO\t0\tPyrrolysine\tO\tPyl\n# 3\tU\t148\tU\t218\tSelenocysteine\tU\tSec\n# just substitute O,U by *\n\nlist_seqs_corrected = [ str(seq).replace('O','*').replace('U','*') for seq in list_seqs ]\nlist_seqs_corrected2 = [ str(seq).replace('O','*').replace('U','*') for seq in list_seqs2]\n\n\n# %%time\n# --------------------------------------------------------------------------------------------------------------\n# ############################################# Plot histograms ################################################\n# --------------------------------------------------------------------------------------------------------------\nl1 = np.array(list_lens)\nl2 = np.array(list_lens2)\n\nd = pd.concat( (pd.Series(l1).describe(), pd.Series(l2).describe() ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.subplot(121)\nplt.hist(l1[l1<1500], bins=150)\nplt.title('Train Lengths')\nplt.subplot(122)\nplt.hist(l2[l2<1500], bins=150)\nplt.title('Test Lenghts')\nplt.show()\nfor t in [10,50,100,200,300,350, 400,410,420,430,440, 450]:\n    print('Threshold:',t, 'Lefter than threshold - ', 'Train:', (l1<t).sum(), 'Test:', (l2<t).sum(), )","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:23.463683Z","iopub.execute_input":"2023-06-07T13:03:23.464106Z","iopub.status.idle":"2023-06-07T13:03:40.257478Z","shell.execute_reply.started":"2023-06-07T13:03:23.464075Z","shell.execute_reply":"2023-06-07T13:03:40.256018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Gene Ontology terms for the proteins ","metadata":{}},{"cell_type":"code","source":"%%time\n# fn =  str(path) + '/Train/train_terms.tsv'\nfn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv'\ndf = pd.read_csv(fn , sep = '\\t')# , index_col = 0)\ndisplay(df.head(3) )\n\nmap_protein_id_to_GOterms = df.groupby('EntryID')['term'].apply(lambda x: list(x))\nmap_protein_id_to_GOterms.head(2)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:40.259265Z","iopub.execute_input":"2023-06-07T13:03:40.259612Z","iopub.status.idle":"2023-06-07T13:03:51.176389Z","shell.execute_reply.started":"2023-06-07T13:03:40.259584Z","shell.execute_reply":"2023-06-07T13:03:51.175221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set the selected sequence - will align all CAFA5 to it","metadata":{}},{"cell_type":"code","source":"sel_prot_id = 'P12931'\nif  sel_prot_id == 'P12931':\n    prot_name = 'SRC human P12931'\n    sel_seq = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\n    print('Selected protein:', 'SRC - P12931- Proto-oncogene tyrosine-protein kinase Src - https://www.uniprot.org/uniprotkb/P12931/entry' )\n    print('https://en.wikipedia.org/wiki/Proto-oncogene_tyrosine-protein_kinase_Src')\n    \n\n    file_fasta = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'#  '/kaggle/input/biopython-genbank/NC_005816.fna'\n    sequences = SeqIO.parse(file_fasta, \"fasta\")\n\n    for seq in sequences:\n        if seq.id == sel_prot_id:\n            print('Description:', seq.description) \n            break\n            \n    list_keywords_describing_protein = ['kinase', 'protein kinase']\n    window_len_for_keywords_indicator_moving_average = 50 # SRC # 1 # P53\n\n    list_cutoffs_to_show_distribution_Levenshtein = [402, 399 ] # SRC  # the threshold is around 74-75 of length \n    list_cutoffs_to_show_distribution_LocalAlignmentBioPython = [70, 100, 300, 568, 574, 575] # SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = [2200,2300, 2250 ] # SRC \n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:51.178418Z","iopub.execute_input":"2023-06-07T13:03:51.179299Z","iopub.status.idle":"2023-06-07T13:03:51.985454Z","shell.execute_reply.started":"2023-06-07T13:03:51.179257Z","shell.execute_reply":"2023-06-07T13:03:51.983248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nN_seqs_to_take = 10_000  # to make fast test we can take not all CAFA5 data but first N of it, to conser all: see it to big number e.g. 200_000\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:51.986893Z","iopub.execute_input":"2023-06-07T13:03:51.987241Z","iopub.status.idle":"2023-06-07T13:03:51.992828Z","shell.execute_reply.started":"2023-06-07T13:03:51.987211Z","shell.execute_reply":"2023-06-07T13:03:51.991598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_scoring_data = {} # Results of various alignments calcualtions will be stored here ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:51.994604Z","iopub.execute_input":"2023-06-07T13:03:51.994998Z","iopub.status.idle":"2023-06-07T13:03:52.018129Z","shell.execute_reply.started":"2023-06-07T13:03:51.994967Z","shell.execute_reply":"2023-06-07T13:03:52.016473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pd.set_option('display.max_columns', 600)\n# pd.set_option('display.width', 1500)\npd.set_option('max_colwidth', 500)\npd.set_option('display.max_rows', 500)","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:52.020417Z","iopub.execute_input":"2023-06-07T13:03:52.021027Z","iopub.status.idle":"2023-06-07T13:03:52.032313Z","shell.execute_reply.started":"2023-06-07T13:03:52.020972Z","shell.execute_reply":"2023-06-07T13:03:52.030951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ParasailSW alignment CAFA5 to the selected sequence and analysis","metadata":{}},{"cell_type":"code","source":"%%time\nsimilarity_name = 'ParasailSW'\nprint('Similarity:', similarity_name, '. Length selected sequence:', len(sel_seq) , '. Protein:', prot_name )\n\n\n# ---------------------------------------------------------------------------------\n# Computes alignments\n# --------------------------------------------------------------------------------\n\nt0 = time.time()\nl1 = np.array([ parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs[:N_seqs_to_take] ] ) # CAFA5 Train\nprint('%.4f seconds on %d alignments,  first  lenghth %d'%(time.time()-t0, len(l1),  len(sel_seq) ))\nt0 = time.time()\nl2 = np.array([ parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs2[:N_seqs_to_take] ] ) # CAFA5 Test\nprint('%.4f seconds on %d alignments,  first  lenghth %d'%(time.time()-t0, len(l2),  len(sel_seq) ))\n\n\n# ---------------------------------------------------------------------------------\n# Show stastics\n# --------------------------------------------------------------------------------\n\nd = pd.concat( (pd.Series(l1).describe(percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999]), pd.Series(l2).describe( percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999] ) ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 20)\nplt.subplot(121)\nplt.hist(l1[l1<3*len(sel_seq)], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<3*len(sel_seq)], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\nfor t in [80, 100, 200]:\n    print('Threshold:',t, 'Greater than threshold - ', 'Train:', (l1>t).sum(), 'Test:', (l2>t).sum(), )\n    \n# ---------------------------------------------------------------------------------\n# Convert alignment scores to dataframes and save to csv     \n# Concat train and test, \n# Add description information (human readble info on proteins) for train part of data . Save to csv \n# --------------------------------------------------------------------------------\ndi1 = pd.DataFrame(index = list_ids[:len(l1)] , data = np.array( [l1, list_lens[:len(l1)]]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi2 = pd.DataFrame(index = list_ids2[:len(l2)] , data = np.array( [l2, list_lens2[:len(l2)]]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN_loc = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N_loc) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\n\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use \n\n\n# ---------------------------------------------------------------------------------\n# Show top scored proteins with human readable information to check that results are as expected:    \n#\n# --------------------------------------------------------------------------------\nprint()\nprint('Show top scored proteins with human readable information')\nN_to_show_top_scored_proteins_with_desription = 10\nfor N1 in [0,400,4000,40000]:\n    if N1+N_to_show_top_scored_proteins_with_desription < len(di):\n        print('Top ', N1, N1+N_to_show_top_scored_proteins_with_desription)\n        display(di.iloc[N1:(N1+N_to_show_top_scored_proteins_with_desription),:])\n\nprint()    \n    \n# ---------------------------------------------------------------------------------\nprint('Statists of the key description terms like \"kinase\" which are expected to be found in description of the  similar proteins  ' )\nprint('One expects - to find more such words for the top similar and lower for less similar ')\nprint('Better aligners are expected to have such curve higher than the worse ones ')\n# --------------------------------------------------------------------------------\nWL = window_len_for_keywords_indicator_moving_average # 50 # SRC # 1 P53\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # list_keywords_describing_protein = ['kinase', 'protein kinase']\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()        \n\nprint()    \n    \n# ---------------------------------------------------------------------------------\nprint('Compare Gene Ontology similarity (functional similarity) vs  structrual  similarity (alignment or embedding) ' )\nprint('One expects - structural similarity implies functional similarity, thus similar proteins to should have similar GO-terms')\nprint('Better aligners are expected to have such curve higher than the worse ones ')\n# --------------------------------------------------------------------------------\n    \n# Compute dice index between GO terms of the selected protein (sel_prot_id) and i-th similar to it protein.\n# Then create plots for that vector.\n# The code supports simultenous processing of the results from the several alignments stored in dict_scoring_data (calcualted above)\n# dict_scoring_data[k][0] - dataframe ordered by similarity , k-th index gives the protein ID of the k-th similar  \ndf_res = pd.DataFrame()\nlist_similarities_to_compare = list(dict_scoring_data.keys() )\nfor i1,k in enumerate( list_similarities_to_compare  ):\n    dt = dict_scoring_data[k][0]\n    s0 = set( map_protein_id_to_GOterms[ sel_prot_id ] )\n#   print(k, len(s0)); display(dt.head(2));  \n    l = []\n    for nn in range(len(dt)):\n        s = set( map_protein_id_to_GOterms[ dt.index[nn] ] )\n        dice = 2 * len( s & s0 ) / ( len(s) + len(s0) )\n        l.append(dice)\n    df_res[k] = l\ndisplay( df_res.head(2) )\ndisplay( df_res.corr() )\nprint('Proteins count:', len(df_res) )\nfor WL in [50,5]:\n    for N in [100,1000,10000,100_000]:\n        plt.figure(figsize = (20,6))\n        plt.title('Dice between GO-terms for ' + prot_name + ' vs i-th similar protein. Cutoff: '+str(N) + ' WL '+str(WL), fontsize = 18 )\n        for col in df_res.columns:\n            plt.plot(df_res[col].iloc[:N].rolling(window=WL).mean(), label = col)\n        plt.legend()\n        plt.grid()\n        plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:03:52.03415Z","iopub.execute_input":"2023-06-07T13:03:52.034542Z","iopub.status.idle":"2023-06-07T13:04:00.839136Z","shell.execute_reply.started":"2023-06-07T13:03:52.034512Z","shell.execute_reply":"2023-06-07T13:04:00.837721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show human readable description for top scored proteins ","metadata":{}},{"cell_type":"code","source":"%%time \ndisplay(di.iloc[:4,:])\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:04:00.841036Z","iopub.execute_input":"2023-06-07T13:04:00.841528Z","iopub.status.idle":"2023-06-07T13:04:00.860022Z","shell.execute_reply.started":"2023-06-07T13:04:00.841486Z","shell.execute_reply":"2023-06-07T13:04:00.858441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# T5embeddings similarity of CAFA5 to the selected sequence and analysis","metadata":{}},{"cell_type":"code","source":"%%time\nsimilarity_name = 'T5embeds PearCorr'\n\nfn = '/kaggle/input/t5embeds/train_embeds.npy'\nX = np.load(fn)\nprint(X.shape)\nprint(X[:2,:3] )\n\nfn = '/kaggle/input/t5embeds/train_ids.npy'\nvec_train_protein_ids = np.load(fn)\nprint(vec_train_protein_ids.shape)\nvec_train_protein_ids\n\nfn = '/kaggle/input/t5embeds/test_embeds.npy'\nX2 = np.load(fn)\nprint(X2.shape)\nprint(X2[:2,:3] )\n\nfn = '/kaggle/input/t5embeds/test_ids.npy'\nvec_test_protein_ids = np.load(fn)\nprint(vec_test_protein_ids.shape)\nvec_test_protein_ids","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:04:00.862354Z","iopub.execute_input":"2023-06-07T13:04:00.862845Z","iopub.status.idle":"2023-06-07T13:04:01.738482Z","shell.execute_reply.started":"2023-06-07T13:04:00.862799Z","shell.execute_reply":"2023-06-07T13:04:01.737276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.feature_selection import r_regression # Computes Pearson correlation between all columns of X, and one column y\n\nprint('Similarity:', similarity_name, '. Length selected sequence:', len(sel_seq) , '. Protein:', prot_name )\n\n# ---------------------------------------------------------------------------------\n# Compute  similarity\n# --------------------------------------------------------------------------------\n\n# Train part:\nIX = np.where( vec_train_protein_ids == sel_prot_id)[0]; sel_vec1 = X[IX[0], :]; print(sel_vec1.shape) # Get vector of the emebedding for the selected sequence\nt0 = time.time()\n# l1 = np.array([ parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs[:N_seqs_to_take] ] ) # CAFA5 Train\nl1 = r_regression(X.T[:,:N_seqs_to_take] ,sel_vec1 ) # CAFA5 Train # Computes Pearson correlation between all columns of X, and one column y\nprint(similarity_name, 'similarity. %.4f seconds on %d,  first  lenghth %d'%(time.time()-t0, len(l1),  len(sel_seq) ))\n\n# Test part: \nIX = np.where( vec_test_protein_ids == sel_prot_id)[0]; sel_vec2 = X2[IX[0], :]; print(sel_vec2.shape) # # Get vector of the emebedding for the selected sequence\nprint('Check should be zero:', (sel_vec2 != sel_vec1).sum() , ' vector we extracted from train part should be same as from the test part - sanity check')\nt0 = time.time()\n# l2 = np.array([ parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs2[:N_seqs_to_take] ] ) # CAFA5 Test\nl2 = r_regression(X2.T[:,:N_seqs_to_take] ,sel_vec2 ) # CAFA5 Test # Computes Pearson correlation between all columns of X, and one column y\nprint(similarity_name, 'similarity. %.4f seconds on %d,  first  lenghth %d'%(time.time()-t0, len(l2),  len(sel_seq) ))\n\n\n# ---------------------------------------------------------------------------------\n# Show stastics\n# --------------------------------------------------------------------------------\n\nd = pd.concat( (pd.Series(l1).describe(percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999]), pd.Series(l2).describe( percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999] ) ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 20)\nplt.subplot(121)\nplt.hist(l1[l1<3*len(sel_seq)], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<3*len(sel_seq)], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\nfor t in [80, 100, 200]:\n    print('Threshold:',t, 'Greater than threshold - ', 'Train:', (l1>t).sum(), 'Test:', (l2>t).sum(), )\n    \n# ---------------------------------------------------------------------------------\n# Convert alignment scores to dataframes and save to csv     \n# Concat train and test, \n# Add description information (human readble info on proteins) for train part of data . Save to csv \n# --------------------------------------------------------------------------------\ndi1 = pd.DataFrame(index = list_ids[:len(l1)] , data = np.array( [l1, list_lens[:len(l1)]]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi2 = pd.DataFrame(index = list_ids2[:len(l2)] , data = np.array( [l2, list_lens2[:len(l2)]]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN_loc = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N_loc) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\n\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use \n\n\n# ---------------------------------------------------------------------------------\n# Show top scored proteins with human readable information to check that results are as expected:    \n#\n# --------------------------------------------------------------------------------\nprint()\nprint('Show top scored proteins with human readable information')\nN_to_show_top_scored_proteins_with_desription = 10\nfor N1 in [0,400,4000,40000]:\n    if N1+N_to_show_top_scored_proteins_with_desription < len(di):\n        print('Top ', N1, N1+N_to_show_top_scored_proteins_with_desription)\n        display(di.iloc[N1:(N1+N_to_show_top_scored_proteins_with_desription),:])\n\nprint()    \n    \n# ---------------------------------------------------------------------------------\nprint('Statists of the key description terms like \"kinase\" which are expected to be found in description of the  similar proteins  ' )\nprint('One expects - to find more such words for the top similar and lower for less similar ')\nprint('Better aligners are expected to have such curve higher than the worse ones ')\n# --------------------------------------------------------------------------------\nWL = window_len_for_keywords_indicator_moving_average # 50 # SRC # 1 P53\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # list_keywords_describing_protein = ['kinase', 'protein kinase']\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()        \n\nprint()    \n    \n# ---------------------------------------------------------------------------------\nprint('Compare Gene Ontology similarity (functional similarity) vs  structrual  similarity (alignment or embedding) ' )\nprint('One expects - structural similarity implies functional similarity, thus similar proteins to should have similar GO-terms')\nprint('Better aligners are expected to have such curve higher than the worse ones ')\n# --------------------------------------------------------------------------------\n    \n# Compute dice index between GO terms of the selected protein (sel_prot_id) and i-th similar to it protein.\n# Then create plots for that vector.\n# The code supports simultenous processing of the results from the several alignments stored in dict_scoring_data (calcualted above)\n# dict_scoring_data[k][0] - dataframe ordered by similarity , k-th index gives the protein ID of the k-th similar  \ndf_res = pd.DataFrame()\nlist_similarities_to_compare = list(dict_scoring_data.keys() )\nfor i1,k in enumerate( list_similarities_to_compare  ):\n    dt = dict_scoring_data[k][0]\n    s0 = set( map_protein_id_to_GOterms[ sel_prot_id ] )\n#   print(k, len(s0)); display(dt.head(2));  \n    l = []\n    for nn in range(len(dt)):\n        s = set( map_protein_id_to_GOterms[ dt.index[nn] ] )\n        dice = 2 * len( s & s0 ) / ( len(s) + len(s0) )\n        l.append(dice)\n    df_res[k] = l\ndisplay( df_res.head(2) )\ndisplay( df_res.corr() )\nprint('Proteins count:', len(df_res) )\n\nfor WL in [50,5]:\n    for N in [100,1000,10000,100_000]:\n        plt.figure(figsize = (20,6))\n        plt.title('Dice between GO-terms for ' + prot_name + ' vs i-th similar protein. Cutoff: '+str(N) + ' WL '+str(WL), fontsize = 18 )\n        for col in df_res.columns:\n            plt.plot(df_res[col].iloc[:N].rolling(window=WL).mean(), label = col)\n        plt.legend()\n        plt.grid()\n        plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:04:01.740397Z","iopub.execute_input":"2023-06-07T13:04:01.741148Z","iopub.status.idle":"2023-06-07T13:04:08.555967Z","shell.execute_reply.started":"2023-06-07T13:04:01.741107Z","shell.execute_reply":"2023-06-07T13:04:08.554731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# BioPython alignment","metadata":{}},{"cell_type":"code","source":"%%time \nsimilarity_name = 'BioPython Local Blossum62'\nprint( similarity_name )\n\nfrom Bio import Align\nfrom Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\naligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\naligner.open_gap_score = -11\naligner.extend_gap_score = -1\naligner.target_end_gap_score = -11\naligner.query_end_gap_score = -1\n# Parameters: Blossum62 and gap-score: 11,1 are used in NCBI web server blastP (\"P\" - for protein) : https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome\n# hopefully the role of params in BioPython and NCBI is the same,\n# at least we get reasonable results \n\nprint(aligner)\n\n\n\n# ---------------------------------------------------------------------------------\n# Compute  similarity\n# --------------------------------------------------------------------------------\n\n# Train part:\nt0 = time.time()\n# l1 = np.array([ parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs[:N_seqs_to_take] ] ) # CAFA5 Train\n# l1 = r_regression(X.T[:,:N_seqs_to_take] ,sel_vec1 ) # CAFA5 Train # Computes Pearson correlation between all columns of X, and one column y\n# N = 10_000 # 4.65 secs  for 1000 , so 40 secs for 10_000, 400 secs for 100_000, 1200 secs( 20 minutes) for 300 k \nlist_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs_corrected[:N_seqs_to_take] ]\nprint(similarity_name, 'similarity. %.4f seconds on %d,  first  lenghth %d'%(time.time()-t0, len(l1),  len(sel_seq) ))\n\n# Test part: \nt0 = time.time()\n# l2 = np.array([ parasail.sw_scan_16(sel_seq, seq, 11, 1, parasail.blosum62).score for seq in list_seqs2[:N_seqs_to_take] ] ) # CAFA5 Test\n# l2 = r_regression(X2.T[:,:N_seqs_to_take] ,sel_vec2 ) # CAFA5 Test # Computes Pearson correlation between all columns of X, and one column y\nlist_scores2 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs_corrected2[:N_seqs_to_take] ]\nprint(similarity_name, 'similarity. %.4f seconds on %d,  first  lenghth %d'%(time.time()-t0, len(l2),  len(sel_seq) ))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:06:56.433925Z","iopub.execute_input":"2023-06-07T13:06:56.434426Z","iopub.status.idle":"2023-06-07T13:10:14.261318Z","shell.execute_reply.started":"2023-06-07T13:06:56.434395Z","shell.execute_reply":"2023-06-07T13:10:14.259664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scores_analysis(list_scores1,list_scores2, similarity_name,  prot_name, sel_seq,    )","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:36:53.805217Z","iopub.execute_input":"2023-06-07T13:36:53.805675Z","iopub.status.idle":"2023-06-07T13:36:53.907197Z","shell.execute_reply.started":"2023-06-07T13:36:53.805641Z","shell.execute_reply":"2023-06-07T13:36:53.905233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def scores_analysis(list_scores1,list_scores2, similarity_name,  prot_name, sel_seq,    ):\n    \n    l1 = np.asarray(list_scores1)\n    l2 = np.asarray(list_scores2)\n    print('Similarity:', similarity_name, 'selected sequence:', prot_name, 'length:', len(sel_seq), 'n_scores:', len(l1))\n    # ---------------------------------------------------------------------------------\n    # Show stastics\n    # --------------------------------------------------------------------------------\n    \n    d = pd.concat( (pd.Series(l1).describe(percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999]), pd.Series(l2).describe( percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999] ) ) , axis = 1 )\n    d.columns = ['train', 'test'] \n    display(d)\n\n    plt.figure(figsize=(20,4))\n    plt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 20)\n    plt.subplot(121)\n    plt.hist(l1[l1<3*len(sel_seq)], bins=150)\n    plt.title('Train', fontsize = 20, loc='left' )\n    plt.subplot(122)\n    plt.hist(l2[l2<3*len(sel_seq)], bins=150)\n    plt.title('Test', fontsize = 20, loc='right')\n    plt.show()\n\n    if 'embed' in similarity_name:\n        for t in [0.7, 0.8, 0.9]:\n            print('Threshold:',t, 'Greater than threshold - ', 'Train:', (l1>t).sum(), 'Test:', (l2>t).sum(), )\n    elif 'Levenshtein' in similarity_name:\n        for t in np.array([0.72, 0.73, 0.74])*len(sel_seq):\n            t = int(t)\n            print('Threshold:',t, 'Lefter than threshold - ', 'Train:', (l1>t).sum(), 'Test:', (l2>t).sum(), )\n    else:\n        for t in [70, 80, 100, 200]:\n            print('Threshold:',t, 'Greater than threshold - ', 'Train:', (l1>t).sum(), 'Test:', (l2>t).sum(), )\n\n            \nscores_analysis(list_scores1,list_scores2, similarity_name,  prot_name, sel_seq,    )            ","metadata":{"execution":{"iopub.status.busy":"2023-06-07T13:41:37.93942Z","iopub.execute_input":"2023-06-07T13:41:37.940164Z","iopub.status.idle":"2023-06-07T13:41:38.955383Z","shell.execute_reply.started":"2023-06-07T13:41:37.940126Z","shell.execute_reply":"2023-06-07T13:41:38.954269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ---------------------------------------------------------------------------------\n# Show stastics\n# --------------------------------------------------------------------------------\n\nd = pd.concat( (pd.Series(l1).describe(percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999]), pd.Series(l2).describe( percentiles = [0.25, 0.5, 0.75, 0.9,0.95,0.99,0.999] ) ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 20)\nplt.subplot(121)\nplt.hist(l1[l1<3*len(sel_seq)], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<3*len(sel_seq)], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\nfor t in [80, 100, 200]:\n    print('Threshold:',t, 'Greater than threshold - ', 'Train:', (l1>t).sum(), 'Test:', (l2>t).sum(), )\n    \n# ---------------------------------------------------------------------------------\n# Convert alignment scores to dataframes and save to csv     \n# Concat train and test, \n# Add description information (human readble info on proteins) for train part of data . Save to csv \n# --------------------------------------------------------------------------------\ndi1 = pd.DataFrame(index = list_ids[:len(l1)] , data = np.array( [l1, list_lens[:len(l1)]]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi2 = pd.DataFrame(index = list_ids2[:len(l2)] , data = np.array( [l2, list_lens2[:len(l2)]]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN_loc = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N_loc) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\n\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use \n\n\n# ---------------------------------------------------------------------------------\n# Show top scored proteins with human readable information to check that results are as expected:    \n#\n# --------------------------------------------------------------------------------\nprint()\nprint('Show top scored proteins with human readable information')\nN_to_show_top_scored_proteins_with_desription = 10\nfor N1 in [0,400,4000,40000]:\n    if N1+N_to_show_top_scored_proteins_with_desription < len(di):\n        print('Top ', N1, N1+N_to_show_top_scored_proteins_with_desription)\n        display(di.iloc[N1:(N1+N_to_show_top_scored_proteins_with_desription),:])\n\nprint()    \n    \n# ---------------------------------------------------------------------------------\nprint('Statists of the key description terms like \"kinase\" which are expected to be found in description of the  similar proteins  ' )\nprint('One expects - to find more such words for the top similar and lower for less similar ')\nprint('Better aligners are expected to have such curve higher than the worse ones ')\n# --------------------------------------------------------------------------------\nWL = window_len_for_keywords_indicator_moving_average # 50 # SRC # 1 P53\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # list_keywords_describing_protein = ['kinase', 'protein kinase']\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()        \n\nprint()    \n    \n# ---------------------------------------------------------------------------------\nprint('Compare Gene Ontology similarity (functional similarity) vs  structrual  similarity (alignment or embedding) ' )\nprint('One expects - structural similarity implies functional similarity, thus similar proteins to should have similar GO-terms')\nprint('Better aligners are expected to have such curve higher than the worse ones ')\n# --------------------------------------------------------------------------------\n    \n# Compute dice index between GO terms of the selected protein (sel_prot_id) and i-th similar to it protein.\n# Then create plots for that vector.\n# The code supports simultenous processing of the results from the several alignments stored in dict_scoring_data (calcualted above)\n# dict_scoring_data[k][0] - dataframe ordered by similarity , k-th index gives the protein ID of the k-th similar  \ndf_res = pd.DataFrame()\nlist_similarities_to_compare = list(dict_scoring_data.keys() )\nfor i1,k in enumerate( list_similarities_to_compare  ):\n    dt = dict_scoring_data[k][0]\n    s0 = set( map_protein_id_to_GOterms[ sel_prot_id ] )\n#   print(k, len(s0)); display(dt.head(2));  \n    l = []\n    for nn in range(len(dt)):\n        s = set( map_protein_id_to_GOterms[ dt.index[nn] ] )\n        dice = 2 * len( s & s0 ) / ( len(s) + len(s0) )\n        l.append(dice)\n    df_res[k] = l\ndisplay( df_res.head(2) )\ndisplay( df_res.corr() )\nprint('Proteins count:', len(df_res) )\n\nfor WL in [50,5]:\n    for N in [100,1000,10000,100_000]:\n        plt.figure(figsize = (20,6))\n        plt.title('Dice between GO-terms for ' + prot_name + ' vs i-th similar protein. Cutoff: '+str(N) + ' WL '+str(WL), fontsize = 18 )\n        for col in df_res.columns:\n            plt.plot(df_res[col].iloc[:N].rolling(window=WL).mean(), label = col)\n        plt.legend()\n        plt.grid()\n        plt.show()","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":"","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-07T13:04:08.557859Z","iopub.execute_input":"2023-06-07T13:04:08.558239Z","iopub.status.idle":"2023-06-07T13:04:08.564473Z","shell.execute_reply.started":"2023-06-07T13:04:08.558208Z","shell.execute_reply":"2023-06-07T13:04:08.563067Z"},"trusted":true},"execution_count":null,"outputs":[]}]}