{ "cells": [ { "cell_type": "code", "execution_count": 2, "id": "402442c7", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import h5py\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import numpy as np\n", "pd.options.mode.chained_assignment = None\n", "\n", "from sklearn.cluster import Birch\n", "from sklearn.cluster import DBSCAN\n", "from sklearn.cluster import OPTICS\n", "from sklearn import mixture\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "8727f5da", "metadata": {}, "outputs": [], "source": [ "# function for work\n", "\n", "def local_max(data, eps, n_rep):\n", " data_sort = data.sort_values('intensity', ascending=False)\n", " data_sort['peak'] = 1\n", " data_sort['region'] = 0\n", " \n", " n_reg = 0\n", " \n", " for n in range(n_rep):\n", " change = 0\n", " print(f\"repetition: {n}\")\n", " for idx in data_sort[data_sort['peak']==1].index:\n", " row = data_sort.loc[idx]\n", " if row['peak'] == 1:\n", " reg = data_sort[(data_sort['rt'] > (row['rt']-eps[0])) & (data_sort['rt'] < (row['rt']+eps[0])) &\n", " (data_sort['mz'] > (row['mz']-eps[1])) & (data_sort['mz'] < (row['mz']+eps[1]))]\n", " if len(reg) > 20:\n", " reg_max = reg['intensity'].idxmax()\n", " data_sort['peak'].loc[reg.index] = 0\n", " data_sort['peak'].loc[reg_max] = 1\n", " \n", " peak_reg = data_sort['region'].loc[reg_max]\n", " if peak_reg == 0:\n", " n_reg = n_reg + 1 \n", " peak_reg = n_reg\n", " \n", " data_sort['region'].loc[reg.index] = peak_reg\n", " \n", " if idx != reg_max:\n", " change = 1\n", " else:\n", " row['peak'] = 0\n", " print(f'found {n_reg} regions')\n", " if change == 0:\n", " print('-----')\n", " print('no change')\n", " print(f'found {n_reg} regions')\n", " return data_sort\n", " #data_sort = data_sort.sort_values(['intensity', 'peak'], ascending=False)\n", " print('-----')\n", " print(f'found {n_reg} regions')\n", " return data_sort\n", "\n", "\n", "def clustering_Birch(data, offset=0):\n", " clusters = Birch(threshold=0.01, n_clusters=None).fit_predict(data[['rtax','mz']])\n", " clusters[clusters>=0] += offset\n", " \n", " return clusters\n", "\n", "\n", "def clustering_DBScan(data, offset=0):\n", " clusters = DBSCAN(eps=0.01, min_samples=20, n_jobs=6).fit_predict(reg[['rtax','mz']])\n", " clusters[clusters>=0] += offset\n", " \n", " return clusters\n", "\n", "\n", "def clustering_Optics(data, offset=0):\n", " clusters = Optics(eps=0.01, min_samples=20, n_jobs=6).fit_predict(reg[['rtax','mz']])\n", " clusters[clusters>=0] += offset\n", " \n", " return clusters\n", "\n", "\n", "def clustering_GMM(data, n, offset=0):\n", " clusters = mixture.BayesianGaussianMixture(n_components=n, covariance_type=\"full\").fit_predict(reg[['rtax','mz']])\n", " clusters[clusters>=0] += offset\n", " \n", " return clusters" ] }, { "cell_type": "code", "execution_count": 4, "id": "b7dc766c", "metadata": {}, "outputs": [], "source": [ "with pd.HDFStore('NEWSTD_500.h5') as store:\n", " raw = store['spectra']\n", "raw['rt'] = raw['rt']*3600\n", "\n", "d_mean = list()\n", "d_median = list()\n", "for i in np.random.choice(raw['index'].unique(), 100):\n", " roi = raw[raw['index']==i]\n", " vector = np.diff(roi['mz'].sort_values())\n", " d_mean.append(np.mean(np.diff(roi['mz'].sort_values())))\n", " d_median.append(np.median(np.diff(roi['mz'].sort_values())))\n", " \n", "data = raw\n", "data = data[data[\"intensity\"] != 0]\n", "data['rtax'] = data['index']*np.median(d_median)\n", "data.rename({'index': 'scan'}, axis=1, inplace=1)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "85c82bdf", "metadata": {}, "outputs": [], "source": [ "y, x, _ = plt.hist(np.log(data['intensity']), bins=100)\n", "t = np.exp(x[np.where(y == y.max())])[0]\n", "data = data[data['intensity'] > t]\n", "data.describe()" ] }, { "cell_type": "code", "execution_count": null, "id": "bafff7d9", "metadata": {}, "outputs": [], "source": [ "# Space division\n", "n_bins = (10, 20)\n", "bins_overlap =(1.2, 1.2)\n", "eps = (10,1)\n", "rep = 10\n", "\n", "eps2 = (10, 0.01)\n", "rep2 = 5\n", "\n", "rt_bins, rt_step = np.linspace(data['rt'].min(), data['rt'].max(), n_bins[0]+1, retstep=True)\n", "mz_bins, mz_step = np.linspace(data['mz'].min(), data['mz'].max(), n_bins[1]+1, retstep=True)\n", "\n", "rt_step = rt_step * bins_overlap[0]\n", "mz_step = mz_step * bins_overlap[1]\n", "\n", "data['clusters_Birch'] = 0\n", "data['clusters_DBScan'] = 0\n", "data['clusters_Optics'] = 0\n", "data['clusters_GMM'] = 0\n", "\n", "for rt_start in rt_bins[:-1]:\n", " rt_bin = (rt_start, rt_start+rt_step)\n", " for mz_start in mz_bins[:-1]:\n", " mz_bin = (mz_start, mz_start+mz_step)\n", " \n", " roi = data[(data['rt']>rt_bin[0]) & (data['rt']mz_bin[0]) & (data['mz']