{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# What is about ? \n\nSequences alignments are extremely important algorithms used in modern bioinformatics: https://en.wikipedia.org/wiki/Sequence_alignment\nSee e.g. very nice BioPython Kaggle tutorial notebook by ANDREY SHTRAUSS : https://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#5-|-PAIRWISE-SEQUENCE-ALIGNMENT (and other notebooks by the same author).\n\nHere we explore local alignment   from another package - \"skbio\":\nhttp://scikit-bio.org/docs/0.4.1/index.html  which works faster than BioPyton, but less documented and there are some other cavets.\n\n**Timing:** For example the simplest version: 142246 alignments  in 13.7 secs (mean len 533). For more complicated version (with use of blosum62 substitution matrix) - 40 secs . \nThe first variant might require 4.5 hours in BioPython - instead of 13.7 seconds with skbio. \n\nWe   make some simple comparisons of the local alignments with different params.\nWe take some proteins from the CAFA5 competition and align them.\n\n**Outcomes:** we see that incorporation  substitution matrix is quite important for the protein alignments (at least for skbio).\nThe results without substitution matrix differs from the one with it. \nWe see the difference in simple examples and well as we see different distributions of scores and similarities on CAFA5 data. \n\nIt seems that some params for BioPython and skbio local alignments are different. \nBecause in the simplest example for seq1 = \"EVSAW\" , seq2 = \"KEVLA\" the skbio gives alignment consisting only from one letter - common \"A\" for both words - that is quite different from others. \n\n\n\nPS \n\nFor the competition purposes one can use alignments in a several ways:\nfind most similar proteins to the given one and try to transfer labels from them to the it. \nOr 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\nhttps://www.kaggle.com/code/alexandervc/cafa5-levenshtein-distance-features). \nOr one can use \"similarities\" to obtain groups for groupwise validation. \nThough the main cavet is that alignment are quite slow.\n\nPS\n\nThanks to : https://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#5-|-PAIRWISE-SEQUENCE-ALIGNMENT\n\n\"Diamond\" package works much faster and allows extremely fast alignment many x many :\nSee Liza Geraseva notebook: https://www.kaggle.com/code/geraseva/diamond\n\nBiopython works quite slower see examples: https://www.kaggle.com/code/alexandervc/cafa5-18-alignments-biopython-compare\n\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# 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":{"execution":{"iopub.status.busy":"2023-05-31T15:04:45.369376Z","iopub.execute_input":"2023-05-31T15:04:45.370185Z","iopub.status.idle":"2023-05-31T15:04:45.379796Z","shell.execute_reply.started":"2023-05-31T15:04:45.37013Z","shell.execute_reply":"2023-05-31T15:04:45.377887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install scikit-bio","metadata":{"execution":{"iopub.status.busy":"2023-05-31T15:04:45.837924Z","iopub.execute_input":"2023-05-31T15:04:45.838504Z","iopub.status.idle":"2023-05-31T15:06:10.176406Z","shell.execute_reply.started":"2023-05-31T15:04:45.838472Z","shell.execute_reply":"2023-05-31T15:06:10.175227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Substitution matrices\n\nhttps://en.wikipedia.org/wiki/Substitution_matrix\n\nSubstitution matrix describes the frequency at which a character in a nucleotide sequence or a protein sequence changes to other character states over evolutionary time\n","metadata":{}},{"cell_type":"markdown","source":"## Identity substitution matrix from skbio - format double dict\n\nskibio has(depricated) make_identity_substitution_matrix function - which returns simple example of the  substitution matrix ","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(sm['A'])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T15:09:07.66178Z","iopub.execute_input":"2023-05-31T15:09:07.662182Z","iopub.status.idle":"2023-05-31T15:09:07.67844Z","shell.execute_reply.started":"2023-05-31T15:09:07.662153Z","shell.execute_reply":"2023-05-31T15:09:07.677054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Take blosum62 substitution matrix for protein's amino-acids from BioPython and convert it to skibio format (double dictionary)","metadata":{}},{"cell_type":"code","source":"from 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-05-31T15:47:08.457773Z","iopub.execute_input":"2023-05-31T15:47:08.458133Z","iopub.status.idle":"2023-05-31T15:47:08.470952Z","shell.execute_reply.started":"2023-05-31T15:47:08.4581Z","shell.execute_reply":"2023-05-31T15:47:08.469458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Protein alignments toy examples  ","metadata":{}},{"cell_type":"markdown","source":"## Strange results for skbio local alignment without substitution matrix - we get just one letter \"A\" as an alignment for \"EVSAW\" and \"KEVLA\"\n","metadata":{}},{"cell_type":"code","source":"from skbio.alignment import StripedSmithWaterman\nseq1 = \"EVSAW\"\nseq2 = \"KEVLA\"\n\nquery = StripedSmithWaterman(seq1)# \"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'] )\ndir(alignment)\ndisplay(alignment )\n\nprint(); \nprint('---------------- alignment with identity substitution matrix ----------------------')\nprint();\n\nquery = StripedSmithWaterman(seq1, substitution_matrix = sm2)# \"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'] )\ndir(alignment)\ndisplay(alignment )\n\nprint(); \nprint('---------------- alignment with blosum62 substitution matrix ----------------------')\nprint();\n\nquery = StripedSmithWaterman(seq1, substitution_matrix = dict_sm_blosum62)# \"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'] )\ndir(alignment)\ndisplay(alignment )\n\n\nprint(); \nprint('---------------- global alignment from BioPython ----------------------')\nprint();\n\nfrom Bio import Align\naligner = Align.PairwiseAligner()\n#alignments = aligner.align(\"TACCG\", \"ACG\")\nalignments = aligner.align(seq1, seq2 )\nprint(alignments[0])\n\nprint(); \nprint('---------------- local alignment from BioPython ----------------------')\nprint();\n\nfrom Bio import Align\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\nalignments = aligner.align(seq1, seq2 )\nprint(alignments[0])\n\n\n\nprint(); \nprint('---------------- local alignment from BioPython with BLOSUM62 ----------------------')\nprint();\n\nfrom Bio import Align\nfrom Bio.Align import substitution_matrices\naligner = Align.PairwiseAligner()\naligner.mode = 'local'\naligner.substitution_matrix = substitution_matrices.load(\"BLOSUM62\")\nalignments = aligner.align(seq1, seq2 )\nprint(alignments[0])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-31T16:48:53.567579Z","iopub.execute_input":"2023-05-31T16:48:53.567934Z","iopub.status.idle":"2023-05-31T16:48:53.598251Z","shell.execute_reply.started":"2023-05-31T16:48:53.567905Z","shell.execute_reply":"2023-05-31T16:48:53.596785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load CAFA5 train data examples ","metadata":{"execution":{"iopub.status.busy":"2023-05-31T12:50:32.380554Z","iopub.execute_input":"2023-05-31T12:50:32.381491Z","iopub.status.idle":"2023-05-31T12:50:32.401601Z","shell.execute_reply.started":"2023-05-31T12:50:32.381455Z","shell.execute_reply":"2023-05-31T12:50:32.40029Z"}}},{"cell_type":"code","source":"from Bio import SeqIO\nfn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\nprint(\"Sequence example:\\n\\n\", next(iter(SeqIO.parse(fn, \"fasta\"))))\nsequences = SeqIO.parse(fn, \"fasta\")\nnum_sequences = sum(1 for seq in sequences)\nprint()\nprint(\"Number of sequences:\", num_sequences)\n\nsequences = SeqIO.parse(fn, \"fasta\")\nlist_seqs = [str(seq.seq) for seq in sequences]\nlist_lens = [len(t) for t in list_seqs]\nprint('Mean len: %.1f'%(np.mean(list_lens)))\nprint(len(list_seqs),  list_seqs[:3])\nsequences = SeqIO.parse(fn, \"fasta\")\nlist_ids = [str(seq.id) for seq in sequences]\nprint(len(list_ids),  list_ids[:3])","metadata":{"execution":{"iopub.status.busy":"2023-05-31T15:31:22.044562Z","iopub.execute_input":"2023-05-31T15:31:22.044953Z","iopub.status.idle":"2023-05-31T15:31:25.981228Z","shell.execute_reply.started":"2023-05-31T15:31:22.044925Z","shell.execute_reply":"2023-05-31T15:31:25.979796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time \nfrom skbio.alignment import StripedSmithWaterman\ns1 = list_seqs[0] # Protein(\"HEAGAWGHEE\")\ns2 = list_seqs[0] # Protein(\"HEAGAWGHEE\")\nquery = StripedSmithWaterman(s1)# \"HEAGAWGHEE\") # \"ACTAAGGCTCTCTACCCCTCTCAGAGA\")\nalignment = query(s2)# \"AAAAAACTCTCTAAACTCACTAAGGCTCTCTACCCCTCTTCAGAGAAGTCGA\")\nprint(alignment)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T15:29:42.640159Z","iopub.execute_input":"2023-05-31T15:29:42.64052Z","iopub.status.idle":"2023-05-31T15:29:42.647389Z","shell.execute_reply.started":"2023-05-31T15:29:42.64049Z","shell.execute_reply":"2023-05-31T15:29:42.646148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Alignments on CAFA5 examples  142K proteins ","metadata":{}},{"cell_type":"code","source":"%%time\nimport time\nquery_sequence =  list_seqs[0]\nprint('len of the first alignment sequence:', len(query_sequence))\n\ndict_df = {}\nfor prm in ['no subs matrix', 'blosum62', 'identity subs matrix' ]:\n    print('Alignment:', prm)\n    t0 = time.time()\n    if prm == 'no subs matrix':\n        query = StripedSmithWaterman(query_sequence) #  , substitution_matrix = sm2)\n    elif prm == 'blosum62':\n        query = StripedSmithWaterman(query_sequence, substitution_matrix = dict_sm_blosum62)# \"ACTAAGGCTCTCTACCCCTCTCAGAGA\")\n    elif prm == 'identity subs matrix':\n        query = StripedSmithWaterman(query_sequence, substitution_matrix = sm2)# \"ACTAAGGCTCTCTACCCCTCTCAGAGA\")\n\n    \n    l = [ query(s2) for s2 in list_seqs]# Alignment made here\n    \n    print(len(l), 'alignments done in %.1f secs'%(time.time()-t0))\n\n    # -------------------------------------------------------------------\n    # Collect statistcs:\n    # -------------------------------------------------------------------\n    \n    df = pd.DataFrame()\n    ls = [t['optimal_alignment_score'] for t in l]\n    l_l = [ (1+t['target_end_optimal'] - t['target_begin']) for t in l] \n    df['Id'] = list_ids\n    df['score'] = ls\n    df['similarity'] = np.array(ls)/np.array(l_l)\n    df['alignment len'] = l_l\n    df['len'] = [len(t) for t in list_seqs]\n    display(df.describe())\n#     display( df.sort_values('similarity', ascending = False).head(5) )\n#     display( df.sort_values('score', ascending = False).head(5) )\n#     display( df.sort_values('score', ascending = False).tail(5) )\n    d1 = ( df.sort_values('similarity', ascending = False).head(5) ).reset_index()\n    d2 = ( df.sort_values('score', ascending = False).head(5) ).reset_index()\n    d3 = ( df.sort_values('score', ascending = False).tail(5) ).reset_index()\n    print('Top5 similarity, score and tail score data:')\n    display( pd.concat( (d1,d2,d3), axis = 1) )\n            \n    ll = df.columns[1:]\n    plt.figure(figsize = (20,5))\n    for i,col in enumerate(ll):\n        plt.subplot(1,len(ll),i+1 )\n        plt.hist(df[col], bins = 100 )\n        plt.title(col ,fontsize = 15 )\n        plt.suptitle(str(prm),fontsize = 15)\n    plt.show()\n    \n    dict_df[prm] = df\n    df.to_csv('df_stat'+str(prm).replace(' ','_')+'.csv' )\n    print()\n    print('--------------------------------------------------------------------------------------')\n    print()\n    \n    ","metadata":{"execution":{"iopub.status.busy":"2023-05-31T16:07:51.190224Z","iopub.execute_input":"2023-05-31T16:07:51.190621Z","iopub.status.idle":"2023-05-31T16:09:09.50684Z","shell.execute_reply.started":"2023-05-31T16:07:51.190589Z","shell.execute_reply":"2023-05-31T16:09:09.505295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis of the agreement methods without substitution matrix - not agree with others, while others agree 30-50%. ","metadata":{}},{"cell_type":"code","source":"%%time\nfrom itertools import combinations\n\nN = 100\nfor col in ['similarity', 'score']:\n    print('Compare by:', col)\n    for prm1,prm2 in combinations(list(dict_df.keys()), 2 ):\n        print(prm1,' vs  ', prm2)\n        df1 = dict_df[prm1]\n        df2 = dict_df[prm2]\n        s1 = set(df1.sort_values(col, ascending = False)['Id'].iloc[:N]) \n        s2 = set(df2.sort_values(col, ascending = False)['Id'].iloc[:N]) \n        s = s1&s2\n        print('Intersection size %d out of %d'%(len(s), N))\n        print()\n    print()\n    print()\n        ","metadata":{"execution":{"iopub.status.busy":"2023-05-31T16:15:11.417465Z","iopub.execute_input":"2023-05-31T16:15:11.41784Z","iopub.status.idle":"2023-05-31T16:15:11.786888Z","shell.execute_reply.started":"2023-05-31T16:15:11.417812Z","shell.execute_reply":"2023-05-31T16:15:11.785499Z"},"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":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{},"execution_count":null,"outputs":[]}]}