{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":87793,"databundleVersionId":11403143,"isSourceIdPinned":false,"sourceType":"competition"},{"sourceId":227203561,"sourceType":"kernelVersion"}],"dockerImageVersionId":30918,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Finding torsion angles from the distribution with multiple peaks","metadata":{}},{"cell_type":"markdown","source":"The purpose of this notebook is to find torsion angles, when there are multiple peaks in the distribution.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport datetime\nfrom scipy import signal\nfrom numpy import cross, eye, dot\nfrom scipy.linalg import expm, norm\nfrom scipy.signal import find_peaks\nimport scipy.stats as stats\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport warnings\nwarnings.simplefilter('ignore')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-14T15:49:44.302955Z","iopub.execute_input":"2025-03-14T15:49:44.303277Z","iopub.status.idle":"2025-03-14T15:49:44.309071Z","shell.execute_reply.started":"2025-03-14T15:49:44.303252Z","shell.execute_reply":"2025-03-14T15:49:44.307851Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"'n_1' is the torsion angle column for the 'res4' units.","metadata":{}},{"cell_type":"code","source":"train_labels = pd.read_csv('/kaggle/input/matrix-lmn/train_labels_perfect.csv')\ntrain_labels[['ID', 'resname', 'res4', 'n_1']]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-14T15:49:44.310274Z","iopub.execute_input":"2025-03-14T15:49:44.310701Z","iopub.status.idle":"2025-03-14T15:49:44.915799Z","shell.execute_reply.started":"2025-03-14T15:49:44.310667Z","shell.execute_reply":"2025-03-14T15:49:44.914668Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"residue_quartet = ['GGGU', 'GGUG', 'GUGC', 'UGCU', 'GCUC', 'CUCA', 'UCAG', 'CAGU', 'AGUA', 'GUAC', 'UACG', 'ACGA', 'CGAG', 'GAGA', 'AGAG', 'GAGG', 'AGGA', 'GGAA', 'GAAC', 'AACC', 'ACCG', 'CCGC', 'CGCA', 'GCAC', 'CACC', 'ACCC', 'GGCG', 'GCGC', 'GCAG', 'AGUG', 'GUGG', 'UGGG', 'GGGC', 'GGCU', 'GCUA', 'CUAG', 'UAGC', 'AGCG', 'CGCC', 'GCCA', 'CCAC', 'CACU', 'ACUC', 'UCAA', 'CAAA', 'AAAA', 'AAAG', 'AAGG', 'AGGC', 'GGCC', 'GCCC', 'CCCA', 'CCAU', 'GGGA', 'GGAC', 'GACU', 'ACUG', 'CUGA', 'UGAC', 'GACG', 'CGAU', 'GAUC', 'AUCA', 'UCAC', 'CACG', 'ACGC', 'AGUC', 'GUCU', 'UCUA', 'CUAU', 'GGAU', 'GAUA', 'AUAA', 'UAAC', 'AACU', 'ACUU', 'CUUC', 'UUCG', 'UCGG', 'CGGU', 'GGUU', 'GUUG', 'UUGU', 'UGUC', 'GUCC', 'UCCC', 'GCGA', 'CGAC', 'GACC', 'CCCU', 'CCUG', 'UGAU', 'GAUG', 'AUGA', 'UGAG', 'GCCG', 'CCGA', 'CGAA', 'GAAA', 'AAAC', 'CCGU', 'CGCU', 'GCUU', 'CUUG', 'UUGC', 'UGCG', 'GCGU', 'CGUC', 'CUCG', 'UCGU', 'CGUA', 'GUAA', 'UAAG', 'AAGA', 'GAGU', 'GUCA', 'ACCA', 'AAGC', 'AGCC', 'CCCG', 'UUAC', 'UACC', 'CCAA', 'CAAG', 'AAGU', 'AGUU', 'GUUU', 'UUUG', 'UUGA', 'AGGU', 'GGUA', 'CGUG', 'GUGU', 'UGUA', 'GUAG', 'AGCU', 'UCAU', 'CAUU', 'AUUA', 'UUAG', 'CUCC', 'UCCG', 'GAGC', 'GGCA', 'CAGA', 'AGAU', 'AUCU', 'UCUG', 'GCCU', 'CUGG', 'GGAG', 'CUCU', 'UCUC', 'CUGC', 'UGCC', 'GCAA', 'GGUC', 'CAGC', 'GCUG', 'ACGG', 'UACA', 'ACAG', 'CAGG', 'GGGG', 'UCUU', 'CGGA', 'UCCA', 'UGUG', 'GUGA', 'UGAA', 'AACA', 'ACAC', 'CGGC', 'GCGG', 'UGGA', 'UACU', 'AGAA', 'CUGU', 'UGUU', 'GUUC', 'UUCC', 'CCAG', 'AGAC', 'GACA', 'ACCU', 'CCUC', 'UCCU', 'UCGC', 'CGCG', 'CCUA', 'CUAA', 'GUUA', 'UUAU', 'UAUG', 'AUGG', 'UGGC', 'UUCA', 'CAAC', 'UUGG', 'GAAG', 'ACGU', 'CGUU', 'UUUC', 'CCUU', 'CGGG', 'ACAU', 'AUUG', 'UGCA', 'ACAA', 'CCCC', 'CUUU', 'UUUU', 'AGGG', 'CAUC', 'AUCG', 'UGGU', 'UAGU', 'CAUG', 'AUGC', 'UAGG', 'GUCG', 'UCGA', 'CCGG', 'AUAU', 'UAUC', 'ACUA', 'GUAU', 'UAUU', 'CAUA', 'AUAC', 'AUAG', 'AGCA', 'AUGU', 'UAAA', 'AAAU', 'AAUC', 'UAUA', 'GAUU', 'AUUC', 'AUCC', 'AACG', 'CAAU', 'AAUG', 'UUUA', 'UUAA', 'UAAU', 'CACA', 'UAGA', 'GCAU', 'UUCU', 'GAAU', 'CUUA', 'AAUA', 'CUAC', 'AAUU', 'AUUU']\n\n# for i in range(len(residue_quartet)):\ni = 8\ndf = train_labels[train_labels.res4==residue_quartet[i]]\nplt.figure(figsize=(10, 5))\nsns.histplot(df['n_1'], bins=50, kde=True)\nplt.xlabel(\"Torsion angle [degree]\")\nplt.ylabel(\"Count\")\nplt.title(f\"Torsion Angle Distribution for {residue_quartet[i]}\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-14T15:49:44.924828Z","iopub.execute_input":"2025-03-14T15:49:44.925187Z","iopub.status.idle":"2025-03-14T15:49:45.363696Z","shell.execute_reply.started":"2025-03-14T15:49:44.925160Z","shell.execute_reply":"2025-03-14T15:49:45.362518Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"https://qiita.com/shokishimada/items/f630a20099e8e4bdc2f7","metadata":{}},{"cell_type":"code","source":"#sampleはKDEの対象データ\ndef kde(x,sample,band_width,kernel):\n    n=len(sample)\n    return np.sum([1/(n*band_width)*kernel((x-sample[i])/band_width) for i in range(n)])\n\n#正規分布の定義\ndef normal(x,mean,sigma):\n    return 1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x-mean)/sigma)**2)\n\n#ガウスカーネル関数の定義\ndef normal_kernel(x):\n    return 1/np.sqrt(2*np.pi)*np.exp(-x**2/2)\n\ndef peaks(df):\n    sample = list(df)     # df['n_1']\n    sample_size = len(sample)\n\n    x = np.linspace(np.min(sample),np.max(sample),360)\n    band_width = np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)\n    y = [kde(x[i],sample,band_width,normal_kernel) for i in range(len(x))]\n    peaks, _ = find_peaks(y)\n\n    peak_x_list = [x[i] for i in peaks[::-1]]\n    peak_y_list = [y[i] for i in peaks[::-1]]\n    ratio = peak_y_list/sum(peak_y_list)\n\n    # print(peak_x_list)\n    # print(ratio)\n    return peak_x_list, list(ratio)\n\nprint(f'torsion angles [degree] for the {residue_quartet[i]} residue unit')\nprint(residue_quartet[i], peaks(df['n_1']))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-14T16:01:53.297812Z","iopub.execute_input":"2025-03-14T16:01:53.298334Z","iopub.status.idle":"2025-03-14T16:01:53.699544Z","shell.execute_reply.started":"2025-03-14T16:01:53.298295Z","shell.execute_reply":"2025-03-14T16:01:53.698530Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"sample = list(df['n_1'])\nsample_size = len(sample)\n\nfig,ax=plt.subplots(nrows=1,figsize=(15,7))\n\nax1=ax#[0]\nx = np.linspace(np.min(sample),np.max(sample),360)\nband_width = np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)\ny = [kde(x[i],sample,band_width,normal_kernel) for i in range(len(x))]\nax1.plot(x,y,label='kde')\n\n# ax2.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2),label='source')\nax1.legend(fontsize=15)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-14T15:49:45.773452Z","iopub.execute_input":"2025-03-14T15:49:45.773860Z","iopub.status.idle":"2025-03-14T15:49:46.433111Z","shell.execute_reply.started":"2025-03-14T15:49:45.773825Z","shell.execute_reply":"2025-03-14T15:49:46.431900Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}