{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"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":59094,"databundleVersionId":7010844,"sourceType":"competition"}],"dockerImageVersionId":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"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\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport scipy.cluster.hierarchy as sch\nfrom scipy.cluster.hierarchy import fcluster\nfrom scipy.spatial import distance\nfrom scipy.cluster import hierarchy\n\n\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-11-20T17:20:29.600250Z","iopub.execute_input":"2023-11-20T17:20:29.600672Z","iopub.status.idle":"2023-11-20T17:20:30.702413Z","shell.execute_reply.started":"2023-11-20T17:20:29.600637Z","shell.execute_reply":"2023-11-20T17:20:30.701315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\ndf_de_train = pd.read_parquet(fn)# , index_col = 0)\nprint(df_de_train.shape)\ndf_de_train","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:20:30.704936Z","iopub.execute_input":"2023-11-20T17:20:30.705554Z","iopub.status.idle":"2023-11-20T17:20:33.773090Z","shell.execute_reply.started":"2023-11-20T17:20:30.705512Z","shell.execute_reply":"2023-11-20T17:20:33.771906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#1: if in list, 0: if not in list\n\npredetermined_list = ['Alvocidib', 'Belinostat', 'Foretinib', 'LDN 193189', 'Linagliptin', 'O-Demethylated Adapa',\n 'Dabrafenib', 'Dactolisib', 'Idelalisib', 'MLN 2238', 'Palbociclib', 'Porcn Inhibitor III',\n'CHIR-99021', 'Crizotinib', 'Oprozomib (ONX 0912)', 'Penfluridol',  'R428']","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:20:33.774615Z","iopub.execute_input":"2023-11-20T17:20:33.774979Z","iopub.status.idle":"2023-11-20T17:20:33.780823Z","shell.execute_reply.started":"2023-11-20T17:20:33.774950Z","shell.execute_reply":"2023-11-20T17:20:33.779625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"public_ids = [\n    'LSM-43216', 'LSM-1050', 'LSM-45849', 'LSM-42800', 'LSM-1131', 'LSM-6335', 'LSM-1211',\n    'LSM-45239', 'LSM-1130', 'LSM-45786', 'LSM-5199', 'LSM-45281',\n    'LSM-6324', # 'ACY-1215' -> 'Ricolinostat'\n    'LSM-3309', 'LSM-1056', 'LSM-45591', 'LSM-46203', 'LSM-5662',\n    'LSM-47134',  # 'SB-2342' -> '5-(9-Isopropyl-8-methyl-2-morpholino-9H-purin-6-yl)pyrimidin-2-amine  '\n    'LSM-45637', 'LSM-1127', 'LSM-46971', 'LSM-1172', 'LSM-46042', 'LSM-1101', 'LSM-45758',\n    'LSM-5218', 'LSM-2287', 'LSM-1014',\n    'LSM-1040', #  'fostamatinib' -> 'Tamatinib'\n    'LSM-1476;LSM-5290',\n    'LSM-45680',  # 'basimglurant' -> 'RG7090'\n    'LSM-4349',  # '5-iodotubercidin' -> 'IN1451'\n    'LSM-3425', 'LSM-45806',\n    'LSM-45616',  # 'SB-683698' -> 'TR-14035'\n    'LSM-1055',\n    'LSM-43281',  # 'C-646' -> 'STK219801'\n    'LSM-5690', 'LSM-1155', 'LSM-2499',\n    'LSM-2382',  # 'JTC-801' -> 'UNII-BXU45ZH6LI'\n    'LSM-45220', 'LSM-1037', 'LSM-1005', 'LSM-1180', 'LSM-36812',\n    'LSM-45924',  # 'filgotinib' -> 'GLPG0634'\n    'LSM-2013',  # 'TL-HRAS-61' -> TL_HRAS26'\n    'LSM-4738'\n]","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:20:33.784853Z","iopub.execute_input":"2023-11-20T17:20:33.785838Z","iopub.status.idle":"2023-11-20T17:20:33.794603Z","shell.execute_reply.started":"2023-11-20T17:20:33.785795Z","shell.execute_reply":"2023-11-20T17:20:33.793475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_de_train_with_ids = df_de_train[df_de_train[\"sm_lincs_id\"].isin(public_ids)]","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:20:33.796287Z","iopub.execute_input":"2023-11-20T17:20:33.796611Z","iopub.status.idle":"2023-11-20T17:20:33.836532Z","shell.execute_reply.started":"2023-11-20T17:20:33.796584Z","shell.execute_reply":"2023-11-20T17:20:33.835314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_de_train_without_ids = df_de_train[~df_de_train[\"sm_lincs_id\"].isin(public_ids)]","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:20:33.838689Z","iopub.execute_input":"2023-11-20T17:20:33.839007Z","iopub.status.idle":"2023-11-20T17:20:33.883471Z","shell.execute_reply.started":"2023-11-20T17:20:33.838979Z","shell.execute_reply":"2023-11-20T17:20:33.882265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# df_de_train","metadata":{}},{"cell_type":"code","source":"N = df_de_train.shape[1] # 5000\nprint(N)\nX = df_de_train[['sm_name'] + list(df_de_train.columns[5:N])].groupby('sm_name').median()\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\nl = list(X.index)\nl = [t[:20] for t in l] # cut long names\nl_dict = {i: 1 if i in predetermined_list else 0 for i in l}\nl_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\ncm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\nprint(cm.shape)\nclustergrid = sns.clustermap(cm,cmap=\"coolwarm\")\nplt.show()\n\nreordered_columns = clustergrid.dendrogram_col.reordered_ind\nreordered_rows = clustergrid.dendrogram_row.reordered_ind\nprint(len(reordered_rows), len(reordered_columns))\nprint(list(cm.index[reordered_rows]))\n\nsns.clustermap(cm, annot=True, fmt=\".2f\", cmap=\"coolwarm\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:20:33.884952Z","iopub.execute_input":"2023-11-20T17:20:33.885371Z","iopub.status.idle":"2023-11-20T17:21:22.978923Z","shell.execute_reply.started":"2023-11-20T17:20:33.885338Z","shell.execute_reply":"2023-11-20T17:21:22.977759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i0,corr_method in enumerate(['pearson']):\n    N = df_de_train.shape[1] # 5000\n    print(N)\n    X = df_de_train[['sm_name'] + list(df_de_train.columns[5:N])].groupby('sm_name').median()\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    l = list(X.index)\n    l = [t[:20] for t in l] # cut long names\n    l_dict = {i: 1 if i in predetermined_list else 0 for i in l}\n    l_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\n    cm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\n    v = cm.values[np.triu_indices(len(cm),k=1)]\n    print(len(v), len(cm)*(len(cm)-1)/2 ) # Check two numbers coincide - number of upper diagonal elements\n    plt.hist(v, bins = 1000)#  , color='k')\n    plt.title(corr_method)\n    plt.show()\n    sr = pd.Series(v).describe()\n    sr.name = corr_method\n    sr = sr.to_frame()\n    \n    display(sr )    \n    if i0 == 0:\n        df_methods_stat = sr\n    else:\n        df_methods_stat = df_methods_stat.join(sr)\n        \n        \n    threshold4top_corr = 0.4\n    list_columns_ids = list( cm.columns)\n    z =    np.where( np.triu(cm.abs().values,1) >= threshold4top_corr )\n    df_top_corrs = pd.DataFrame()\n    IX = 0\n    for (i,j) in zip( z[0], z[1]):\n        if i>=j : continue \n        #print(i,j, list_genes_ids[i], list_genes_ids[j], cm.iloc[i,j])\n        IX += 1\n#         df_top_corrs.loc[IX, 'I1'] = i\n#         df_top_corrs.loc[IX, 'I2'] = j\n        df_top_corrs.loc[IX, 'Id1'] = list_columns_ids[i]\n        df_top_corrs.loc[IX, 'Id2'] = list_columns_ids[j]\n        df_top_corrs.loc[IX, 'Correlation'] = cm.iloc[i,j]\n\n#     df_top_corrs['I1'] = df_top_corrs['I1'].astype(int) \n#     df_top_corrs['I2'] = df_top_corrs['I2'].astype(int) \n\n    print(df_top_corrs.shape)\n\n    display ( df_top_corrs.sort_values('Correlation' , ascending = False, key = abs ).head(5) )\n\n    n_top = 20\n    if i0 == 0:\n        df_top_corrs_join = df_top_corrs.sort_values('Correlation' , ascending = False, key = abs).head(n_top)\n        df_top_corrs_join.index = range(len(df_top_corrs_join))\n    else:\n        df_tmp = df_top_corrs.sort_values('Correlation' , ascending = False, key = abs ).head(n_top)\n        df_tmp.index = range(len(df_tmp))\n        df_top_corrs_join = df_top_corrs_join.join( df_tmp, rsuffix = '_'+corr_method) \n    \ndisplay( df_methods_stat        )\ndisplay( df_top_corrs_join )","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:22.980363Z","iopub.execute_input":"2023-11-20T17:21:22.981120Z","iopub.status.idle":"2023-11-20T17:21:25.749797Z","shell.execute_reply.started":"2023-11-20T17:21:22.981081Z","shell.execute_reply":"2023-11-20T17:21:25.748661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters_list = []\n\nfor corr_method in ['pearson']:\n    print(corr_method)\n    fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\n    df_de_train = pd.read_parquet(fn)# , index_col = 0)\n    \n    N = df_de_train.shape[1] # 5000\n    print(N)\n    X = df_de_train[['sm_name'] + list(df_de_train.columns[5:N])].groupby('sm_name').median()\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    l = list(X.index)\n    l = [t[:20] for t in l] # cut long names\n    l_dict = {i: 1 if i in predetermined_list else 0 for i in l}\n    l_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\n    cm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\n\n        \n    correlations = np.abs(cm)\n    correlations_array = np.abs(cm)\n\n    m_linkage = hierarchy.linkage(\n        distance.pdist(correlations_array), method='average')  \n    \n    \n\n    clusters = hierarchy.fcluster(m_linkage, 0.5*distance.pdist(correlations_array).max(), 'distance')\n    \n        \n    index_clusters = []    \n    val, cnt = np.unique(clusters, return_counts=True)\n    \n    for i,cluster in enumerate(clusters):        \n        clusters_list.append([cm.index[i], cluster, corr_method])\n        \n        index_cluster = f\"{cm.index[i]}_cl{cluster}_sz{int(cnt[val == cluster])}\"\n        index_clusters.append(index_cluster)\n    \n    \n    \n    index_dict = dict(zip(cm.index, index_clusters))\n    \n       \n    cm.index = cm.index.map(index_dict.get)    \n    \n    \n    sns.clustermap(correlations, row_linkage=m_linkage, col_linkage=m_linkage, method=\"average\",cmap='vlag', figsize=(25, 25), \n                   xticklabels=cm.index, yticklabels=cm.index).fig.suptitle(f'{corr_method}', fontsize=20, ha='center'); plt.show()  \n","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:25.751128Z","iopub.execute_input":"2023-11-20T17:21:25.751551Z","iopub.status.idle":"2023-11-20T17:21:31.545613Z","shell.execute_reply.started":"2023-11-20T17:21:25.751522Z","shell.execute_reply":"2023-11-20T17:21:31.544512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_clusters = pd.DataFrame(clusters_list, columns=['cd', 'cluster', 'corr_method'])","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.550045Z","iopub.execute_input":"2023-11-20T17:21:31.550505Z","iopub.status.idle":"2023-11-20T17:21:31.559556Z","shell.execute_reply.started":"2023-11-20T17:21:31.550473Z","shell.execute_reply":"2023-11-20T17:21:31.558531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_clusters.shape","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.560786Z","iopub.execute_input":"2023-11-20T17:21:31.561572Z","iopub.status.idle":"2023-11-20T17:21:31.576428Z","shell.execute_reply.started":"2023-11-20T17:21:31.561541Z","shell.execute_reply":"2023-11-20T17:21:31.575348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_clusters.head()","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.577916Z","iopub.execute_input":"2023-11-20T17:21:31.578255Z","iopub.status.idle":"2023-11-20T17:21:31.591246Z","shell.execute_reply.started":"2023-11-20T17:21:31.578215Z","shell.execute_reply":"2023-11-20T17:21:31.590220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_clusters = df_clusters[df_clusters['corr_method'] == 'pearson']","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.592624Z","iopub.execute_input":"2023-11-20T17:21:31.592970Z","iopub.status.idle":"2023-11-20T17:21:31.599909Z","shell.execute_reply.started":"2023-11-20T17:21:31.592941Z","shell.execute_reply":"2023-11-20T17:21:31.598544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters_df_list = [pearson_clusters]","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.601710Z","iopub.execute_input":"2023-11-20T17:21:31.602072Z","iopub.status.idle":"2023-11-20T17:21:31.610649Z","shell.execute_reply.started":"2023-11-20T17:21:31.602044Z","shell.execute_reply":"2023-11-20T17:21:31.609559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for s in clusters_df_list:\n    min_value = len(s['cluster'].value_counts())","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.611595Z","iopub.execute_input":"2023-11-20T17:21:31.611898Z","iopub.status.idle":"2023-11-20T17:21:31.623582Z","shell.execute_reply.started":"2023-11-20T17:21:31.611872Z","shell.execute_reply":"2023-11-20T17:21:31.622508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_v_counts = pearson_clusters['cluster'].value_counts().head(min_value)\n\npearson_df = pearson_clusters.merge(pearson_v_counts.to_frame(),\n                                left_on='cluster',\n                                right_index=True)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.624739Z","iopub.execute_input":"2023-11-20T17:21:31.625074Z","iopub.status.idle":"2023-11-20T17:21:31.648955Z","shell.execute_reply.started":"2023-11-20T17:21:31.625045Z","shell.execute_reply":"2023-11-20T17:21:31.648099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_group = pearson_df.groupby('cluster')\n\npearson_df2 = pearson_group.apply(lambda x: x['cd'].unique())\ns = pearson_df2.str.len().sort_values(ascending=False).index\npearson_df2 = pearson_df2.reindex(s)\npearson = pearson_df2.reset_index(drop=True)\n\nfor x in pearson:\n    print('число элементов в кластере: ', len(x), ', кластер: ', x)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.650245Z","iopub.execute_input":"2023-11-20T17:21:31.650770Z","iopub.status.idle":"2023-11-20T17:21:31.666246Z","shell.execute_reply.started":"2023-11-20T17:21:31.650733Z","shell.execute_reply":"2023-11-20T17:21:31.665349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# df_de_train_with_ids","metadata":{}},{"cell_type":"code","source":"N = df_de_train_with_ids.shape[1] # 5000\nprint(N)\nX = df_de_train_with_ids[['sm_name'] + list(df_de_train_with_ids.columns[5:N])].groupby('sm_name').median()\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\nl = list(X.index)\nl = [t[:20] for t in l] # cut long names\nl_dict = {i: 1 if i in predetermined_list else 0 for i in l}\nl_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\ncm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\nprint(cm.shape)\nclustergrid = sns.clustermap(cm,cmap=\"coolwarm\")\nplt.show()\n\nreordered_columns = clustergrid.dendrogram_col.reordered_ind\nreordered_rows = clustergrid.dendrogram_row.reordered_ind\nprint(len(reordered_rows), len(reordered_columns))\nprint(list(cm.index[reordered_rows]))\n\nsns.clustermap(cm, annot=True, fmt=\".2f\", cmap=\"coolwarm\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:31.667288Z","iopub.execute_input":"2023-11-20T17:21:31.668171Z","iopub.status.idle":"2023-11-20T17:21:38.717702Z","shell.execute_reply.started":"2023-11-20T17:21:31.668125Z","shell.execute_reply":"2023-11-20T17:21:38.716380Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i0,corr_method in enumerate(['pearson']):\n    N = df_de_train_with_ids.shape[1] # 5000\n    print(N)\n    X = df_de_train_with_ids[['sm_name'] + list(df_de_train_with_ids.columns[5:N])].groupby('sm_name').median()\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    l = list(X.index)\n    l = [t[:20] for t in l] # cut long names\n    l_dict = {i: 1 if i in predetermined_list else 0 for i in l}\n    l_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\n    cm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\n    v = cm.values[np.triu_indices(len(cm),k=1)]\n    print(len(v), len(cm)*(len(cm)-1)/2 ) # Check two numbers coincide - number of upper diagonal elements\n    plt.hist(v, bins = 1000)#  , color='k')\n    plt.title(corr_method)\n    plt.show()\n    sr = pd.Series(v).describe()\n    sr.name = corr_method\n    sr = sr.to_frame()\n    \n    display(sr )    \n    if i0 == 0:\n        df_methods_stat = sr\n    else:\n        df_methods_stat = df_methods_stat.join(sr)\n        \n        \n    threshold4top_corr = 0.4\n    list_columns_ids = list( cm.columns)\n    z =    np.where( np.triu(cm.abs().values,1) >= threshold4top_corr )\n    df_top_corrs = pd.DataFrame()\n    IX = 0\n    for (i,j) in zip( z[0], z[1]):\n        if i>=j : continue \n        #print(i,j, list_genes_ids[i], list_genes_ids[j], cm.iloc[i,j])\n        IX += 1\n#         df_top_corrs.loc[IX, 'I1'] = i\n#         df_top_corrs.loc[IX, 'I2'] = j\n        df_top_corrs.loc[IX, 'Id1'] = list_columns_ids[i]\n        df_top_corrs.loc[IX, 'Id2'] = list_columns_ids[j]\n        df_top_corrs.loc[IX, 'Correlation'] = cm.iloc[i,j]\n\n#     df_top_corrs['I1'] = df_top_corrs['I1'].astype(int) \n#     df_top_corrs['I2'] = df_top_corrs['I2'].astype(int) \n\n    print(df_top_corrs.shape)\n\n    display ( df_top_corrs.sort_values('Correlation' , ascending = False, key = abs ).head(5) )\n\n    n_top = 20\n    if i0 == 0:\n        df_top_corrs_join = df_top_corrs.sort_values('Correlation' , ascending = False, key = abs).head(n_top)\n        df_top_corrs_join.index = range(len(df_top_corrs_join))\n    else:\n        df_tmp = df_top_corrs.sort_values('Correlation' , ascending = False, key = abs ).head(n_top)\n        df_tmp.index = range(len(df_tmp))\n        df_top_corrs_join = df_top_corrs_join.join( df_tmp, rsuffix = '_'+corr_method) \n    \ndisplay( df_methods_stat        )\ndisplay( df_top_corrs_join )","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:38.719573Z","iopub.execute_input":"2023-11-20T17:21:38.720002Z","iopub.status.idle":"2023-11-20T17:21:40.689555Z","shell.execute_reply.started":"2023-11-20T17:21:38.719965Z","shell.execute_reply":"2023-11-20T17:21:40.688433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters_list = []\n\nfor corr_method in ['pearson']:\n    print(corr_method)\n    fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\n    df_de_train = pd.read_parquet(fn)# , index_col = 0)\n    df_de_train_with_ids = df_de_train[df_de_train[\"sm_lincs_id\"].isin(public_ids)]\n    N = df_de_train_with_ids.shape[1] # 5000\n    print(N)\n    X = df_de_train_with_ids[['sm_name'] + list(df_de_train_with_ids.columns[5:N])].groupby('sm_name').median()\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    l = list(X.index)\n    l = [t[:20] for t in l] # cut long names\n    l_dict = {i: 1 if i in predetermined_list else 0 for i in l}\n    l_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\n    cm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\n\n        \n    correlations = np.abs(cm)\n    correlations_array = np.abs(cm)\n\n    m_linkage = hierarchy.linkage(\n        distance.pdist(correlations_array), method='average')  \n    \n    \n\n    clusters = hierarchy.fcluster(m_linkage, 0.5*distance.pdist(correlations_array).max(), 'distance')\n    \n        \n    index_clusters = []    \n    val, cnt = np.unique(clusters, return_counts=True)\n    \n    for i,cluster in enumerate(clusters):        \n        clusters_list.append([cm.index[i], cluster, corr_method])\n        \n        index_cluster = f\"{cm.index[i]}_cl{cluster}_sz{int(cnt[val == cluster])}\"\n        index_clusters.append(index_cluster)\n    \n    \n    \n    index_dict = dict(zip(cm.index, index_clusters))\n    \n       \n    cm.index = cm.index.map(index_dict.get)    \n    \n    \n    sns.clustermap(correlations, row_linkage=m_linkage, col_linkage=m_linkage, method=\"average\",cmap='vlag', figsize=(25, 25), \n                   xticklabels=cm.index, yticklabels=cm.index).fig.suptitle(f'{corr_method}', fontsize=20, ha='center'); plt.show()  \n","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:40.690858Z","iopub.execute_input":"2023-11-20T17:21:40.691312Z","iopub.status.idle":"2023-11-20T17:21:45.375409Z","shell.execute_reply.started":"2023-11-20T17:21:40.691272Z","shell.execute_reply":"2023-11-20T17:21:45.374292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_clusters = pd.DataFrame(clusters_list, columns=['cd', 'cluster', 'corr_method'])","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.376658Z","iopub.execute_input":"2023-11-20T17:21:45.376977Z","iopub.status.idle":"2023-11-20T17:21:45.383103Z","shell.execute_reply.started":"2023-11-20T17:21:45.376948Z","shell.execute_reply":"2023-11-20T17:21:45.382358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_clusters = df_clusters[df_clusters['corr_method'] == 'pearson']","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.384375Z","iopub.execute_input":"2023-11-20T17:21:45.384740Z","iopub.status.idle":"2023-11-20T17:21:45.394443Z","shell.execute_reply.started":"2023-11-20T17:21:45.384711Z","shell.execute_reply":"2023-11-20T17:21:45.393548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters_df_list = [pearson_clusters]","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.395636Z","iopub.execute_input":"2023-11-20T17:21:45.396360Z","iopub.status.idle":"2023-11-20T17:21:45.406367Z","shell.execute_reply.started":"2023-11-20T17:21:45.396331Z","shell.execute_reply":"2023-11-20T17:21:45.405504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for s in clusters_df_list:\n    min_value = len(s['cluster'].value_counts())","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.407397Z","iopub.execute_input":"2023-11-20T17:21:45.408323Z","iopub.status.idle":"2023-11-20T17:21:45.422877Z","shell.execute_reply.started":"2023-11-20T17:21:45.408293Z","shell.execute_reply":"2023-11-20T17:21:45.422093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_v_counts = pearson_clusters['cluster'].value_counts().head(min_value)\n\npearson_df = pearson_clusters.merge(pearson_v_counts.to_frame(),\n                                left_on='cluster',\n                                right_index=True)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.424013Z","iopub.execute_input":"2023-11-20T17:21:45.424891Z","iopub.status.idle":"2023-11-20T17:21:45.437814Z","shell.execute_reply.started":"2023-11-20T17:21:45.424861Z","shell.execute_reply":"2023-11-20T17:21:45.436851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_group = pearson_df.groupby('cluster')\n\npearson_df2 = pearson_group.apply(lambda x: x['cd'].unique())\ns = pearson_df2.str.len().sort_values(ascending=False).index\npearson_df2 = pearson_df2.reindex(s)\npearson = pearson_df2.reset_index(drop=True)\n\nfor x in pearson:\n    print('число элементов в кластере: ', len(x), ', кластер: ', x)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.439242Z","iopub.execute_input":"2023-11-20T17:21:45.440230Z","iopub.status.idle":"2023-11-20T17:21:45.454612Z","shell.execute_reply.started":"2023-11-20T17:21:45.440190Z","shell.execute_reply":"2023-11-20T17:21:45.453293Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# df_de_train_without_ids","metadata":{}},{"cell_type":"code","source":"N = df_de_train_without_ids.shape[1] # 5000\nprint(N)\nX = df_de_train_without_ids[['sm_name'] + list(df_de_train_without_ids.columns[5:N])].groupby('sm_name').median()\nprint(X.shape)\ncm = np.corrcoef(X)\nprint(cm[:3,:2])\ncm = np.abs(cm)\nl = list(X.index)\nl = [t[:20] for t in l] # cut long names\nl_dict = {i: 1 if i in predetermined_list else 0 for i in l}\nl_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\ncm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\nprint(cm.shape)\nclustergrid = sns.clustermap(cm,cmap=\"coolwarm\")\nplt.show()\n\nreordered_columns = clustergrid.dendrogram_col.reordered_ind\nreordered_rows = clustergrid.dendrogram_row.reordered_ind\nprint(len(reordered_rows), len(reordered_columns))\nprint(list(cm.index[reordered_rows]))\n\nsns.clustermap(cm, annot=True, fmt=\".2f\", cmap=\"coolwarm\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:21:45.456290Z","iopub.execute_input":"2023-11-20T17:21:45.456634Z","iopub.status.idle":"2023-11-20T17:22:07.733068Z","shell.execute_reply.started":"2023-11-20T17:21:45.456586Z","shell.execute_reply":"2023-11-20T17:22:07.731760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i0,corr_method in enumerate(['pearson']):\n    N = df_de_train_without_ids.shape[1] # 5000\n    print(N)\n    X = df_de_train_without_ids[['sm_name'] + list(df_de_train_without_ids.columns[5:N])].groupby('sm_name').median()\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    l = list(X.index)\n    l = [t[:20] for t in l] # cut long names\n    l_dict = {i: 1 if i in predetermined_list else 0 for i in l}\n    l_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\n    cm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\n    v = cm.values[np.triu_indices(len(cm),k=1)]\n    print(len(v), len(cm)*(len(cm)-1)/2 ) # Check two numbers coincide - number of upper diagonal elements\n    plt.hist(v, bins = 1000)#  , color='k')\n    plt.title(corr_method)\n    plt.show()\n    sr = pd.Series(v).describe()\n    sr.name = corr_method\n    sr = sr.to_frame()\n    \n    display(sr )    \n    if i0 == 0:\n        df_methods_stat = sr\n    else:\n        df_methods_stat = df_methods_stat.join(sr)\n        \n        \n    threshold4top_corr = 0.4\n    list_columns_ids = list( cm.columns)\n    z =    np.where( np.triu(cm.abs().values,1) >= threshold4top_corr )\n    df_top_corrs = pd.DataFrame()\n    IX = 0\n    for (i,j) in zip( z[0], z[1]):\n        if i>=j : continue \n        #print(i,j, list_genes_ids[i], list_genes_ids[j], cm.iloc[i,j])\n        IX += 1\n#         df_top_corrs.loc[IX, 'I1'] = i\n#         df_top_corrs.loc[IX, 'I2'] = j\n        df_top_corrs.loc[IX, 'Id1'] = list_columns_ids[i]\n        df_top_corrs.loc[IX, 'Id2'] = list_columns_ids[j]\n        df_top_corrs.loc[IX, 'Correlation'] = cm.iloc[i,j]\n\n#     df_top_corrs['I1'] = df_top_corrs['I1'].astype(int) \n#     df_top_corrs['I2'] = df_top_corrs['I2'].astype(int) \n\n    print(df_top_corrs.shape)\n\n    display ( df_top_corrs.sort_values('Correlation' , ascending = False, key = abs ).head(5) )\n\n    n_top = 20\n    if i0 == 0:\n        df_top_corrs_join = df_top_corrs.sort_values('Correlation' , ascending = False, key = abs).head(n_top)\n        df_top_corrs_join.index = range(len(df_top_corrs_join))\n    else:\n        df_tmp = df_top_corrs.sort_values('Correlation' , ascending = False, key = abs ).head(n_top)\n        df_tmp.index = range(len(df_tmp))\n        df_top_corrs_join = df_top_corrs_join.join( df_tmp, rsuffix = '_'+corr_method) \n    \ndisplay( df_methods_stat        )\ndisplay( df_top_corrs_join )","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:07.740534Z","iopub.execute_input":"2023-11-20T17:22:07.740926Z","iopub.status.idle":"2023-11-20T17:22:10.478748Z","shell.execute_reply.started":"2023-11-20T17:22:07.740894Z","shell.execute_reply":"2023-11-20T17:22:10.477616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters_list = []\n\nfor corr_method in ['pearson']:\n    print(corr_method)\n    fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\n    df_de_train = pd.read_parquet(fn)# , index_col = 0)\n    df_de_train_without_ids = df_de_train[~df_de_train[\"sm_lincs_id\"].isin(public_ids)]\n    N = df_de_train_without_ids.shape[1] # 5000\n    print(N)\n    X = df_de_train_without_ids[['sm_name'] + list(df_de_train_without_ids.columns[5:N])].groupby('sm_name').median()\n    print(X.shape)\n    cm = np.corrcoef(X)\n    print(cm[:3,:2])\n    cm = np.abs(cm)\n    l = list(X.index)\n    l = [t[:20] for t in l] # cut long names\n    l_dict = {i: 1 if i in predetermined_list else 0 for i in l}\n    l_combined = [f\"{name} ({num})\" for name, num in l_dict.items()]\n    cm = pd.DataFrame(cm, index=l_combined, columns=l_combined)\n\n        \n    correlations = np.abs(cm)\n    correlations_array = np.abs(cm)\n\n    m_linkage = hierarchy.linkage(\n        distance.pdist(correlations_array), method='average')  \n    \n    \n\n    clusters = hierarchy.fcluster(m_linkage, 0.5*distance.pdist(correlations_array).max(), 'distance')\n    \n        \n    index_clusters = []    \n    val, cnt = np.unique(clusters, return_counts=True)\n    \n    for i,cluster in enumerate(clusters):        \n        clusters_list.append([cm.index[i], cluster, corr_method])\n        \n        index_cluster = f\"{cm.index[i]}_cl{cluster}_sz{int(cnt[val == cluster])}\"\n        index_clusters.append(index_cluster)\n    \n    \n    \n    index_dict = dict(zip(cm.index, index_clusters))\n    \n       \n    cm.index = cm.index.map(index_dict.get)    \n    \n    \n    sns.clustermap(correlations, row_linkage=m_linkage, col_linkage=m_linkage, method=\"average\",cmap='vlag', figsize=(25, 25), \n                   xticklabels=cm.index, yticklabels=cm.index).fig.suptitle(f'{corr_method}', fontsize=20, ha='center'); plt.show()  \n","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:10.480438Z","iopub.execute_input":"2023-11-20T17:22:10.480794Z","iopub.status.idle":"2023-11-20T17:22:14.762730Z","shell.execute_reply.started":"2023-11-20T17:22:10.480764Z","shell.execute_reply":"2023-11-20T17:22:14.761916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_clusters = pd.DataFrame(clusters_list, columns=['cd', 'cluster', 'corr_method'])","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:14.763976Z","iopub.execute_input":"2023-11-20T17:22:14.764819Z","iopub.status.idle":"2023-11-20T17:22:14.772055Z","shell.execute_reply.started":"2023-11-20T17:22:14.764787Z","shell.execute_reply":"2023-11-20T17:22:14.770829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_clusters = df_clusters[df_clusters['corr_method'] == 'pearson']","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:14.773686Z","iopub.execute_input":"2023-11-20T17:22:14.774334Z","iopub.status.idle":"2023-11-20T17:22:14.786476Z","shell.execute_reply.started":"2023-11-20T17:22:14.774297Z","shell.execute_reply":"2023-11-20T17:22:14.784991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters_df_list = [pearson_clusters]","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:14.788530Z","iopub.execute_input":"2023-11-20T17:22:14.789003Z","iopub.status.idle":"2023-11-20T17:22:14.802556Z","shell.execute_reply.started":"2023-11-20T17:22:14.788964Z","shell.execute_reply":"2023-11-20T17:22:14.801228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for s in clusters_df_list:\n    min_value = len(s['cluster'].value_counts())","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:14.803764Z","iopub.execute_input":"2023-11-20T17:22:14.804080Z","iopub.status.idle":"2023-11-20T17:22:14.814874Z","shell.execute_reply.started":"2023-11-20T17:22:14.804050Z","shell.execute_reply":"2023-11-20T17:22:14.813962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_v_counts = pearson_clusters['cluster'].value_counts().head(min_value)\n\npearson_df = pearson_clusters.merge(pearson_v_counts.to_frame(),\n                                left_on='cluster',\n                                right_index=True)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:14.816773Z","iopub.execute_input":"2023-11-20T17:22:14.817789Z","iopub.status.idle":"2023-11-20T17:22:14.829607Z","shell.execute_reply.started":"2023-11-20T17:22:14.817756Z","shell.execute_reply":"2023-11-20T17:22:14.828647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pearson_group = pearson_df.groupby('cluster')\n\npearson_df2 = pearson_group.apply(lambda x: x['cd'].unique())\ns = pearson_df2.str.len().sort_values(ascending=False).index\npearson_df2 = pearson_df2.reindex(s)\npearson = pearson_df2.reset_index(drop=True)\n\nfor x in pearson:\n    print('число элементов в кластере: ', len(x), ', кластер: ', x)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:22:14.831048Z","iopub.execute_input":"2023-11-20T17:22:14.831386Z","iopub.status.idle":"2023-11-20T17:22:14.846356Z","shell.execute_reply.started":"2023-11-20T17:22:14.831358Z","shell.execute_reply":"2023-11-20T17:22:14.845174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Top correlations with IDs from list - NK Cells","metadata":{}},{"cell_type":"code","source":"N = df_de_train.shape[1] # 5000\nX = df_de_train[['sm_name'] + list(df_de_train.columns[5:N])].groupby('sm_name').median()\ncm = np.corrcoef(X)\ncm = np.abs(cm)\nl = list(X.index)\nl = [t[:20] for t in l] # cut long names\ncm = pd.DataFrame(cm, index=l, columns=l)\n\n# Create a subset of df_de_train\nsubset = df_de_train[(df_de_train[\"cell_type\"] == \"NK cells\") | df_de_train[\"sm_name\"].isin(predetermined_list)]\n\n# Filter the index values of the cm dataframe\ncm_filtered = cm[cm.index.isin(subset['sm_name'])]\n\n# Create a new predetermined_list_filtered\npredetermined_list_filtered = [item for item in predetermined_list if item in cm_filtered.index]\n\n# Print missing items\nif len(predetermined_list_filtered) < len(predetermined_list):\n    missing_items = set(predetermined_list) - set(predetermined_list_filtered)\n    print(f\"The following items were not found in the correlation matrix and will be excluded: {missing_items}\")\n\n# Find the most correlated items\nfrom scipy.stats import pearsonr\n\n# Find the most correlated items\nmost_correlated = {}\nfor item in predetermined_list_filtered:\n    correlations = cm_filtered.loc[item, :]\n    correlations[item] = 0  # exclude self-correlation\n    most_correlated_item = correlations.idxmax()\n    correlation_coefficient, p_value = pearsonr(cm_filtered[item], cm_filtered[most_correlated_item])\n    correlation_coefficient = round(correlation_coefficient, 3)  # Round to 3 decimal places\n    most_correlated[item] = (most_correlated_item, correlation_coefficient, p_value)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:24:19.890209Z","iopub.execute_input":"2023-11-20T17:24:19.890623Z","iopub.status.idle":"2023-11-20T17:24:20.290553Z","shell.execute_reply.started":"2023-11-20T17:24:19.890589Z","shell.execute_reply":"2023-11-20T17:24:20.289466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Sort items by correlation in descending order\nmost_correlated_sorted = sorted(most_correlated.items(), key=lambda x: x[1][1], reverse=True)","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:24:22.898599Z","iopub.execute_input":"2023-11-20T17:24:22.898975Z","iopub.status.idle":"2023-11-20T17:24:22.904968Z","shell.execute_reply.started":"2023-11-20T17:24:22.898946Z","shell.execute_reply":"2023-11-20T17:24:22.903691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Print items with correlation coefficient and p-value\nfor item, (most_correlated_item, correlation, p_value) in most_correlated_sorted:\n    print(f\"Item: {item}, Most Correlated Item: {most_correlated_item}, Correlation: {correlation}, P-Value: {p_value}\")","metadata":{"execution":{"iopub.status.busy":"2023-11-20T17:24:24.998986Z","iopub.execute_input":"2023-11-20T17:24:25.000027Z","iopub.status.idle":"2023-11-20T17:24:25.007021Z","shell.execute_reply.started":"2023-11-20T17:24:24.999988Z","shell.execute_reply":"2023-11-20T17:24:25.005739Z"},"trusted":true},"execution_count":null,"outputs":[]}]}