{"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\nHere we give some examples how to work with these packages and make speed-tests.\nAlignments can be used for CAFA5 like tasks as discussed below. \n\nCurrently we consider Parasail, Diamond, BioPython, Skbio, Pyalign and for the comparison we also add Levenshtein distance. \nWe hope to add other aligners to the benchmark - please suggest ones which can work on Kaggle. \nWe benchmark on the set of random sequences as well as on the real protein sequences from CAFA5 challenge. \nThe number of sequences is a parameter in the notebook, we consider cases with 10_000, 100_000, 142246 sequences (the entire train of CAFA5). \n\n### Current outcome : use Parasail aligner for \"one vs many\", and Diamond for \"many vs many\" (with quality sacrifice)\n\n#### Parasail \n\n Parasail is 40-70 times faster than BioPython alignment (depending on sequences) and gives exactly the same results (for considered cases - Smith-Waterman (local), use blossum62, etc). It is about 7 faster than skbio and only 3-4 times slower than Levenshtein distance (not an alignment, but similar and simpler). \n\nLinks: https://github.com/jeffdaily/parasail#table-of-contents , Daily, Jeff. (2016). Parasail: SIMD C library for global, semi-global, and local pairwise sequence alignments. BMC Bioinformatics, 17(1), 1-11. http://dx.doi.org/10.1186/s12859-016-0930-z\n\n#### Parasail vs. Diamond\n\n For the task \"one vs many\" considered here - Parasail seems to be the best choice, because\n despite Diamond \"--fast\" is faster than Parasail about 2-2.5 times, \n but \"diamond--fast\" outputs at least 20-50 times less similar proteins than it is expected\n (in particular 20-30 less then the same diamond, but with the option \"----ultra-sensitive\" - see tests below).\n \n So, to conclude, since speed-up 2 times, with loss of quality 20-50 times, does not seems to be a good choice.\n We would prefer Parasail for the task \"one vs many\".\n \n As we discuss below \"diamond\" is much better for the case \"many vs many\". \n\n Note: the speed of diamond with the next after \"--fast\" option: \"--mid-sensitive\" is already about 1.5-2 slower than Parasail (for \"one vs many\").\n Remark: strangely enough the option \" --very-sensitive\" is better in many cases both in speed and quality. But still it is 1.2-1.5 slower than Parasail with less quality. \n\n Remark: The fact that diamond has worse quality is not surprsing since blast (used in diamond) is a suboptimal algorithm, while Parasail uses optimal Smith-Waterman algorithm (which is more expensive computationally). \n\n\n#### Diamond - tradeoff speed/quality \n\n Diamond reported to be sometimes up to 20_000 times faster than blast (which is already fast by itself, and faster that Smith Waterman local alignment (however is suboptimal - quality is sacrificed)). \n It can align \"many x many\" e.g. (100_000 x 100_000) proteins in 10 minutes (fastest mode). It is nearly impossible to do such things with other packages. Also it supports clustering mode.\n However one should note that as any blast-like algorithm it is suboptimal, comparing to optimal Smith-Waterman local alignment algorithm.\n Diamond always returns NOT the full alignment results, but only part which is \"well-aligned\" from diamond's point of view.\n Thus there is possible lose of quality , if diamond's \"well-aligned\" does not covery full \"well-alignable\", which is indeed the case.\n To control that issue there are several parameters, most important is \"sensitivity\".\n It varies from \"--fast\" to \"--ultra-sensitive\".\n Pay attention that it affects the execution time quite strongly in our experiments we see 20 times difference between \"--fast\" and \"--ultra-sensitive\".\n\n In our tests (one vs many):\n Parasail is SLOWER than Diamond--fast-mode about 2-2.5 times. For the price of quality drop 20-50 times. \n Parasail is FASTER than the next mode: Diamond--mid-sensitive about 1.5 times. The quality drop is much less: 2-5 times. \n However even most slow mode: Diamond--ultra-sensitive (4 times slower than Parasail) returns LESS similar proteins then really are and can be found by Parasail, according to some our estimates which we will discuss in the next notebook. (Not suprising, since blast is suboptimal, while Smith-Waterman is optimal). \n \n Diamond's tradeoff Time/Quality (on 40000 CAFA5):\n \n| Mode | Time | Count Total Returned | Count Returned e-value < 0.001 |\n| --- | --- | --- | --- |\n|diamond --fast | 1.99| 45.0 |\t\t\t\t\t45.0|\n|diamond --mid-sensitive | 7.31|\t491.0|\t\t\t\t\t472.0|\n|diamond --sensitive |13.46|721.0 |\t\t\t\t\t690.0|\n|diamond --more-sensitive |13.00\t|825.0|\t\t\t\t\t793.0|\n|diamond --very-sensitiv |5.14\t|1193.0|\t\t\t\t\t1060.0|\n|diamond --ultra-sensitive | 19.93\t|1341.0|\t\t\t\t\t1140.0|\n\n Note: strangely enough the mode \" --very-sensitive\" seems to work very fast and produce good results, so it is better in both respects than \"--mid-sensitive\", \" --sensitive\", \"--more-sensitive\". \n \n Another point - what about diamons's precision in opposite to sensitivity: as we see above - diamond is not perfectly sensitive. Our tests suggests that it is quite precise, i.e. for completely random words - it does not return results or return them with very high e-value - more than 1, and so much more than standard threshold 0.001. In other words we do not see false positive in that experiment even in \"ultra-sensitive\"-mode, so we conclude that precision of diamond is quite goood. \n \n Let us emphasize again that the main advantage of Diamond is ability to process \"many by many\" alignments in short time, which would be nearly impossible with other aligners. We plan to discuss it in a separate notebook. \n\nThanks to Liza Geraseva's notebook on diamond: https://www.kaggle.com/code/geraseva/diamond - please upvote !\n \nLinks: https://github.com/bbuchfink/diamondhttps://github.com/bbuchfink/diamond , Papers: Buchfink B, Reuter K, Drost HG, \"Sensitive protein alignments at tree-of-life scale using DIAMOND\", Nature Methods 18, 366–368 (2021). doi:10.1038/s41592-021-01101-x . Clustering: Buchfink B, Ashkenazy H, Reuter K, Kennedy JA, Drost HG, \"Sensitive clustering of protein sequences at tree-of-life scale using DIAMOND DeepClust\", bioRxiv 2023.01.24.525373; doi: https://doi.org/10.1101/2023.01.24.525373 . Original: Buchfink B, Xie C, Huson DH, \"Fast and sensitive protein alignment using DIAMOND\", Nature Methods 12, 59-60 (2015). doi:10.1038/nmeth.3176\n\n\n#### Skbio, Levenshtein, Pyalign \n\n Skbio local SW aligner is about 7 times slower than Parasail and gives scores DIFFERENT from Parasail and BioPython - some checks suggest they are less reasonable: histogram of scores has strange shape, complete disagreement with Levenshtein distances, top scored proteins are less biologically reasonable. It might be the reason - that I have incorrectly formated Blossum62 matrix for skbio aligner - will investigate that in future. \n \n Levenshtein distance is somewhat similar to alignment scores (though not exactly the same). But its calculation is much faster.\n It is about 3-4 times faster than the Parasail alignment. However similarity obtained by it - is not that much good. We observe only about 10% agreement between topN for Levenshtein and Parasail lists in some cases. \n \n Pyalign - in our tests it is quite slower than Parasail, more or less similar to BioPython - sometimes 2 times faster, sometimes slower. \n Also, at the moment we were unable to find how to use substitution matrices like blossum62 and open/extention penalties - what we need for the protein alignments. Accodring to author adavantages of the package that it can work with large alphabets.\n \"Pyalign\" is new package. See the notebook from the package author how to install it on Kaggle: https://www.kaggle.com/code/lieblb/pyalign-example/notebook - please upvote. (The direct ways to install fails: https://www.kaggle.com/code/alexandervc/pyalign-install-problem-resolved). The github of the package is quite nice: https://github.com/poke1024/pyalign - contains lots of useful material, see also discussions e.g. here: https://www.biostars.org/p/9505903/ - that discussion mentions slower results from Parasail comparing to BioPython - which is very different from our results. \n\n\nIn the benchmark we mostly focus on use of local alignment (Smith-Waterman algorithm, or blast) with blossum62 substituition matrix,\nand gap open penalty 11, gap extension penalty 1. (Parameters used by blastP NCBI web server: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome ). \n\n#### Some further details on Parasail \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 Parameters used: blossum62 substitution matrix which is quite standard choice for the protein alignments; gap-open penalty 11, gap extension penalty 1 - all these params are used for the blastP at NCBI web-site: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome . Despite blast and Smith Waterman are not exactly the same algoithms, these parameters seems to produce quite resonable results for Smith Waterman local alignment we use here. \n\nLinks: \"Parasail\": https://github.com/jeffdaily/parasail-python Daily, 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\nSome discussion and comparison can be found at: https://www.biostars.org/p/9505903/\n\n\n### Welcome to suggest another packages (useable on Kaggle) for benchmark !\n\nThere are many alignment packages however not all of them possible/easy to install on Kaggle. Nice collection of links can be found: \nhttps://github.com/poke1024/pyalign . See also: https://github.com/danielecook/Awesome-Bioinformatics#pairwise and https://en.wikipedia.org/wiki/List_of_sequence_alignment_software The area is quite big and some aligners are specific to certain tasks like alignment of short reads used for NGS sequencing (not proteins) e.g. https://en.wikipedia.org/wiki/TopHat_(bioinformatics) , https://en.wikipedia.org/wiki/Bowtie_(sequence_analysis), https://bio-bwa.sourceforge.net/ (used for ATAC-seq, ChIP-seq, etc...). \n\n\nThere are some with GPU support : \n\nCUDASW++ - GPU accelerated Smith Waterman algorithm \nhttps://cudasw.sourceforge.net/homepage.htm#latest\n\nGASAL2: a GPU accelerated sequence alignment library for high-throughput NGS data - PMC\nhttps://www.ncbi.nlm.nih.gov/pmc/articles/PMC6815017/\n\nAnd many others. \n\nThere is very interesting project from Stanford NLP team: \"String2String\" https://github.com/stanfordnlp/string2string ,\nhowever it is kind of more NLP focused, and it seems it is not easy to insert e.g. blossum62 matrix in alignment which is important for bioinformatics protein alignment tasks.\n\n### What will happen if one uses Levenshtein distance or other method known to be NOT so good for protein alignment ? \n\nThe answer - loss of quality would be quite visible, it is not small. \n\nLocal protein alignment with substitution matrices blossum62, or pam, with penalty parameters 11 for gap open, 1 for gap extension - are typical choices for the problem. \n\nHowever it might be interesting to check what if we use other method for example Levenshtein distance, or discard blossum62, or set other penalties. \n\nTo get a first idea what happens we analyze intersection percent of sets of topN scored proteins obtained by one method and by another. \nIntersection nearly 100% would mean - near complete ageement. However we see that for e.g. Local/Global it is only 50%, for Local/Levenshtein goes down to 2% for large N - that means agreement is not that much good. \n\n\n### Biologically motivated visual inspection of the results. \n\nIn the present notebook we take SRC human protein as an example which would be aligned to CAFA5 proteins. \nSo we have biological expectation what should be on top of the list of the most similar proteins: \n\n 1) Most similar - homologs from very similar organisms e.g. for human - it is mouse, rat, etc...\n 2) homologs for more distant orgnisms - like birds, fishes (chicken, Dani Rario = zebrafish)\n 3) Members of the known family: e.g. for SRC: Yes, Fyn, and Fgr - SrcA subfamily, Lck, Hck, Blk, and Lyn - the SrcB subfamily, and Frk in its own subfamily. \n 4) protein tyrosine kinases in general - other tyrosine protein kinases \n 5) protein kinases in general - other kinases (not necessarily tyrosine, but Serine/threonine-protein kinases)\n\nWhen we take large enough number of CAFA5 proteins we will indeed see that. Set param: N_seqs_to_take_CAFA5 = 100_000 \n\nIn the next notebook we will give more details on that. \n\n### PS\n\nFor the competition purposes one can use \"similarity\" measures (alignment scores) in a several ways: find the most similar proteins to the given one and try to transfer labels from them to it (See e.g. Liza Geraseva notebook: https://www.kaggle.com/code/geraseva/diamond ). Or one can use \"similarities\" obained by alignment as features - i.e. a kind of feature engineering (similar to Levenshtein distance features 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\nFor gentle introduction to alignments with BioPython one can see: https://www.kaggle.com/code/shtrausslearning/biological-sequence-alignment \nand other Kaggle notebooks by the same author. \n\nOne can use alignments from BioPython: https://www.kaggle.com/code/alexandervc/cafa5-18-alignments-biopython-compare , some faster version from skbio: https://www.kaggle.com/code/alexandervc/cafa5-19-alignments-skbio , or very slow NCBI web-server way: https://www.kaggle.com/code/alexandervc/cafa5-20-ncbiwww-blast-biopython or very fast \"Diamond\": https://www.kaggle.com/code/geraseva/diamond However for any option there are certain concerns - either speed or quality.\n \n \n ### Versions\n \n V3 - sizes of the tests sets for alignment - 142_246 sequences Execution time: 3 hours \n V2- sizes of the tests sets for alignment - 100_000 sequences. Execution time: 2 hours\n V1 - sizes of the tests sets for alignment - 10_000 sequences. Execution time: 20 minutes\n","metadata":{}},{"cell_type":"markdown","source":"# Key params","metadata":{}},{"cell_type":"code","source":"# Params to control the number of sequences to be aligned.\n# If we take too many - notebook can execute quite a long time.\n# For 1000 - it takes just about 1.5 minute . The bottleneck is the BioPython protein alignment - its executes 10-100 times slowly than otehr methods. \n# For 10_000 - about 6 minutes, for 100_000 - about 35 minutes. See version of the notebooks - different params explored. \n\nN_seqs_to_take_for_random_sequences = 142246# 100_000 # How many random sequences to generate. \n\nN_seqs_to_take_CAFA5 = N_seqs_to_take_for_random_sequences # 142246# 100_000 # How many proteins to take from the CAFA5 train data. \n# 142246 - full size of train CAFAT5 \n\nlist_aligners = [ 'BioPython Local Blossum62 11_1', 'Parasail SW Local Blossum62 11_1', 'skbio Local Blossum62 11_1' , 'Levenshtein' , 'diamond', \n 'pyalign Local', 'pyalign Global', 'BioPython Global Blossum62 11_1', 'BioPython Global Default', 'Parasail NW Global Blossum62 11_1',]\n\n# Choose list parameters to test for diamond sensitivity \nlist_diamond_sensititivity_modes = ['--fast','--mid-sensitive','--sensitive','--more-sensitive','--very-sensitive','--ultra-sensitive' ]\n#Full list: ['--fast','--mid-sensitive','--sensitive','--more-sensitive','--very-sensitive','--ultra-sensitive' ]\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:15.541424Z","iopub.execute_input":"2023-06-09T15:08:15.541858Z","iopub.status.idle":"2023-06-09T15:08:15.578652Z","shell.execute_reply.started":"2023-06-09T15:08:15.541823Z","shell.execute_reply":"2023-06-09T15:08:15.577728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Results of all benchamrks will be save here: \n\nimport pandas as pd \n\ndict_align_results = {}\ndf_stat = pd.DataFrame()\ndf_stat_diamond = pd.DataFrame()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:15.580326Z","iopub.execute_input":"2023-06-09T15:08:15.580842Z","iopub.status.idle":"2023-06-09T15:08:15.591053Z","shell.execute_reply.started":"2023-06-09T15:08:15.580811Z","shell.execute_reply":"2023-06-09T15:08:15.590051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show previously computed benchmarks results","metadata":{}},{"cell_type":"code","source":"# Let show an example of the the benchmark generated by the present notebook - precomputed and stored. \n# Parameters might be different from the current version. \n\nimport pandas as pd\npd.set_option('max_colwidth', 500)\npd.set_option('display.max_rows', 500)\n\nf1 = '/kaggle/input/cafa5-features-etc/df_aligners_benchmark.csv'\ndf_stat_previous = pd.read_csv(f1, index_col = 0 )\ndisplay( df_stat_previous )\n\nf2 = '/kaggle/input/cafa5-features-etc/df_diamond_stat.csv'\ndf_stat_diamond_previous = pd.read_csv(f2, index_col = 0 )\ndisplay( df_stat_diamond_previous )\n\nf3 = '/kaggle/input/cafa5-features-etc/df_scores_with_descriptions.csv'\ndf_scores_ext_previous = pd.read_csv(f3, index_col = 0 )\ndisplay( df_scores_ext_previous.head(5) )\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:15.592609Z","iopub.execute_input":"2023-06-09T15:08:15.593048Z","iopub.status.idle":"2023-06-09T15:08:15.779607Z","shell.execute_reply.started":"2023-06-09T15:08:15.593007Z","shell.execute_reply":"2023-06-09T15:08:15.778163Z"},"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\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-09T15:08:15.782228Z","iopub.execute_input":"2023-06-09T15:08:15.782681Z","iopub.status.idle":"2023-06-09T15:08:16.735126Z","shell.execute_reply.started":"2023-06-09T15:08:15.782649Z","shell.execute_reply":"2023-06-09T15:08:16.733982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Parasail aligner toy example \n\nLinks: \"Parasail\": https://github.com/jeffdaily/parasail-python Daily, 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\nSome discussion and comparison can be found at: https://www.biostars.org/p/9505903/\n","metadata":{}},{"cell_type":"code","source":"%%time \nif ('Parasail SW Local Blossum62 11_1' in list_aligners) or ( 'Parasail NW Local Blossum62 11_1' in list_aligners): # \n !pip install parasail\n import parasail\n\n print('Local score:')\n result = parasail.sw_scan_16(\"ARNDCQEGHILKMFPSTWYVBZX*\", \"AAAA\" , 11, 1, parasail.blosum62)\n print(result.score)\n print()\n \n print('Global score:')\n result = parasail.nw_scan_16(\"ARNDCQEGHILKMFPSTWYVBZX*\", \"AAAA\" , 11, 1, parasail.blosum62)\n print(result.score)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:16.736708Z","iopub.execute_input":"2023-06-09T15:08:16.737497Z","iopub.status.idle":"2023-06-09T15:08:33.284839Z","shell.execute_reply.started":"2023-06-09T15:08:16.737464Z","shell.execute_reply":"2023-06-09T15:08:33.283333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Warning - be careful (do not use) \"_8\" mode of parasail - it truncates scores to 255 = 2^8 or to 0 . \n\nSimilar problem might appear for very long proteins in \"_16\" mode - if scores increase 2^16 = 65536 - not clear how to resolve, but it seems not relevant for most of the proteins, since they are not that much long ","metadata":{}},{"cell_type":"markdown","source":"### In \"_stats\" mode truncation goes to 255","metadata":{}},{"cell_type":"code","source":"%%time \nif ('Parasail SW Local Blossum62 11_1' in list_aligners) or ( 'Parasail NW Local Blossum62 11_1' in list_aligners): # \n import parasail\n print(\"_16 mode works okay: \")\n result = parasail.sw_scan_16(\"asdf\"*100, \"asdf\"*100, 11, 1, parasail.blosum62)\n print(result.score)\n result = parasail.sw_scan_16(\"eklmn\"*100, \"klmnxyz\"*100, 11, 1, parasail.pam100)\n print(result.score)\n print()\n\n print('WARNING(!) the score truncation for \"_8\" mode - different results are truncated to the same 255 ')\n result = parasail.sw_stats_striped_8(\"asdf\"*100, \"asdf\"*100, 11, 1, parasail.pam100)\n print('WARNING(!) the score truncation for \"_8\" mode: ')\n print(result.score)\n result = parasail.sw_stats_striped_8(\"eklmn\"*100, \"klmnxyz\"*100, 11, 1, parasail.pam100)\n print('WARNING(!) the score truncation for \"_8\" mode: ')\n print(result.score)","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:33.286747Z","iopub.execute_input":"2023-06-09T15:08:33.287224Z","iopub.status.idle":"2023-06-09T15:08:33.299644Z","shell.execute_reply.started":"2023-06-09T15:08:33.287175Z","shell.execute_reply":"2023-06-09T15:08:33.298318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### In \"_scan\" mode truncation goes to 0","metadata":{}},{"cell_type":"code","source":"%%time\nif ('Parasail SW Local Blossum62 11_1' in list_aligners) or ( 'Parasail NW Local Blossum62 11_1' in list_aligners): # \n print('WARNING(!) the score truncation for \"_8\" mode - different results are truncated to the same 255 ')\n result = parasail.sw_scan_8(\"asdf\"*100, \"asdf\"*100, 11, 1, parasail.pam100)\n print('WARNING(!) the score truncation for \"_8\" mode: ')\n print(result.score)\n result = parasail.sw_scan_8(\"eklmn\"*100, \"klmnxyz\"*100, 11, 1, parasail.pam100)\n print('WARNING(!) the score truncation for \"_8\" mode: ')\n print(result.score)","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:33.301406Z","iopub.execute_input":"2023-06-09T15:08:33.30178Z","iopub.status.idle":"2023-06-09T15:08:33.316525Z","shell.execute_reply.started":"2023-06-09T15:08:33.30175Z","shell.execute_reply":"2023-06-09T15:08:33.31512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ByoPython aligner toy example \n\nSee more examples in https://www.kaggle.com/code/alexandervc/cafa5-18-alignments-biopython-compare\n\nAnd very nice tutorial(s) by ANDREY SHTRAUSS : https://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#5-|-PAIRWISE-SEQUENCE-ALIGNMENT (and other notebooks by the same author). (However, pay attention that syntax used there refers to older-depricated and slower BioPython aligners).\n","metadata":{}},{"cell_type":"code","source":"%%time\nif ('BioPython Local Blossum62 11_1' in list_aligners) or ('BioPython Global Blossum62 11_1' in list_aligners):\n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n print('Local ')\n aligner = Align.PairwiseAligner()\n aligner.mode = 'local'\n aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n aligner.open_gap_score = -11\n aligner.extend_gap_score = -1\n aligner.target_end_gap_score = -11\n aligner.query_end_gap_score = -1\n\n\n result2 = aligner.align(\"ARNDCQEGHILKMFPSTWYVBZX*\", \"AAAA\" )\n print(result2[0])\n print(result2[0].score)\n print()\n \n print('Global ')\n aligner = Align.PairwiseAligner()\n aligner.mode = 'global'\n aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n aligner.open_gap_score = -11\n aligner.extend_gap_score = -1\n aligner.target_end_gap_score = -11\n aligner.query_end_gap_score = -1\n\n\n result2 = aligner.align(\"ARNDCQEGHILKMFPSTWYVBZX*\", \"AAAA\" )\n print(result2[0])\n print(result2[0].score) ","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:33.318158Z","iopub.execute_input":"2023-06-09T15:08:33.318606Z","iopub.status.idle":"2023-06-09T15:08:33.335476Z","shell.execute_reply.started":"2023-06-09T15:08:33.318568Z","shell.execute_reply":"2023-06-09T15:08:33.33415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nif 'BioPython Global Default' in list_aligners: # \n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n print('Local ')\n aligner = Align.PairwiseAligner()\n aligner.mode = 'global'\n print(aligner)\n \n print()\n result2 = aligner.align(\"ARNDCQEGHILKMFPSTWYVBZX*\", \"AAAA\" )\n print(result2[0])\n print(result2[0].score) ","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:33.336656Z","iopub.execute_input":"2023-06-09T15:08:33.336967Z","iopub.status.idle":"2023-06-09T15:08:33.454998Z","shell.execute_reply.started":"2023-06-09T15:08:33.336941Z","shell.execute_reply":"2023-06-09T15:08:33.453782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# skbio aligner toy example (and install)","metadata":{}},{"cell_type":"code","source":"%%time\nif 'skbio Local Blossum62 11_1' in list_aligners: # \n \n !pip install scikit-bio\n\n #-----------------------------------------------------------------------------------------------\n print('Take blossum62 substitution matrix from the BioPython and convert it to format acceptable by skbio')\n print('unfortunately we do not see built-in blossum62 matrix in skbio, so need to prepare it')\n # see also https://www.kaggle.com/code/alexandervc/cafa5-19-alignments-skbio?scriptVersionId=131755045&cellId=6\n #-----------------------------------------------------------------------------------------------\n\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n names = substitution_matrices.load()\n print(list(names) )\n sm = substitution_matrices.load(\"BLOSUM62\") \n # sm = np.array(sm)\n # print(type(sm), sm.shape )\n # print(np.array(sm)[:3,:5])\n print( sm.alphabet )\n list(str(sm.alphabet))\n ll = list(str(sm.alphabet))\n dict_sm_blosum62 = {}\n for i in range(len(ll)):\n dict_tmp = {}\n for j in range(len(ll)):\n dict_tmp[ll[j]] = np.array(sm)[i,j]\n #print(dict_tmp)\n dict_sm_blosum62[ll[i]] = dict_tmp.copy()\n print( dict_sm_blosum62['A'] )\n\n #-----------------------------------------------------------------------------------------------\n # Compute example \n #-----------------------------------------------------------------------------------------------\n\n print()\n print()\n\n from skbio.alignment import StripedSmithWaterman\n seq1 = \"ARNDCQEGHILKMFPSTWYVBZX*\"\n seq2 = \"AAAA\" \n\n query = StripedSmithWaterman(seq1, substitution_matrix = dict_sm_blosum62, gap_open_penalty = 11, gap_extend_penalty = 1 )\n alignment = query(seq2)\n\n print('score:', alignment.optimal_alignment_score)\n print()\n print(alignment)\n print()\n print('aligned_query_sequence:', alignment.aligned_query_sequence, len(alignment.aligned_query_sequence))\n print('aligned_target_sequence:',alignment.aligned_target_sequence, len(alignment.aligned_target_sequence))\n t = alignment\n print( t['target_end_optimal'] - t['target_begin'], t['query_end'] - t['query_begin'] )\n display(alignment )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:08:33.459425Z","iopub.execute_input":"2023-06-09T15:08:33.459777Z","iopub.status.idle":"2023-06-09T15:10:20.034736Z","shell.execute_reply.started":"2023-06-09T15:08:33.459748Z","shell.execute_reply":"2023-06-09T15:10:20.033294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Levenshtein distance - toy example (and install) \n\nIt is NOT an aligner , still quite close. The point is that it can be calculated much faster than any alignment. \nAnd can serve as kind of weak analogue of global (not local) alignment score. So we include it as benchmark. \n\nhttps://en.wikipedia.org/wiki/Levenshtein_distance\n\nFor the GPU way - see Chris Deotte : https://www.kaggle.com/code/cdeotte/train-data-contains-mutations-like-test-data\n\nhttps://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","metadata":{}},{"cell_type":"code","source":"%%time\nif 'Levenshtein' in list_aligners: # \n !pip install python-Levenshtein\n from Levenshtein import distance\n edit_dist = distance(\"ah\", \"aho\")\n edit_dist","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:10:20.036571Z","iopub.execute_input":"2023-06-09T15:10:20.036925Z","iopub.status.idle":"2023-06-09T15:10:33.595133Z","shell.execute_reply.started":"2023-06-09T15:10:20.03689Z","shell.execute_reply":"2023-06-09T15:10:33.59359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Diamond toy example (and install)\n\nDiamond claimed to be sometimes 20_000 faster than blastP (which is already faster (but less precise) than Smith Waterman local alignment). \n\nIt can align \"many by many\" proteins. It supports clustering via special command. \nIt returns NOT all alignment results, but only those which somewhat \"well\"-aligned. \n\nWe follow approach of Liza Geraseva: https://www.kaggle.com/code/geraseva/diamond - please upvote ! \n\nIts interface is not Python simple (at least I do not know one). So even for the simple alignment of two sequences, we need to several steps:\n\n 1) we first need to create two fasta files - one contains just the first sequence and anoter just the second\n 2) then one fast is converted into internal format of diamon (\"database\" creation)\n 3) only then we launch aligner\n 4) need to load csv file with results\n \nPay attention on parameters:\n\n \"k\" - maximal possible number of returned results for each input protein \n \"evalue\" - admissible e-values - greater than threshold will not be returned. \n Sensitivity - # --fast # --mid-sensitive # '--sensitive' # '--more-sensitive' # --very-sensitive # '--ultra-sensitive' \n Sensitivity - quite affects the execution time from fast to very senstive it might grow up to 20 times \n while \"k\" and e-value does not seem to affect.\n ","metadata":{}},{"cell_type":"code","source":"%%time \nif 'diamond' in list_aligners:\n \n print('Diamond')\n !wget http://github.com/bbuchfink/diamond/releases/download/v2.1.6/diamond-linux64.tar.gz\n !tar xzf diamond-linux64.tar.gz\n !rm diamond-linux64.tar.gz\n\n #------------------------------------------------------------------------------------\n print()\n print('create two FASTA files for both data: what is aligned (query(ies)), and on-what is aligned (subject(s))')\n print()\n #------------------------------------------------------------------------------------\n\n # Examples should contain only protein letters and should NOT be very short - otherwise diamond will drop them out \n seq1 = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\n seq2 = seq1# 'MGSN'#'KS' \n\n fn1_query = \"fasta_for_diamond_test_1_query.txt\"\n ofile = open(fn1_query, \"w\")\n for i in range(1):\n ofile.write(\">\" + 'test_sequence1' + \"\\n\" +seq1 + \"\\n\")\n ofile.close()\n fn2_subject = \"fasta_for_diamond_test_2_subject.txt\"\n ofile = open(fn2_subject, \"w\")\n for i in range(1):\n ofile.write(\">\" + 'test_sequence2' + \"\\n\" +seq2 + \"\\n\")\n ofile.close()\n\n #------------------------------------------------------------------------------------\n print()\n print('---------------------------------------------------- ')\n print('Create a database file - on-what we align (subjects) ')\n print('---------------------------------------------------- ')\n print()\n #------------------------------------------------------------------------------------\n import time\n from subprocess import Popen, PIPE\n db_name='train_db'\n\n t0 = time.time()\n p = Popen(['./diamond', 'makedb', \n '--in',fn2_subject ,\n '-d', db_name], stdin=PIPE, stdout=PIPE)\n stdout, stderr = p.communicate()\n print('Database created in %.1f seconds'%(time.time() - t0 ))\n\n #------------------------------------------------------------------------------------------------\n print()\n print('---------------------------------------------------- ')\n print('Run a blastp-like search')\n print('---------------------------------------------------- ')\n print()\n #------------------------------------------------------------------------------------------------\n k = 15 # 142050 # maximal number of \"similar\" proteins for each query protein to be returned \n evalue_threshold = 100_000_000 # default=0.001 # but we want to get as much as possible # controls how similar sequences will be returned \n str_sensitivity = '--ultra-sensitive' # --fast # --mid-sensitive # '--sensitive' # '--more-sensitive' # --very-sensitive # '--ultra-sensitive'\n\n outfile_name = 'matches.tsv'\n time0 = time.time() \n p = Popen(['./diamond', 'blastp', '-d', db_name, '--evalue',str(evalue_threshold),\n '-q', fn1_query,\n '-o', outfile_name, '--max-target-seqs', str(k), # '--quiet', str_sensitivity\n ], stdin=PIPE, stdout=PIPE) # '--quiet' # # \n stdout, stderr = p.communicate()\n print(f'Execution time: {time.time()-time0}s')\n\n # Cannot make work additional output \n # '-outfmt 6 qseqid sseqid',\n # diamond blastp --query query.fasta --db database.dmnd --outfmt 6 qseqid sseqid qlen slen\n # '--outfmt 6 qseqid sseqid qlen slen pident length mismatch gapopen qstart qend sstart send evalue bitscore' \n\n #------------------------------------------------------------------------------------------------\n print()\n print('---------------------------------------------------- ')\n print('Load results of alignment in pandas dataframe')\n print('---------------------------------------------------- '); \n print()\n #------------------------------------------------------------------------------------------------\n\n matches=pd.read_csv(outfile_name, sep='\\t', header=None, \n names=['qseqid', 'sseqid', 'pident', 'length', 'mismatch', \n 'gapopen', 'qstart', 'qend', 'sstart','send', 'evalue', 'bitscore'])\n matches['qseqid']=matches['qseqid'].apply(lambda x: x.split('\\\\t')[0])\n print(matches.shape)\n print('Bitscore:', matches['bitscore'].iat[0])\n matches.head(10)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:10:33.597539Z","iopub.execute_input":"2023-06-09T15:10:33.598047Z","iopub.status.idle":"2023-06-09T15:10:39.313906Z","shell.execute_reply.started":"2023-06-09T15:10:33.597997Z","shell.execute_reply":"2023-06-09T15:10:39.312338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pyalign toy example (and install)\n\n\"Pyalign\" is a new package. \nPyalign - in our tests it is quite slower than Parasail, more or less similar to BioPython - sometimes 2 times faster, sometimes slower. \nAlso, at the moment we were unable to find how to use substitution matrices like blossum62 and open/extention penalties - what we need for the protein alignments. Accodring to author adavantages of the package that it can work with large alphabets.\nSee the notebook from the package author how to install it on Kaggle: https://www.kaggle.com/code/lieblb/pyalign-example/notebook - please upvote. (The direct ways to install fails: https://www.kaggle.com/code/alexandervc/pyalign-install-problem-resolved). \n\nThe github of the package is quite nice: https://github.com/poke1024/pyalign , see also discussions e.g. here: https://www.biostars.org/p/9505903/ - that discussion mentions slower results from Parasail comparing to BioPython - which is very different to our results.\n\n\nNote that installing and compilation may take up to 8 minutes on Kaggle. As the author exlpained: \"installing via pip seems to fail because the xtensor libraries are missing\"\n\n","metadata":{}},{"cell_type":"code","source":"%%time \nif ('pyalign Local' in list_aligners) or ('pyalign Global' in list_aligners ):\n \n !git clone https://github.com/xtensor-stack/xtl && \\\n cd xtl && cmake . && make install && cd .. && rm -r xtl\n\n !git clone https://github.com/xtensor-stack/xtensor && \\\n cd xtensor && cmake . && make install && cd .. && rm -r xtensor \n\n !git clone https://github.com/xtensor-stack/xsimd && \\\n cd xsimd && cmake . && make install && cd .. && rm -r xsimd \n\n !git clone https://github.com/xtensor-stack/xtensor-python && \\\n cd xtensor-python && cmake -Dpybind11_DIR=`pybind11-config --cmakedir` . && \\\n make install && cd .. && rm -r xtensor-python \n\n !pip install pyalign \n\n import pyalign\n\n alignment = pyalign.global_alignment(\"INDUSTRY\", \"INTEREST\", gap_cost=0, eq=1, ne=-1)\n print(alignment.score)\n display(alignment)\n alignment = pyalign.local_alignment(\"INDUSTRY\", \"INTEREST\", gap_cost=0, eq=1, ne=-1)\n print(alignment.score)\n display(alignment) \n ","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:10:39.31621Z","iopub.execute_input":"2023-06-09T15:10:39.316734Z","iopub.status.idle":"2023-06-09T15:18:42.662162Z","shell.execute_reply.started":"2023-06-09T15:10:39.316682Z","shell.execute_reply":"2023-06-09T15:18:42.660793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Benchmark on random data\n\nHere we take one sequence (\"selected sequence\") and generate many (N_seqs_to_take_for_random_sequences) random shuffles of it. \nThen make alignments of the selected sequence to all randomly shuffled. \n\nCurrently as a selected sequence we take SRC human gene - length 536. But probably it is not important what to take. ","metadata":{}},{"cell_type":"markdown","source":"## Create data by random permutation of a given string","metadata":{}},{"cell_type":"code","source":"%%time\n\nsel_seq = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\n# Sequence for SRC human protein \n# https://en.wikipedia.org/wiki/Proto-oncogene_tyrosine-protein_kinase_Src\n# Uniprot: https://www.uniprot.org/uniprotkb/P12931/entry\n \nprint(len(sel_seq))\nlist_seqs_random = []\nfor i in range(N_seqs_to_take_for_random_sequences):\n characters = list(sel_seq)\n np.random.shuffle(characters)\n s = ''.join(characters ) \n list_seqs_random.append( s )\n \nprint(len(list_seqs_random), list_seqs_random[0][:10])\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:18:42.664188Z","iopub.execute_input":"2023-06-09T15:18:42.665075Z","iopub.status.idle":"2023-06-09T15:18:43.133391Z","shell.execute_reply.started":"2023-06-09T15:18:42.665033Z","shell.execute_reply":"2023-06-09T15:18:43.132437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Aligners speed test","metadata":{}},{"cell_type":"code","source":"%%time \nimport time \nstr_info = 'random len'+str(len(sel_seq))\n\nverbose = 1\nN_seqs_to_take_loc = len( list_seqs_random )\nlist_seqs = list_seqs_random\n\n\n# ---------------------------------------------------------------------------------\n# BioPython aligner(s) \n# --------------------------------------------------------------------------------\nsimilarity_name = 'BioPython Local Blossum62 11_1'\nif similarity_name in list_aligners: # \n print( similarity_name )\n\n\n # ---------------------------------------------------------------------------------\n # Prepare BioPython aligner \n # --------------------------------------------------------------------------------\n\n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n aligner = Align.PairwiseAligner()\n aligner.mode = 'local'\n aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n aligner.open_gap_score = -11\n aligner.extend_gap_score = -1\n aligner.target_end_gap_score = -11\n aligner.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\n if verbose >= 100:\n print(aligner)\n\n # ---------------------------------------------------------------------------------\n # Compute similarity BioPython \n # ---------------------------------------------------------------------------------\n\n t0 = time.time()\n list_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs ]\n tm = np.round(time.time()-t0,2)\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(time.time()-t0, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\nsimilarity_name = 'BioPython Global Blossum62 11_1'\nif similarity_name in list_aligners: # \n print( similarity_name )\n\n # ---------------------------------------------------------------------------------\n # Prepare BioPython aligner \n # --------------------------------------------------------------------------------\n\n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n aligner = Align.PairwiseAligner()\n aligner.mode = 'global'\n aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n aligner.open_gap_score = -11\n aligner.extend_gap_score = -1\n aligner.target_end_gap_score = -11\n aligner.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\n if verbose >= 100:\n print(aligner)\n\n # ---------------------------------------------------------------------------------\n # Compute similarity BioPython \n # ---------------------------------------------------------------------------------\n\n t0 = time.time()\n list_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs ]\n tm = np.round(time.time()-t0,2)\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(time.time()-t0, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\nsimilarity_name = 'BioPython Global Default'\nif similarity_name in list_aligners: # \n print( similarity_name )\n\n # ---------------------------------------------------------------------------------\n # Prepare BioPython aligner \n # --------------------------------------------------------------------------------\n\n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n aligner = Align.PairwiseAligner()\n aligner.mode = 'global'\n# aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n# aligner.open_gap_score = -11\n# aligner.extend_gap_score = -1\n# aligner.target_end_gap_score = -11\n# aligner.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\n if verbose >= 100:\n print(aligner)\n\n # ---------------------------------------------------------------------------------\n # Compute similarity BioPython \n # ---------------------------------------------------------------------------------\n\n t0 = time.time()\n list_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs ]\n tm = np.round(time.time()-t0,2)\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(time.time()-t0, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n \n# ---------------------------------------------------------------------------------\n# Parasail Aligner \n# ---------------------------------------------------------------------------------\nsimilarity_name = 'Parasail SW Local Blossum62 11_1'\nif similarity_name in list_aligners: # \n print()\n print( similarity_name )\n t0 = time.time()\n list_scores1 = [parasail.sw_scan_16(sel_seq, seq , 11, 1, parasail.blosum62).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n \nsimilarity_name = 'Parasail NW Global Blossum62 11_1'\nif similarity_name in list_aligners: # \n print()\n print( similarity_name )\n t0 = time.time()\n list_scores1 = [parasail.nw_scan_16(sel_seq, seq , 11, 1, parasail.blosum62).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n\n# ---------------------------------------------------------------------------------\n# skbio Local Blossum62\n# ---------------------------------------------------------------------------------\nsimilarity_name = 'skbio Local Blossum62 11_1'\nif similarity_name in list_aligners: # \n print()\n print( similarity_name )\n query = StripedSmithWaterman(sel_seq, substitution_matrix = dict_sm_blosum62, gap_open_penalty = 11, gap_extend_penalty = 1 )\n t0 = time.time()\n list_scores1 = [query(seq)['optimal_alignment_score'] for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n\n# ---------------------------------------------------------------------------------\n# Levenshtein\n# --------------------------------------------------------------------------------\nif 'Levenshtein' in list_aligners: # \n print()\n similarity_name = 'Levenshtein'\n print( similarity_name )\n t0 = time.time()\n list_scores1 = [ distance(sel_seq, seq) for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n# ---------------------------------------------------------------------------------\n# Diamond\n# --------------------------------------------------------------------------------\nif 'diamond' in list_aligners:\n similarity_name = 'diamond'\n print(); print('diamond')\n #------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamnod. Create two FASTA files for both data: what is aligned (query(ies)), and on-what is aligned (subject(s))')\n print()\n #------------------------------------------------------------------------------------\n\n fn1_query = \"fasta_for_diamond_test_1_query.txt\"\n ofile = open(fn1_query, \"w\")\n for i in range(1):\n ofile.write(\">\" + 'test_sequence1' + \"\\n\" +sel_seq + \"\\n\")\n ofile.close()\n fn2_subject = \"fasta_for_diamond_test_2_subject.txt\"\n ofile = open(fn2_subject, \"w\")\n list_ids_for_random = ['test_sequence2_'+str(i) for i in range(len(list_seqs))]\n for i in range(len(list_seqs)):\n ofile.write(\">\" + list_ids_for_random[i] + \"\\n\" +list_seqs[i] + \"\\n\")\n ofile.close()\n\n #------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamond. Create a database file - on-what we align (subjects) ')\n print()\n #------------------------------------------------------------------------------------\n import time\n from subprocess import Popen, PIPE\n db_name='train_db'\n\n t0 = time.time()\n p = Popen(['./diamond', 'makedb', \n '--in',fn2_subject ,\n '-d', db_name], stdin=PIPE, stdout=PIPE)\n stdout, stderr = p.communicate()\n print('Database created in %.1f seconds'%(time.time() - t0 ))\n\n #------------------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamond. Run a blastp-like search')\n print()\n #------------------------------------------------------------------------------------------------\n for str_sensitivity in list_diamond_sensititivity_modes:\n #str_sensitivity = '--ultra-sensitive' # --fast # --mid-sensitive # '--sensitive' # '--more-sensitive' # --very-sensitive # '--ultra-sensitive'\n k = len(list_seqs) # maximal number of \"similar\" proteins for each query protein to be returned \n evalue_threshold = 100_000_000 # default=0.001 # but we want to get as much as possible # controls how similar sequences will be returned \n\n outfile_name = 'matches.tsv'\n t0 = time.time() \n p = Popen(['./diamond', 'blastp', '-d', db_name, '--evalue',str(evalue_threshold),\n '-q', fn1_query,\n '-o', outfile_name, '--max-target-seqs', str(k), '--quiet', str_sensitivity,\n ], stdin=PIPE, stdout=PIPE) # '--quiet' # # \n stdout, stderr = p.communicate()\n tm2 = np.round(time.time()-t0,2)\n print(f'Execution time: {time.time()-t0}s')\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_seqs), len(sel_seq) ))\n\n\n # Cannot make work additional output \n # '-outfmt 6 qseqid sseqid',\n # diamond blastp --query query.fasta --db database.dmnd --outfmt 6 qseqid sseqid qlen slen\n # '--outfmt 6 qseqid sseqid qlen slen pident length mismatch gapopen qstart qend sstart send evalue bitscore' \n\n #------------------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamond. Load results of alignment in pandas dataframe')\n print()\n #------------------------------------------------------------------------------------------------\n\n matches=pd.read_csv(outfile_name, sep='\\t', header=None, \n names=['qseqid', 'sseqid', 'pident', 'length', 'mismatch', \n 'gapopen', 'qstart', 'qend', 'sstart','send', 'evalue', 'bitscore'])\n matches['qseqid']=matches['qseqid'].apply(lambda x: x.split('\\\\t')[0])\n print(str_sensitivity , 'Diamond returned: ', len(matches), 'results. Others will be filled by zero score')\n print(str_sensitivity , 'Diamond returned with e-value < 0.001: ', (matches['evalue']<0.001).sum(), 'results. Others will be filled by zero score')\n str_tmp = ' ' + str(len(list_seqs)) + ' ' + str_info\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity, 'Time on'+str_tmp] = tm\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity, 'Count Returned'+str_tmp] = len(matches)\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity,'Count Returned e-value < 0.001'+str_tmp]=(matches['evalue']<0.001).sum()\n if len(matches) > 0:\n print('Min evalue:', matches['evalue'].min(), 'Max evalue:', matches['evalue'].max() , '1percent evalue:', np.percentile( matches['evalue'], 0.01))\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity,'Min e-value'+str_tmp]=matches['evalue'].min()\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity,'1%% Percentile e-value'+str_tmp]=np.percentile( matches['evalue'], 0.01)\n \n\n # \n if 1: # Diamond returns only aligned sequences, we need to get list_scores in the original order with 0 for scores for those missed\n df_tmp = pd.DataFrame(index = list_ids_for_random[:len(list_seqs)] , data = range(len(list_seqs)) )\n df_tmp = df_tmp.join( matches[['sseqid', 'bitscore']].set_index('sseqid'), how = 'left')\n df_tmp = df_tmp.sort_values(0) # Just in case sorting was lost during join \n df_tmp = df_tmp.fillna(0)\n list_scores1 = df_tmp['bitscore'].values # Diamond returns only aligned sequences, we need to get list_scores in the original order with 0 for scores for those missed\n\n df_stat.loc[similarity_name + ' ' + str_sensitivity, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_sensitivity + ' ' + str_info] = np.asarray(list_scores1.copy() )\n \n if verbose >= 100:\n print()\n print('---------------------------------------------------- ')\n print(' diamond finished ')\n print('---------------------------------------------------- '); \n print()\n print()\n print()\n \n# ---------------------------------------------------------------------------------\n# pyalign Local\n# --------------------------------------------------------------------------------\nsimilarity_name = 'pyalign Local'\nif similarity_name in list_aligners:\n import pyalign\n # alignment = pyalign.local_alignment(\"INDUSTRY\", \"INTEREST\", gap_cost=0, eq=1, ne=-1)\n print()\n \n print( similarity_name )\n t0 = time.time()\n list_scores1 = [ pyalign.local_alignment(sel_seq,seq, gap_cost=0, eq=1, ne=-1).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n# ---------------------------------------------------------------------------------\n# pyalign Global\n# --------------------------------------------------------------------------------\nsimilarity_name = 'pyalign Global'\nif similarity_name in list_aligners:\n import pyalign\n # alignment = pyalign.local_alignment(\"INDUSTRY\", \"INTEREST\", gap_cost=0, eq=1, ne=-1)\n print()\n \n print( similarity_name )\n t0 = time.time()\n list_scores1 = [ pyalign.global_alignment(sel_seq,seq, gap_cost=0, eq=1, ne=-1).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n \n \n \n \n \n \n \n \n \n \n# ---------------------------------------------------------------------------------\n# Compare results \n# --------------------------------------------------------------------------------\nprint(); print();\ndf_scores = pd.DataFrame( dict_align_results )\nprint('Correlation matrix between obtained scores (expected to be all 1, but it is not: )')\ndisplay(df_scores.corr() )\nprint('Spearman correlations')\ndisplay(df_scores.corr(method = 'spearman') )\n\nprint()\ndisplay(df_stat)","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:18:43.135078Z","iopub.execute_input":"2023-06-09T15:18:43.135733Z","iopub.status.idle":"2023-06-09T15:24:27.242719Z","shell.execute_reply.started":"2023-06-09T15:18:43.1357Z","shell.execute_reply":"2023-06-09T15:24:27.241451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('statistics for diamond')\ndf_stat_diamond.to_csv('df_diamond_stat.csv')\ndisplay(df_stat_diamond) \n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:24:27.244474Z","iopub.execute_input":"2023-06-09T15:24:27.244859Z","iopub.status.idle":"2023-06-09T15:24:27.271126Z","shell.execute_reply.started":"2023-06-09T15:24:27.244825Z","shell.execute_reply":"2023-06-09T15:24:27.269625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Speed test on random data:')\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:24:27.273035Z","iopub.execute_input":"2023-06-09T15:24:27.273426Z","iopub.status.idle":"2023-06-09T15:24:27.286999Z","shell.execute_reply.started":"2023-06-09T15:24:27.273393Z","shell.execute_reply":"2023-06-09T15:24:27.285462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.to_csv('df_aligners_benchmark.csv')","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:24:27.289299Z","iopub.execute_input":"2023-06-09T15:24:27.289776Z","iopub.status.idle":"2023-06-09T15:24:27.303045Z","shell.execute_reply.started":"2023-06-09T15:24:27.289732Z","shell.execute_reply":"2023-06-09T15:24:27.302109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some analysis ","metadata":{}},{"cell_type":"code","source":"display(df_scores.describe() )\nprint()\nfor col in df_scores.columns:\n plt.figure(figsize = (15,3))\n plt.hist( df_scores[col], bins = 50 )\n plt.title('Scores '+ col + ' count '+str(len(df_scores)) , fontsize = 20)\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:24:27.304786Z","iopub.execute_input":"2023-06-09T15:24:27.305895Z","iopub.status.idle":"2023-06-09T15:24:33.311519Z","shell.execute_reply.started":"2023-06-09T15:24:27.30586Z","shell.execute_reply":"2023-06-09T15:24:33.310328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load real protein data CAFA5 challenge ","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= 100:\n print(aligner)\n\n # ---------------------------------------------------------------------------------\n # Compute similarity BioPython \n # ---------------------------------------------------------------------------------\n\n t0 = time.time()\n list_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs ]\n tm = np.round(time.time()-t0,2)\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(time.time()-t0, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\nsimilarity_name = 'BioPython Global Blossum62 11_1'\nif similarity_name in list_aligners: # \n print( similarity_name )\n\n # ---------------------------------------------------------------------------------\n # Prepare BioPython aligner \n # --------------------------------------------------------------------------------\n\n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n aligner = Align.PairwiseAligner()\n aligner.mode = 'global'\n aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n aligner.open_gap_score = -11\n aligner.extend_gap_score = -1\n aligner.target_end_gap_score = -11\n aligner.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\n if verbose >= 100:\n print(aligner)\n\n # ---------------------------------------------------------------------------------\n # Compute similarity BioPython \n # ---------------------------------------------------------------------------------\n\n t0 = time.time()\n list_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs ]\n tm = np.round(time.time()-t0,2)\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(time.time()-t0, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\nsimilarity_name = 'BioPython Global Default'\nif similarity_name in list_aligners: # \n print( similarity_name )\n\n # ---------------------------------------------------------------------------------\n # Prepare BioPython aligner \n # --------------------------------------------------------------------------------\n\n from Bio import Align\n from Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\n aligner = Align.PairwiseAligner()\n aligner.mode = 'global'\n# aligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\n# aligner.open_gap_score = -11\n# aligner.extend_gap_score = -1\n# aligner.target_end_gap_score = -11\n# aligner.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\n if verbose >= 100:\n print(aligner)\n\n # ---------------------------------------------------------------------------------\n # Compute similarity BioPython \n # ---------------------------------------------------------------------------------\n\n t0 = time.time()\n list_scores1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs ]\n tm = np.round(time.time()-t0,2)\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(time.time()-t0, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n \n# ---------------------------------------------------------------------------------\n# Parasail Aligner \n# ---------------------------------------------------------------------------------\nsimilarity_name = 'Parasail SW Local Blossum62 11_1'\nif similarity_name in list_aligners: # \n print()\n print( similarity_name )\n t0 = time.time()\n list_scores1 = [parasail.sw_scan_16(sel_seq, seq , 11, 1, parasail.blosum62).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n \nsimilarity_name = 'Parasail NW Global Blossum62 11_1'\nif similarity_name in list_aligners: # \n print()\n print( similarity_name )\n t0 = time.time()\n list_scores1 = [parasail.nw_scan_16(sel_seq, seq , 11, 1, parasail.blosum62).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n\n# ---------------------------------------------------------------------------------\n# skbio Local Blossum62\n# ---------------------------------------------------------------------------------\nsimilarity_name = 'skbio Local Blossum62 11_1'\nif similarity_name in list_aligners: # \n print()\n print( similarity_name )\n query = StripedSmithWaterman(sel_seq, substitution_matrix = dict_sm_blosum62, gap_open_penalty = 11, gap_extend_penalty = 1 )\n t0 = time.time()\n list_scores1 = [query(seq)['optimal_alignment_score'] for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, first lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n\n# ---------------------------------------------------------------------------------\n# Levenshtein\n# --------------------------------------------------------------------------------\nif 'Levenshtein' in list_aligners: # \n print()\n similarity_name = 'Levenshtein'\n print( similarity_name )\n t0 = time.time()\n list_scores1 = [ distance(sel_seq, seq) for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, first lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n\n# ---------------------------------------------------------------------------------\n# Diamond\n# --------------------------------------------------------------------------------\nif 'diamond' in list_aligners:\n similarity_name = 'diamond'\n print(); print('diamond')\n #------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamnod. Create two FASTA files for both data: what is aligned (query(ies)), and on-what is aligned (subject(s))')\n print()\n #------------------------------------------------------------------------------------\n\n fn1_query = \"fasta_for_diamond_test_1_query.txt\"\n ofile = open(fn1_query, \"w\")\n for i in range(1):\n ofile.write(\">\" + 'test_sequence1' + \"\\n\" +sel_seq + \"\\n\")\n ofile.close()\n fn2_subject = \"fasta_for_diamond_test_2_subject.txt\"\n ofile = open(fn2_subject, \"w\")\n for i in range(len(list_seqs)):\n ofile.write(\">\" + list_ids[i] + \"\\n\" +list_seqs[i] + \"\\n\")\n ofile.close()\n\n #------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamond. Create a database file - on-what we align (subjects) ')\n print()\n #------------------------------------------------------------------------------------\n import time\n from subprocess import Popen, PIPE\n db_name='train_db'\n\n t0 = time.time()\n p = Popen(['./diamond', 'makedb', \n '--in',fn2_subject ,\n '-d', db_name], stdin=PIPE, stdout=PIPE)\n stdout, stderr = p.communicate()\n print('Database created in %.1f seconds'%(time.time() - t0 ))\n\n #------------------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamond. Run a blastp-like search')\n print()\n #------------------------------------------------------------------------------------------------\n for str_sensitivity in list_diamond_sensititivity_modes:\n #str_sensitivity = '--ultra-sensitive' # --fast # --mid-sensitive # '--sensitive' # '--more-sensitive' # --very-sensitive # '--ultra-sensitive'\n k = len(list_seqs) # maximal number of \"similar\" proteins for each query protein to be returned \n evalue_threshold = 100_000_000 # default=0.001 # but we want to get as much as possible # controls how similar sequences will be returned \n\n outfile_name = 'matches.tsv'\n t0 = time.time() \n p = Popen(['./diamond', 'blastp', '-d', db_name, '--evalue',str(evalue_threshold),\n '-q', fn1_query,\n '-o', outfile_name, '--max-target-seqs', str(k), '--quiet', str_sensitivity,\n ], stdin=PIPE, stdout=PIPE) # '--quiet' # # \n stdout, stderr = p.communicate()\n tm2 = np.round(time.time()-t0,2)\n print(f'Execution time: {time.time()-t0}s')\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_seqs), len(sel_seq) ))\n\n\n # Cannot make work additional output \n # '-outfmt 6 qseqid sseqid',\n # diamond blastp --query query.fasta --db database.dmnd --outfmt 6 qseqid sseqid qlen slen\n # '--outfmt 6 qseqid sseqid qlen slen pident length mismatch gapopen qstart qend sstart send evalue bitscore' \n\n #------------------------------------------------------------------------------------------------\n if verbose >= 100:\n print()\n print('Diamond. Load results of alignment in pandas dataframe')\n print()\n #------------------------------------------------------------------------------------------------\n\n matches=pd.read_csv(outfile_name, sep='\\t', header=None, \n names=['qseqid', 'sseqid', 'pident', 'length', 'mismatch', \n 'gapopen', 'qstart', 'qend', 'sstart','send', 'evalue', 'bitscore'])\n matches['qseqid']=matches['qseqid'].apply(lambda x: x.split('\\\\t')[0])\n print(str_sensitivity , 'Diamond returned: ', len(matches), 'results. Others will be filled by zero score')\n print(str_sensitivity , 'Diamond returned with e-value < 0.001: ', (matches['evalue']<0.001).sum(), 'results. Others will be filled by zero score')\n if len(matches) > 0:\n print('Min evalue:', matches['evalue'].min(), 'Max evalue:', matches['evalue'].max() , '1percent evalue:', np.percentile( matches['evalue'], 0.01))\n str_tmp = ' ' + str(len(list_seqs)) + ' ' + str_info\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity, 'Time on'+str_tmp] = tm\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity, 'Count Returned'+str_tmp] = len(matches)\n df_stat_diamond.loc[similarity_name + ' ' + str_sensitivity,'Count Returned e-value < 0.001'+str_tmp]=(matches['evalue']<0.001).sum()\n \n\n # \n if 1: # Diamond returns only aligned sequences, we need to get list_scores in the original order with 0 for scores for those missed\n df_tmp = pd.DataFrame(index = list_ids[:len(list_seqs)] , data = range(len(list_seqs)) )\n df_tmp = df_tmp.join( matches[['sseqid', 'bitscore']].set_index('sseqid'), how = 'left')\n df_tmp = df_tmp.sort_values(0) # Just in case sorting was lost during join \n df_tmp = df_tmp.fillna(0)\n list_scores1 = df_tmp['bitscore'].values # Diamond returns only aligned sequences, we need to get list_scores in the original order with 0 for scores for those missed\n\n df_stat.loc[similarity_name + ' ' + str_sensitivity, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_sensitivity + ' ' + str_info] = np.asarray(list_scores1.copy() )\n \n if verbose >= 100:\n print()\n print('---------------------------------------------------- ')\n print(' diamond finished ')\n print('---------------------------------------------------- '); \n print()\n print()\n print()\n \n# ---------------------------------------------------------------------------------\n# pyalign Local\n# --------------------------------------------------------------------------------\nsimilarity_name = 'pyalign Local'\nif similarity_name in list_aligners:\n import pyalign\n # alignment = pyalign.local_alignment(\"INDUSTRY\", \"INTEREST\", gap_cost=0, eq=1, ne=-1)\n print()\n \n print( similarity_name )\n t0 = time.time()\n list_scores1 = [ pyalign.local_alignment(sel_seq,seq, gap_cost=0, eq=1, ne=-1).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n\n# ---------------------------------------------------------------------------------\n# pyalign Global\n# --------------------------------------------------------------------------------\nsimilarity_name = 'pyalign Global'\nif similarity_name in list_aligners:\n import pyalign\n # alignment = pyalign.local_alignment(\"INDUSTRY\", \"INTEREST\", gap_cost=0, eq=1, ne=-1)\n print()\n \n print( similarity_name )\n t0 = time.time()\n list_scores1 = [ pyalign.global_alignment(sel_seq,seq, gap_cost=0, eq=1, ne=-1).score for seq in list_seqs ]\n tm2 = np.round(time.time()-t0,2)\n print( '%.3f times faster than previous'%(tm/tm2))\n tm = tm2\n print(similarity_name, 'similarity. %.4f seconds on %d, lenghth %d'%(tm, len(list_scores1), len(sel_seq) ))\n\n df_stat.loc[similarity_name, 'Time on '+str(len(list_scores1)) + ' ' + str_info ] = tm\n dict_align_results[similarity_name + ' ' + str_info] = np.array(list_scores1.copy() )\n \n\n \n \n \n \n \n \n \n# ---------------------------------------------------------------------------------\n# Compare results \n# --------------------------------------------------------------------------------\n\ndf_scores = pd.DataFrame( dict_align_results )\nprint('Correlation matrix between obtained scores (expected to be all 1, but it is not: )')\ndisplay(df_scores.corr() )\nprint('Spearman correlations')\ndisplay(df_scores.corr(method = 'spearman') )\n\nprint()\ndisplay(df_stat)","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:24:50.732825Z","iopub.execute_input":"2023-06-09T15:24:50.733175Z","iopub.status.idle":"2023-06-09T15:30:33.250656Z","shell.execute_reply.started":"2023-06-09T15:24:50.733145Z","shell.execute_reply":"2023-06-09T15:30:33.248952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('statistics for diamond')\ndf_stat_diamond.to_csv('df_diamond_stat.csv')\ndisplay(df_stat_diamond) \n","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:33.253556Z","iopub.execute_input":"2023-06-09T15:30:33.254001Z","iopub.status.idle":"2023-06-09T15:30:33.279563Z","shell.execute_reply.started":"2023-06-09T15:30:33.25396Z","shell.execute_reply":"2023-06-09T15:30:33.278338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Speed tests on random (%d) and real (%d) data:'%(N_seqs_to_take_for_random_sequences, N_seqs_to_take_CAFA5 ))\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:33.281186Z","iopub.execute_input":"2023-06-09T15:30:33.281572Z","iopub.status.idle":"2023-06-09T15:30:33.298758Z","shell.execute_reply.started":"2023-06-09T15:30:33.28154Z","shell.execute_reply":"2023-06-09T15:30:33.297292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.to_csv('df_aligners_benchmark.csv')","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:33.300942Z","iopub.execute_input":"2023-06-09T15:30:33.301528Z","iopub.status.idle":"2023-06-09T15:30:33.316539Z","shell.execute_reply.started":"2023-06-09T15:30:33.301484Z","shell.execute_reply":"2023-06-09T15:30:33.315173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis: histograms for scores on CAFA5","metadata":{}},{"cell_type":"code","source":"display(df_scores.describe() )\nprint()\nfor col in df_scores.columns:\n plt.figure(figsize = (15,3))\n plt.hist( df_scores[col], bins = 50 )\n plt.title('Scores '+ col + ' count '+str(len(df_scores)) , fontsize = 20)\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:33.318474Z","iopub.execute_input":"2023-06-09T15:30:33.318831Z","iopub.status.idle":"2023-06-09T15:30:39.59365Z","shell.execute_reply.started":"2023-06-09T15:30:33.318801Z","shell.execute_reply":"2023-06-09T15:30:39.59212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Check: human readable descriptions of top scored proteins are as expected - tyrosine kinases \n\nThe selected sequence is SRC human protein - a tirosine kinase. \nhttps://en.wikipedia.org/wiki/Proto-oncogene_tyrosine-protein_kinase_Src\nUniprot: https://www.uniprot.org/uniprotkb/P12931/entry\n\nSo the most similar should in particular be other tyrosine kinases.\n\nThat is indeed the case. See outcomes below. \n\n\nWhat exactly we get depends on the size of the considered part of CAFA5 data - if we would take all - we will see lots of very close relatives to SRC. If we take only 1000 samples we see only some tirosine kinases. How many proteins is taken is controlled by variable: N_seqs_to_take_CAFA5\n\nIf take rather big number of proteins like N_seqs_to_take_CAFA5 = 100_000 (see versions 2-3 of the notebook and output csv file),\nwe will see on top of the most alingments methods - the following expected proteins and in the expected order:\n\n Visual manual inspection of the top similar lists, have the following expectations.\n 1) Most similar - homologs from very similar organisms e.g. for human - it is mouse, rat, etc...\n 2) homologs for more distant orgnisms - like birds, fishes (chicken, Dani Rario = zebrafish)\n 3) Members of the known family: e.g. for SRC: Yes, Fyn, and Fgr - SrcA subfamily, Lck, Hck, Blk, and Lyn - the SrcB subfamily, and Frk in its own subfamily. For P53 - it is P63, P73\n 4) For SRC and protein tyrosine kinases in general - other tyrosine protein kinases \n 5) For SRC and protein kinases in general - other kinases (not necessarily tyrosine, but Serine/threonine-protein kinases)\n ","metadata":{}},{"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)\n\ndf_scores.index = list_ids[:len(df_scores)]\ndf_scores_ext = pd.DataFrame(index = df_scores.index) \ndf_scores_ext['description'] = [dict_ids2description[t] for t in list_ids[:len(df_scores)] ]\ndf_scores_ext = pd.concat( (df_scores_ext, df_scores), axis = 1 )\ndisplay(df_scores_ext.head(5))\ndf_scores_ext.sort_values(df_scores_ext.columns[0], ascending = False ).to_csv('df_scores_with_descriptions.csv')","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:39.595196Z","iopub.execute_input":"2023-06-09T15:30:39.595584Z","iopub.status.idle":"2023-06-09T15:30:39.847913Z","shell.execute_reply.started":"2023-06-09T15:30:39.595554Z","shell.execute_reply":"2023-06-09T15:30:39.84684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N_to_show = 15\nfor col in df_scores.columns:\n if col == 'description': continue\n print('Sort by', col)\n if 'Levenshtein' not in col:\n display(df_scores_ext.sort_values(col, ascending = False).head(N_to_show) )\n else:\n display(df_scores_ext.sort_values(col, ascending = True).head(N_to_show) )\n print()","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:39.85007Z","iopub.execute_input":"2023-06-09T15:30:39.850947Z","iopub.status.idle":"2023-06-09T15:30:40.485856Z","shell.execute_reply.started":"2023-06-09T15:30:39.850897Z","shell.execute_reply":"2023-06-09T15:30:40.484554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Check agreement between score ranking \n\nThe scores themselves might be different by the ranking (ordering) they generate might be quite similar. Here we explore it. \n\n","metadata":{}},{"cell_type":"code","source":"%%time\nfrom itertools import combinations\n\ndf_intersect_stat = pd.DataFrame(); IX_loc = 0\nfor N_loc in [100,200]:\n IX_loc = 0\n for col1,col2 in combinations(df_scores.columns, 2 ):\n if col1 == 'description': continue \n if col2 == 'description': continue \n\n if 'Levenshtein' not in col1:\n s1 = set(df_scores.sort_values(col1, ascending = False).index[:N_loc] )\n else:\n s1 = set(df_scores.sort_values(col1, ascending = True).index[:N_loc] )\n if 'Levenshtein' not in col2:\n s2 = set(df_scores.sort_values(col2, ascending = False).index[:N_loc] )\n else:\n s2 = set(df_scores.sort_values(col2, ascending = True).index[:N_loc] )\n s = s1&s2\n # print('Top%d scored by '%N_loc, col1,' vs ', ' by ', col2 , '. Intersection size is: %d'%(len(s)))\n # print()\n df_intersect_stat.loc[IX_loc,'Top scored by' ] = col1\n df_intersect_stat.loc[IX_loc,'vs top scored by' ] = col2\n df_intersect_stat.loc[IX_loc,'Intersection size of top'+str(N_loc) ] = len(s)\n IX_loc +=1 \n # print()\n # print()\ndf_intersect_stat.to_csv('df_intersect_stat.csv')\ndisplay(df_intersect_stat)","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:40.492266Z","iopub.execute_input":"2023-06-09T15:30:40.492696Z","iopub.status.idle":"2023-06-09T15:30:41.288137Z","shell.execute_reply.started":"2023-06-09T15:30:40.492661Z","shell.execute_reply":"2023-06-09T15:30:41.286847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Advanced agreement analysis\n\nMethodology discussed in : https://www.kaggle.com/code/alexandervc/curve-defined-by-permutations\nand many other notebooks of the project: https://www.kaggle.com/datasets/alexandervc/research-project-01-around-multimodal-singlecell\n\nThe idea is: when the curve below becomes parabola - that is starting point of complete disconcordance between the rankings\n\nWhile slope of the linear approximation on the initial period - characterizes percent of concordance between two rankings.\n\nRoughly speaking taking topN from one ranking and topN from another ranking - how many we will have common elements - f(N) ?\n\nThe point is when rankigs well agree we that f(N) - is linear, while full disagreement (random) f(N) is parabola.\n\nPS\n\nWhen we have more than two rakings - cubic, quartic, etc polinom appears.","metadata":{}},{"cell_type":"code","source":"%%time\nl_x = np.arange(1,1001)* (len(df_scores)/1000) # To speed-up calculate not for every point , but 1000 points \nl_x = l_x.astype(int)\n\nprint('For brevity - compare only the first alinger with other ones.')\n# print('Full option - compare all pairs - that would not easy to grasp')\ncol1 = df_scores.columns[0]\nfor col2 in df_scores.columns[1:]:\n# for col1,col2 in combinations(list_similarities_to_compare, 2 ):\n print(); print(); print();\n print(col1, ' vs ', col2 )\n l_y = []\n for i0,N_loc in enumerate(l_x):\n \n if 'Levenshtein' not in col1:\n s1 = set(df_scores.sort_values(col1, ascending = False).index[:N_loc] )\n else:\n s1 = set(df_scores.sort_values(col1, ascending = True).index[:N_loc] )\n if 'Levenshtein' not in col2:\n s2 = set(df_scores.sort_values(col2, ascending = False).index[:N_loc] )\n else:\n s2 = set(df_scores.sort_values(col2, ascending = True).index[:N_loc] )\n s = s1&s2\n l_y.append(len(s))\n l_y = np.array(l_y)\n\n \n deg = 2 # len( list_similarities_to_compare )\n pol = np.polyfit(l_x,l_y, deg)\n\n N1 = 10\n pol2 = np.polyfit(l_x[:N1],l_y[:N1], 1) \n print('Percent of agreement:', np.round(100*pol2[0],2) )\n# print(); print(); \n\n print('Linear approximating polynom %+.5f x %+.1f :'%(pol2[0], pol2[1] ) )\n print('Parabolic approximating polynom %-.10f x^2 %+.1f x %+.1f :'%(pol[0], pol[1], pol[2] ) )\n\n plt.figure(figsize = (20,6))\n plt.suptitle(col1 + ' vs ' +col2, fontsize = 18)\n plt.subplot(121)\n plt.plot(l_x,l_y , label = 'Data')\n plt.plot(l_x, np.polyval(pol,l_x ) , label = 'Approximation Parabola' )\n plt.plot(l_x, np.polyval(pol2,l_x ) , label = 'Approximation Linear' )\n plt.legend()\n plt.grid()\n\n N2 = 20\n plt.subplot(122)\n plt.title( 'Cutoff: '+str(N2) )\n plt.plot(l_x[:N2],l_y[:N2] , label = 'Data' )\n plt.plot(l_x[:N2], np.polyval(pol,l_x[:N2]) , label = 'Approximation Parabola' )\n plt.plot(l_x[:N2], np.polyval(pol2,l_x[:N2] ) , label = 'Approximation Linear' )\n plt.legend()\n plt.grid()\n plt.show()\n\n\n plt.figure(figsize = (20,6))\n plt.suptitle('Plots showing descrepancy between data and approximations:\\n y/x - should be constant if y is linear\\n or linear if y is x^2' , fontsize = 16)\n plt.subplot(121)\n plt.plot(l_x,l_y/l_x , label = 'Data Normalized')\n plt.legend()\n plt.grid()\n\n N2 = 20\n plt.subplot(122)\n plt.title( 'Cutoff: '+str(N2) )\n plt.plot(l_x[:N2],l_y[:N2]/l_x[:N2] , label = 'Data Normalized' )\n plt.legend()\n plt.grid()\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:30:41.289968Z","iopub.execute_input":"2023-06-09T15:30:41.290369Z","iopub.status.idle":"2023-06-09T15:32:15.481516Z","shell.execute_reply.started":"2023-06-09T15:30:41.290334Z","shell.execute_reply":"2023-06-09T15:32:15.480307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Show tests results ","metadata":{}},{"cell_type":"code","source":"print('statistics for diamond')\n# df_stat_diamond.to_csv('df_diamond_stat.csv')\ndisplay(df_stat_diamond) \n\nprint('Speed tests on random (%d) and real (%d) data:'%(N_seqs_to_take_for_random_sequences, N_seqs_to_take_CAFA5 ))\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2023-06-09T15:32:15.483462Z","iopub.execute_input":"2023-06-09T15:32:15.484614Z","iopub.status.idle":"2023-06-09T15:32:15.518858Z","shell.execute_reply.started":"2023-06-09T15:32:15.484566Z","shell.execute_reply":"2023-06-09T15:32:15.5175Z"},"trusted":true},"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-09T15:32:15.520663Z","iopub.execute_input":"2023-06-09T15:32:15.521537Z","iopub.status.idle":"2023-06-09T15:32:15.531564Z","shell.execute_reply.started":"2023-06-09T15:32:15.521488Z","shell.execute_reply":"2023-06-09T15:32:15.530319Z"},"trusted":true},"execution_count":null,"outputs":[]}],"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"}}