{"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\nHere we split protein sequences from the train part of the CAFA5 competition to the similarity groups, and  create groupwise CV folds.\nThe key ingredient is precomputed local alignment similarity by DIAMOND package. The computation is done in the separate notebook: [to-do link to notebook] (Thanks to Liza Geraseva for highlighting the DIAMOND: https://www.kaggle.com/code/geraseva/diamond ). \n\nIt is the standard approach to group similar samples into the same folds, otherwise the CV-score of the model might not reflect the performance on some  unseen samples. Some experiments showed that random CV for the train does not perfectly fit LB, it might be the reason is that random CV is not perfect, thus we need to try groupwise CV.\n\n**Setup:** We consider several similarities between samples (protein sequences) - the key one - similarity by DIAMOND blast-like local alignment package (which is faster than blastP). \n\n**Key points:** rougly speaking if proteins are similar (with respect to any of the considered similarities) we want to put them into one group. To achieve that - it is more or less the same as to cluster the vertices of the following graph:  vertices - samples (proteins), connected by an edge if there is  similarity relation between them. Connected components of that graph are exactly the desired groups. More relaxed version of connected compenents - clusters obtained by some graph clustering algorithm. We use Louvain clustering (we do not expect much differences from other choices). That are the main points. \n\n**Outline:** Having several similarity relations - for  each one we create list of pairs of samples (proteins) which are similar with respect to that relation. We join these lists. Create graph and cluster it. \n\n**Results/output:** we save csv files with group markers for each sample (protein), as well as another csv file with several  options for the folds markers - from 2 folds to 20 folds.  To give a flavor what groups we obtain:\nsome characteristics: n_clusters 5178, Sizes of top clusters: [6274, 5290, 5255, 5129, 4554..],\nmost frequent size clusters: {1: 3107, 2: 833, 3: 367, 4: 216... }. See details below. \n\n**Why not alignment similarity alone?** why do we consider several relations ? Alignment similarity is the key one, and most probably it alone is enough. However we are not 100% sure, and prefer to have an option to add additional relations to strengthen the main one. The reasons to think that it alone is not the perfect choice are the following. We observed that DIAMOND does not work well for e.g. very short peptides like \"Tachykinin peptides\" (about 10 amino acids) - they were not found to be similar by Diamond - despite expected to be so; and in general it is sub-optimial blast-like algorithm not so sensitive like optimal Smith-Waterman local alignment (which is however much slower - would take weeks - instead of 3 hours by DIAMOND).  Thus we prefer to have an option to strengthen the main similarity relation.\n(These steps are optional - can be easily dropped out - see \"key params\" section - we can choose any subpart of similarity relations).\n\n\n**Similarity relations:** \n\n    duplicated sequences\n    some special groups like \"Tachykinin peptides\"\n    proteins with the same gene name\n    proteins with the similar gene name\n    proteins with the similar description\n    KEY: proteins which are similar by DIAMOND blastP local alignment (we use \"ultra-senstity\" mode - the most sensitive mode of DIAMOND - still not as senstive as Smith-Waterman (SW) local alignment, but quite sensitive and works in reasonable time - about 3 hours, while SW will take weeks). \n    \n    As mentioned above: \n    After that we create graph with vertices - samples (proteins), connected by edge - if any of similarities exist   (no weights are used in the current version - just connected or not.)\n    After that we cluster the graph by e.g. Louvain algorithm and get the groups. \n    \n   \n**Why not the standard way?** Briefly - it seems it is not sensitive enough. \nWe  obtain larger groups  comparing to more standard ways to cluster sequences by CD-HIT, MMseq, DIAMOND-clustering-mode etc.\nOur analysis showed [to-do link to notebook] that groups obtained by these methods are not so big, i.e. they group only very similar proteins, but not capturing middle-strength similarity. Packages do it intentionally for the sake of speed.\nSo an idea -  a bit sacrifice speed to get more sensitive similarity relations (i.e. larger groups). To add more details - we see that some proteins like SRC protein tyrosine kinase have about 4000 proteins in train which are statistically significantly similar to it, but clustering methods CD-HIT, MMseq or DIAMOND-clustering-mode give groups sized at most hundrends -  [to-do link to notebook] - thus order of magnitude lower than really existing similarity. \n\nAny way we have parameters in the notebook - like e-value-threshold for DIAMOND matches - so we can increase or decrease the sensitivity of the imposed relations. See section \"Key params\".\n\n**Two ways to get groups from graph:** To obtain groups (clusters) from graph there are two ways: the simple and most strong - just to take connected components (it is enough when we do not have many similarity relations), the second way - to take clusters by some graph clustering algorithm - e.g. Louvain (quite good and fast algorithm - standard choice for many tasks). The problem with the first way is that we can obtain \"giant connected component(s)\" i.e. some of group(s) would be extremely large and that is not good for creating folds, that indeed happens even for our key relation - DIAMOND similarity alone. Thus we prefer Louvain clustering. So, in the other words, taking connected components we obtain too strong similarity relation - nearly to \"everything similar to everything\" and we need to relax it by taking Louvain clusters. \n\n\n**Is our sensitivity enough?** How to improve ? Sizes of our largest groups  - 7k,5k, etc - that is quite enough since our estimate about 4000 (it was obtained by sensitive Smith-Waterman on part of train). However even imposing all our similirity relations we see that some (thousands) proteins are \"orphans\" i.e. their group size is 1. Something like that is expected,  but still  it is not clear - either it is really the truth, or it might be that DIAMOND  is not sensitive enough. There are ways to improve, but we do not currently pursue them, since we do not expect huge improvements. These ways are the following: one may try to run more sensitive local Smith-Waterman (e.g. Parasail package: https://www.kaggle.com/code/alexandervc/protein-aligners-benchmark-parasail-diamond-etc ) for these particular proteins and try to find proteins similar to them and join them into the groups. Or to use the embeddings to get similarities by these methods. In general it might be better to have more uniformly distributed sizes of the groups, nevetheless the current approach migth be a good for start. \n\nIn general short proteins seems to be \"underclustered\" , while long - \"over-clustered\". It might be because it is more easy to get similarity by random chance for short proteins - and thus DIAMOND estimating e-value or making other test inside will drop them out. Despite we run DIAMOND with extremely large e-value threshold  1_000_000 , it seems still does not prevent DIAMOND to drop out some short proteins. \n\n**Why not to cluster embeddings?** Well, it might be even better and more simple way ! However it is not clear at the moment. Protein embeddings are new technology, and before using them it might be more natural to have a baseline by more standard technique of protein sequence analysis. It would be interesting to compare the methods. \n\n**How to improve further:** It might be interesting to think towards iterative stratified folds, somehow combined with approach here.  \n\nPS\n\nBackground: Sequences alignments are extremely important algorithms used in modern bioinformatics: https://en.wikipedia.org/wiki/Sequence_alignment . For gentle introduction to alignments with BioPython one can see: https://www.kaggle.com/code/shtrausslearning/biological-sequence-alignment and other Kaggle notebooks by the same author.\nMore aligners and benchmarks: https://www.kaggle.com/code/alexandervc/protein-aligners-benchmark-parasail-diamond-etc\n\nDIAMOND: 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\nThanks to Liza Geraseva's notebook on diamond: https://www.kaggle.com/code/geraseva/diamond - please upvote !","metadata":{}},{"cell_type":"markdown","source":"## Versions\n\n### 9-... relaxing folds\n    9 evalue_threshold = 0.0001 \n    10 evalue_threshold = 0.00001 \n    11 evalue_threshold = 0.00001 \n    12 evalue_threshold = 0.000001 \n    13 evalue_threshold = 0.000001 , less similarities\n    14 evalue_threshold = 0.000001 , even more less similarities\n    15 evalue_threshold = 1e-10 , 3 sims\n    16 evalue_threshold = 1e-20 , 3 sims\n    17 evalue_threshold = 1e-30 , 3 sims\n    18 evalue_threshold = 1e-50 , 3 sims\n    \n    19 evalue_threshold = 1e-6 , 5 sims (except: genes names cutted digits) \n\n\n### 1-8 basic version with \"strict\" folds \n\n### Detailed results on how relaxed groups we have depending on e-value and similarities: \n\n    3sims = 'seq duplicate', 'special groups',   'diamond local alignment'\n\n    V18:  18 evalue_threshold = 1e-50 , 3 sims\n    n_clusters: 47353 Sizes of clusters: [881, 816, 649, ... \n\n    V17: evalue_threshold = 1e-30 , 3 sims\n    n_clusters: 37157 Sizes of clusters: [1502, 1429, 941, ...\n    Value counts of component sizes:  1: 23782, 2: 4660,  3: 2168, ... \n\n    V16:  evalue_threshold = 1e-20 , 3 sims\n    n_clusters: 31191 Sizes of clusters: 2570, 1966, 1556, ...\n    Value counts of component sizes: 1: 19345, 2: 3924, 3: 1876,\n\n    V15: evalue_threshold = 1e-10 , 3 sims\n    n_clusters: 23702 Sizes of clusters: 4301, 2927, 1540, ...\n    Value counts of component sizes: 1: 14260, 2: 2984, 3: 1434, ...  \n\n    V14: evalue_threshold = 1e-7 , 3 sims\n    n_clusters: 20776 Sizes of clusters: 4426, 4117, 3075, ... \n    Value counts of component sizes: 1: 12443, 2: 2628, 3: 1230,\n\n    V13: evalue_threshold = 1e-7, 4sims: ['seq duplicate', 'special groups', 'gene name',  'diamond local alignment'\n    n_clusters: 17252 Sizes of clusters: 4455, 3658, 3083, ... \n    Value counts of component sizes: 1: 10483, 2: 2279, 3: 1029, ...\n\n    THAT WAS TESTED: \n    V12: evalue_threshold = 1e-7, 6sims (full list)\n    n_clusters: 6464  Sizes of clusters: [5515, 4684, 4482,...\n    Value counts of component sizes: 1: 3867, 2: 1024, 3: 467, 4: 296,\n\n    V11: evalue_threshold = 1e-6, 6sims (full list)\n    n_clusters: 6179 Sizes of clusters: 6878, 4801, 4682,...\n    Value counts of component sizes: 1: 3692, 2: 989, 3: 448,...\n\n    V10: evalue_threshold = 1e-5, 6sims (full list)\n    n_clusters: 5857 Sizes of clusters: 6818, 5518, 4911,...\n    Value counts of component sizes: 1: 3501, 2: 935, 3: 416,... \n\n    V9: evalue_threshold = 1e-4, 6sims (full list)\n    n_clusters: 5518  Sizes of clusters: 5353, 5093, 5012,...\n    Value counts of component sizes: 1: 3317, 2: 872, 3: 388,...\n\n    V8: evalue_threshold = 1e-3, 6sims (full list)\n    n_clusters: 5186 Sizes of clusters: 7346, 6870, 5611,... \n    Value counts of component sizes: 1: 3108, 2: 833, 3: 367,...\n","metadata":{}},{"cell_type":"markdown","source":"# Key params \n\n","metadata":{}},{"cell_type":"code","source":"evalue_threshold = 1e-6 # 0.001 is default threshold for DIAMOND \n# smaller evalue_threshold - give smaller groups, higher - bigger (but we do not expect much change from 0.001 - somehow Diamond by itself is restrictive on the output - and will not return too many hits despite all our efforts )\n# So higher values - increase sensitivity, while lower - decrease \n\n\nlist_allowed_similarity_relations = ['seq duplicate', 'special groups',  \n 'gene name',  'description',\n 'diamond local alignment' ]# 'gene name', 'gene name end digits cutted', 'description',\n\n# Full list currently: \n#  ['seq duplicate', 'special groups', 'gene name', \n# 'gene name end digits cutted', 'description', 'diamond local alignment' ]\n\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:49:40.956823Z","iopub.execute_input":"2023-07-10T11:49:40.957331Z","iopub.status.idle":"2023-07-10T11:49:40.962435Z","shell.execute_reply.started":"2023-07-10T11:49:40.957293Z","shell.execute_reply":"2023-07-10T11:49:40.9616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparations and data load","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport time\nt0start = time.time()\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n\n# 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-07-10T11:49:39.547834Z","iopub.execute_input":"2023-07-10T11:49:39.548561Z","iopub.status.idle":"2023-07-10T11:49:40.954872Z","shell.execute_reply.started":"2023-07-10T11:49:39.548527Z","shell.execute_reply":"2023-07-10T11:49:40.95373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load CAFA5 train data","metadata":{}},{"cell_type":"code","source":"%%time \n# --------------------------------------------------------------------------------------------------------------\n# ############################################# Load CAFA5 train data ##################################\n# --------------------------------------------------------------------------------------------------------------\nfrom Bio import SeqIO\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()\nprint(len(list_ids), len(set(list_ids))) # check no duplicates for ids  - okya for train \n\nmap_id2num =  {seq_id:i for  i, seq_id in enumerate(list_ids) }\nprint( list(map_id2num.keys())[:2], list(map_id2num.values())[:2] )\n\ndf_seqs = pd.DataFrame(data = [list_ids, list_lens, list_seqs]).T\ndf_seqs.columns = ['id','len','seq']\ndf_seqs.head(3)","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:49:49.655267Z","iopub.execute_input":"2023-07-10T11:49:49.655645Z","iopub.status.idle":"2023-07-10T11:49:59.893195Z","shell.execute_reply.started":"2023-07-10T11:49:49.655617Z","shell.execute_reply":"2023-07-10T11:49:59.892189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load parsed info from fasta train","metadata":{}},{"cell_type":"code","source":"%%time\ndf_eda = pd.read_csv('/kaggle/input/cafa5-features-etc/df_train_eda.csv')\nprint(df_eda.shape)\ndf_eda","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:53:00.704834Z","iopub.execute_input":"2023-07-10T11:53:00.706225Z","iopub.status.idle":"2023-07-10T11:53:02.679579Z","shell.execute_reply.started":"2023-07-10T11:53:00.706172Z","shell.execute_reply":"2023-07-10T11:53:02.678442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Check - should be zero:')\n(df_eda['id'].values != df_seqs['id'].values).sum()\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:53:08.689723Z","iopub.execute_input":"2023-07-10T11:53:08.690096Z","iopub.status.idle":"2023-07-10T11:53:08.703515Z","shell.execute_reply.started":"2023-07-10T11:53:08.690068Z","shell.execute_reply":"2023-07-10T11:53:08.702213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Storages for key results ","metadata":{}},{"cell_type":"code","source":"%%time\ndict_list_edges_for_graph = {} # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \nlist_edges_for_graph = [] # here is union of all above ","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:53:13.391976Z","iopub.execute_input":"2023-07-10T11:53:13.392812Z","iopub.status.idle":"2023-07-10T11:53:13.39934Z","shell.execute_reply.started":"2023-07-10T11:53:13.39277Z","shell.execute_reply":"2023-07-10T11:53:13.398102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity - sequence duplicates ","metadata":{}},{"cell_type":"code","source":"%%time\nprint( pd.Series(list_seqs).value_counts().head(10) )\ndf_eda[df_eda['sequence'] == 'MTQSNPNEQNVELNRTSLYWGLLLIFVLAVLFSNYFFN' ].head(10)","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:53:18.030714Z","iopub.execute_input":"2023-07-10T11:53:18.031128Z","iopub.status.idle":"2023-07-10T11:53:18.281819Z","shell.execute_reply.started":"2023-07-10T11:53:18.031088Z","shell.execute_reply":"2023-07-10T11:53:18.281027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom itertools import combinations\nimport time \n\n# Settings:\nstr_similarity_id = 'seq duplicate'\nif str_similarity_id in list_allowed_similarity_relations:  \n    series_selected = df_eda['sequence'] \n\n    verbose = 100\n    if verbose >= 1:\n        print('Start processing similarity:', str_similarity_id, '. Number of groups: ', (series_selected.value_counts() > 1 ).sum()  )\n        print()\n\n    # Local results storage:\n    list_edges_for_graph_loc = []  # Local results will be here - pairs of similar samples (edges of graph)\n\n    # Start processing: \n    v_duplication = series_selected.value_counts()\n    print('Top dubplicated:' )\n    print(v_duplication.head(8) )\n    print()\n    list_with_duplicate = list(v_duplication[ v_duplication > 1 ].index) # one object for each similarity class \n    m = series_selected.isin(list_with_duplicate)\n    print('Count groups of duplicates/similarity:', (v_duplication>1).sum() )\n    print('Total samples having dubplicate/similar:', m.sum())\n    print()\n    print('Start process. Create list of pairs - all pairs from each group will be listed.')\n    t0 = time.time()\n    for i0, uv in enumerate(list_with_duplicate): # loop of over similarity groups\n\n        # Technical checks, needed for some similarities , but not the others:  \n        if len(str(uv)) <= 2 : continue # too short genes names \n        if str(uv).isnumeric() : continue # numbers are not names \n        flag_stop = 0\n        for stop_kw in ['uncharacterized', 'putative', 'probable']:\n            if stop_kw in str(uv).lower(): flag_stop = 1\n        if flag_stop > 0: continue \n\n        # \n        m2 = series_selected[m] == uv # Select samples within current similarity group\n        s1 = set(series_selected[m][m2].index) # make list all of them\n        # Main part - append to results list  -  each pair of samples within current similarity group\n        list_edges_for_graph_loc += list( combinations(s1,2) )\n\n\n        # Optional output to control process/show progress \n        if verbose >= 10:\n            if i0%1000 == 0:\n                print('Progress:',i0,' out of ', len(list_with_duplicate) , '. Current pairs count:', len(list_edges_for_graph_loc), len(set(list_edges_for_graph_loc)), '%.1f secs passed'%(time.time() - t0) )\n        if verbose >= 100:\n            if i0<3:\n                print('Detailed information.', i0, m2.sum(), uv)\n                print(len(s1), 'Tail10 indices of current group:', list(s1)[-10:])\n                print('Tail20 pairs:', list_edges_for_graph_loc[-20:])\n                print()\n\n    print('Finished. %.1f secs passed'%(time.time() - t0))\n\n\n    dict_list_edges_for_graph[str_similarity_id] = list_edges_for_graph_loc.copy() # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \n    list_edges_for_graph += list_edges_for_graph_loc # here is union of all above \n\n    print()\n    print('Results info:')\n    print('Similarity:',str_similarity_id ,'. Number of pairs obtained:', len(list_edges_for_graph_loc), '. Number of non-identical pairs:',  len(set(list_edges_for_graph_loc)))\n    print('Total:', len(list_edges_for_graph), len(set(list_edges_for_graph)))\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:13.146144Z","iopub.execute_input":"2023-07-10T11:56:13.147545Z","iopub.status.idle":"2023-07-10T11:56:20.240695Z","shell.execute_reply.started":"2023-07-10T11:56:13.147494Z","shell.execute_reply":"2023-07-10T11:56:20.239112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create graph and сluster it ","metadata":{}},{"cell_type":"code","source":"%%time \nimport time\nimport igraph\nt0 = time.time()\ng = igraph.Graph()\ng.add_vertices(len(list_ids))\ng.es['label'] = list_ids # list(matches['qseqid'])  # \ng.add_edges(set(list_edges_for_graph_loc) )\nprint('%.1f graph created'%(time.time() - t0 ) )\nprint() \n\n\n\nt0 = time.time()\nlist_cc = g.connected_components(mode='WEAK')\nl3 = [len(t) for t in list_cc]\nprint('%.1f got connected_components'%(time.time() - t0 ) )\nprint('n_components:', len(list_cc),'. Max size:', np.max(l3), '. Min size:' , np.min(l3), '. Count of size one: ', \n      (np.array(l3)==1).sum(), '. Sizes (top20):',  np.sort(l3)[::-1][:20],  )# ,ascending = False)\nprint('Value counts of component sizes:')\nprint(dict( pd.Series(l3).value_counts().head(30) )  )\n\n########################################################################\n# Cluster by Louvain algorithm \n# https://igraph.org/python/doc/igraph.Graph-class.html#community_multilevel\n# See examples e.g. here: https://www.kaggle.com/code/alexandervc/igraph-cheatsheet?scriptVersionId=59716528&cellId=42\n########################################################################\nprint()\nt0 = time.time()\nlouvain_partition = g.community_multilevel()# weights=graph.es['weight'], return_levels=False)\nprint('%.1f got louvain_partition'%(time.time() - t0 ) )\nprint('n_clusters:', pd.Series(louvain_partition.membership).nunique() )\nprint('Sizes of clusters:')\nprint(list( pd.Series(louvain_partition.membership).value_counts().head(30) )  )\nprint('Value counts of component sizes:')\nprint(dict( pd.Series(louvain_partition.membership).value_counts().value_counts().head(30) )  )\nprint(); print();\n\nfor k_ in ['connected components', 'Louvain clusters']:\n    N1=300\n    if k_ == 'connected components':\n        sr_loc1 = np.sort(l3)[::-1][:N1]\n        sr_loc2 = pd.Series(l3)\n    else:\n        sr_loc1 = pd.Series(louvain_partition.membership).value_counts().head(N1).values\n        sr_loc2 = pd.Series(louvain_partition.membership).value_counts()\n        \n    plt.figure(figsize = (20,8))\n    plt.suptitle(str_similarity_id,fontsize = 18)\n    plt.subplot(121)\n    \n    plt.plot(sr_loc1)\n    plt.title('Top'+str(N1)+' '+ k_ + ' sizes',fontsize = 16)\n    plt.grid()\n    plt.subplot(122)\n    N2 = 300\n    d_ = sr_loc2.value_counts().head(N2).sort_index()\n    plt.plot( d_.index , np.log10(1+d_.values) )\n    plt.title('Top'+str(N2) + ' log10 count ' + k_ + ' of given size',fontsize = 16)\n    plt.xlabel('Size', fontsize = 16)\n    plt.ylabel('Count Log10')\n    plt.grid()\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:26.377419Z","iopub.execute_input":"2023-07-10T11:56:26.378008Z","iopub.status.idle":"2023-07-10T11:56:28.561364Z","shell.execute_reply.started":"2023-07-10T11:56:26.377977Z","shell.execute_reply":"2023-07-10T11:56:28.560294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Graph - same code as above wrapped into function","metadata":{}},{"cell_type":"code","source":"%%time \nimport time\nimport igraph\n\ndef get_groups_by_graph_clustering(list_edges_for_graph_loc, str_similarity_id, plot_mode = 'default' ):\n    if len(list_edges_for_graph_loc) == 0: return np.array([]) # exceptional case \n    t0 = time.time()\n    g = igraph.Graph()\n    g.add_vertices(len(list_ids))\n    g.es['label'] = list_ids # list(matches['qseqid'])  # \n    g.add_edges(set(list_edges_for_graph_loc) )\n    print('%.1f graph created'%(time.time() - t0 ) )\n    print() \n\n\n\n    t0 = time.time()\n    list_cc = g.connected_components(mode='WEAK')\n    l3 = [len(t) for t in list_cc]\n    print('%.1f got connected_components'%(time.time() - t0 ) )\n    print('n_components:', len(list_cc),'. Max size:', np.max(l3), '. Min size:' , np.min(l3), '. Count of size one: ', \n          (np.array(l3)==1).sum(), '. Sizes (top20):',  np.sort(l3)[::-1][:20],  )# ,ascending = False)\n    print('Value counts of component sizes:')\n    print(dict( pd.Series(l3).value_counts().head(30) )  )\n\n    ########################################################################\n    # Cluster by Louvain algorithm \n    # https://igraph.org/python/doc/igraph.Graph-class.html#community_multilevel\n    # See examples e.g. here: https://www.kaggle.com/code/alexandervc/igraph-cheatsheet?scriptVersionId=59716528&cellId=42\n    ########################################################################\n    print()\n    t0 = time.time()\n    louvain_partition = g.community_multilevel()# weights=graph.es['weight'], return_levels=False)\n    print('%.1f got louvain_partition'%(time.time() - t0 ) )\n    print('n_clusters:', pd.Series(louvain_partition.membership).nunique() )\n    print('Sizes of clusters:')\n    print(list( pd.Series(louvain_partition.membership).value_counts().head(30) )  )\n    print('Value counts of component sizes:')\n    print(dict( pd.Series(louvain_partition.membership).value_counts().value_counts().head(30) )  )\n    print(); print();\n\n    if plot_mode == 'default':\n        for k_ in ['connected components', 'Louvain clusters']:\n            N1=300\n            if k_ == 'connected components':\n                sr_loc1 = np.sort(l3)[::-1][:N1]\n                sr_loc2 = pd.Series(l3)\n            else:\n                sr_loc1 = pd.Series(louvain_partition.membership).value_counts().head(N1).values\n                sr_loc2 = pd.Series(louvain_partition.membership).value_counts()\n\n            plt.figure(figsize = (20,8))\n            plt.suptitle(str_similarity_id,fontsize = 18)\n            plt.subplot(121)\n\n            plt.plot(sr_loc1)\n            plt.title('Top'+str(N1)+' '+ k_ + ' sizes',fontsize = 16)\n            plt.grid()\n            plt.subplot(122)\n            N2 = 300\n            d_ = sr_loc2.value_counts().head(N2).sort_index()\n            plt.plot( d_.index , np.log10(1+d_.values) )\n            plt.title('Top'+str(N2) + ' log10 count ' + k_ + ' of given size',fontsize = 16)\n            plt.xlabel('Size', fontsize = 16)\n            plt.ylabel('Count Log10')\n            plt.grid()\n            plt.show()\n\n    return np.asarray(louvain_partition.membership ) \n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:28.562998Z","iopub.execute_input":"2023-07-10T11:56:28.563346Z","iopub.status.idle":"2023-07-10T11:56:28.580403Z","shell.execute_reply.started":"2023-07-10T11:56:28.563318Z","shell.execute_reply":"2023-07-10T11:56:28.579162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"get_groups_by_graph_clustering(list_edges_for_graph, str_similarity_id,  plot_mode = 'default' )","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:29.29018Z","iopub.execute_input":"2023-07-10T11:56:29.290822Z","iopub.status.idle":"2023-07-10T11:56:31.128958Z","shell.execute_reply.started":"2023-07-10T11:56:29.290792Z","shell.execute_reply":"2023-07-10T11:56:31.127891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity - special proteins groups (e.g.  'Tachykinin-related peptide'  )\n\nPrevious analysis showed that many short 'Tachykinin-related peptide' were not aligned into one group by DIAMOND (even in ultra-sensitive mode). We need to correct it manually","metadata":{"execution":{"iopub.status.busy":"2023-06-12T07:23:45.762339Z","iopub.execute_input":"2023-06-12T07:23:45.762712Z","iopub.status.idle":"2023-06-12T07:23:45.779408Z","shell.execute_reply.started":"2023-06-12T07:23:45.762683Z","shell.execute_reply":"2023-06-12T07:23:45.778463Z"}}},{"cell_type":"code","source":"str_kw = 'Tachykinin-related peptide'\nmsk1 = df_eda['description'].apply( lambda x: str_kw in x )\nprint(msk1.sum(), 'count proteins containing', str_kw)\ndf_eda[msk1].head(10)","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:31.130826Z","iopub.execute_input":"2023-07-10T11:56:31.131832Z","iopub.status.idle":"2023-07-10T11:56:31.204568Z","shell.execute_reply.started":"2023-07-10T11:56:31.131787Z","shell.execute_reply":"2023-07-10T11:56:31.203525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nstr_similarity_id = 'special groups'\nif str_similarity_id in list_allowed_similarity_relations:  \n\n    str_kw = 'Tachykinin-related peptide'\n    msk1 = df_eda['description'].apply( lambda x: str_kw in x )\n    print(msk1.sum(), 'count proteins containing', str_kw)\n\n    list_edges_for_graph_loc = []\n\n    s1 = set(df_eda[msk1].index)\n    print(s1)\n    list_edges_for_graph_loc += list( combinations(s1,2) ) \n\n    dict_list_edges_for_graph[str_similarity_id] = list_edges_for_graph_loc.copy() # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \n    list_edges_for_graph += list_edges_for_graph_loc # here is union of all above \n\n    print()\n    print('Results info:')\n    print('Similarity:',str_similarity_id ,'. Number of pairs obtained:', len(list_edges_for_graph_loc), '. Number of non-identical pairs:',  len(set(list_edges_for_graph_loc)))\n    print('Total:', len(list_edges_for_graph), len(set(list_edges_for_graph)))\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:31.576414Z","iopub.execute_input":"2023-07-10T11:56:31.57738Z","iopub.status.idle":"2023-07-10T11:56:31.632893Z","shell.execute_reply.started":"2023-07-10T11:56:31.57734Z","shell.execute_reply":"2023-07-10T11:56:31.632121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nget_groups_by_graph_clustering(list_edges_for_graph_loc, str_similarity_id,  plot_mode = 'default' )\n\nget_groups_by_graph_clustering(list_edges_for_graph, str(list(dict_list_edges_for_graph.keys())),  plot_mode = 'default' )\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:32.896529Z","iopub.execute_input":"2023-07-10T11:56:32.897269Z","iopub.status.idle":"2023-07-10T11:56:36.629031Z","shell.execute_reply.started":"2023-07-10T11:56:32.897213Z","shell.execute_reply":"2023-07-10T11:56:36.627966Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity - same gene name ","metadata":{}},{"cell_type":"code","source":"print(len(list_edges_for_graph), len(set(list_edges_for_graph)))\ndf_eda[df_eda['gene name lower'] == 'fur'].head(10) \n# \"fur\" seems to be really gene: Transcriptional regulation by Ferric Uptake Regulator (Fur) in pathogenic bacteria \n# https://www.frontiersin.org/articles/10.3389/fcimb.2013.00059/full\n\n# Another strange named example:\n# # Trypanosoma brucei Tb927.2.6100 Is an Essential Protein Associated with Kinetoplast DNA\n# # https://www.frontiersin.org/articles/10.3389/fcimb.2013.00059/full\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:36.630729Z","iopub.execute_input":"2023-07-10T11:56:36.631051Z","iopub.status.idle":"2023-07-10T11:56:36.668798Z","shell.execute_reply.started":"2023-07-10T11:56:36.631022Z","shell.execute_reply":"2023-07-10T11:56:36.668019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom itertools import combinations\nimport time \n\n# Settings:\nstr_similarity_id = 'gene name'\nif str_similarity_id in list_allowed_similarity_relations:  \n    series_selected = df_eda['gene name lower'] \n\n    verbose = 10\n    if verbose >= 1:\n        print('Start processing similarity:', str_similarity_id, '. Number of groups: ', (series_selected.value_counts() > 1 ).sum()  )\n        print()\n\n    # Local results storage:\n    list_edges_for_graph_loc = []  # Local results will be here - pairs of similar samples (edges of graph)\n\n    # Start processing: \n    v_duplication = series_selected.value_counts()\n    print('Top dubplicated:' )\n    print(v_duplication.head(8) )\n    print()\n    list_with_duplicate = list(v_duplication[ v_duplication > 1 ].index) # one object for each similarity class \n    m = series_selected.isin(list_with_duplicate)\n    print('Count groups of duplicates/similarity:', (v_duplication>1).sum() )\n    print('Total samples having dubplicate/similar:', m.sum())\n    print()\n    print('Start process. Create list of pairs - all pairs from each group will be listed.')\n    t0 = time.time()\n    for i0, uv in enumerate(list_with_duplicate): # loop of over similarity groups\n\n        # Technical checks, needed for some similarities , but not the others:  \n        if len(str(uv)) <= 2 : continue # too short genes names \n        if str(uv).isnumeric() : continue # numbers are not names \n        flag_stop = 0\n        for stop_kw in ['uncharacterized', 'putative', 'probable']:\n            if stop_kw in str(uv).lower(): flag_stop = 1\n        if flag_stop > 0: continue \n\n        # \n        m2 = series_selected[m] == uv # Select samples within current similarity group\n        s1 = set(series_selected[m][m2].index) # make list all of them\n        # Main part - append to results list  -  each pair of samples within current similarity group\n        list_edges_for_graph_loc += list( combinations(s1,2) )\n\n        # Optional output to control process/show progress \n        if verbose >= 10:\n            if i0%1000 == 0:\n                print('Progress:',i0,' out of ', len(list_with_duplicate) , '. Current pairs count:', len(list_edges_for_graph_loc), len(set(list_edges_for_graph_loc)), '%.1f secs passed'%(time.time() - t0) )\n        if verbose >= 100:\n            if i0<3:\n                print('Detailed information.', i0, m2.sum(), uv)\n                print(len(s1), 'Tail10 indices of current group:', list(s1)[-10:])\n                print('Tail20 pairs:', list_edges_for_graph_loc[-20:])\n                print()\n\n    print('Finished. %.1f secs passed'%(time.time() - t0))\n\n\n    dict_list_edges_for_graph[str_similarity_id] = list_edges_for_graph_loc.copy() # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \n    list_edges_for_graph += list_edges_for_graph_loc # here is union of all above \n\n    print()\n    print('Results info:')\n    print('Similarity:',str_similarity_id ,'. Number of pairs obtained:', len(list_edges_for_graph_loc), '. Number of non-identical pairs:',  len(set(list_edges_for_graph_loc)))\n    print('Total:', len(list_edges_for_graph), len(set(list_edges_for_graph)))\n","metadata":{"execution":{"iopub.status.busy":"2023-07-10T11:56:36.670314Z","iopub.execute_input":"2023-07-10T11:56:36.671029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nget_groups_by_graph_clustering(list_edges_for_graph_loc, str_similarity_id,  plot_mode = 'default' )\n\nget_groups_by_graph_clustering(list_edges_for_graph, str(list(dict_list_edges_for_graph.keys())),  plot_mode = 'default' )\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity -  gene name - same after deleting digits at the end \n\nMany genes are named like CDK1,CDK2 ... even naming shows their similarity (though , it sometimes can be misleading, but for us it is not a problem to have bigger groups, unless theirs sizes exceed limits like 10_000, so for the moment we include that option.). ","metadata":{}},{"cell_type":"code","source":"import re\ndef delete_digits_at_end(string):\n    string = str(string)\n    pattern = r'\\d+$'  # Matches one or more digits at the end of the string\n    modified_string = re.sub(pattern, '', string)\n    return modified_string\n\nv_gene_no_last_digits = df_eda['gene name lower'].apply( delete_digits_at_end )\nseries_selected = v_gene_no_last_digits \nprint(series_selected.value_counts().head(5) )\ndf_eda[series_selected == 'tb11.01.'].head(10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom itertools import combinations\nimport time \n\n# Settings:\nstr_similarity_id = 'gene name end digits cutted'\nif str_similarity_id in list_allowed_similarity_relations:  \n    series_selected = v_gene_no_last_digits\n\n    verbose = 10\n    if verbose >= 1:\n        print('Start processing similarity:', str_similarity_id, '. Number of groups: ', (series_selected.value_counts() > 1 ).sum()  )\n        print()\n\n    # Local results storage:\n    list_edges_for_graph_loc = []  # Local results will be here - pairs of similar samples (edges of graph)\n\n    # Start processing: \n    v_duplication = series_selected.value_counts()\n    print('Top dubplicated:' )\n    print(v_duplication.head(8) )\n    print()\n    list_with_duplicate = list(v_duplication[ v_duplication > 1 ].index) # one object for each similarity class \n    m = series_selected.isin(list_with_duplicate)\n    print('Count groups of duplicates/similarity:', (v_duplication>1).sum() )\n    print('Total samples having dubplicate/similar:', m.sum())\n    print()\n    print('Start process. Create list of pairs -  all pairs from each group will be listed.')\n    t0 = time.time()\n    for i0, uv in enumerate(list_with_duplicate): # loop of over similarity groups\n\n        # Technical checks, needed for some similarities , but not the others:  \n        if len(str(uv)) <= 2 : continue # too short genes names \n        if str(uv).isnumeric() : continue # numbers are not names \n        flag_stop = 0\n        for stop_kw in ['uncharacterized', 'putative', 'probable']:\n            if stop_kw in str(uv).lower(): flag_stop = 1\n        if flag_stop > 0: continue \n\n        # \n        m2 = series_selected[m] == uv # Select samples within current similarity group\n        s1 = set(series_selected[m][m2].index) # make list all of them\n        # Main part - append to results list  -  each pair of samples within current similarity group\n        list_edges_for_graph_loc += list( combinations(s1,2) )\n\n        # Optional output to control process/show progress \n        if verbose >= 10:\n            if i0%1000 == 0:\n                print('Progress:',i0,' out of ', len(list_with_duplicate) , '. Current pairs count:', len(list_edges_for_graph_loc), len(set(list_edges_for_graph_loc)), '%.1f secs passed'%(time.time() - t0) )\n        if verbose >= 100:\n            if i0<3:\n                print('Detailed information.', i0, m2.sum(), uv)\n                print(len(s1), 'Tail10 indices of current group:', list(s1)[-10:])\n                print('Tail20 pairs:', list_edges_for_graph_loc[-20:])\n                print()\n\n    print('Finished. %.1f secs passed'%(time.time() - t0))\n\n\n    dict_list_edges_for_graph[str_similarity_id] = list_edges_for_graph_loc.copy() # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \n    list_edges_for_graph += list_edges_for_graph_loc # here is union of all above \n\n    print()\n    print('Results info:')\n    print('Similarity:',str_similarity_id ,'. Number of pairs obtained:', len(list_edges_for_graph_loc), '. Number of non-identical pairs:',  len(set(list_edges_for_graph_loc)))\n    print('Total:', len(list_edges_for_graph), len(set(list_edges_for_graph)))\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nget_groups_by_graph_clustering(list_edges_for_graph_loc, str_similarity_id,  plot_mode = 'default' )\n\nget_groups_by_graph_clustering(list_edges_for_graph, str(list(dict_list_edges_for_graph.keys())),  plot_mode = 'default' )\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity - same description \n\nGenes have many aliases and so it is not always easy to match - becasue in one sample it can be called like AAA , in the other like BBB - there are databases with aliases for genes names, but will not use that way now.\n\nSo we will additionally check that if descriptions of genes is the same - they will assigned to the same group.  With some exceptions like \"Uncharacterized protein\"\n","metadata":{"execution":{"iopub.status.busy":"2023-06-12T09:43:18.947292Z","iopub.execute_input":"2023-06-12T09:43:18.947686Z","iopub.status.idle":"2023-06-12T09:43:18.95516Z","shell.execute_reply.started":"2023-06-12T09:43:18.947654Z","shell.execute_reply":"2023-06-12T09:43:18.954059Z"}}},{"cell_type":"code","source":"print( df_eda['Description Cleaned'].value_counts().head(30) )\ndf_eda[df_eda['Description Cleaned'] == 'Conserved protein ' ].head(15)\ndf_eda[df_eda['Description Cleaned'] == 'Tyrosine-protein kinase ' ].head(10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom itertools import combinations\nimport time \n\n# Settings:\nstr_similarity_id = 'description'\nif str_similarity_id in list_allowed_similarity_relations:  \n    series_selected = df_eda['Description Cleaned']\n\n    verbose = 10\n    if verbose >= 1:\n        print('Start processing similarity:', str_similarity_id, '. Number of groups: ', (series_selected.value_counts() > 1 ).sum()  )\n        print()\n\n    # Local results storage:\n    list_edges_for_graph_loc = []  # Local results will be here - pairs of similar samples (edges of graph)\n\n    # Start processing: \n    v_duplication = series_selected.value_counts()\n    print('Top dubplicated:' )\n    print(v_duplication.head(8) )\n    print()\n    list_with_duplicate = list(v_duplication[ v_duplication > 1 ].index) # one object for each similarity class \n    m = series_selected.isin(list_with_duplicate)\n    print('Count groups of duplicates/similarity:', (v_duplication>1).sum() )\n    print('Total samples having dubplicate/similar:', m.sum())\n    print()\n    print('Start process. Create list of pairs -  all pairs from each group will be listed.')\n    t0 = time.time()\n    for i0, uv in enumerate(list_with_duplicate): # loop of over similarity groups\n\n        # Technical checks, needed for some similarities , but not the others:  \n        if len(str(uv)) <= 2 : continue # too short genes names \n        if str(uv).isnumeric() : continue # numbers are not names \n        flag_stop = 0\n        for stop_kw in ['uncharacterized', 'putative', 'probable']:\n            if stop_kw in str(uv).lower(): flag_stop = 1\n        if flag_stop > 0: continue \n\n        # \n        m2 = series_selected[m] == uv # Select samples within current similarity group\n        s1 = set(series_selected[m][m2].index) # make list all of them\n        # Main part - append to results list  -  each pair of samples within current similarity group\n        list_edges_for_graph_loc += list( combinations(s1,2) )\n\n        # Optional output to control process/show progress \n        if verbose >= 10:\n            if i0%1000 == 0:\n                print('Progress:',i0,' out of ', len(list_with_duplicate) , '. Current pairs count:', len(list_edges_for_graph_loc), len(set(list_edges_for_graph_loc)), '%.1f secs passed'%(time.time() - t0) )\n        if verbose >= 100:\n            if i0<3:\n                print('Detailed information.', i0, m2.sum(), uv)\n                print(len(s1), 'Tail10 indices of current group:', list(s1)[-10:])\n                print('Tail20 pairs:', list_edges_for_graph_loc[-20:])\n                print()\n\n    print('Finished. %.1f secs passed'%(time.time() - t0))\n\n\n    dict_list_edges_for_graph[str_similarity_id] = list_edges_for_graph_loc.copy() # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \n    list_edges_for_graph += list_edges_for_graph_loc # here is union of all above \n\n    print()\n    print('Results info:')\n    print('Similarity:',str_similarity_id ,'. Number of pairs obtained:', len(list_edges_for_graph_loc), '. Number of non-identical pairs:',  len(set(list_edges_for_graph_loc)))\n    print('Total:', len(list_edges_for_graph), len(set(list_edges_for_graph)))\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nget_groups_by_graph_clustering(list_edges_for_graph_loc, str_similarity_id,  plot_mode = 'default' )\n\nget_groups_by_graph_clustering(list_edges_for_graph, str(list(dict_list_edges_for_graph.keys())),  plot_mode = 'default' )\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity by DIAMOND local alignment\n\nLoad data with alignment \"train x train\" calculated by DIAMOND package.\n\n\nThanks to Liza Geraseva: https://www.kaggle.com/code/geraseva/diamond\n\n\nDIAMOND as any blast-like algorithm is sub-optimal, and thus  not that much sensitive as Smith-Waterman local alignment. \nBut for the case of obtaining groups we do need very strong sensitivity, otherwise we can get too big groups - like \"everything is similar to everything\". Also running Smith-Waterman for train-x-train even with the fast Parasail package would take weeks on Kaggle.\nSo DIAMOND seems to be the reasonable choice. \n","metadata":{}},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/cafa5-features-etc/matches_train_diamond_ultra_sensitive/matches_train_diamond_ultra_sensitive.tsv'\nmatches=pd.read_csv(fn, sep='\\t', header=None,  # 'qlen', 'slen',\n                    names=['qseqid', 'sseqid',  'pident', 'length', 'mismatch', \n                           'gapopen', 'qstart', 'qend', 'sstart','send', 'evalue', 'bitscore'])\nmatches['qseqid']=matches['qseqid'].apply(lambda x: x.split('\\\\t')[0])\nprint(matches.shape)\nprint(matches.evalue.max())\nmatches.head(5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nstr_similarity_id = 'diamond local alignment'\nif str_similarity_id in list_allowed_similarity_relations:  \n    msk = matches['evalue'] < evalue_threshold \n    print(msk.sum() )\n    lt1 = matches[msk]['qseqid'].map(map_id2num)#  [ dict_id2num[matches[]]]\n    lt2 = matches[msk]['sseqid'].map(map_id2num)#\n    list_edges_for_graph_loc = zip(lt1,lt2)\n    \n    list_edges_for_graph_loc = list(list_edges_for_graph_loc)\n    \n    dict_list_edges_for_graph[str_similarity_id] = list_edges_for_graph_loc.copy() # For each similarity type here will be a list  of  pairs of similar samples by that similarity relation  \n    list_edges_for_graph += list_edges_for_graph_loc # here is union of all above \n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nget_groups_by_graph_clustering(list_edges_for_graph_loc, str_similarity_id,  plot_mode = 'default' )\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Final groups","metadata":{}},{"cell_type":"code","source":"%%time\nvec_groups_final = get_groups_by_graph_clustering(list_edges_for_graph, str(list(dict_list_edges_for_graph.keys())),  plot_mode = 'default' )\nprint(type(vec_groups_final), vec_groups_final.shape )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_groups = pd.DataFrame( index = list_ids, data = vec_groups_final, columns = ['group'] )\nn_ = len( list(dict_list_edges_for_graph.keys()) )\nif 'diamond local alignment' in dict_list_edges_for_graph.keys():\n    str_evalue = '_evalue' + str( evalue_threshold).replace('.','d')\nelse:\n    str_evalue = ''\nfs = \"df_groups_Nsims\"+str(n_)+ str_evalue +\".csv\" \nprint(fs)\ndf_groups.to_csv(fs)\nprint(df_groups.iloc[:,0].nunique())\ndf_groups.head(10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Folds ","metadata":{}},{"cell_type":"code","source":"%%time\ndf_folds = pd.DataFrame(index = list_ids) \nfrom sklearn.model_selection import GroupKFold\ngroups = df_groups.iloc[:,0]\nX = np.arange(len(list_ids)).reshape(-1,1)\ny = np.arange(len(list_ids))\ncc = 0\nfor n_splits in range(2,21):\n    group_kfold = GroupKFold(n_splits=n_splits)\n    col = 'cv_folds_'+str(n_splits)\n    df_folds[col] = 0\n    for i, (train_index, test_index) in enumerate(group_kfold.split(X, y, groups)):\n        df_folds[col].iloc[test_index] = i\n        if cc <= 5:\n            print(n_splits, i, len(test_index), test_index[:15])\n        cc+=1\n","metadata":{"execution":{"iopub.status.busy":"2023-06-16T15:51:27.246126Z","iopub.execute_input":"2023-06-16T15:51:27.246514Z","iopub.status.idle":"2023-06-16T15:51:28.697776Z","shell.execute_reply.started":"2023-06-16T15:51:27.246485Z","shell.execute_reply":"2023-06-16T15:51:28.696698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# save df_folds \n\nn_ = len( list(dict_list_edges_for_graph.keys()) )\nif 'diamond local alignment' in dict_list_edges_for_graph.keys():\n    str_evalue = '_evalue' + str( evalue_threshold).replace('.','d')\nelse:\n    str_evalue = ''\nfs = \"df_folds_Nsims\"+str(n_)+ str_evalue +\".csv\" \nprint(fs)\n\ndf_folds.to_csv(fs)\nprint(df_folds.iloc[:,0].nunique())\ndf_folds.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-06-16T15:51:28.698917Z","iopub.execute_input":"2023-06-16T15:51:28.69959Z","iopub.status.idle":"2023-06-16T15:51:29.506466Z","shell.execute_reply.started":"2023-06-16T15:51:28.699561Z","shell.execute_reply":"2023-06-16T15:51:29.505301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time        \nfor col in df_folds.columns:\n    print(col, df_folds[col].nunique(), dict(df_folds[col].value_counts()) )","metadata":{"execution":{"iopub.status.busy":"2023-06-16T15:51:29.50845Z","iopub.execute_input":"2023-06-16T15:51:29.508893Z","iopub.status.idle":"2023-06-16T15:51:29.571037Z","shell.execute_reply.started":"2023-06-16T15:51:29.508855Z","shell.execute_reply":"2023-06-16T15:51:29.569919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Visual inspection that folds are randomly scattered along the index \nplt.figure(figsize = (20,5))\nfor col in df_folds.columns[:5]:\n    v = df_folds[col]\n    plt.plot(v.values[:1000])\nplt.grid()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-16T15:51:29.572323Z","iopub.execute_input":"2023-06-16T15:51:29.572623Z","iopub.status.idle":"2023-06-16T15:51:30.08977Z","shell.execute_reply.started":"2023-06-16T15:51:29.572598Z","shell.execute_reply":"2023-06-16T15:51:30.088293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# one more check that folds are not correlated with index - thus randomly shuffled\nd_ = df_folds.copy()\nd_['index'] = range(len(d_))\nd_.corr().round(2)","metadata":{"execution":{"iopub.status.busy":"2023-06-16T15:51:30.091319Z","iopub.execute_input":"2023-06-16T15:51:30.091645Z","iopub.status.idle":"2023-06-16T15:51:30.376926Z","shell.execute_reply.started":"2023-06-16T15:51:30.091617Z","shell.execute_reply":"2023-06-16T15:51:30.37595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Models examples ","metadata":{}},{"cell_type":"code","source":"import gc\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-16T17:09:01.942888Z","iopub.execute_input":"2023-06-16T17:09:01.943303Z","iopub.status.idle":"2023-06-16T17:09:03.453533Z","shell.execute_reply.started":"2023-06-16T17:09:01.943272Z","shell.execute_reply":"2023-06-16T17:09:03.452682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n######################################################################################3\n###############  Load embeddings features \n######################################################################################3\nfn = '/kaggle/input/t5embeds/train_embeds.npy'\nprint(fn)\nX = np.load(fn)\nprint(X.shape)\nprint(X[:2,:3])\n\nif 1:\n    fn = '/kaggle/input/t5embeds/test_embeds.npy'\n    print(fn)\n    X_submit = np.load(fn)\n    print(X_submit.shape)\n    print(X_submit[:2,:3])\n\n######################################################################################3\n###############  Load targets (processed)\n######################################################################################3\nfn = '/kaggle/input/cafa5-features-etc/Y_1499.npy'# Y_1499_labels.npy'\nY = np.load(fn)\nY_save = Y.copy()\nprint( Y.shape )\nY[:2,:3]\n","metadata":{"execution":{"iopub.status.busy":"2023-06-16T17:00:04.745656Z","iopub.execute_input":"2023-06-16T17:00:04.746127Z","iopub.status.idle":"2023-06-16T17:00:36.52402Z","shell.execute_reply.started":"2023-06-16T17:00:04.74609Z","shell.execute_reply":"2023-06-16T17:00:36.522355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_df_models_stat = {}\ndict_Y_pred_oof = {} # double dict: access: d[str_model_id][CV_id]\nY = Y_save[:,0:500] # \n","metadata":{"execution":{"iopub.status.busy":"2023-06-16T17:01:24.595613Z","iopub.execute_input":"2023-06-16T17:01:24.596025Z","iopub.status.idle":"2023-06-16T17:01:24.754608Z","shell.execute_reply.started":"2023-06-16T17:01:24.595995Z","shell.execute_reply":"2023-06-16T17:01:24.753091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.metrics import f1_score\nfrom sklearn.linear_model import Ridge\nimport gc\ngc.collect()\n\nmodel = Ridge(alpha=1.0)\nstr_model_id = 'Ridge1'\ndict_Y_pred_oof[str_model_id] = {} \n\n# col4folds = 'cv_folds_3'\n\nfor col4folds in  ['cv_folds_5','cv_folds_10', 'cv_folds_15']: \n    df_scores = pd.DataFrame()\n    vf = df_folds[col4folds]\n    print(str_model_id, col4folds, 'Y.shape:', Y.shape )\n    t00 = time.time()\n    Y_pred_oof = np.zeros( Y.shape )\n    Y_pred_submit = np.zeros( (X_submit.shape[0], Y.shape[1]) )\n    for i_fold,uv in enumerate(vf.unique()):\n        m = vf == uv\n        IX_train = np.where( ~m)[0]\n        IX_test  = np.where(  m)[0]\n        print(i_fold, 'Fold:',uv, 'n_samples in test:', m.sum(), 'Indexes and their lens:', len(IX_train), len(IX_test), IX_train[:10], IX_test[:10] )\n\n        t0 = time.time()\n        model.fit(X[IX_train,:], Y[IX_train,:] )\n        print('%.1f secs passed on fit'%(time.time() - t0 ))\n\n        t0 = time.time()\n        Y_pred_submit = (Y_pred_submit * i_fold  + model.predict(X_submit) )/ (i_fold + 1)\n        print('%.1f secs passed on submit predict'%(time.time() - t0 ))\n\n        ################################ predicts and scores calculations  #################################\n\n        IX_loc = IX_test\n        Y_pred = model.predict(X[IX_loc,:])\n        t0 = time.time()\n        l_ = []\n        for i_ in range(Y.shape[1]):\n            s = roc_auc_score(Y[IX_loc,i_].ravel(),  Y_pred[:,i_].ravel() ); l_.append(s)\n        print('%.1f secs passed on test roc-auc'%(time.time() - t0 ), 'mean rocauc: %.3f'%(np.mean(l_)))\n        df_scores[str_model_id + ' Test'+str(uv)] = l_\n        Y_pred_oof[IX_loc,:] = Y_pred\n\n        IX_loc = IX_train \n        Y_pred = model.predict(X[IX_loc,:])\n        t0 = time.time()\n        l_ = []\n        for i_ in range(Y.shape[1]):\n            s = roc_auc_score(Y[IX_loc,i_].ravel(),  Y_pred[:,i_].ravel() ); l_.append(s)\n        print('%.1f secs passed on train roc-auc'%(time.time() - t0 ), 'mean rocauc: %.3f'%(np.mean(l_)))\n        df_scores[str_model_id + ' Train'+str(uv)] = l_\n        \n        gc.collect()\n        \n    np.save('Y_pred_submit|'+str_model_id+'|' + col4folds , Y_pred_submit )         \n    dict_Y_pred_oof[str_model_id][col4folds] = Y_pred_oof\n    dict_df_models_stat[str_model_id + ' '+col4folds] =   df_scores.copy()\n\n    print('%.1f secs passed on model'%(time.time() - t00 ), str_model_id , 'Folds:', col4folds)\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        print(kw,'mean rocauc score:',df_scores[l_col].mean().mean().round(3), str_model_id + ' '+col4folds)\n    print()\n    \n    ######################################### Plots ######################################### \n    score_name = 'rocauc'\n    plt.figure(figsize  = (20, 6) )    \n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        v = df_scores[l_col].mean(axis = 1)\n        plt.plot(v, label = kw)\n        print(kw,'mean rocauc score:',v.mean(), str_model_id + ' '+col4folds)\n    plt.grid()\n    plt.legend()\n    plt.title(str_model_id + ' '+col4folds + ' ' + score_name, fontsize = 20 )\n    plt.ylabel(score_name, fontsize = 16)\n    plt.xlabel('index of target (ordered by 1 frequency)', fontsize = 16)\n    plt.show()    \n    \n    plt.figure(figsize  = (20, 6) )    \n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        v = df_scores[l_col].std(axis = 1)\n        plt.plot(v, label = kw)\n        print(kw,'std rocauc score:',v.mean(), str_model_id + ' '+col4folds)\n    plt.grid()\n    plt.legend()\n    plt.title('STD on folds '+str_model_id + ' '+col4folds + ' ' + score_name, fontsize = 20 )\n    plt.ylabel('std '+score_name, fontsize = 16)\n    plt.xlabel('index of target (ordered by 1 frequency)', fontsize = 16)\n    plt.show()    \n    \n    \n    plt.figure(figsize  = (20, 12) )    \n    for col in df_scores.columns:\n        v = df_scores[col]\n        plt.plot( v , label = col)\n    plt.grid()\n    plt.legend()\n    plt.title(str_model_id + ' '+col4folds+ ' ' + score_name, fontsize = 20 )\n    plt.ylabel(score_name, fontsize = 16)\n    plt.xlabel('index of target (ordered by 1 frequency)', fontsize = 16)\n    plt.show()    \n\n    \nfor k in dict_df_models_stat:\n    df_scores = dict_df_models_stat[k]\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        print(kw,'mean score:',df_scores[l_col].mean().mean().round(3), k)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-16T17:10:11.49162Z","iopub.execute_input":"2023-06-16T17:10:11.492063Z","iopub.status.idle":"2023-06-16T17:33:09.982892Z","shell.execute_reply.started":"2023-06-16T17:10:11.492032Z","shell.execute_reply":"2023-06-16T17:33:09.98141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfor k in dict_df_models_stat:\n    df_scores = dict_df_models_stat[k]\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        print(kw,'mean score:',df_scores[l_col].mean().mean().round(3),'std:',df_scores[l_col].std(axis=1).median().round(3) ,k)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:13:04.230202Z","iopub.execute_input":"2023-06-16T16:13:04.230668Z","iopub.status.idle":"2023-06-16T16:13:04.261699Z","shell.execute_reply.started":"2023-06-16T16:13:04.230636Z","shell.execute_reply":"2023-06-16T16:13:04.260595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Anti-correlation between train-std and test-mean might be natural - big std on train - poor predictability on test\n# zero correlation between test-std and test-mean is unclear\n# here std,mean are calculated across folds \nfor k in dict_df_models_stat:\n    print(k)\n    df_scores = dict_df_models_stat[k]\n    d_ = pd.DataFrame()\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        v = df_scores[l_col].std(axis = 1)\n        d_[kw+' std'] = v\n        d_[kw+' mean'] = df_scores[l_col].mean(axis = 1)\n    d_['index'] = range(len(d_))\n    display( d_.corr().round(2) )","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:13:04.263219Z","iopub.execute_input":"2023-06-16T16:13:04.263695Z","iopub.status.idle":"2023-06-16T16:13:04.339037Z","shell.execute_reply.started":"2023-06-16T16:13:04.263657Z","shell.execute_reply":"2023-06-16T16:13:04.337892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfor k in dict_df_models_stat:\n    plt.figure(figsize = (20,7))\n    df_scores = dict_df_models_stat[k]\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        print(l_col)\n        v = df_scores[l_col].std(axis = 1)\n        plt.plot(v, label = kw)\n        print(kw,'mean score:',df_scores[l_col].mean().mean().round(3), k)\n        print(kw,'average std %.4f'%v.mean(), k )\n\n    plt.title('STD Folds scores (rocauc) ' + k, fontsize = 20 )  \n    plt.legend(fontsize = 16)\n    plt.grid()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:13:04.340506Z","iopub.execute_input":"2023-06-16T16:13:04.340933Z","iopub.status.idle":"2023-06-16T16:13:05.699648Z","shell.execute_reply.started":"2023-06-16T16:13:04.340904Z","shell.execute_reply":"2023-06-16T16:13:05.69832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Second model \n\nThanks to SIMON VEITNER \"Simple MLP\": https://www.kaggle.com/code/simonveitner/simple-mlp - please upvote ! ","metadata":{}},{"cell_type":"code","source":"%%time\nimport gc\ngc.collect()\n\nfrom tqdm import tqdm\ntqdm.pandas()\n\nfrom keras.models import Sequential\nfrom keras.layers import Dense\n# measure roc auc score metric \nfrom tensorflow.keras.metrics import AUC\n\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.metrics import f1_score\nfrom sklearn.linear_model import Ridge\n\nepochs = 15; batch_size = 128\nstr_model_id = 'MLPkeras'\ndict_Y_pred_oof[str_model_id] = {} \n\n\nnfeats = X.shape[1]\nnlabels = Y.shape[1]\nmodel = Sequential()\nmodel.add(Dense(256, activation='relu', input_dim=nfeats))\nmodel.add(Dense(128, activation='relu'))\nmodel.add(Dense(nlabels, activation='sigmoid'))\nmodel.compile(loss='binary_crossentropy',\n                optimizer='adam',\n                metrics=[AUC()])\nmodel.summary()\n\n# col4folds = 'cv_folds_3'\n\nfor col4folds in  ['cv_folds_5']: # ,'cv_folds_10', 'cv_folds_15']: \n    df_scores = pd.DataFrame()\n    vf = df_folds[col4folds]\n    print(str_model_id, col4folds, 'Y.shape:', Y.shape )\n    t00 = time.time()\n    Y_pred_oof = np.zeros( Y.shape )\n    for uv in vf.unique():\n\n        # Set Keras MLP model\n        nfeats = X.shape[1]\n        nlabels = Y.shape[1]\n        model = Sequential()\n        model.add(Dense(256, activation='relu', input_dim=nfeats))\n        model.add(Dense(128, activation='relu'))\n        model.add(Dense(nlabels, activation='sigmoid'))\n        model.compile(loss='binary_crossentropy',\n                        optimizer='adam',\n                        metrics=[AUC()])\n        \n        \n        m = vf == uv\n        IX_train = np.where( ~m)[0]\n        IX_test  = np.where(  m)[0]\n        print('Fold:',uv, 'n_samples in test:', m.sum(), 'Indexes and their lens:', len(IX_train), len(IX_test), IX_train[:10], IX_test[:10] )\n\n        t0 = time.time()\n#         model.fit(X[IX_train,:], Y[IX_train,:] )\n        model.fit(X[IX_train,:],Y[IX_train,:], epochs=epochs, batch_size=batch_size) # \n        print('%.1f secs passed on fit'%(time.time() - t0 ))\n\n        t0 = time.time()\n        Y_pred = model.predict(X[IX_train,:])\n        print('%.1f secs passed on predict'%(time.time() - t0 ))\n\n        ################################ predicts and scores calculations  #################################\n        IX_loc = IX_test\n        Y_pred = model.predict(X[IX_loc,:])\n        t0 = time.time()\n        l_ = []\n        for i_ in range(Y.shape[1]):\n            s = roc_auc_score(Y[IX_loc,i_].ravel(),  Y_pred[:,i_].ravel() ); l_.append(s)\n        print('%.1f secs passed on test roc-auc'%(time.time() - t0 ), 'mean rocauc: %.3f'%(np.mean(l_)))\n        df_scores[str_model_id + ' Test'+str(uv)] = l_\n        Y_pred_oof[IX_loc,:] = Y_pred\n\n        IX_loc = IX_train \n        Y_pred = model.predict(X[IX_loc,:])\n        t0 = time.time()\n        l_ = []\n        for i_ in range(Y.shape[1]):\n            s = roc_auc_score(Y[IX_loc,i_].ravel(),  Y_pred[:,i_].ravel() ); l_.append(s)\n        print('%.1f secs passed on train roc-auc'%(time.time() - t0 ), 'mean rocauc: %.3f'%(np.mean(l_)))\n        df_scores[str_model_id + ' Train'+str(uv)] = l_\n        \n        gc.collect()\n\n        \n    dict_Y_pred_oof[str_model_id][col4folds] = Y_pred_oof\n    dict_df_models_stat[str_model_id + ' '+col4folds] =   df_scores.copy()\n\n    print('%.1f secs passed on model'%(time.time() - t00 ), str_model_id , 'Folds:', col4folds)\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        print(kw,'mean rocauc score:',df_scores[l_col].mean().mean().round(3), str_model_id + ' '+col4folds)\n    print()\n    \n    ######################################### Plots ######################################### \n    score_name = 'rocauc'\n    plt.figure(figsize  = (20, 6) )    \n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        v = df_scores[l_col].mean(axis = 1)\n        plt.plot(v, label = kw)\n        print(kw,'mean rocauc score:',v.mean(), str_model_id + ' '+col4folds)\n    plt.grid()\n    plt.legend()\n    plt.title(str_model_id + ' '+col4folds + ' ' + score_name, fontsize = 20 )\n    plt.ylabel(score_name, fontsize = 16)\n    plt.xlabel('index of target (ordered by 1 frequency)', fontsize = 16)\n    plt.show()    \n    \n    plt.figure(figsize  = (20, 6) )    \n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        v = df_scores[l_col].std(axis = 1)\n        plt.plot(v, label = kw)\n        print(kw,'std rocauc score:',v.mean(), str_model_id + ' '+col4folds)\n    plt.grid()\n    plt.legend()\n    plt.title('STD on folds '+str_model_id + ' '+col4folds + ' ' + score_name, fontsize = 20 )\n    plt.ylabel(score_name, fontsize = 16)\n    plt.xlabel('index of target (ordered by 1 frequency)', fontsize = 16)\n    plt.show()    \n    \n    \n    plt.figure(figsize  = (20, 12) )    \n    for col in df_scores.columns:\n        v = df_scores[col]\n        plt.plot( v , label = col)\n    plt.grid()\n    plt.legend()\n    plt.title(str_model_id + ' '+col4folds+ ' ' + score_name, fontsize = 20 )\n    plt.ylabel(score_name, fontsize = 16)\n    plt.xlabel('index of target (ordered by 1 frequency)', fontsize = 16)\n    plt.show()    \n\n    \nfor k in dict_df_models_stat:\n    df_scores = dict_df_models_stat[k]\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        print(kw,'mean score:',df_scores[l_col].mean().mean().round(3), k)\n    ","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:13:05.701431Z","iopub.execute_input":"2023-06-16T16:13:05.701965Z","iopub.status.idle":"2023-06-16T16:32:23.633218Z","shell.execute_reply.started":"2023-06-16T16:13:05.701932Z","shell.execute_reply":"2023-06-16T16:32:23.631792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# Anti-correlation between train-std and test-mean might be natural - big std on train - poor predictability on test\n# zero correlation between test-std and test-mean is unclear\n# here std,mean are calculated across folds \nfor k in dict_df_models_stat:\n    print(k)\n    df_scores = dict_df_models_stat[k]\n    d_ = pd.DataFrame()\n    for kw in ['Test', 'Train']:\n        l_col = [col for col in df_scores.columns if kw in col ]\n        v = df_scores[l_col].std(axis = 1)\n        d_[kw+' std'] = v\n        d_[kw+' mean'] = df_scores[l_col].mean(axis = 1)\n    d_['index'] = range(len(d_))\n    display( d_.corr().round(2) )","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:32:23.635198Z","iopub.execute_input":"2023-06-16T16:32:23.635732Z","iopub.status.idle":"2023-06-16T16:32:23.736717Z","shell.execute_reply.started":"2023-06-16T16:32:23.635651Z","shell.execute_reply":"2023-06-16T16:32:23.735595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for k in dict_df_models_stat:\n    print(k,dict_df_models_stat[k].keys() )","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:32:23.73812Z","iopub.execute_input":"2023-06-16T16:32:23.738566Z","iopub.status.idle":"2023-06-16T16:32:23.745918Z","shell.execute_reply.started":"2023-06-16T16:32:23.738531Z","shell.execute_reply":"2023-06-16T16:32:23.744753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2023-06-16T16:32:23.747528Z","iopub.execute_input":"2023-06-16T16:32:23.747892Z","iopub.status.idle":"2023-06-16T16:32:23.759504Z","shell.execute_reply.started":"2023-06-16T16:32:23.747864Z","shell.execute_reply":"2023-06-16T16:32:23.758089Z"},"trusted":true},"execution_count":null,"outputs":[]}]}