{"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\nWe benchmark \"similarity\"-measures for proteins  obtained by different methods: in particular by classical sequence alignment related, vs modern\nembedding (deep-learning related). The outcome is that embeddings produce similar quality results to the most slow-hard alignment methods, but 1000+ times faster. Some faster sequence alignment like methods give much lower quality, while being still 100+ times slower than embeddings. \n\nThe benchmarks is performed on a set of some well studied proteins like SRC oncogene, P53 \"genome guardian\" onco-supressor (top number of publications ever), etc. And organized as follows:\n\n    Compute \"similarity\" measures for selected protein and 140 000+ proteins from the CAFA5 challenge.\n    Order the proteins by the obtained similarity scores.\n    \n    Test 1. Compare to Gene-Ontology-similarity.\n    \n    Compute the difference between gene-ontology terms for the selected protein and i-th similar protein. (Dice coefficient for measuring the \"distance\" between the two sets of gene-ontology terms. Also smooth by moving average to get easy to analyse results).\n    One expects better similarity scores is - the better it is agreed with the gene ontology score, i.e. \n    Outcomes: \n        The embedding similarity is almost the same well agreed with Gene-Ontology-similarity as the top slow-hard (Smith Waterman) local alignment (being 1000 times faster).\n        The similarities obtained by faster alignments (stripped banded Smith-Waterman from Skbio), as well as simple Levenshtein-distinance-similarity much worse agreed with Gene-Ontology-similarity \n        \n        The embedding similarity even for 1000-th (and even more 4000-th) ranked protein still indicates certain non-trivial relation between proteins and in particular certain Gene-Ontology-similarity. Thus the similarity is quite \"long-range\" - similar to top slowe-hard local alignment (Smith Waterman).\n        \n    Test 2. Compare to \"keyword\"-similarity.\n    \n    Proteins comes with short textual human readable desription. E.g. for SRC - \"Proto-oncogene tyrosine-protein kinase Src\"\n    For each protein create and indicator - 0,1:  does its description contains a selected (manually) keyword(s) - like \"kinase\", \"protein kinase\", \"tumor\"... Optionally compute a moving average of that indicator (moving along index of ordering of the proteins by similarity).\n    One expectes that the obtained curves start high values (near 1) and goes down. \n    We see that embedding the curve goes down the same slow for slow-hard (Smith Waterman) local alignment.\n    That means e.g. that similar proteins \n    \n    Test 3. Manual checks.\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-necesarily tyrosine, but  Serine/threonine-protein kinases)\n    \n\nSome other outcomes:\n\n    Hundreds vs Thousands seen similar proteins. \n    Simpler-faster sequence methods like Levenshtein distance or skibio alignment can see HUNDREDS similar proteins, while emebeddings and slow-hard alignment can see THOUSANDs similar proteins. That means we can assign meaningful threshold on  similarity score which can distinguish hundrends proteins in one case and thousands in the other. That is real uplift of performance for emebeddings comparing to these alignment methods. \n    \n\n\n\n#### More on SRC :\n\nhttps://en.wikipedia.org/wiki/Proto-oncogene_tyrosine-protein_kinase_Src \n\nUniprot: https://www.uniprot.org/uniprotkb/P12931/entry\n\nSrc is in the Src family of kinase - https://en.wikipedia.org/wiki/Src_family_kinase : \n\n    Src kinase family is a family of non-receptor tyrosine kinases that includes nine members: Src, Yes, Fyn, and Fgr, forming the SrcA subfamily, Lck, Hck, Blk, and Lyn in the SrcB subfamily, and Frk in its own subfamily. Frk has homologs in invertebrates such as flies and worms, and Src homologs exist in organisms as diverse as unicellular choanoflagellates, but the SrcA and SrcB subfamilies are specific to vertebrates. Src family kinases contain six conserved domains: a N-terminal myristoylated segment, a SH2 domain, a SH3 domain, a linker region, a tyrosine kinase domain, and C-terminal tail.[1]\n\n    Src family kinases interact with many cellular cytosolic, nuclear and membrane proteins, modifying these proteins by phosphorylation of tyrosine residues. A number of substrates have been discovered for these enzymes.[2][3][4] Deregulation, including constitutive activation or over expression, may contribute to the progression of cellular transformation and oncogenic activity.[5]\n\n### SRC in cancer, in particular in Medulloblastoma\n\nSRC is famous oncogene, its disregulation appears in many cancers.\nThere is  certain research on drugs  - SRC inhibitors, which might be helpful in some cancer cases.\nAs usually with cancers drugs do not work in all cases, and better understanding when they can be sucessfully applied is a typical research task.\n    \nAt the moment I am trying to   understand better the role of SRC and other oncogenes in Medulloblastoma,\nsee some datasets and notebooks on Kaggle:\n\nhttps://www.kaggle.com/datasets/alexandervc/medulloblastoma-omics-data\n\nhttps://www.kaggle.com/datasets/alexandervc/medulloblastoma-rna-seq-phosphoproteomics-ptms\n\nhttps://www.kaggle.com/datasets/alexandervc/scrnaseq-medulloblastoma-nature2019\n\nhttps://www.kaggle.com/datasets/alexandervc/scrna-seq-mouse-model-shh-driven-medulloblastoma\n\nhttps://www.kaggle.com/search?q=medulloblastoma+in%3Adatasets\n\nWould be happy to discuss with someone interested in similar topics.\nHopefully some CAFA like models can be used to understand better these questions.\n\n\n### PS\n\nFor the competition purposes one can use \"similarity\" measures in a several ways: find most similar proteins to the given one and try to transfer labels from them to it. Or one can use \"similarities\" obained by alignment as features, for examples selecting several sequences and \"similarities\" to them - are features (similar to Levenshtein distance feature here https://www.kaggle.com/code/alexandervc/cafa5-levenshtein-distance-features). Or one can use \"similarities\" to obtain groups for groupwise validation. Though the main cavet is that alignments are quite slow.\n\nOne can use alignments from BioPython: 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","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time()\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-06-06T14:33:21.231058Z","iopub.execute_input":"2023-06-06T14:33:21.233381Z","iopub.status.idle":"2023-06-06T14:33:21.257831Z","shell.execute_reply.started":"2023-06-06T14:33:21.233203Z","shell.execute_reply":"2023-06-06T14:33:21.256157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Imports","metadata":{}},{"cell_type":"code","source":"from Bio.pairwise2 import format_alignment\n# from Bio.SubsMat import MatrixInfo \nfrom Bio import pairwise2\nfrom Bio import SeqIO, SearchIO\nfrom Bio.Seq import Seq\nfrom Bio.SeqRecord import SeqRecord\nfrom Bio.Blast import NCBIWWW\nfrom Bio.Blast import NCBIXML\n\nfrom Bio.Phylo.TreeConstruction import DistanceTreeConstructor\nfrom Bio.Phylo.TreeConstruction import DistanceCalculator\nfrom Bio.Phylo.PhyloXML import Phylogeny\nfrom Bio import Phylo\nfrom pprint import pprint","metadata":{"execution":{"iopub.status.busy":"2023-06-06T14:33:21.979726Z","iopub.execute_input":"2023-06-06T14:33:21.980826Z","iopub.status.idle":"2023-06-06T14:33:21.989738Z","shell.execute_reply.started":"2023-06-06T14:33:21.980783Z","shell.execute_reply":"2023-06-06T14:33:21.988435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load/Prepare Data","metadata":{}},{"cell_type":"code","source":"%%time\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]","metadata":{"execution":{"iopub.status.busy":"2023-06-06T14:33:23.396765Z","iopub.execute_input":"2023-06-06T14:33:23.397232Z","iopub.status.idle":"2023-06-06T14:33:33.937529Z","shell.execute_reply.started":"2023-06-06T14:33:23.397198Z","shell.execute_reply":"2023-06-06T14:33:33.936268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\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(t, (l1<t).sum(),(l2<t).sum(), )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T14:33:33.941542Z","iopub.execute_input":"2023-06-06T14:33:33.941922Z","iopub.status.idle":"2023-06-06T14:33:34.972368Z","shell.execute_reply.started":"2023-06-06T14:33:33.941889Z","shell.execute_reply":"2023-06-06T14:33:34.971051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set the selected protein sequence (e.g. SRC)","metadata":{}},{"cell_type":"code","source":"sel_prot_id = 'Q99592' # ZBTB18 'Q9H4A3' # WNK1 'Q9BYP7'# WNK3 'Q96J92' # WNK4 'Q9Y3S1'# WNK2 #  'P06493'# CDK1 '# P04637' # P53  # 'P12931' # SRC    \n\n# ['P04637', 'O15350', 'P06493', 'Q9Y3S1', 'Q96J92', 'Q9BYP7' , 'Q9H4A3', 'Q99592']: # TP53, TP73 , CDK1, WNK2 WNK4 WNK3 WNK1 ZBTB18\n\nprint('Selected protein Id:', sel_prot_id )\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:14:31.523451Z","iopub.execute_input":"2023-06-06T15:14:31.524037Z","iopub.status.idle":"2023-06-06T15:14:31.534494Z","shell.execute_reply.started":"2023-06-06T15:14:31.523999Z","shell.execute_reply":"2023-06-06T15:14:31.532011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Get sequence for the protin and configure some params\n","metadata":{}},{"cell_type":"code","source":"%%time\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\")\n\n\nlist_keywords_describing_protein = []\nwindow_len_for_keywords_indicator_moving_average = 50# 50 # SRC # 1 # P53\nlist_cutoffs_to_show_distribution_Levenshtein = [] # [402, 399 ] # SRC  # [296, 295 ] # P53\nlist_cutoffs_to_show_distribution_LocalAlignmentBioPython = []#  [70, 100, 300, 568, 574, 575] # SRC \nlist_cutoffs_to_show_distribution_LocalAlignmentSkbio = [] # \nif  sel_prot_id == 'P12931':\n    prot_name = 'SRC human P12931'\n    sel_seq = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\n    print('Selected protein:', 'SRC - P12931- Proto-oncogene tyrosine-protein kinase Src - https://www.uniprot.org/uniprotkb/P12931/entry' )\n    print('https://en.wikipedia.org/wiki/Proto-oncogene_tyrosine-protein_kinase_Src')\n    for seq in sequences:\n        if seq.id == sel_prot_id:\n            print('Description:', seq.description) \n            break\n            \n    list_keywords_describing_protein = ['kinase', 'protein kinase']\n    window_len_for_keywords_indicator_moving_average = 50 # SRC # 1 # P53\n\n    list_cutoffs_to_show_distribution_Levenshtein = [402, 399 ] # SRC  # the threshold is around 74-75 of length \n    list_cutoffs_to_show_distribution_LocalAlignmentBioPython = [70, 100, 300, 568, 574, 575] # SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = [2200,2300, 2250 ] # SRC \n    \nelif sel_prot_id == 'P04637':\n    prot_name = 'P53 human ' + sel_prot_id\n    print('P53 - \"Guardian of the genome\" - the main onco-supressor - most studied gene ever')\n    print('https://www.uniprot.org/uniprotkb/P04637/entry')\n    print('https://en.wikipedia.org/wiki/P53')\n    print(\"The TP53 gene is the most frequently mutated gene (>50%) in human cancer, indicating that the TP53 gene plays a crucial role in preventing cancer formation.[5] TP53 gene encodes proteins that bind to DNA and regulate gene expression to prevent mutations of the genome.[12] In addition to the full-length protein, the human TP53 gene encodes at least 15 protein isoforms.\")\n    for seq in sequences:\n        if seq.id == sel_prot_id:\n            sel_seq = str( seq.seq )\n            print('Description:', seq.description) \n            break\n    list_keywords_describing_protein = ['tumor']\n    window_len_for_keywords_indicator_moving_average = 1 # P53 50 # SRC # \n    list_cutoffs_to_show_distribution_Levenshtein = [296, 295 ] # P53  # the threshold is around 74-75 of length \n    list_cutoffs_to_show_distribution_LocalAlignmentBioPython = np.array([70, 100, 300, 568, 574, 575])*0.75 #   rescale params for SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = [1750, 1800, 1810,  1820] # P53\n\nelif sel_prot_id == 'P06493':\n    prot_name = 'CDK1 human ' + sel_prot_id\n    print('CDK1 - key role in cell proliferation cycle regulation for eukaryots from yeasts to humans')\n    print('https://www.uniprot.org/uniprotkb/P06493/entry')\n    print('https://en.wikipedia.org/wiki/Cyclin-dependent_kinase_1')\n    print(\"Cyclin-dependent kinase 1 also known as CDK1 or cell division cycle protein 2 homolog is a highly conserved protein that functions as a serine/threonine protein kinase, and is a key player in cell cycle regulation.[5] It has been highly studied in the budding yeast S. cerevisiae, and the fission yeast S. pombe, where it is encoded by genes cdc28 and cdc2, respectively.[6] With its cyclin partners, Cdk1 forms complexes that phosphorylate a variety of target substrates (over 75 have been identified in budding yeast); phosphorylation of these proteins leads to cell cycle progression.[7]\")\n    for seq in sequences:\n        if seq.id == sel_prot_id:\n            sel_seq = str( seq.seq )\n            print('Description:', seq.description) \n            break\n    list_keywords_describing_protein = ['kinase', 'cyclin-dependent kinase']\n    window_len_for_keywords_indicator_moving_average = 5 # P53 50 # SRC # \n    list_cutoffs_to_show_distribution_Levenshtein = np.array([0.72, 0.73, 0.74, 0.75, 0.76 ]) * len(sel_seq)  # the threshold is around 74-75 of length \n    list_cutoffs_to_show_distribution_Levenshtein = list_cutoffs_to_show_distribution_Levenshtein.astype(int)\n    list_cutoffs_to_show_distribution_LocalAlignmentBioPython = np.array([70, 100, 300, 568, 574, 575])*len(sel_seq) / 536 #   rescale params for SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = np.array([  2500, 2600, 2700])*len(sel_seq)/536 # rescale params for SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = list_cutoffs_to_show_distribution_LocalAlignmentSkbio.astype(int)\n    \nelse:\n    prot_name = sel_prot_id \n    for seq in sequences:\n        if seq.id == sel_prot_id:\n            sel_seq = str( seq.seq )\n            print('Description:', seq.description) \n            break\n    list_keywords_describing_protein = []#  ['kinase', 'cyclin-dependent kinase']\n    for kw in ['kinase', 'cyclin-dependent kinase',  'protein kinase','cyclin-dependent kinase', 'Serine/threonine-protein kinase' ]:\n        if kw.lower() in str(seq.description).lower(): \n            list_keywords_describing_protein += [ kw ]#  ['kinase', 'cyclin-dependent kinase']\n    s = str( seq.description )\n    if '|' in s:\n        l = s.split('|')\n        if len(l)>=2:\n            prot_name = l[2].split(' ')[0] +' ' +  sel_prot_id\n    print(prot_name )\n        \n    window_len_for_keywords_indicator_moving_average = 5 # P53 50 # SRC # \n    list_cutoffs_to_show_distribution_Levenshtein = np.array([0.72, 0.73, 0.74, 0.75, 0.76 ]) * len(sel_seq)  # the threshold is around 74-75 of length \n    list_cutoffs_to_show_distribution_Levenshtein = list_cutoffs_to_show_distribution_Levenshtein.astype(int)\n    list_cutoffs_to_show_distribution_LocalAlignmentBioPython = np.array([70, 100, 300, 568, 574, 575])*len(sel_seq) / 536 #   rescale params for SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = np.array([1800, 2000,  2500, 2600, 2700])*len(sel_seq)/536 # rescale params for SRC \n    list_cutoffs_to_show_distribution_LocalAlignmentSkbio = list_cutoffs_to_show_distribution_LocalAlignmentSkbio.astype(int)\n    \nprint('Sequence length:', len(sel_seq) )     \nprint(sel_seq )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:58:35.60589Z","iopub.execute_input":"2023-06-06T15:58:35.606316Z","iopub.status.idle":"2023-06-06T15:58:36.092955Z","shell.execute_reply.started":"2023-06-06T15:58:35.606284Z","shell.execute_reply":"2023-06-06T15:58:36.091972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s = 'Q9Y3S1 sp|Q9Y3S1|WNK2_HUMAN Serine/threonine-protein kinase WNK2 OS=Homo sapiens OX=9606 GN=WNK2 PE=1 SV=4'\nif '|' in s:\n    l = s.split('|')\n    if len(l)>=2:\n        print(l[2].split(' ')[0] )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:56:25.09892Z","iopub.execute_input":"2023-06-06T15:56:25.099441Z","iopub.status.idle":"2023-06-06T15:56:25.10762Z","shell.execute_reply.started":"2023-06-06T15:56:25.099403Z","shell.execute_reply":"2023-06-06T15:56:25.106Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_scoring_data = {} # Results will be saved here ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:18:09.609567Z","iopub.execute_input":"2023-06-06T15:18:09.610025Z","iopub.status.idle":"2023-06-06T15:18:09.615734Z","shell.execute_reply.started":"2023-06-06T15:18:09.609985Z","shell.execute_reply":"2023-06-06T15:18:09.614401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pd.set_option('display.max_columns', 600)\n# pd.set_option('display.width', 1500)\npd.set_option('max_colwidth', 500)\npd.set_option('display.max_rows', 500)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T14:43:30.479136Z","iopub.execute_input":"2023-06-06T14:43:30.479536Z","iopub.status.idle":"2023-06-06T14:43:30.486174Z","shell.execute_reply.started":"2023-06-06T14:43:30.479507Z","shell.execute_reply":"2023-06-06T14:43:30.484858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Levenshtein","metadata":{}},{"cell_type":"code","source":"!pip install python-Levenshtein\nfrom Levenshtein import distance\nedit_dist = distance(\"ah\", \"aho\")\nedit_dist","metadata":{"execution":{"iopub.status.busy":"2023-06-06T14:42:14.039453Z","iopub.execute_input":"2023-06-06T14:42:14.040603Z","iopub.status.idle":"2023-06-06T14:42:25.552237Z","shell.execute_reply.started":"2023-06-06T14:42:14.040552Z","shell.execute_reply":"2023-06-06T14:42:25.550442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Compute distances. Look for threshold which indicates proteins are outliers of the distribution - that suggests have biologically meaningful similarity to the selected one ","metadata":{}},{"cell_type":"code","source":"%%time\nsimilarity_name = 'Levenshtein'\nprint('Length selected sequence:', len(sel_seq) , 'protein:', prot_name )\n\nl1 = np.array([ distance(sel_seq, seq) for seq in list_seqs ] )\nl2 = np.array([ distance(sel_seq, seq) for seq in list_seqs2 ] )\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.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 20)\nplt.subplot(121)\nplt.hist(l1[l1<3*len(sel_seq)], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<3*len(sel_seq)], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\n# list_cutoffs_to_show_distribution_Levenshtein = [402, 399 ] # SRC \n# list_cutoffs_to_show_distribution_Levenshtein = [296, 295 ] # P53\nprint('Attempt to show \"outlier\" part of the score distribution - which corresponds to biologically significant similarity. And start of the main part od the distribution - which is probably mostly random')\nfor N in list_cutoffs_to_show_distribution_Levenshtein: \n    print('Cutoff:', N , 'count:', (l1<N).sum()  )\n    plt.figure(figsize=(20,4))\n    plt.suptitle(similarity_name +' similarity to ' + prot_name + ' cutoff: '+str(N), fontsize = 20)\n    plt.subplot(121)\n    plt.hist(l1[l1<N], bins=150)\n    plt.title('Train', fontsize = 15)\n    plt.subplot(122)\n    plt.hist(l2[l2<N], bins=150)\n    plt.title('Test', fontsize = 15)\n    plt.show()\n\nfor t in [10,50,100,200,300,350, 397,398,399, 400,401,402,403,404,405,406,407,408,409,410,420,430,440, 450]:\n    print(t, 'cumulative:', (l1<t).sum(),(l2<t).sum(), 'equals to', (l1==t).sum(),(l2==t).sum(), )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:19:13.746266Z","iopub.execute_input":"2023-06-06T15:19:13.74704Z","iopub.status.idle":"2023-06-06T15:19:44.765979Z","shell.execute_reply.started":"2023-06-06T15:19:13.746989Z","shell.execute_reply":"2023-06-06T15:19:44.764573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check protein annotations - really correspond to expected biological similarity to the selected protein","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","metadata":{"execution":{"iopub.status.busy":"2023-06-06T14:49:14.602767Z","iopub.execute_input":"2023-06-06T14:49:14.603282Z","iopub.status.idle":"2023-06-06T14:49:14.610253Z","shell.execute_reply.started":"2023-06-06T14:49:14.603248Z","shell.execute_reply":"2023-06-06T14:49:14.608677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"di1 = pd.DataFrame(index = list_ids , data = np.array( [l1, list_lens]).T , columns = ['distance', 'len']).sort_values('distance', ascending = True)\ndi2 = pd.DataFrame(index = list_ids2 , data = np.array( [l2, list_lens2]).T , columns = ['distance', 'len']).sort_values('distance', ascending = True)\n\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use \n\nprint('Top250')\ndisplay(di.iloc[:250,:])\nprint('400-420')\ndisplay(di.iloc[400:420,:])\nprint('4000-4020')\ndisplay(di.iloc[4000:4020,:])\nprint('40000-40020')\ndisplay(di.iloc[40000:40020,:])\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:20:34.783069Z","iopub.execute_input":"2023-06-06T15:20:34.783565Z","iopub.status.idle":"2023-06-06T15:20:35.766518Z","shell.execute_reply.started":"2023-06-06T15:20:34.783527Z","shell.execute_reply":"2023-06-06T15:20:35.765301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nWL = window_len_for_keywords_indicator_moving_average # 50 # SRC # 1 P53\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # list_keywords_describing_protein = ['kinase', 'protein kinase']\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:24:33.122317Z","iopub.execute_input":"2023-06-06T15:24:33.122808Z","iopub.status.idle":"2023-06-06T15:24:34.84929Z","shell.execute_reply.started":"2023-06-06T15:24:33.122762Z","shell.execute_reply":"2023-06-06T15:24:34.847922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Local Alignment with BioPython ","metadata":{}},{"cell_type":"code","source":"%%time\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\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","metadata":{"execution":{"iopub.status.busy":"2023-06-05T17:27:41.238807Z","iopub.execute_input":"2023-06-05T17:27:41.239573Z","iopub.status.idle":"2023-06-05T17:27:41.374819Z","shell.execute_reply.started":"2023-06-05T17:27:41.23954Z","shell.execute_reply":"2023-06-05T17:27:41.373077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#dir(list_align1[0][0])\n# for k in ['aligned', 'coordinates',  'counts',  'format',  'indices',  'infer_coordinates',\n#  'inverse_indices',  'map', 'path',  'query',  'score',  'sequences',  'shape',  'sort',  'substitutions', 'target']:\n#     print(k, getattr( list_align1[0][0], k) )\n# print( list_align1[0][0] )\n","metadata":{"execution":{"iopub.status.busy":"2023-06-05T17:27:41.376222Z","iopub.execute_input":"2023-06-05T17:27:41.37657Z","iopub.status.idle":"2023-06-05T17:27:41.381093Z","shell.execute_reply.started":"2023-06-05T17:27:41.376542Z","shell.execute_reply":"2023-06-05T17:27:41.380073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom Bio import Align\nfrom Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \n\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\naligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\naligner.open_gap_score = -11\naligner.extend_gap_score = -1\naligner.target_end_gap_score = -11\naligner.query_end_gap_score = -1\n# Parameters: Blossum62 and gap-score: 11,1 are used in NCBI web server blastP (\"P\" - for protein) : https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome\n# hopefully the role of params in BioPython and NCBI is the same,\n# at least we get reasonable results \n\nprint(aligner)\n\nsimilarity_name = 'BioPython Local Blossum62'\n\n# N = 10_000 # 4.65 secs  for 1000 , so 40 secs for 10_000, 400 secs for 100_000, 1200 secs( 20 minutes) for 300 k \nlist_align1 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs_corrected ]\nlist_align2 = [aligner.align(sel_seq, seq )[0].score for seq in list_seqs_corrected2 ]\n","metadata":{"execution":{"iopub.status.busy":"2023-06-05T17:27:41.382516Z","iopub.execute_input":"2023-06-05T17:27:41.38291Z","iopub.status.idle":"2023-06-05T18:10:53.686843Z","shell.execute_reply.started":"2023-06-05T17:27:41.382875Z","shell.execute_reply":"2023-06-05T18:10:53.685876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Look for threshold which indicates proteins are outliers of the distribution - that suggests they have biologically meaningful similarity to the selected one ","metadata":{}},{"cell_type":"code","source":"%%time \n# l1 = np.array([ t[0].score for t in list_align1 ] )\n# l2 = np.array([ t[0].score for t in list_align2 ] )\nl1 = np.array(list_align1)\nl2 = np.array(list_align2 )\n\nd = pd.concat( (pd.Series(l1).describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]), pd.Series(l2).describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]) ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 17)\nplt.subplot(121)\nplt.hist(l1[l1<300], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<300], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\n# list_cutoffs_to_show_distribution_LocalAlignmentBioPython = [70, 100, 300, 568, 574, 575] # SRC \nprint('Attempt to show \"outlier\" part of the score distribution - which corresponds to biologically significant similarity. And start of the main part od the distribution - which is probably mostly random')\nfor N in list_cutoffs_to_show_distribution_LocalAlignmentBioPython: \n    plt.figure(figsize=(20,4))\n    plt.suptitle(similarity_name +' similarity to ' + prot_name + ' cutoff: '+str(N), fontsize = 17)\n    plt.subplot(121)\n    plt.hist(l1[l1>N], bins=150)\n    plt.title('Train', fontsize = 15)\n    plt.subplot(122)\n    plt.hist(l2[l2>N], bins=150)\n    plt.title('Test', fontsize = 15)\n    plt.show()\n    print('count train',(l1>N).sum(),'count test',(l2>N).sum() )\n\nprint(); print();\nfor t in range(560,580): # [10,50,100, 450]:\n    print(t, 'cumulative:', (l1>t).sum(),(l2>t).sum(), 'equals to', (l1==t).sum(),(l2==t).sum(), )","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:10:53.688319Z","iopub.execute_input":"2023-06-05T18:10:53.688732Z","iopub.status.idle":"2023-06-05T18:11:00.128005Z","shell.execute_reply.started":"2023-06-05T18:10:53.688705Z","shell.execute_reply":"2023-06-05T18:11:00.126819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check protein annotations - really correspond to expected similarity to the selected protein","metadata":{}},{"cell_type":"code","source":"di1 = pd.DataFrame(index = list_ids , data = np.array( [l1, list_lens]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi2 = pd.DataFrame(index = list_ids2 , data = np.array( [l2, list_lens2]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\n\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndisplay(di.iloc[:20,:])\nprint('400-420')\ndisplay(di.iloc[400:420,:])\nprint('4000-4020')\ndisplay(di.iloc[4000:4020,:])\nprint('40000-40020')\ndisplay(di.iloc[40000:40020,:])\n\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\n\n\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use ","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:11:00.129492Z","iopub.execute_input":"2023-06-05T18:11:00.129806Z","iopub.status.idle":"2023-06-05T18:11:01.264145Z","shell.execute_reply.started":"2023-06-05T18:11:00.129779Z","shell.execute_reply":"2023-06-05T18:11:01.262971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nWL = window_len_for_keywords_indicator_moving_average # 50 # SRC # 1 # P53\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # ['kinase', 'protein kinase']:\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compare similarity rankings\n\nMethodology discussed in : https://www.kaggle.com/code/alexandervc/curve-defined-by-permutations\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 perid - 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. \n","metadata":{}},{"cell_type":"code","source":"%%time\nl = np.arange(100,50000,100)\nl2 = []\nfor i0,N in enumerate(l):\n    for i1,k in enumerate(dict_scoring_data.keys()):\n        dt = dict_scoring_data[k][0]\n        if i1 == 0:\n            s = set(dt.index[:N])\n        else:\n            s = s & set(dt.index[:N])\n        if i0 == 0:\n            print(k)\n            display(dt.head(2))\n            \n    l2.append(len(s))\n    \nprint(len(l2), l2[:10])","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:11:01.265621Z","iopub.execute_input":"2023-06-05T18:11:01.266284Z","iopub.status.idle":"2023-06-05T18:11:07.20192Z","shell.execute_reply.started":"2023-06-05T18:11:01.266247Z","shell.execute_reply":"2023-06-05T18:11:07.20093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndeg = len( dict_scoring_data )\npol = np.polyfit(l,l2, deg)\n\nN = 50\npol2 = np.polyfit(l[:N],l2[:N], 1) \nprint('Percent of agreement:', np.round(100*pol2[0],2) )\nprint(); print(); \n\nl = np.asarray(l)\nl2 = np.asarray(l2)\n\nstr_inf = 'Compare: ' + str(list(dict_scoring_data.keys()))\n\n\nprint('Linear approximating polynom  %+.5f x %+.1f :'%(pol2[0], pol2[1] ) )\nprint('Parabolic approximating polynom %-.4f x^2 %+.1f x %+.1f :'%(pol[0], pol[1], pol[2] ) )\n\nplt.figure(figsize = (20,6))\nplt.suptitle(str_inf, fontsize = 20 )\n\nplt.subplot(121)\nplt.plot(l,l2 , label = 'Data')\nplt.plot(l, np.polyval(pol,l ) , label = 'Approximation Parabola' )\nplt.plot(l, np.polyval(pol2,l ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N] , label = 'Data' )\nplt.plot(l[:N], np.polyval(pol,l[:N]) , label = 'Approximation Parabola' )\nplt.plot(l[:N], np.polyval(pol2,l[:N] ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\nplt.show()\n\n\nplt.figure(figsize = (20,6))\nplt.suptitle(str_inf, fontsize = 20 )\nplt.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')\nplt.subplot(121)\nplt.plot(l,l2/l , label = 'Data Normalized')\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N]/l[:N] , label = 'Data Normalized' )\nplt.legend()\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:11:07.204105Z","iopub.execute_input":"2023-06-05T18:11:07.204533Z","iopub.status.idle":"2023-06-05T18:11:08.237277Z","shell.execute_reply.started":"2023-06-05T18:11:07.204495Z","shell.execute_reply":"2023-06-05T18:11:08.23613Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Skbio local alignment\n\nSee also: \nhttps://www.kaggle.com/code/alexandervc/cafa5-19-alignments-skbio","metadata":{}},{"cell_type":"code","source":"!pip install scikit-bio","metadata":{"execution":{"iopub.status.busy":"2023-06-06T11:46:33.826514Z","iopub.execute_input":"2023-06-06T11:46:33.826964Z","iopub.status.idle":"2023-06-06T11:48:11.300674Z","shell.execute_reply.started":"2023-06-06T11:46:33.82693Z","shell.execute_reply":"2023-06-06T11:48:11.29891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Prepare substitution matrix for skbio - format double dict\n\nTake blosum62 substitution matrix for protein's amino-acids from BioPython and convert it to skibio format (double dictionary)\n\nTo understand what format skbio is using for substitution matrices - we can look on:\nskibio has(depricated) make_identity_substitution_matrix function - which returns simple example of the substitution matrix\n\n","metadata":{}},{"cell_type":"code","source":"from skbio.alignment import make_identity_substitution_matrix\nsm2 = make_identity_substitution_matrix(1,-1)# , list(str(sm.alphabet)) )\nprint('substitution matrix for DNA:')\ndisplay(sm2 )\nprint()\n\nprint('substitution matrix for Protein:')\nfrom Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \nsm = substitution_matrices.load(\"BLOSUM62\") \nprint('blosum62[:3,:5]:')\nprint(np.array(sm)[:3,:5])\nsm2 = make_identity_substitution_matrix(1,-1 , list(str(sm.alphabet)) )\nprint(len(sm.alphabet), sm.alphabet)\nprint(list(sm['A'])[:3] )\n\n\nfrom Bio.Align import substitution_matrices # works for the new versions of BioPython instead of from Bio.SubsMat import MatrixInfo \nnames = substitution_matrices.load()\nprint(list(names) )\nsm = substitution_matrices.load(\"BLOSUM62\") \n# sm = np.array(sm)\n# print(type(sm), sm.shape )\n# print(np.array(sm)[:3,:5])\nprint( sm.alphabet )\nlist(str(sm.alphabet))\nll = list(str(sm.alphabet))\ndict_sm_blosum62 = {}\nfor 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()\nprint( dict_sm_blosum62['A'] )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T11:58:49.147876Z","iopub.execute_input":"2023-06-06T11:58:49.148264Z","iopub.status.idle":"2023-06-06T11:58:49.170922Z","shell.execute_reply.started":"2023-06-06T11:58:49.148236Z","shell.execute_reply":"2023-06-06T11:58:49.169638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Toy example\n\nSee also: https://www.kaggle.com/code/alexandervc/cafa5-19-alignments-skbio?scriptVersionId=131755045&cellId=11\n","metadata":{}},{"cell_type":"code","source":"%%time \n\nfrom skbio.alignment import StripedSmithWaterman\n\nseq1 = \"EVSAW\"\nseq2 = \"KEVLA\"\n\nquery = StripedSmithWaterman(seq1, substitution_matrix = dict_sm_blosum62, gap_open_penalty = 11, gap_extend_penalty = 1 )# \"ACTAAGGCTCTCTACCCCTCTCAGAGA\")\n# gap_open_penalty : int, optional The penalty applied to creating a gap in the alignment. This CANNOT be 0. Default is 5.\n# gap_extend_penalty : int, optional The penalty applied to extending a gap in the alignment. This CANNOT be 0. Default is 2.\n\nalignment = query(seq2) # \"AAAAAACTCTCTAAACTCACTAAGGCTCTCTACCCCTCTTCAGAGAAGTCGA\"\nprint(alignment)\nprint('aligned_query_sequence:', alignment.aligned_query_sequence, len(alignment.aligned_query_sequence))\nprint('aligned_target_sequence:',alignment.aligned_target_sequence, len(alignment.aligned_target_sequence))\nt = alignment\nprint( t['target_end_optimal'] - t['target_begin'], t['query_end'] - t['query_begin'] )\ndisplay(alignment )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T11:58:51.454909Z","iopub.execute_input":"2023-06-06T11:58:51.455295Z","iopub.status.idle":"2023-06-06T11:58:51.465317Z","shell.execute_reply.started":"2023-06-06T11:58:51.455258Z","shell.execute_reply":"2023-06-06T11:58:51.464352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom skbio.alignment import StripedSmithWaterman\n\nsimilarity_name = 'skbio Local Blossum62'\n\nquery = StripedSmithWaterman(sel_seq, substitution_matrix = dict_sm_blosum62, gap_open_penalty = 11, gap_extend_penalty = 1 )# \"ACTAAGGCTCTCTACCCCTCTCAGAGA\")\n# gap_open_penalty : int, optional The penalty applied to creating a gap in the alignment. This CANNOT be 0. Default is 5.\n# gap_extend_penalty : int, optional The penalty applied to extending a gap in the alignment. This CANNOT be 0. Default is 2.\n\nt0 = time.time()\nlist_align1 = [ query(s2) for s2 in list_seqs_corrected]# Alignment made here\nprint(len(list_align1), 'alignments for train done in %.1f secs'%(time.time()-t0))\nt0 = time.time()\nlist_align2 = [ query(s2) for s2 in list_seqs_corrected2]# Alignment made here\nprint(len(list_align2), 'alignments for test done in %.1f secs'%(time.time()-t0))\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:25:21.386794Z","iopub.execute_input":"2023-06-06T15:25:21.387286Z","iopub.status.idle":"2023-06-06T15:39:19.638267Z","shell.execute_reply.started":"2023-06-06T15:25:21.387253Z","shell.execute_reply":"2023-06-06T15:39:19.636008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nl1 = np.array([t['optimal_alignment_score'] for t in list_align1])\nl2 = np.array([t['optimal_alignment_score'] for t in list_align2] )\n\nd = pd.concat( (pd.Series(l1).describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]), pd.Series(l2).describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]) ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 17)\nplt.subplot(121)\nplt.hist(l1[l1<len(sel_seq)*5], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<len(sel_seq)*5], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\nprint('Attempt to show \"outlier\" part of the score distribution - which corresponds to biologically significant similarity. And start of the main part od the distribution - which is probably mostly random')\n# list_cutoffs_to_show_distribution_LocalAlignmentSkbio = [2200,2300, 2250 ] # SRC \n# list_cutoffs_to_show_distribution_LocalAlignmentSkbio = [1750, 1800, 1810,  1820] # P53\nfor N in list_cutoffs_to_show_distribution_LocalAlignmentSkbio: #  = [2200,2300, 2250 ] # SRC \n    plt.figure(figsize=(20,4))\n    plt.suptitle(similarity_name +' similarity to ' + prot_name + ' cutoff: '+str(N), fontsize = 17)\n    plt.subplot(121)\n    plt.hist(l1[(l1>N)&(l1<10_000)], bins=150)\n    plt.title('Train', fontsize = 15)\n    plt.subplot(122)\n    plt.hist(l2[(l2>N)&(l2<10_000)], bins=150)\n    plt.title('Test', fontsize = 15)\n    plt.show()\n    print('count train',(l1>N).sum(),'count test',(l2>N).sum() )\n\n\nprint(); print();\nfor t in [2200,2250, 2300]:#  range(560,580): # [10,50,100, 450]:\n    print(t, 'cumulative:', (l1>t).sum(),(l2>t).sum(), 'equals to', (l1==t).sum(),(l2==t).sum(), )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:47:26.55946Z","iopub.execute_input":"2023-06-06T15:47:26.560043Z","iopub.status.idle":"2023-06-06T15:47:33.446428Z","shell.execute_reply.started":"2023-06-06T15:47:26.560005Z","shell.execute_reply":"2023-06-06T15:47:33.444407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"di1 = pd.DataFrame(index = list_ids , data = np.array( [l1, list_lens]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi2 = pd.DataFrame(index = list_ids2 , data = np.array( [l2, list_lens2]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\n\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndisplay(di.iloc[:250,:])\nprint('400-420')\ndisplay(di.iloc[400:420,:])\nprint('4000-4020')\ndisplay(di.iloc[4000:4020,:])\nprint('40000-40020')\ndisplay(di.iloc[40000:40020,:])\n\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\n\n\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:47:52.110202Z","iopub.execute_input":"2023-06-06T15:47:52.110636Z","iopub.status.idle":"2023-06-06T15:47:53.192314Z","shell.execute_reply.started":"2023-06-06T15:47:52.110603Z","shell.execute_reply":"2023-06-06T15:47:53.191261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# window_len_for_keywords_indicator_moving_average = 1 # P53 #50 # SRC \nWL = window_len_for_keywords_indicator_moving_average\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # ['kinase', 'protein kinase']:\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:48:08.498656Z","iopub.execute_input":"2023-06-06T15:48:08.50017Z","iopub.status.idle":"2023-06-06T15:48:10.505168Z","shell.execute_reply.started":"2023-06-06T15:48:08.500116Z","shell.execute_reply":"2023-06-06T15:48:10.503641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compare similarity rankings - 2","metadata":{}},{"cell_type":"code","source":"dict_scoring_data.keys()","metadata":{"execution":{"iopub.status.busy":"2023-06-05T19:58:10.529Z","iopub.execute_input":"2023-06-05T19:58:10.530091Z","iopub.status.idle":"2023-06-05T19:58:10.540677Z","shell.execute_reply.started":"2023-06-05T19:58:10.530042Z","shell.execute_reply":"2023-06-05T19:58:10.538917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Levenshtein vs skbio local blossum62 ","metadata":{}},{"cell_type":"code","source":"%%time\nlist_similarities_to_compare = ['Levenshtein', 'skbio Local Blossum62'] # list( dict_scoring_data.keys() )\nl = np.arange(100,50000,100)\nl2 = []\nfor i0,N in enumerate(l):\n    for i1,k in enumerate( list_similarities_to_compare ):\n        dt = dict_scoring_data[k][0]\n        if i1 == 0:\n            s = set(dt.index[:N])\n        else:\n            s = s & set(dt.index[:N])\n        if i0 == 0:\n            print(k)\n            display(dt.head(2))\n            \n    l2.append(len(s))\n    \nprint(len(l2), l2[:10])\n\n\ndeg = len( list_similarities_to_compare )\npol = np.polyfit(l,l2, deg)\n\nN = 50\npol2 = np.polyfit(l[:N],l2[:N], 1) \nprint('Percent of agreement:', np.round(100*pol2[0],2) )\nprint(); print(); \n\nl = np.asarray(l)\nl2 = np.asarray(l2)\n\nprint('Linear approximating polynom  %+.5f x %+.1f :'%(pol2[0], pol2[1] ) )\nprint('Parabolic approximating polynom %-.10f x^2 %+.1f x %+.1f :'%(pol[0], pol[1], pol[2] ) )\n\nplt.figure(figsize = (20,6))\nplt.subplot(121)\nplt.plot(l,l2 , label = 'Data')\nplt.plot(l, np.polyval(pol,l ) , label = 'Approximation Parabola' )\nplt.plot(l, np.polyval(pol2,l ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N] , label = 'Data' )\nplt.plot(l[:N], np.polyval(pol,l[:N]) , label = 'Approximation Parabola' )\nplt.plot(l[:N], np.polyval(pol2,l[:N] ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\nplt.show()\n\n\nplt.figure(figsize = (20,6))\nplt.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')\nplt.subplot(121)\nplt.plot(l,l2/l , label = 'Data Normalized')\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N]/l[:N] , label = 'Data Normalized' )\nplt.legend()\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:48:26.506871Z","iopub.execute_input":"2023-06-06T15:48:26.507303Z","iopub.status.idle":"2023-06-06T15:48:38.1815Z","shell.execute_reply.started":"2023-06-06T15:48:26.507269Z","shell.execute_reply":"2023-06-06T15:48:38.180051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## BioPython local blossum62 vs skbio local blossum62 ","metadata":{}},{"cell_type":"code","source":"%%time\nlist_similarities_to_compare = ['BioPython Local Blossum62', 'skbio Local Blossum62'] # list( dict_scoring_data.keys() )\nl = np.arange(100,50000,100)\nl2 = []\nfor i0,N in enumerate(l):\n    for i1,k in enumerate( list_similarities_to_compare ):\n        dt = dict_scoring_data[k][0]\n        if i1 == 0:\n            s = set(dt.index[:N])\n        else:\n            s = s & set(dt.index[:N])\n        if i0 == 0:\n            print(k)\n            display(dt.head(2))\n            \n    l2.append(len(s))\n    \nprint(len(l2), l2[:10])","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:15:35.603699Z","iopub.execute_input":"2023-06-05T18:15:35.604622Z","iopub.status.idle":"2023-06-05T18:15:41.350331Z","shell.execute_reply.started":"2023-06-05T18:15:35.604594Z","shell.execute_reply":"2023-06-05T18:15:41.349254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndeg = len( list_similarities_to_compare )\n\nstr_inf = 'Compare: ' + str(list_similarities_to_compare)\n\npol = np.polyfit(l,l2, deg)\n\nN = 50\npol2 = np.polyfit(l[:N],l2[:N], 1) \nprint('Percent of agreement:', np.round(100*pol2[0],2) )\nprint(); print(); \n\nl = np.asarray(l)\nl2 = np.asarray(l2)\n\nprint('Linear approximating polynom  %+.5f x %+.1f :'%(pol2[0], pol2[1] ) )\nprint('Parabolic approximating polynom %-.10f x^2 %+.1f x %+.1f :'%(pol[0], pol[1], pol[2] ) )\n\nplt.figure(figsize = (20,6))\nplt.suptitle(str_inf, fontsize = 17 )\nplt.subplot(121)\nplt.plot(l,l2 , label = 'Data')\nplt.plot(l, np.polyval(pol,l ) , label = 'Approximation Parabola' )\nplt.plot(l, np.polyval(pol2,l ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N] , label = 'Data' )\nplt.plot(l[:N], np.polyval(pol,l[:N]) , label = 'Approximation Parabola' )\nplt.plot(l[:N], np.polyval(pol2,l[:N] ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\nplt.show()\n\n\nplt.figure(figsize = (20,6))\nplt.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')\nplt.subplot(121)\nplt.plot(l,l2/l , label = 'Data Normalized')\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N]/l[:N] , label = 'Data Normalized' )\nplt.legend()\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:15:41.351572Z","iopub.execute_input":"2023-06-05T18:15:41.351887Z","iopub.status.idle":"2023-06-05T18:15:42.280095Z","shell.execute_reply.started":"2023-06-05T18:15:41.351861Z","shell.execute_reply":"2023-06-05T18:15:42.27901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# GO-Terms similarity vs sequence similarity  analysis ","metadata":{}},{"cell_type":"code","source":"%%time\n# fn =  str(path) + '/Train/train_terms.tsv'\nfn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv'\ndf = pd.read_csv(fn , sep = '\\t')# , index_col = 0)\ndf","metadata":{"execution":{"iopub.status.busy":"2023-06-05T19:59:38.652389Z","iopub.execute_input":"2023-06-05T19:59:38.65283Z","iopub.status.idle":"2023-06-05T19:59:42.765394Z","shell.execute_reply.started":"2023-06-05T19:59:38.652798Z","shell.execute_reply":"2023-06-05T19:59:42.763999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmap_id2GOterms = df.groupby('EntryID')['term'].apply(lambda x: list(x))\nmap_id2GOterms.head(2)","metadata":{"execution":{"iopub.status.busy":"2023-06-05T19:59:42.769324Z","iopub.execute_input":"2023-06-05T19:59:42.769762Z","iopub.status.idle":"2023-06-05T19:59:50.507472Z","shell.execute_reply.started":"2023-06-05T19:59:42.769721Z","shell.execute_reply":"2023-06-05T19:59:50.505977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nN = 140_000\ndf_res = pd.DataFrame()\nlist_similarities_to_compare = list(dict_scoring_data.keys() )\nfor i1,k in enumerate( list_similarities_to_compare  ):\n    dt = dict_scoring_data[k][0]\n    print(k)\n    display(dt.head(2))\n    s0 = set( map_id2GOterms[ dt.index[0] ] )\n    print(len(s0))\n    l = []\n    for nn in range(N):\n        s = set( map_id2GOterms[ dt.index[nn] ] )\n        dice = 2 * len( s & s0 ) / ( len(s) + len(s0) )\n        l.append(dice)\n    df_res[k] = l\n    \ndf_res    ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:49:36.381239Z","iopub.execute_input":"2023-06-06T15:49:36.381741Z","iopub.status.idle":"2023-06-06T15:49:40.449556Z","shell.execute_reply.started":"2023-06-06T15:49:36.381706Z","shell.execute_reply":"2023-06-06T15:49:40.447716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_res.corr()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:49:49.047464Z","iopub.execute_input":"2023-06-06T15:49:49.048301Z","iopub.status.idle":"2023-06-06T15:49:49.066525Z","shell.execute_reply.started":"2023-06-06T15:49:49.04826Z","shell.execute_reply":"2023-06-06T15:49:49.065151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfor N in [100,1000,10000,100_000]:\n    plt.figure(figsize = (20,6))\n    plt.title('Dice between GO-terms for '+prot_name+' vs i-th similar protein. Cutoff: '+str(N), fontsize = 20 )\n    for col in df_res.columns:\n        plt.plot(df_res[col].iloc[:N].rolling(window=50).mean(), label = col)\n    plt.legend()\n    plt.grid()\n    plt.show()\n\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:49:49.556379Z","iopub.execute_input":"2023-06-06T15:49:49.556915Z","iopub.status.idle":"2023-06-06T15:49:51.693078Z","shell.execute_reply.started":"2023-06-06T15:49:49.556873Z","shell.execute_reply":"2023-06-06T15:49:51.691821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfor N in [100,1000,10000,100_000]:\n    plt.figure(figsize = (20,6))\n    plt.title('Dice between GO-terms for '+prot_name+'  vs i-th similar protein. Cutoff: '+str(N), fontsize = 20 )\n    for col in df_res.columns:\n        plt.plot(df_res[col].iloc[:N].rolling(window=5).mean(), label = col)\n    plt.legend()\n    plt.grid()\n    plt.show()\n\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:10.665539Z","iopub.execute_input":"2023-06-06T15:50:10.666665Z","iopub.status.idle":"2023-06-06T15:50:13.71364Z","shell.execute_reply.started":"2023-06-06T15:50:10.666615Z","shell.execute_reply":"2023-06-06T15:50:13.711933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity by T5 embeddings ","metadata":{}},{"cell_type":"code","source":"%%time\nsimilarity_name = 'T5embeds PearCorr'\n\nfn = '/kaggle/input/t5embeds/train_embeds.npy'\nX = np.load(fn)\nprint(X.shape)\nX\n\nfn = '/kaggle/input/t5embeds/train_ids.npy'\nvec_train_protein_ids = np.load(fn)\nprint(vec_train_protein_ids.shape)\nvec_train_protein_ids\n\nfn = '/kaggle/input/t5embeds/test_embeds.npy'\nX2 = np.load(fn)\nprint(X2.shape)\nX2\n\nfn = '/kaggle/input/t5embeds/test_ids.npy'\nvec_test_protein_ids = np.load(fn)\nprint(vec_test_protein_ids.shape)\nvec_test_protein_ids\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:15.05659Z","iopub.execute_input":"2023-06-06T15:50:15.057708Z","iopub.status.idle":"2023-06-06T15:50:18.202258Z","shell.execute_reply.started":"2023-06-06T15:50:15.057648Z","shell.execute_reply":"2023-06-06T15:50:18.200906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nfrom sklearn.feature_selection import r_regression\n\nIX = np.where( vec_train_protein_ids == sel_prot_id)[0]\nsel_vec1 = X[IX[0], :]\nprint(sel_vec1.shape)\nIX = np.where( vec_test_protein_ids == sel_prot_id)[0]\nsel_vec2 = X2[IX[0], :]\nprint(sel_vec2.shape)\nprint('Check should be zero:', (sel_vec2 != sel_vec1).sum()  )\n\nl1 = r_regression(X.T ,sel_vec1 )\nl2 = r_regression(X2.T ,sel_vec2 )\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:18.203911Z","iopub.execute_input":"2023-06-06T15:50:18.20482Z","iopub.status.idle":"2023-06-06T15:50:19.175178Z","shell.execute_reply.started":"2023-06-06T15:50:18.204781Z","shell.execute_reply":"2023-06-06T15:50:19.173339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nd = pd.concat( (pd.Series(l1).describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]), pd.Series(l2).describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]) ) , axis = 1 )\nd.columns = ['train', 'test'] \ndisplay(d)\n\nplt.figure(figsize=(20,4))\nplt.suptitle(similarity_name +' similarity to ' + prot_name, fontsize = 17)\nplt.subplot(121)\nplt.hist(l1[l1<2500], bins=150)\nplt.title('Train', fontsize = 20)\nplt.subplot(122)\nplt.hist(l2[l2<2500], bins=150)\nplt.title('Test', fontsize = 20)\nplt.show()\n\nfor N in [0.7, 0.8]:\n    plt.figure(figsize=(20,4))\n    plt.suptitle(similarity_name +' similarity to ' + prot_name + ' cutoff: '+str(N), fontsize = 17)\n    plt.subplot(121)\n    plt.hist(l1[(l1>N)&(l1<10_000)], bins=150)\n    plt.title('Train', fontsize = 15)\n    plt.subplot(122)\n    plt.hist(l2[(l2>N)&(l2<10_000)], bins=150)\n    plt.title('Test', fontsize = 15)\n    plt.show()\n    print('count train',(l1>N).sum(),'count test',(l2>N).sum() )\n\n\nprint(); print();\nfor t in [0.7, 0.75, 0.76, 0.8,0.9,0.95,0.99,0.999]:#  range(560,580): # [10,50,100, 450]:\n    print(t, 'cumulative:', (l1>t).sum(),(l2>t).sum(), 'equals to', (l1==t).sum(),(l2==t).sum(), )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:19.496859Z","iopub.execute_input":"2023-06-06T15:50:19.497246Z","iopub.status.idle":"2023-06-06T15:50:22.226967Z","shell.execute_reply.started":"2023-06-06T15:50:19.497217Z","shell.execute_reply":"2023-06-06T15:50:22.225002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"di1 = pd.DataFrame(index = list_ids , data = np.array( [l1, list_lens]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\ndi2 = pd.DataFrame(index = list_ids2 , data = np.array( [l2, list_lens2]).T , columns = ['score', 'len']).sort_values('score', ascending = False)\n\ndi1.to_csv('df_train_'+similarity_name.replace(' ','_') + '.csv')\ndi2.to_csv('df_test_'+similarity_name.replace(' ','_') + '.csv')\n\nN = 50000\ndi = di1.reset_index().head(N).join(   di2.reset_index().head(N) , lsuffix='_train', rsuffix='_test' )\nl = [dict_ids2description[t] for t in di['index_train']]\ndi['train description'] = l\ndisplay(di.iloc[:250,:])\nprint('400-420')\ndisplay(di.iloc[400:420,:])\nprint('4000-4020')\ndisplay(di.iloc[4000:4020,:])\nprint('40000-40020')\ndisplay(di.iloc[40000:40020,:])\n\ndi.to_csv('df_train_test_'+similarity_name.replace(' ','_') + '.csv')\n\n\ndict_scoring_data[similarity_name] = (di1,di2,di) # Save for possible future use ","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:22.381936Z","iopub.execute_input":"2023-06-06T15:50:22.383122Z","iopub.status.idle":"2023-06-06T15:50:24.074156Z","shell.execute_reply.started":"2023-06-06T15:50:22.383073Z","shell.execute_reply":"2023-06-06T15:50:24.07252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nWL = window_len_for_keywords_indicator_moving_average # 50 # SRC # 1 P53\n\nstr_inf = similarity_name +' similarity to ' + prot_name  + '\\n'\nfor kw in list_keywords_describing_protein: # ['kinase', 'protein kinase']:\n    l_ = [kw.lower() in str(t).lower() for t in di['train description'].values ]\n    v = pd.Series(l_).rolling(window=WL).mean()\n    plt.figure(figsize = (20,5))\n    plt.suptitle(str_inf + 'presence of \"'+kw + '\" in desription. Moving average '+str(WL), fontsize = 16 )\n    plt.subplot(121)\n    N=200\n    plt.plot(v[:N])\n    plt.title('Cutoff '+str(N)+'        ', fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N=5000\n    plt.plot(v[:N])\n    plt.title(' Cutoff '+str(N) +'', fontsize = 16 )\n    plt.grid()\n    \n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:31.182161Z","iopub.execute_input":"2023-06-06T15:50:31.182574Z","iopub.status.idle":"2023-06-06T15:50:32.927589Z","shell.execute_reply.started":"2023-06-06T15:50:31.182543Z","shell.execute_reply":"2023-06-06T15:50:32.926088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# GO-Terms similarity vs sequence/embedding similarity analysis - 2 ","metadata":{}},{"cell_type":"code","source":"%%time\n# fn =  str(path) + '/Train/train_terms.tsv'\nfn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv'\ndf = pd.read_csv(fn , sep = '\\t')# , index_col = 0)\ndisplay(df.head(2))\nmap_id2GOterms = df.groupby('EntryID')['term'].apply(lambda x: list(x))\nmap_id2GOterms.head(2)","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:50:44.001364Z","iopub.execute_input":"2023-06-06T15:50:44.001862Z","iopub.status.idle":"2023-06-06T15:50:52.835576Z","shell.execute_reply.started":"2023-06-06T15:50:44.001805Z","shell.execute_reply":"2023-06-06T15:50:52.834663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print( list(dict_scoring_data.keys() ) )","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:06:12.068085Z","iopub.execute_input":"2023-06-06T15:06:12.068506Z","iopub.status.idle":"2023-06-06T15:06:12.074066Z","shell.execute_reply.started":"2023-06-06T15:06:12.06847Z","shell.execute_reply":"2023-06-06T15:06:12.072958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nN = 140_000\ndf_res = pd.DataFrame()\nlist_similarities_to_compare = list(dict_scoring_data.keys() )\nfor i1,k in enumerate( list_similarities_to_compare  ):\n    dt = dict_scoring_data[k][0]\n#     print(k)\n#     display(dt.head(2))\n    s0 = set( map_id2GOterms[ dt.index[0] ] )\n    print(k, len(s0))\n    l = []\n    for nn in range(N):\n        s = set( map_id2GOterms[ dt.index[nn] ] )\n        dice = 2 * len( s & s0 ) / ( len(s) + len(s0) )\n        l.append(dice)\n    df_res[k] = l\n    \ndisplay( df_res.head(2) )\n\ndisplay( df_res.corr() )\n\nfor WL in [50,5]:\n    for N in [100,1000,10000,100_000]:\n        plt.figure(figsize = (20,6))\n        plt.title('Dice between GO-terms for selected protein vs i-th similar protein. Cutoff: '+str(N) + ' WL '+str(WL), fontsize = 20 )\n        for col in df_res.columns:\n            plt.plot(df_res[col].iloc[:N].rolling(window=WL).mean(), label = col)\n        plt.legend()\n        plt.grid()\n        plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:51:00.296942Z","iopub.execute_input":"2023-06-06T15:51:00.298242Z","iopub.status.idle":"2023-06-06T15:51:12.233033Z","shell.execute_reply.started":"2023-06-06T15:51:00.298183Z","shell.execute_reply":"2023-06-06T15:51:12.231639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compare similarity rankings - 3","metadata":{}},{"cell_type":"code","source":"dict_scoring_data.keys()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:51:55.160698Z","iopub.execute_input":"2023-06-06T15:51:55.161555Z","iopub.status.idle":"2023-06-06T15:51:55.170431Z","shell.execute_reply.started":"2023-06-06T15:51:55.161511Z","shell.execute_reply":"2023-06-06T15:51:55.16891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Levenshtein vs T5 embeds","metadata":{}},{"cell_type":"code","source":"%%time\nlist_similarities_to_compare = ['Levenshtein','T5embeds PearCorr']#  'skbio Local Blossum62'] # list( dict_scoring_data.keys() )\nl = np.arange(100,50000,100)\nl2 = []\nfor i0,N in enumerate(l):\n    for i1,k in enumerate( list_similarities_to_compare ):\n        dt = dict_scoring_data[k][0]\n        if i1 == 0:\n            s = set(dt.index[:N])\n        else:\n            s = s & set(dt.index[:N])\n        if i0 == 0:\n            print(k)\n            display(dt.head(2))\n            \n    l2.append(len(s))\n    \nprint(len(l2), l2[:10])\n\n\n\ndeg = len( list_similarities_to_compare )\npol = np.polyfit(l,l2, deg)\n\nN = 50\npol2 = np.polyfit(l[:N],l2[:N], 1) \nprint('Percent of agreement:', np.round(100*pol2[0],2) )\nprint(); print(); \n\nl = np.asarray(l)\nl2 = np.asarray(l2)\n\nprint('Linear approximating polynom  %+.5f x %+.1f :'%(pol2[0], pol2[1] ) )\nprint('Parabolic approximating polynom %-.10f x^2 %+.1f x %+.1f :'%(pol[0], pol[1], pol[2] ) )\n\nplt.figure(figsize = (20,6))\nplt.subplot(121)\nplt.plot(l,l2 , label = 'Data')\nplt.plot(l, np.polyval(pol,l ) , label = 'Approximation Parabola' )\nplt.plot(l, np.polyval(pol2,l ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N] , label = 'Data' )\nplt.plot(l[:N], np.polyval(pol,l[:N]) , label = 'Approximation Parabola' )\nplt.plot(l[:N], np.polyval(pol2,l[:N] ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\nplt.show()\n\n\nplt.figure(figsize = (20,6))\nplt.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')\nplt.subplot(121)\nplt.plot(l,l2/l , label = 'Data Normalized')\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N]/l[:N] , label = 'Data Normalized' )\nplt.legend()\nplt.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-06T15:51:56.879096Z","iopub.execute_input":"2023-06-06T15:51:56.879542Z","iopub.status.idle":"2023-06-06T15:52:09.193399Z","shell.execute_reply.started":"2023-06-06T15:51:56.879505Z","shell.execute_reply":"2023-06-06T15:52:09.192165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 'BioPython Local Blossum62' vs 'T5embeds PearCorr'","metadata":{}},{"cell_type":"code","source":"%%time\nlist_similarities_to_compare = ['BioPython Local Blossum62','T5embeds PearCorr']#  'skbio Local Blossum62'] # list( dict_scoring_data.keys() )\nl = np.arange(100,50000,100)\nl2 = []\nfor i0,N in enumerate(l):\n    for i1,k in enumerate( list_similarities_to_compare ):\n        dt = dict_scoring_data[k][0]\n        if i1 == 0:\n            s = set(dt.index[:N])\n        else:\n            s = s & set(dt.index[:N])\n        if i0 == 0:\n            print(k)\n            display(dt.head(2))\n            \n    l2.append(len(s))\n    \nprint(len(l2), l2[:10])\n\n\n\ndeg = len( list_similarities_to_compare )\npol = np.polyfit(l,l2, deg)\n\nN = 50\npol2 = np.polyfit(l[:N],l2[:N], 1) \nprint('Percent of agreement:', np.round(100*pol2[0],2) )\nprint(); print(); \n\nl = np.asarray(l)\nl2 = np.asarray(l2)\n\nprint('Linear approximating polynom  %+.5f x %+.1f :'%(pol2[0], pol2[1] ) )\nprint('Parabolic approximating polynom %-.10f x^2 %+.1f x %+.1f :'%(pol[0], pol[1], pol[2] ) )\n\nplt.figure(figsize = (20,6))\nplt.subplot(121)\nplt.plot(l,l2 , label = 'Data')\nplt.plot(l, np.polyval(pol,l ) , label = 'Approximation Parabola' )\nplt.plot(l, np.polyval(pol2,l ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N] , label = 'Data' )\nplt.plot(l[:N], np.polyval(pol,l[:N]) , label = 'Approximation Parabola' )\nplt.plot(l[:N], np.polyval(pol2,l[:N] ) , label = 'Approximation Linear' )\nplt.legend()\nplt.grid()\nplt.show()\n\n\nplt.figure(figsize = (20,6))\nplt.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')\nplt.subplot(121)\nplt.plot(l,l2/l , label = 'Data Normalized')\nplt.legend()\nplt.grid()\n\nN = 100\nplt.subplot(122)\nplt.title( 'Cutoff: '+str(N) )\nplt.plot(l[:N],l2[:N]/l[:N] , label = 'Data Normalized' )\nplt.legend()\nplt.grid()\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2023-06-05T18:15:59.313845Z","iopub.execute_input":"2023-06-05T18:15:59.314161Z","iopub.status.idle":"2023-06-05T18:15:59.32018Z","shell.execute_reply.started":"2023-06-05T18:15:59.314134Z","shell.execute_reply":"2023-06-05T18:15:59.3181Z"},"trusted":true},"execution_count":null,"outputs":[]}]}