{"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 run \"Diamond\" package to get similirity relations for the proteins in CAFA5 Kaggle challenge.\n\nWe consider several options - from \"fast\" to \"ultra-senstive\", etc to get various quality of similarity - see different versions of the notebook - different parameters.\n\nThe results are used for example to create groups and groupwise validation for CAFA5: \nhttps://www.kaggle.com/code/alexandervc/cafa5-23-groups-and-folds-diamond-igraph\n\nThanks to L. Gereseva for bringing Diamond to Kaggle: \nhttps://www.kaggle.com/code/geraseva/diamond","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-10T12:48:00.104813Z","iopub.execute_input":"2023-06-10T12:48:00.105175Z","iopub.status.idle":"2023-06-10T12:48:00.114596Z","shell.execute_reply.started":"2023-06-10T12:48:00.105148Z","shell.execute_reply":"2023-06-10T12:48:00.113584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!wget http://github.com/bbuchfink/diamond/releases/download/v2.1.6/diamond-linux64.tar.gz\n!tar xzf diamond-linux64.tar.gz\n!rm diamond-linux64.tar.gz","metadata":{"execution":{"iopub.status.busy":"2023-06-10T12:48:06.02148Z","iopub.execute_input":"2023-06-10T12:48:06.021831Z","iopub.status.idle":"2023-06-10T12:48:10.786649Z","shell.execute_reply.started":"2023-06-10T12:48:06.021805Z","shell.execute_reply":"2023-06-10T12:48:10.785522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create a database file from training data - 5 seconds","metadata":{}},{"cell_type":"code","source":"%%time \nimport time\nfrom subprocess import Popen, PIPE\ndb_name='train_db'\n\n#------------------------------------------------------------------------------------------------\nprint('Create a database file from training data - 5 seconds')\n#------------------------------------------------------------------------------------------------\nt0 = time.time()\np = Popen(['./diamond', 'makedb', \n           '--in', '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta',\n            '-d', db_name], stdin=PIPE, stdout=PIPE)\nstdout, stderr = p.communicate()\nprint('Database created in %.1f seconds'%(time.time() - t0 ))\n","metadata":{"execution":{"iopub.status.busy":"2023-06-10T12:48:27.24975Z","iopub.execute_input":"2023-06-10T12:48:27.250205Z","iopub.status.idle":"2023-06-10T12:48:32.833641Z","shell.execute_reply.started":"2023-06-10T12:48:27.25017Z","shell.execute_reply":"2023-06-10T12:48:32.832641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n!./diamond deepclust -d \"train_db.dmnd\" -o \"clusters2.csv\" --approx-id 30 -M 32G","metadata":{"execution":{"iopub.status.busy":"2023-06-10T13:03:09.041814Z","iopub.execute_input":"2023-06-10T13:03:09.042868Z","iopub.status.idle":"2023-06-10T13:07:20.21619Z","shell.execute_reply.started":"2023-06-10T13:03:09.042832Z","shell.execute_reply":"2023-06-10T13:07:20.214685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(\"clusters2.csv\",header = None, sep = '\\t')\ndf[0].value_counts().value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-06-10T13:12:53.937694Z","iopub.execute_input":"2023-06-10T13:12:53.938192Z","iopub.status.idle":"2023-06-10T13:12:54.116803Z","shell.execute_reply.started":"2023-06-10T13:12:53.938151Z","shell.execute_reply":"2023-06-10T13:12:54.115707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!./diamond cluster -d \"train_db.dmnd\" -o \"clusters.csv\" --approx-id 30 -M 64G","metadata":{"execution":{"iopub.status.busy":"2023-06-10T12:50:24.964787Z","iopub.execute_input":"2023-06-10T12:50:24.965324Z","iopub.status.idle":"2023-06-10T12:54:05.739609Z","shell.execute_reply.started":"2023-06-10T12:50:24.965294Z","shell.execute_reply":"2023-06-10T12:54:05.738092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df[0].value_counts().value_counts()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-10T13:00:02.303573Z","iopub.execute_input":"2023-06-10T13:00:02.303984Z","iopub.status.idle":"2023-06-10T13:00:02.35618Z","shell.execute_reply.started":"2023-06-10T13:00:02.303952Z","shell.execute_reply":"2023-06-10T13:00:02.355242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!head -n 100 \"clusters.csv\"\ndf = pd.read_csv(\"clusters.csv\",header = None, sep = '\\t')\ndf[1].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-06-10T12:57:12.140758Z","iopub.execute_input":"2023-06-10T12:57:12.141617Z","iopub.status.idle":"2023-06-10T12:57:12.358615Z","shell.execute_reply.started":"2023-06-10T12:57:12.141583Z","shell.execute_reply":"2023-06-10T12:57:12.357939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Run a blastp-like search ","metadata":{}},{"cell_type":"code","source":"%%time \n# k=50 , train:\n# Execution time: 550.8059980869293s\n# CPU times: user 30.5 ms, sys: 6.54 ms, total: 37.1 ms\n# Wall time: 9min 10s\n# k=16 , test (original)\n# Execution time: 240.35088419914246s\n\noutfile_name='matches_train.tsv'\nk = 142000\n\nif 1:\n    #------------------------------------------------------------------------------------------------\n    print('Run a blastp-like search for test set')\n    #------------------------------------------------------------------------------------------------\n    time0 = time.time() \n    #\n    fn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta' # Train\n    # fn = '/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta' # Test \n    p = Popen(['./diamond', 'blastp', '-d', db_name, '--evalue',str(1_000_000),\n               '-q', fn,\n                '-o', outfile_name, '--max-target-seqs', str(k),'--ultra-sensitive'], stdin=PIPE, stdout=PIPE) # , '--very-sensitive' #  '--quiet' #  '--sensitive' # '--more-sensitive' # \n    stdout, stderr = p.communicate()\n    print(f'Execution time: {time.time()-time0}s')\n\n    matches=pd.read_csv(outfile_name, sep='\\t', header=None, \n                        names=['qseqid', 'sseqid', 'pident', 'length', 'mismatch', \n                               'gapopen', 'qstart', 'qend', 'sstart','send', 'evalue', 'bitscore'])\n    matches['qseqid']=matches['qseqid'].apply(lambda x: x.split('\\\\t')[0])\n    matches.head(10)    \n    \n    print(matches.shape)\n    display( matches.describe() )\n    \n    v = matches['qseqid'].value_counts()\n    plt.hist(v, bins = 100)\n    plt.show()\n    v.describe()    ","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:23:00.252183Z","iopub.execute_input":"2023-06-08T14:23:00.252929Z","iopub.status.idle":"2023-06-08T14:23:00.261869Z","shell.execute_reply.started":"2023-06-08T14:23:00.252888Z","shell.execute_reply":"2023-06-08T14:23:00.26077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2023-06-08T09:57:28.591316Z","iopub.execute_input":"2023-06-08T09:57:28.592027Z","iopub.status.idle":"2023-06-08T09:57:35.152763Z","shell.execute_reply.started":"2023-06-08T09:57:28.591972Z","shell.execute_reply":"2023-06-08T09:57:35.151337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2023-06-08T10:02:20.238342Z","iopub.execute_input":"2023-06-08T10:02:20.238796Z","iopub.status.idle":"2023-06-08T10:02:21.314249Z","shell.execute_reply.started":"2023-06-08T10:02:20.238758Z","shell.execute_reply":"2023-06-08T10:02:21.313176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sel_seq = 'MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAEPKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGDWWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRESETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGLCHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTLKPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKYLRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYTARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVERGYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL'\n# Sequence for SRC human protein \n# https://en.wikipedia.org/wiki/Proto-oncogene_tyrosine-protein_kinase_Src\n# Uniprot: https://www.uniprot.org/uniprotkb/P12931/entry\nsel_prot_id = 'P12931'\nprot_name = 'SRC human P12931'        ","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:10:59.384213Z","iopub.execute_input":"2023-06-08T14:10:59.384621Z","iopub.status.idle":"2023-06-08T14:10:59.390394Z","shell.execute_reply.started":"2023-06-08T14:10:59.384591Z","shell.execute_reply":"2023-06-08T14:10:59.389152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ofile = open(\"small_fasta.txt\", \"w\")\n\nfor i in range(1):\n    ofile.write(\">\" + prot_name + \"\\n\" +sel_seq + \"\\n\")\n\n#do not forget to close it\nofile.close()\nprint(os.listdir())","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:11:00.825311Z","iopub.execute_input":"2023-06-08T14:11:00.825751Z","iopub.status.idle":"2023-06-08T14:11:00.83268Z","shell.execute_reply.started":"2023-06-08T14:11:00.825714Z","shell.execute_reply":"2023-06-08T14:11:00.8315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"outfile_name='matches_SRC.tsv'\nk = 142050 # Execution time: 72.21724462509155s # 4657 \n#------------------------------------------------------------------------------------------------\nprint('Run a blastp-like search')\n#------------------------------------------------------------------------------------------------\ntime0 = time.time() \n#\n# fn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta' # Train\nfn = 'small_fasta.txt'\n# fn = '/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta' # Test \np = Popen(['./diamond', 'blastp', '-d', db_name,  '--evalue',str(1000000000000000.1),\n           '-q', fn,\n            '-o', outfile_name, '--max-target-seqs', str(k),'--quiet','--ultra-sensitive' \n          ], stdin=PIPE, stdout=PIPE) # '--quiet' # #   '--sensitive' # '--more-sensitive' # \nstdout, stderr = p.communicate()\nprint(f'Execution time: {time.time()-time0}s')\n\n# Cannot make work additional output \n#  '-outfmt 6 qseqid sseqid',\n# diamond blastp --query query.fasta --db database.dmnd --outfmt 6 qseqid sseqid qlen slen\n#           '--outfmt 6 qseqid sseqid qlen slen pident length mismatch gapopen qstart qend sstart send evalue bitscore' \n","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:28:10.44191Z","iopub.execute_input":"2023-06-08T14:28:10.442366Z","iopub.status.idle":"2023-06-08T14:29:22.455458Z","shell.execute_reply.started":"2023-06-08T14:28:10.44233Z","shell.execute_reply":"2023-06-08T14:29:22.454316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"matches=pd.read_csv(outfile_name, 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(50)","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:30:06.51455Z","iopub.execute_input":"2023-06-08T14:30:06.514925Z","iopub.status.idle":"2023-06-08T14:30:06.565259Z","shell.execute_reply.started":"2023-06-08T14:30:06.514896Z","shell.execute_reply":"2023-06-08T14:30:06.564159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(matches['bitscore'], bins = 1000)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:12:32.602271Z","iopub.execute_input":"2023-06-08T14:12:32.602693Z","iopub.status.idle":"2023-06-08T14:12:34.246322Z","shell.execute_reply.started":"2023-06-08T14:12:32.602659Z","shell.execute_reply":"2023-06-08T14:12:34.2441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display( matches.describe(percentiles = [0.25,0.5,0.75,0.9,0.95,0.96,0.97,0.99]) )\ndisplay(matches.corr())\ndisplay(matches.corr(method = 'spearman'))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-06-08T14:12:34.248292Z","iopub.execute_input":"2023-06-08T14:12:34.248793Z","iopub.status.idle":"2023-06-08T14:12:34.354771Z","shell.execute_reply.started":"2023-06-08T14:12:34.248754Z","shell.execute_reply":"2023-06-08T14:12:34.353797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}