{"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":"# The notebook is about comparison of Lasso top features and correlation results","metadata":{}},{"cell_type":"markdown","source":"## Path definition & Imports","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 numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\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-01-23T20:05:58.968169Z","iopub.execute_input":"2023-01-23T20:05:58.968542Z","iopub.status.idle":"2023-01-23T20:05:59.008264Z","shell.execute_reply.started":"2023-01-23T20:05:58.968472Z","shell.execute_reply":"2023-01-23T20:05:59.007606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\n\nl = [] \ni = 0\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    l = l + [ dirname ] \n    for filename in filenames:\n        i += 1\n        if i<=10:\n            print('One of the first 10 files: ', os.path.join(dirname, filename))\nprint()\nprint('Dirs:', set(l))        \nprint('\\n n_files:', i)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:05:59.009727Z","iopub.execute_input":"2023-01-23T20:05:59.010147Z","iopub.status.idle":"2023-01-23T20:05:59.018610Z","shell.execute_reply.started":"2023-01-23T20:05:59.010121Z","shell.execute_reply":"2023-01-23T20:05:59.017571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:05:59.020186Z","iopub.execute_input":"2023-01-23T20:05:59.020735Z","iopub.status.idle":"2023-01-23T20:05:59.914500Z","shell.execute_reply.started":"2023-01-23T20:05:59.020700Z","shell.execute_reply":"2023-01-23T20:05:59.913804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dataset engineering","metadata":{}},{"cell_type":"code","source":"df_lasso = pd.read_csv('/kaggle/input/cd44-lasso/CD44_model_on_full_sample_Lasso  (1).csv')\ndf_lasso.rename(columns = {'Unnamed: 0': 'Gene'}, inplace=True)\ndisplay(df_lasso.head(10))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:05:59.917102Z","iopub.execute_input":"2023-01-23T20:05:59.917478Z","iopub.status.idle":"2023-01-23T20:05:59.957554Z","shell.execute_reply.started":"2023-01-23T20:05:59.917445Z","shell.execute_reply":"2023-01-23T20:05:59.956744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"RNA = pd.read_csv('../input/citeseqgse148127/GSE148127_SCT.normalized.RNA.counts.csv', index_col=0)\ndisplay(RNA.head(3))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:05:59.958596Z","iopub.execute_input":"2023-01-23T20:05:59.959048Z","iopub.status.idle":"2023-01-23T20:06:49.446179Z","shell.execute_reply.started":"2023-01-23T20:05:59.959021Z","shell.execute_reply":"2023-01-23T20:06:49.445175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"count = pd.read_csv('/kaggle/input/citeseqgse148127/GSE148127_ADT.counts.csv', index_col=0)\ndisplay(count.head(3))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:06:49.447374Z","iopub.execute_input":"2023-01-23T20:06:49.447708Z","iopub.status.idle":"2023-01-23T20:06:49.858386Z","shell.execute_reply.started":"2023-01-23T20:06:49.447683Z","shell.execute_reply":"2023-01-23T20:06:49.857154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_44 = count.loc['CITE-CD44']\ndisplay(df_44)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:06:49.859606Z","iopub.execute_input":"2023-01-23T20:06:49.860495Z","iopub.status.idle":"2023-01-23T20:06:49.868347Z","shell.execute_reply.started":"2023-01-23T20:06:49.860467Z","shell.execute_reply":"2023-01-23T20:06:49.867607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Source: https://www.kaggle.com/code/yoshifumimiya/simple-models-for-cds-all-proteins-gse148127/notebook?scriptVersionId=114395502\ntable = pd.read_csv('/kaggle/input/gse148127-gene-rna-list/GSE148127_gene_rna_list.csv')\ngene_id = table[table[\"In RNA dataset\"]=='Yes'][\"ID in RNA dataset\"]\nantibody_id = table[table[\"In RNA dataset\"]=='Yes'][\"Description\"]\ndict_list = dict(zip(gene_id, antibody_id))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:06:49.869277Z","iopub.execute_input":"2023-01-23T20:06:49.869949Z","iopub.status.idle":"2023-01-23T20:06:49.887072Z","shell.execute_reply.started":"2023-01-23T20:06:49.869922Z","shell.execute_reply":"2023-01-23T20:06:49.886032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dict_list","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:06:49.887926Z","iopub.execute_input":"2023-01-23T20:06:49.888152Z","iopub.status.idle":"2023-01-23T20:06:49.896374Z","shell.execute_reply.started":"2023-01-23T20:06:49.888130Z","shell.execute_reply":"2023-01-23T20:06:49.895802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### First try on dataset","metadata":{}},{"cell_type":"code","source":"# Spearman and Pearson correlation for each pair\nnames = []\ncorrelations = []\npvalues = []\ncorr_pear = []\nfor index, row in RNA.iterrows():\n    names.append(index)\n    \n    # Pearson\n    r = np.corrcoef(row, df_44)\n    corr_pear.append(r[0,1])\n    \n    #Spearman\n    rho, pval = stats.spearmanr(row, df_44)\n    correlations.append(rho)\n    pvalues.append(pval)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:06:49.899082Z","iopub.execute_input":"2023-01-23T20:06:49.899700Z","iopub.status.idle":"2023-01-23T20:07:28.993877Z","shell.execute_reply.started":"2023-01-23T20:06:49.899672Z","shell.execute_reply":"2023-01-23T20:07:28.992915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = {'Gene':names, 'Spearman correlation':correlations, 'Pval': pvalues, \n     'Pearson correlation': corr_pear}\ndf_stat = pd.DataFrame(d)\ndisplay(df_stat.head(5))\ndf_stat_sort = df_stat.sort_values(by=['Spearman correlation'], ascending=False)\ndisplay(df_stat_sort.head(5))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:28.994993Z","iopub.execute_input":"2023-01-23T20:07:28.995233Z","iopub.status.idle":"2023-01-23T20:07:29.028798Z","shell.execute_reply.started":"2023-01-23T20:07:28.995210Z","shell.execute_reply":"2023-01-23T20:07:29.027568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat.to_csv('Statistics_cd44.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:29.030189Z","iopub.execute_input":"2023-01-23T20:07:29.031136Z","iopub.status.idle":"2023-01-23T20:07:29.118825Z","shell.execute_reply.started":"2023-01-23T20:07:29.031096Z","shell.execute_reply":"2023-01-23T20:07:29.117798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Analysis","metadata":{}},{"cell_type":"code","source":"df_spearm = df_stat_sort.head(100)\ndf_pears = df_stat.sort_values(by=['Pearson correlation'], ascending=False).head(100)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:29.121924Z","iopub.execute_input":"2023-01-23T20:07:29.122178Z","iopub.status.idle":"2023-01-23T20:07:29.130972Z","shell.execute_reply.started":"2023-01-23T20:07:29.122150Z","shell.execute_reply":"2023-01-23T20:07:29.130067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(df_spearm.head(3))\ndisplay(df_pears.head(3))\ndisplay(df_lasso.head(3))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:29.132034Z","iopub.execute_input":"2023-01-23T20:07:29.132279Z","iopub.status.idle":"2023-01-23T20:07:29.155790Z","shell.execute_reply.started":"2023-01-23T20:07:29.132255Z","shell.execute_reply":"2023-01-23T20:07:29.154815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Spearman, Pearson, Lasso:', df_spearm.shape, df_pears.shape, df_lasso.shape)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:29.157268Z","iopub.execute_input":"2023-01-23T20:07:29.157895Z","iopub.status.idle":"2023-01-23T20:07:29.163633Z","shell.execute_reply.started":"2023-01-23T20:07:29.157849Z","shell.execute_reply":"2023-01-23T20:07:29.162492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"spear = set(df_spearm.head(100)['Gene'])\npears = set(df_pears.head(100)['Gene'])\nlass = set(df_lasso.head(100)['Gene'])\n\n# Spearman & Pearson\nsp_corr = spear & pears\n\n# Lasso & Spearman\nls_corr = lass & spear\n\n# Lasso & Pearson\nlp_corr = lass & pears\n\nprint('S&P:', len(sp_corr), '\\nL&S:', len(ls_corr), '\\nL&P:', len(lp_corr))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:29.164669Z","iopub.execute_input":"2023-01-23T20:07:29.165195Z","iopub.status.idle":"2023-01-23T20:07:29.175581Z","shell.execute_reply.started":"2023-01-23T20:07:29.165168Z","shell.execute_reply":"2023-01-23T20:07:29.174778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('S&P', sp_corr)\nprint('\\n\\n')\nprint('L^S', ls_corr)\nprint('\\n\\n')\nprint('L&P', lp_corr)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:07:29.176549Z","iopub.execute_input":"2023-01-23T20:07:29.176854Z","iopub.status.idle":"2023-01-23T20:07:29.187709Z","shell.execute_reply.started":"2023-01-23T20:07:29.176823Z","shell.execute_reply":"2023-01-23T20:07:29.186761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualisation","metadata":{}},{"cell_type":"code","source":"df_sp = df_stat_sort.dropna()\ndf_ps = df_stat.dropna().sort_values(by=['Pearson correlation'], ascending=False)\nprint(df_sp.shape, df_ps.shape)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:18:07.989218Z","iopub.execute_input":"2023-01-23T20:18:07.990078Z","iopub.status.idle":"2023-01-23T20:18:08.011145Z","shell.execute_reply.started":"2023-01-23T20:18:07.990027Z","shell.execute_reply":"2023-01-23T20:18:08.009923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def permutation_feature_importance(all_intersections):\n    fig = plt.figure(figsize = (20,6))\n    plt.plot(all_intersections ,label = 'original data - counts intersection' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:50],y[:50],1)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-50 points' )\n\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n\n\n    fig = plt.figure(figsize = (20,10))\n    plt.plot(all_intersections ,label = 'original data - counts intersection' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:50],y[:50],1)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'linear approximation based on 0-50 points' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:150],y[:150],1)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-150 points' )\n\n    plt.xlim([0,500])\n    plt.ylim([0,250])\n\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()\n    \n    fig = plt.figure(figsize = (20,10))\n    plt.plot(all_intersections ,label = 'original data - counts intersection' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x,y,2)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'parabolic approximation' )\n\n    y = np.array(all_intersections)\n    x = np.arange(len(y))\n    p = np.polyfit(x[:250],y[:250],1)\n    print(p)\n    plt.plot(np.polyval(p,x),label = 'linear approximation  based on 0-250 points' )\n\n    plt.xlim([0,3000])\n    plt.ylim([0,1500])\n\n    plt.legend(fontsize = 20 )\n    plt.grid()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:24:00.004186Z","iopub.execute_input":"2023-01-23T20:24:00.004529Z","iopub.status.idle":"2023-01-23T20:24:00.018399Z","shell.execute_reply.started":"2023-01-23T20:24:00.004503Z","shell.execute_reply":"2023-01-23T20:24:00.017773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Spearman & Pearson intersection function","metadata":{}},{"cell_type":"code","source":"# Spearman & Pearson\nall_intersections = []\nfor k in range(df_sp.shape[0]):\n    intersection = set(df_sp.head(k)['Gene']) & set(df_ps.head(k)['Gene'])\n    all_intersections.append(len(intersection))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:47:31.191753Z","iopub.execute_input":"2023-01-23T20:47:31.192073Z","iopub.status.idle":"2023-01-23T20:48:23.825535Z","shell.execute_reply.started":"2023-01-23T20:47:31.192048Z","shell.execute_reply":"2023-01-23T20:48:23.824488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"permutation_feature_importance(all_intersections)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:48:23.827785Z","iopub.execute_input":"2023-01-23T20:48:23.828170Z","iopub.status.idle":"2023-01-23T20:48:31.078106Z","shell.execute_reply.started":"2023-01-23T20:48:23.828136Z","shell.execute_reply":"2023-01-23T20:48:31.077312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,6))\nplt.title(\"Normalized intersection top K\")\nnormalised_intersections = [x / (idx+1) for idx, x in enumerate(all_intersections)]\nplt.plot(normalised_intersections)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:48:31.079320Z","iopub.execute_input":"2023-01-23T20:48:31.079560Z","iopub.status.idle":"2023-01-23T20:48:31.262392Z","shell.execute_reply.started":"2023-01-23T20:48:31.079537Z","shell.execute_reply":"2023-01-23T20:48:31.261321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Spearman & Lasso","metadata":{}},{"cell_type":"code","source":"all_intersections = []\nfor k in range(df_sp.shape[0]):\n    intersection = set(df_sp.head(k)['Gene']) & set(df_lasso.head(k)['Gene'])\n    all_intersections.append(len(intersection))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:48:31.264737Z","iopub.execute_input":"2023-01-23T20:48:31.265419Z","iopub.status.idle":"2023-01-23T20:49:25.317328Z","shell.execute_reply.started":"2023-01-23T20:48:31.265390Z","shell.execute_reply":"2023-01-23T20:49:25.316350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"permutation_feature_importance(all_intersections)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:49:25.319591Z","iopub.execute_input":"2023-01-23T20:49:25.320122Z","iopub.status.idle":"2023-01-23T20:49:26.030284Z","shell.execute_reply.started":"2023-01-23T20:49:25.320087Z","shell.execute_reply":"2023-01-23T20:49:26.029341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,6))\nplt.title(\"Normalized intersection top K\")\nnormalised_intersections = [x / (idx+1) for idx, x in enumerate(all_intersections)]\nplt.plot(normalised_intersections)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:49:26.031402Z","iopub.execute_input":"2023-01-23T20:49:26.031666Z","iopub.status.idle":"2023-01-23T20:49:26.218801Z","shell.execute_reply.started":"2023-01-23T20:49:26.031629Z","shell.execute_reply":"2023-01-23T20:49:26.217973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Pearson & Lasso","metadata":{}},{"cell_type":"code","source":"all_intersections = []\nfor k in range(df_ps.shape[0]):\n    intersection = set(df_ps.head(k)['Gene']) & set(df_lasso.head(k)['Gene'])\n    all_intersections.append(len(intersection))","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:49:26.219894Z","iopub.execute_input":"2023-01-23T20:49:26.220483Z","iopub.status.idle":"2023-01-23T20:50:20.534669Z","shell.execute_reply.started":"2023-01-23T20:49:26.220451Z","shell.execute_reply":"2023-01-23T20:50:20.533777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"permutation_feature_importance(all_intersections)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:50:20.536447Z","iopub.execute_input":"2023-01-23T20:50:20.537391Z","iopub.status.idle":"2023-01-23T20:50:21.260360Z","shell.execute_reply.started":"2023-01-23T20:50:20.537343Z","shell.execute_reply":"2023-01-23T20:50:21.259758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize = (20,6))\nplt.title(\"Normalized intersection top K\")\nnormalised_intersections = [x / (idx+1) for idx, x in enumerate(all_intersections)]\nplt.plot(normalised_intersections)","metadata":{"execution":{"iopub.status.busy":"2023-01-23T20:50:21.261440Z","iopub.execute_input":"2023-01-23T20:50:21.262459Z","iopub.status.idle":"2023-01-23T20:50:21.453313Z","shell.execute_reply.started":"2023-01-23T20:50:21.262406Z","shell.execute_reply":"2023-01-23T20:50:21.452210Z"},"trusted":true},"execution_count":null,"outputs":[]}]}