{"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":"**In the previous notebook we transformed the CiteSeq data into a normal distribution.**\n\n[CiteSeq Target Analysis(Blom function)](https://www.kaggle.com/code/yoshifumimiya/citeseq-target-analysis-blom-function)\n\n**With normally distributed data, ANOVA can be used as feature selection.**","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom sklearn.feature_selection import SelectKBest\nfrom sklearn.feature_selection import f_classif","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:00.373533Z","iopub.execute_input":"2022-11-04T13:52:00.374462Z","iopub.status.idle":"2022-11-04T13:52:01.591544Z","shell.execute_reply.started":"2022-11-04T13:52:00.374343Z","shell.execute_reply":"2022-11-04T13:52:01.590215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"raw = pd.read_csv('../input/citeseq-blom/CiteSeq_blom.csv', index_col=0)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:01.593764Z","iopub.execute_input":"2022-11-04T13:52:01.594497Z","iopub.status.idle":"2022-11-04T13:52:06.435906Z","shell.execute_reply.started":"2022-11-04T13:52:01.594440Z","shell.execute_reply":"2022-11-04T13:52:06.434639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"raw.head()","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:06.437361Z","iopub.execute_input":"2022-11-04T13:52:06.438483Z","iopub.status.idle":"2022-11-04T13:52:06.486942Z","shell.execute_reply.started":"2022-11-04T13:52:06.438445Z","shell.execute_reply":"2022-11-04T13:52:06.485702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Identification of \"cell_type\" related proteins","metadata":{}},{"cell_type":"code","source":"columns = raw.columns[0:140]\nX = raw.iloc[:,0:140]\ny = raw[\"cell_type\"]","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:06.499339Z","iopub.execute_input":"2022-11-04T13:52:06.499806Z","iopub.status.idle":"2022-11-04T13:52:06.530602Z","shell.execute_reply.started":"2022-11-04T13:52:06.499765Z","shell.execute_reply":"2022-11-04T13:52:06.529562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selector = SelectKBest(f_classif, k=10) # k is the number of features to be selected\nX_new = selector.fit_transform(X, y)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:06.531934Z","iopub.execute_input":"2022-11-04T13:52:06.532484Z","iopub.status.idle":"2022-11-04T13:52:06.830893Z","shell.execute_reply.started":"2022-11-04T13:52:06.532452Z","shell.execute_reply":"2022-11-04T13:52:06.829843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_scores = list(zip(selector.scores_,columns))\nsorted_feature_scores = sorted(feature_scores,reverse=True)\n\nnum_list = []\ncol_list = []\nfor i in range(140):\n   num_list.append((sorted_feature_scores[i])[0])\n   col_list.append((sorted_feature_scores [i])[1])","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:06.832109Z","iopub.execute_input":"2022-11-04T13:52:06.832438Z","iopub.status.idle":"2022-11-04T13:52:06.839809Z","shell.execute_reply.started":"2022-11-04T13:52:06.832409Z","shell.execute_reply":"2022-11-04T13:52:06.838487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.bar(col_list[0:10],num_list[0:10])\nplt.xticks(rotation=90)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:06.841182Z","iopub.execute_input":"2022-11-04T13:52:06.841574Z","iopub.status.idle":"2022-11-04T13:52:07.120780Z","shell.execute_reply.started":"2022-11-04T13:52:06.841542Z","shell.execute_reply":"2022-11-04T13:52:07.119350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.stats import f\nimport matplotlib.pyplot as plt\nimport numpy as np","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:07.122110Z","iopub.execute_input":"2022-11-04T13:52:07.122523Z","iopub.status.idle":"2022-11-04T13:52:07.128210Z","shell.execute_reply.started":"2022-11-04T13:52:07.122488Z","shell.execute_reply":"2022-11-04T13:52:07.126963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfn = 8-1 # Inter-level degrees of freedom (dfA=k-1)\ndfd = 140*8-8 # Intra-level degrees of freedom (dfE=N-k) (N=n1＋n2＋… nk)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:07.132154Z","iopub.execute_input":"2022-11-04T13:52:07.132548Z","iopub.status.idle":"2022-11-04T13:52:07.141021Z","shell.execute_reply.started":"2022-11-04T13:52:07.132495Z","shell.execute_reply":"2022-11-04T13:52:07.139686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 1)\n\nplt.xlim(-1,10)\nplt.ylim(0,1)\nx = np.linspace(f.ppf(0.0000000001, dfn, dfd),f.ppf(0.9999999999, dfn, dfd), 100)\nax.plot(x, f.pdf(x, dfn, dfd), 'r-')\nax.axvline(f.ppf(0.95, dfn, dfd), ls = \"--\", color = \"navy\")\nprint('upper 5%:', f.ppf(0.95, dfn, dfd))","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:07.142677Z","iopub.execute_input":"2022-11-04T13:52:07.143143Z","iopub.status.idle":"2022-11-04T13:52:07.354983Z","shell.execute_reply.started":"2022-11-04T13:52:07.143092Z","shell.execute_reply":"2022-11-04T13:52:07.353698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame(num_list,index=col_list,columns=['importance'])\ndf","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:07.356169Z","iopub.execute_input":"2022-11-04T13:52:07.356507Z","iopub.status.idle":"2022-11-04T13:52:07.371255Z","shell.execute_reply.started":"2022-11-04T13:52:07.356478Z","shell.execute_reply":"2022-11-04T13:52:07.369929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DE = df[df[\"importance\"]>f.ppf(0.95, dfn, dfd)]\nDE","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:52:07.372949Z","iopub.execute_input":"2022-11-04T13:52:07.373663Z","iopub.status.idle":"2022-11-04T13:52:07.389109Z","shell.execute_reply.started":"2022-11-04T13:52:07.373599Z","shell.execute_reply":"2022-11-04T13:52:07.387884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the ANOVA feature selection, only \"TCRVa7.2\" was considered as noise.","metadata":{}},{"cell_type":"markdown","source":"# Identification of \"donor\" related proteins","metadata":{}},{"cell_type":"code","source":"raw[\"donor\"].unique()","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:07:04.850515Z","iopub.execute_input":"2022-11-04T14:07:04.851428Z","iopub.status.idle":"2022-11-04T14:07:04.861543Z","shell.execute_reply.started":"2022-11-04T14:07:04.851378Z","shell.execute_reply":"2022-11-04T14:07:04.860556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_2 = raw.iloc[:,0:140]\ny_2 = raw[\"donor\"]","metadata":{"execution":{"iopub.status.busy":"2022-11-04T13:57:20.342927Z","iopub.execute_input":"2022-11-04T13:57:20.343358Z","iopub.status.idle":"2022-11-04T13:57:20.386213Z","shell.execute_reply.started":"2022-11-04T13:57:20.343327Z","shell.execute_reply":"2022-11-04T13:57:20.384695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selector2 = SelectKBest(f_classif, k=10) # k is the number of features to be selected\nX_new_2 = selector2.fit_transform(X_2, y_2)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:01:22.431340Z","iopub.execute_input":"2022-11-04T14:01:22.431794Z","iopub.status.idle":"2022-11-04T14:01:22.632708Z","shell.execute_reply.started":"2022-11-04T14:01:22.431757Z","shell.execute_reply":"2022-11-04T14:01:22.631450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_scores_2 = list(zip(selector2.scores_,columns))\nsorted_feature_scores_2 = sorted(feature_scores_2,reverse=True)\n\nnum_list_2 = []\ncol_list_2 = []\nfor i in range(140):\n   num_list_2.append((sorted_feature_scores_2[i])[0])\n   col_list_2.append((sorted_feature_scores_2 [i])[1])","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:03:44.788604Z","iopub.execute_input":"2022-11-04T14:03:44.789081Z","iopub.status.idle":"2022-11-04T14:03:44.796859Z","shell.execute_reply.started":"2022-11-04T14:03:44.789045Z","shell.execute_reply":"2022-11-04T14:03:44.795741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfn_2 = 3-1 # Inter-level degrees of freedom (dfA=k-1)\ndfd_2 = 140*3-3 # Intra-level degrees of freedom (dfE=N-k) (N=n1＋n2＋… nk)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:05:21.999188Z","iopub.execute_input":"2022-11-04T14:05:21.999636Z","iopub.status.idle":"2022-11-04T14:05:22.005243Z","shell.execute_reply.started":"2022-11-04T14:05:21.999583Z","shell.execute_reply":"2022-11-04T14:05:22.004076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 1)\n\nplt.xlim(-1,10)\nplt.ylim(0,1)\nx = np.linspace(f.ppf(0.0000000001, dfn_2, dfd_2),f.ppf(0.9999999999, dfn_2, dfd_2), 100)\nax.plot(x, f.pdf(x, dfn_2, dfd_2), 'r-')\nax.axvline(f.ppf(0.95, dfn_2, dfd_2), ls = \"--\", color = \"navy\")\nprint('upper 5%:', f.ppf(0.95, dfn_2, dfd_2))","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:05:57.462980Z","iopub.execute_input":"2022-11-04T14:05:57.463422Z","iopub.status.idle":"2022-11-04T14:05:57.678249Z","shell.execute_reply.started":"2022-11-04T14:05:57.463385Z","shell.execute_reply":"2022-11-04T14:05:57.676993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_2 = pd.DataFrame(num_list_2,index=col_list_2,columns=['importance'])\ndf_2","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:03:45.720153Z","iopub.execute_input":"2022-11-04T14:03:45.721473Z","iopub.status.idle":"2022-11-04T14:03:45.736403Z","shell.execute_reply.started":"2022-11-04T14:03:45.721410Z","shell.execute_reply":"2022-11-04T14:03:45.734975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DE_2 = df_2[df_2[\"importance\"]>f.ppf(0.95, dfn_2, dfd_2)]\nDE_2","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:07:41.276400Z","iopub.execute_input":"2022-11-04T14:07:41.276859Z","iopub.status.idle":"2022-11-04T14:07:41.294553Z","shell.execute_reply.started":"2022-11-04T14:07:41.276823Z","shell.execute_reply":"2022-11-04T14:07:41.293049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All protein expression levels are needed to identify donors.","metadata":{}},{"cell_type":"markdown","source":"# Identification of \"day\" related proteins","metadata":{}},{"cell_type":"code","source":"raw[\"day\"].unique()","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:09:12.863777Z","iopub.execute_input":"2022-11-04T14:09:12.864207Z","iopub.status.idle":"2022-11-04T14:09:12.875796Z","shell.execute_reply.started":"2022-11-04T14:09:12.864175Z","shell.execute_reply":"2022-11-04T14:09:12.873906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_3 = raw.iloc[:,0:140]\ny_3 = raw[\"day\"]\n\nselector3 = SelectKBest(f_classif, k=10) # k is the number of features to be selected\nX_new_3 = selector3.fit_transform(X_3, y_3)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:10:01.385232Z","iopub.execute_input":"2022-11-04T14:10:01.385682Z","iopub.status.idle":"2022-11-04T14:10:01.668556Z","shell.execute_reply.started":"2022-11-04T14:10:01.385647Z","shell.execute_reply":"2022-11-04T14:10:01.666256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_scores_3 = list(zip(selector3.scores_,columns))\nsorted_feature_scores_3 = sorted(feature_scores_3,reverse=True)\n\nnum_list_3 = []\ncol_list_3 = []\nfor i in range(140):\n   num_list_3.append((sorted_feature_scores_3[i])[0])\n   col_list_3.append((sorted_feature_scores_3 [i])[1])","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:10:34.701102Z","iopub.execute_input":"2022-11-04T14:10:34.701534Z","iopub.status.idle":"2022-11-04T14:10:34.710404Z","shell.execute_reply.started":"2022-11-04T14:10:34.701499Z","shell.execute_reply":"2022-11-04T14:10:34.708544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfn_3 = 3-1 # Inter-level degrees of freedom (dfA=k-1)\ndfd_3 = 140*3-3 # Intra-level degrees of freedom (dfE=N-k) (N=n1＋n2＋… nk)","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:10:53.068768Z","iopub.execute_input":"2022-11-04T14:10:53.070172Z","iopub.status.idle":"2022-11-04T14:10:53.076214Z","shell.execute_reply.started":"2022-11-04T14:10:53.070113Z","shell.execute_reply":"2022-11-04T14:10:53.074655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 1)\n\nplt.xlim(-1,10)\nplt.ylim(0,1)\nx = np.linspace(f.ppf(0.0000000001, dfn_3, dfd_3),f.ppf(0.9999999999, dfn_3, dfd_3), 100)\nax.plot(x, f.pdf(x, dfn_3, dfd_3), 'r-')\nax.axvline(f.ppf(0.95, dfn_3, dfd_3), ls = \"--\", color = \"navy\")\nprint('upper 5%:', f.ppf(0.95, dfn_3, dfd_3))","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:11:23.418802Z","iopub.execute_input":"2022-11-04T14:11:23.419202Z","iopub.status.idle":"2022-11-04T14:11:23.628384Z","shell.execute_reply.started":"2022-11-04T14:11:23.419171Z","shell.execute_reply":"2022-11-04T14:11:23.627111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_3 = pd.DataFrame(num_list_3,index=col_list_3,columns=['importance'])\ndf_3","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:11:51.665821Z","iopub.execute_input":"2022-11-04T14:11:51.666894Z","iopub.status.idle":"2022-11-04T14:11:51.679734Z","shell.execute_reply.started":"2022-11-04T14:11:51.666854Z","shell.execute_reply":"2022-11-04T14:11:51.678906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DE_3 = df_3[df_3[\"importance\"]>f.ppf(0.95, dfn_3, dfd_3)]\nDE_3","metadata":{"execution":{"iopub.status.busy":"2022-11-04T14:12:20.692761Z","iopub.execute_input":"2022-11-04T14:12:20.693174Z","iopub.status.idle":"2022-11-04T14:12:20.707760Z","shell.execute_reply.started":"2022-11-04T14:12:20.693139Z","shell.execute_reply":"2022-11-04T14:12:20.706516Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All protein expression levels are needed to identify donors.","metadata":{}},{"cell_type":"markdown","source":"# Conclusion","metadata":{}},{"cell_type":"markdown","source":"**Although the analysis was based on only limited data from input, all proteins were likely to be necessary for the identification of \"donor\", \"day\", and \"cell_type\".**\n\n**However, only \"TCRVa7.2\" may not be necessary for \"cell_type\".**","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}