{"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":"I used the CAFA5 dataset to explore ML with the Naive Bayes Method. Plesae see Akhil Bhargava's really great notebook (link below). With the help of Akhil's notebook I was able to add the protein classifcation information to the CAFA5 dataset and perform ML using protein classification information. \nhttps://www.kaggle.com/code/abharg16/predicting-protein-classification","metadata":{}},{"cell_type":"code","source":"# 1). ----- Import Libraries and Datasets ------\n\nimport pandas as pd\nimport numpy as np\nfrom matplotlib import pyplot as plt\nimport seaborn as sns\nfrom sklearn.feature_extraction.text import CountVectorizer\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import accuracy_score, confusion_matrix, classification_report\nfrom Bio import SeqIO\nimport tensorflow as tf\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.preprocessing import MultiLabelBinarizer\n\n\n","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2023-05-01T00:21:15.947886Z","iopub.execute_input":"2023-05-01T00:21:15.948249Z","iopub.status.idle":"2023-05-01T00:21:15.960105Z","shell.execute_reply.started":"2023-05-01T00:21:15.948187Z","shell.execute_reply":"2023-05-01T00:21:15.959174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Import Datasets\ntrain_set_file_path = '../input/cafa-5-protein-function-prediction/Train/train_terms.tsv'\ntrain_set = pd.read_csv(train_set_file_path, sep='\\t').dropna().drop_duplicates()\n\nUNIQUE_LABELS = train_set['term'].unique()\n\n# Load train_taxonomy.tsv\ntrain_taxonomy = pd.read_csv('/kaggle/input/cafa-5-protein-function-prediction/Train/train_taxonomy.tsv', sep='\\t')\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-30T23:43:49.984077Z","iopub.execute_input":"2023-04-30T23:43:49.984632Z","iopub.status.idle":"2023-04-30T23:43:58.149894Z","shell.execute_reply.started":"2023-04-30T23:43:49.984573Z","shell.execute_reply":"2023-04-30T23:43:58.149232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\nfastaPath = ('/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta')\n\ndef GetDfFromFasta(fastaPath):    \n    input_file = fastaPath\n    fasta_sequences = SeqIO.parse(open(input_file),'fasta')\n    nameList = [] \n    sequenceList = []\n    for fasta in fasta_sequences:\n        name, sequence = fasta.id, str(fasta.seq)\n        nameList.append(name)\n        sequenceList.append(sequence)\n        \n    return pd.DataFrame(list(zip(nameList, sequenceList)), columns={\"EntryID\",\"seq\"})\n","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:43:58.205556Z","iopub.execute_input":"2023-04-30T23:43:58.205845Z","iopub.status.idle":"2023-04-30T23:43:58.218649Z","shell.execute_reply.started":"2023-04-30T23:43:58.205787Z","shell.execute_reply":"2023-04-30T23:43:58.217814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load train_sequences.fasta and parse it\ntrain_sequences = []\nfor record in SeqIO.parse('/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta', 'fasta'):\n    train_sequences.append((str(record.id), str(record.seq)))\n","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:43:58.219861Z","iopub.execute_input":"2023-04-30T23:43:58.220179Z","iopub.status.idle":"2023-04-30T23:44:00.863995Z","shell.execute_reply.started":"2023-04-30T23:43:58.220131Z","shell.execute_reply":"2023-04-30T23:44:00.863045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimport pandas as pd \ndf = pd.DataFrame(train_sequences) \n\n","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:59:19.590808Z","iopub.execute_input":"2023-04-30T23:59:19.591151Z","iopub.status.idle":"2023-04-30T23:59:19.617263Z","shell.execute_reply.started":"2023-04-30T23:59:19.591102Z","shell.execute_reply":"2023-04-30T23:59:19.616445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"   \ndf.set_axis(['EntryID', 'seq'], axis='columns', inplace=True)\n   ","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:44:00.904904Z","iopub.execute_input":"2023-04-30T23:44:00.905201Z","iopub.status.idle":"2023-04-30T23:44:00.910413Z","shell.execute_reply.started":"2023-04-30T23:44:00.905143Z","shell.execute_reply":"2023-04-30T23:44:00.909657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df3 = pd.merge(train_set, train_taxonomy, on=[\"EntryID\"])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-30T23:44:00.911794Z","iopub.execute_input":"2023-04-30T23:44:00.912157Z","iopub.status.idle":"2023-04-30T23:44:02.340615Z","shell.execute_reply.started":"2023-04-30T23:44:00.912089Z","shell.execute_reply":"2023-04-30T23:44:02.339901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df4 = pd.merge(df3, df, on=[\"EntryID\"])","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-30T23:44:02.341611Z","iopub.execute_input":"2023-04-30T23:44:02.341997Z","iopub.status.idle":"2023-04-30T23:44:03.908255Z","shell.execute_reply.started":"2023-04-30T23:44:02.341943Z","shell.execute_reply":"2023-04-30T23:44:03.907305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"\n# Import Datasets\ndf_seq = pd.read_csv('/kaggle/input/protein-data-set/pdb_data_seq.csv')\ndf_char = pd.read_csv('/kaggle/input/protein-data-set/pdb_data_no_dups.csv')\n\n\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-30T23:44:03.931051Z","iopub.execute_input":"2023-04-30T23:44:03.931405Z","iopub.status.idle":"2023-04-30T23:44:08.327752Z","shell.execute_reply.started":"2023-04-30T23:44:03.931347Z","shell.execute_reply":"2023-04-30T23:44:08.326932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 2). ----- Filter and Process Dataset ------\n\n# Filter for only proteins\nprotein_char = df_char[df_char.macromoleculeType == 'Protein']\nprotein_seq = df_seq[df_seq.macromoleculeType == 'Protein']\n\n# Select only necessary variables to join\nprotein_char = protein_char[['structureId','classification']]\nprotein_seq = protein_seq[['structureId','sequence']]\nprotein_seq.head()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:44:08.328904Z","iopub.execute_input":"2023-04-30T23:44:08.329161Z","iopub.status.idle":"2023-04-30T23:44:08.443096Z","shell.execute_reply.started":"2023-04-30T23:44:08.329117Z","shell.execute_reply":"2023-04-30T23:44:08.442371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df5 = pd.merge(protein_char, protein_seq, on=['structureId'])","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:50:21.115536Z","iopub.execute_input":"2023-04-30T23:50:21.115988Z","iopub.status.idle":"2023-04-30T23:50:21.246814Z","shell.execute_reply.started":"2023-04-30T23:50:21.115942Z","shell.execute_reply":"2023-04-30T23:50:21.24612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = {'seq': 'sequence'}\ndf7 = df4.rename(columns=d)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:51:24.221116Z","iopub.execute_input":"2023-04-30T23:51:24.221492Z","iopub.status.idle":"2023-04-30T23:51:25.113614Z","shell.execute_reply.started":"2023-04-30T23:51:24.221432Z","shell.execute_reply":"2023-04-30T23:51:25.112821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df6 = pd.merge(df7, df5, on=['sequence'])","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:51:34.459839Z","iopub.execute_input":"2023-04-30T23:51:34.460154Z","iopub.status.idle":"2023-04-30T23:51:35.809158Z","shell.execute_reply.started":"2023-04-30T23:51:34.460107Z","shell.execute_reply":"2023-04-30T23:51:35.808297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df6)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:51:48.888583Z","iopub.execute_input":"2023-04-30T23:51:48.888888Z","iopub.status.idle":"2023-04-30T23:51:48.929981Z","shell.execute_reply.started":"2023-04-30T23:51:48.888828Z","shell.execute_reply":"2023-04-30T23:51:48.929299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%d is the number of rows in the joined dataset' %df6.shape[0])","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:56:08.26156Z","iopub.execute_input":"2023-04-30T23:56:08.261861Z","iopub.status.idle":"2023-04-30T23:56:08.266044Z","shell.execute_reply.started":"2023-04-30T23:56:08.261821Z","shell.execute_reply":"2023-04-30T23:56:08.265457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look at classification type counts\ncounts = df6.classification.value_counts()\nprint(counts)\n\n#plot counts\nplt.figure()\nsns.distplot(counts, hist = False, color = 'purple')\nplt.title('Count Distribution for Family Types')\nplt.ylabel('% of records')\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:57:30.193252Z","iopub.execute_input":"2023-04-30T23:57:30.19366Z","iopub.status.idle":"2023-04-30T23:57:30.656006Z","shell.execute_reply.started":"2023-04-30T23:57:30.193614Z","shell.execute_reply":"2023-04-30T23:57:30.655202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get classification types where counts are under 1000\ntypes = np.asarray(counts[(counts < 1000)].index)\n\n# Filter dataset's records for classification types < 1000\ndata = df6[df6.classification.isin(types)]\n\nprint(types)\nprint('%d is the number of records in the final filtered dataset' %data.shape[0])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-01T00:05:28.517771Z","iopub.execute_input":"2023-05-01T00:05:28.518441Z","iopub.status.idle":"2023-05-01T00:05:28.639857Z","shell.execute_reply.started":"2023-05-01T00:05:28.518386Z","shell.execute_reply":"2023-05-01T00:05:28.638289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 3). ----- Train Test Split -----\n\n# Split Data\nX_train, X_test,y_train,y_test = train_test_split(data['sequence'], data['classification'], test_size = 0.2, random_state = 1)\n\n# Create a Count Vectorizer to gather the unique elements in sequence\nvect = CountVectorizer(analyzer = 'char_wb', ngram_range = (4,4))\n\n# Fit and Transform CountVectorizer\nvect.fit(X_train)\nX_train_df = vect.transform(X_train)\nX_test_df = vect.transform(X_test)\n\n#Print a few of the features\nprint(vect.get_feature_names()[-20:])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-01T00:06:10.538727Z","iopub.execute_input":"2023-05-01T00:06:10.539301Z","iopub.status.idle":"2023-05-01T00:07:25.4121Z","shell.execute_reply.started":"2023-05-01T00:06:10.539246Z","shell.execute_reply":"2023-05-01T00:07:25.411232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 4). ------ Machine Learning Models ------\n\n# Make a prediction dictionary to store accuracys\nprediction = dict()\n\n# Naive Bayes Model\nfrom sklearn.naive_bayes import MultinomialNB\nmodel = MultinomialNB()\nmodel.fit(X_train_df, y_train)\nNB_pred = model.predict(X_test_df)\nprediction[\"MultinomialNB\"] = accuracy_score(NB_pred, y_test)\nprint( prediction['MultinomialNB'])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-01T00:08:02.911995Z","iopub.execute_input":"2023-05-01T00:08:02.912285Z","iopub.status.idle":"2023-05-01T00:08:47.948072Z","shell.execute_reply.started":"2023-05-01T00:08:02.912237Z","shell.execute_reply":"2023-05-01T00:08:47.946528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 5). ----- Plot Confusion Matrix for NB -----\n\n# Plot confusion matrix\nconf_mat = confusion_matrix(y_test, NB_pred, labels = types)\n\n#Normalize confusion_matrix\nconf_mat = conf_mat.astype('float')/ conf_mat.sum(axis=1)[:, np.newaxis]\n\n# Plot Heat Map\nfig , ax = plt.subplots()\nfig.set_size_inches(13, 8)\nsns.heatmap(conf_mat)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-01T00:10:35.429852Z","iopub.execute_input":"2023-05-01T00:10:35.430204Z","iopub.status.idle":"2023-05-01T00:10:37.266467Z","shell.execute_reply.started":"2023-05-01T00:10:35.430156Z","shell.execute_reply":"2023-05-01T00:10:37.265671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(types[24])\nprint(types[84])\n","metadata":{"execution":{"iopub.status.busy":"2023-05-01T00:33:09.343771Z","iopub.execute_input":"2023-05-01T00:33:09.344181Z","iopub.status.idle":"2023-05-01T00:33:09.349981Z","shell.execute_reply.started":"2023-05-01T00:33:09.344107Z","shell.execute_reply":"2023-05-01T00:33:09.349148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\nThe confusion matrix shows label index 12 and 84 being misclassified quite a bit, these are DNA_REPAIR and ribosomal proetins","metadata":{"_kg_hide-input":true}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"#Print F1 score metrics\nprint(classification_report(y_test, NB_pred, target_names = types))\n","metadata":{"execution":{"iopub.status.busy":"2023-05-01T00:12:04.162325Z","iopub.execute_input":"2023-05-01T00:12:04.162684Z","iopub.status.idle":"2023-05-01T00:12:04.399873Z","shell.execute_reply.started":"2023-05-01T00:12:04.162626Z","shell.execute_reply":"2023-05-01T00:12:04.399058Z"},"trusted":true},"execution_count":null,"outputs":[]}]}