{"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":"# What is about ?\n\nWe take CD protein data. First we transpose it and make various dimensional reductions - thus we visualize PROTEINs (not cells!). Thus we might get an intuation what are the clusters among the proteins. We do not get very clear results, but probably we can see 2-3-4 clusters. \n\nThen we proceed with direct clustering the proteins. \n\nWe use clustring method: from scipy.cluster import hierarchy \n    \n    The key param here: \"PARAM\" is defined in the line:\n    clusters = hierarchy.fcluster(m_linkage, PARAM*distance.pdist(correlations_array).max(), 'distance')\n    For the PARAM = 0.3 we get many small clusters - in particular many 1 element clusters - that probably not what we want.\n    PARAM = 0.5 - gives more reasonable results - 6 clusters - see details below\n  \nComparison with clusters obtained by the other method:\n\n    It appears to be that the biggest (61 protein) cluster with low correlated proteins has been mainly splitted into 2: with about 20 , 30 elements in each.\n\n    The cluster with 40 elements has been mainly splitted into two with 23, 22 elements.\n\n    The cluster with 24 elemtns fits cluster with 17 elements.\n\n    Other clusters are small so we will not discuss them.\n\nThe clusters from here: https://www.kaggle.com/code/alekse1pakhomov/multimodal-singlecell-correlation?scriptVersionId=119682605&cellId=1\n\n\nPS \n\nSimilar analysis for the other dataset can be found in:\nhttps://www.kaggle.com/visualcomments/multi-modal-cite-seq-clusters\n","metadata":{}},{"cell_type":"markdown","source":"## Some outcomes on clustering \n\n### For PARAM = 0.5 we get 6 clusters\n\nnumber of clusters: 6\n\nThe biggest cluster - proteins NOT correlated with the others - so it is does not make sense we should not look at it.\n\nSome small clusters probably should be unified into one. \n\ncluster sizes: [3, 4, 8, 24, 40, 61]\nmean correlations: [0.68 0.6  0.51 0.39 0.32 0.11]\n\n\ncluster size:  40 , mean abs correlation per cluster:  0.322 , genes:  ['CD270', 'CD154', 'CD105', 'CD45RO', 'CD279', 'TIGIT', 'CD335', 'Podoplanin', 'CD103', 'CD152', 'CD107a', 'CD134', 'CD141', 'CD1d', 'CD314', 'CD57', 'CD272', 'CD278', 'CX3CR1', 'CD24', 'integrinB7', 'TCR', 'CD192', 'FceRIa', 'CD137', 'CD83', 'CD124', 'CD226', 'CD28', 'CD71', 'CD26', 'CD115', 'CD158', 'LOX-1', 'CD158b', 'CD142', 'CD85j', 'HLA-E', 'CD82', 'CD88']\n\ncluster size:  24 , mean abs correlation per cluster:  0.393 , genes:  ['CD112', 'CD47', 'CD48', 'CD52', 'HLA-A-B-C', 'CD45RA', 'CD123', 'CD44', 'CD31', 'CD69', 'CD62L', 'CD95', 'HLA-DR', 'CD11a', 'CD244', 'CD54', 'CD13', 'CD49b', 'CD81', 'CD18', 'CD45', 'CD22', 'CD72', 'CD9']\n\ncluster size:  8 , mean abs correlation per cluster:  0.507 , genes:  ['CD155', 'CD49f', 'KLRG1', 'CD58', 'CD119', 'CD29', 'CD63', 'CD224']\n\ncluster size:  4 , mean abs correlation per cluster:  0.6 , genes:  ['CD33', 'CD38', 'CD49d', 'CD162']\n\ncluster size:  3 , mean abs correlation per cluster:  0.681 , genes:  ['CD32', 'CD41', 'CD36']\n\n\nThe biggest cluster - proteins NOT correlated with the others - so it is does not make sense we should not look at it:\ncluster size:  61 , mean abs correlation per cluster:  0.106 , genes:  ['CD86', 'CD274', 'CD40', 'CD3', 'CD8', 'CD56', 'CD19', 'CD11c', 'CD7', 'CD194', 'CD4', 'CD14', 'CD16', 'CD25', 'Mouse-IgG1', 'Mouse-IgG2a', 'Mouse-IgG2b', 'Rat-IgG2b', 'CD20', 'CD146', 'IgM', 'CD5', 'CD195', 'CD196', 'CD185', 'CD161', 'CD223', 'CD27', 'CD1c', 'CD11b', 'CD64', 'CD35', 'CD39', 'CD21', 'CD79b', 'CD169', 'CD268', 'CD42b', 'CD62P', 'Rat-IgG1', 'Rat-IgG2a', 'CD122', 'CD163', 'CD2', 'CD303', 'IgD', 'CD127', 'CD304', 'CD172a', 'CD93', 'CD49a', 'CD73', 'TCRVa7.2', 'TCRVd2', 'CD158e1', 'CD319', 'CD352', 'CD94', 'CD23', 'CD328', 'CD101']\n","metadata":{}},{"cell_type":"markdown","source":"# Load and prepare data","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        if 'open-problems-multimodal' in dirname:\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-02-20T17:13:37.861359Z","iopub.execute_input":"2023-02-20T17:13:37.861889Z","iopub.status.idle":"2023-02-20T17:13:37.930482Z","shell.execute_reply.started":"2023-02-20T17:13:37.861788Z","shell.execute_reply":"2023-02-20T17:13:37.929490Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\nimport matplotlib.pyplot as plt\n","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:37.932258Z","iopub.execute_input":"2023-02-20T17:13:37.933434Z","iopub.status.idle":"2023-02-20T17:13:38.550807Z","shell.execute_reply.started":"2023-02-20T17:13:37.933384Z","shell.execute_reply":"2023-02-20T17:13:38.549946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf_protein = pd.read_hdf('/kaggle/input/open-problems-multimodal/train_cite_targets.h5')\ndf = df_protein\nprint(df.shape)\ndf","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:38.552488Z","iopub.execute_input":"2023-02-20T17:13:38.552922Z","iopub.status.idle":"2023-02-20T17:13:39.717475Z","shell.execute_reply.started":"2023-02-20T17:13:38.552879Z","shell.execute_reply":"2023-02-20T17:13:39.716396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v4color = df.std(axis = 0)\nv4color.name = 'std'\ndisplay( v4color.sort_values(ascending = False).head(20) )\ndisplay( v4color.describe() )\nplt.figure(figsize = (20,5))\nplt.plot(v4color.sort_values(ascending = False).values , '*-')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:39.719975Z","iopub.execute_input":"2023-02-20T17:13:39.720908Z","iopub.status.idle":"2023-02-20T17:13:40.210411Z","shell.execute_reply.started":"2023-02-20T17:13:39.720873Z","shell.execute_reply":"2023-02-20T17:13:40.209154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nX = scaler.fit_transform(df)\nX = X.T\nprint(X.shape)\nX","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:40.211706Z","iopub.execute_input":"2023-02-20T17:13:40.212077Z","iopub.status.idle":"2023-02-20T17:13:40.431649Z","shell.execute_reply.started":"2023-02-20T17:13:40.212029Z","shell.execute_reply":"2023-02-20T17:13:40.430654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA","metadata":{}},{"cell_type":"code","source":"%%time\nfrom sklearn.decomposition import PCA\nreducer = PCA(n_components=20)\nr = reducer.fit_transform(X)\n","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:40.432915Z","iopub.execute_input":"2023-02-20T17:13:40.433302Z","iopub.status.idle":"2023-02-20T17:13:41.208930Z","shell.execute_reply.started":"2023-02-20T17:13:40.433269Z","shell.execute_reply":"2023-02-20T17:13:41.207409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"palette1 = 'rainbow' #  ['blue','red']\nmarker1 = 'o'\nalpha1 = 1\nstr_data_inf = 'NIPS22 CDprots'\n#v4color = df.std(axis = 0)\n#v4color.name = 'std'\n\n\nfig = plt.figure(figsize = (20,10) ); c = 0\nax = sns.scatterplot(x = r[:,0],y=r[:,1], hue = v4color,  palette = palette1,  marker = marker1, alpha = alpha1 )\nif 1:\n    plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n    plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title        \nplt.title(' PCA '+   str_data_inf + ' n_samples='+str(len(r)),fontsize = 20 ) \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:41.210991Z","iopub.execute_input":"2023-02-20T17:13:41.211921Z","iopub.status.idle":"2023-02-20T17:13:41.688083Z","shell.execute_reply.started":"2023-02-20T17:13:41.211886Z","shell.execute_reply":"2023-02-20T17:13:41.687229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_x_subplots = 2\nc = 0\ncc = 0\nfor (i,j) in [(0,1),(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    cc += 1\n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        plt.suptitle(str_data_inf +  ' n_samples='+str(len(df)),fontsize = 20 ) \n\n\n    c += 1; fig.add_subplot(1,n_x_subplots ,c)\n\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1,  marker = marker1 , alpha = alpha1 )\n    if c == 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title    \n        plt.legend( bbox_to_anchor = (1.15, 1.), loc='upper right')\n    else:\n        plt.legend('')\n    plt.title('PCA '+str(i)+' '+str(j), fontsize = 20)\n    \n    \n    \nplt.show()    ","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:41.689390Z","iopub.execute_input":"2023-02-20T17:13:41.690311Z","iopub.status.idle":"2023-02-20T17:13:43.396742Z","shell.execute_reply.started":"2023-02-20T17:13:41.690274Z","shell.execute_reply":"2023-02-20T17:13:43.395210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# UMAP","metadata":{}},{"cell_type":"code","source":"%%time\nimport umap\nreducer = umap.UMAP()\nprint('UMAP params: min_dist, n_neighbors: ', reducer.min_dist, reducer.n_neighbors )\n\nr = reducer.fit_transform(X)","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:13:43.404826Z","iopub.execute_input":"2023-02-20T17:13:43.405390Z","iopub.status.idle":"2023-02-20T17:14:14.101281Z","shell.execute_reply.started":"2023-02-20T17:13:43.405340Z","shell.execute_reply":"2023-02-20T17:14:14.099975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# reducer = umap.UMAP()\n# r = reducer.fit_transform(df)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1, alpha = alpha1  )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=15) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=15) # for legend title        \n    plt.title( ' UMAP '+ str_data_inf +  ' n_samples='+str(len(X)), fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:14:14.106116Z","iopub.execute_input":"2023-02-20T17:14:14.106808Z","iopub.status.idle":"2023-02-20T17:14:14.622452Z","shell.execute_reply.started":"2023-02-20T17:14:14.106773Z","shell.execute_reply":"2023-02-20T17:14:14.621200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport umap\nmetric = 'cosine'\nreducer = umap.UMAP(metric = metric, random_state=42 )\nprint('UMAP params: min_dist, n_neighbors: ', reducer.min_dist, reducer.n_neighbors )\n\nr = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1, alpha = alpha1  )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=15) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=15) # for legend title        \n    plt.title( metric + ' UMAP '+ str_data_inf +  ' n_samples='+str(len(X)), fontsize = 20)\n    plt.show()\n    \n#%%time\nimport umap\nmetric = 'correlation'\nreducer = umap.UMAP(metric = metric , random_state=42)\nprint('UMAP params: min_dist, n_neighbors: ', reducer.min_dist, reducer.n_neighbors )\n\nr = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1, alpha = alpha1  )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=15) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=15) # for legend title        \n    plt.title( metric + ' UMAP '+ str_data_inf +  ' n_samples='+str(len(X)), fontsize = 20)\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:14:14.624397Z","iopub.execute_input":"2023-02-20T17:14:14.624756Z","iopub.status.idle":"2023-02-20T17:14:22.210620Z","shell.execute_reply.started":"2023-02-20T17:14:14.624723Z","shell.execute_reply":"2023-02-20T17:14:22.209524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# 8 min for (70988, 140) - six times \nfig = plt.figure(figsize = (20,16)); c = 0; cc = 0\nplt.suptitle('UMAP ' + str_data_inf+   ' n_samples='+str(len(X)) , fontsize = 20 )\nfor min_dist in [0.1, 0.9]: # 0.1 - defailt min_dist\n    for n_neighbors in [5,15, 100]: # 15 - default n_neighbors,\n        c += 1; fig.add_subplot(2,3,c)\n        str_inf = 'n_neighbors='+str(n_neighbors) + ' min_dist='+str(min_dist) \n        reducer = umap.UMAP(n_neighbors = n_neighbors, min_dist = min_dist, n_components= 2  )# random_state=42 # metric = \n        r2 = reducer.fit_transform(X)\n\n        i,j = 0,1\n        ax = sns.scatterplot(x = r2[:,i],y = r2[:,j], hue = v4color, palette = palette1 , marker =  marker1 , alpha = alpha1 )\n        if 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title        \n        plt.title(str_inf , fontsize = 20)\n        plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n        plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n        \n#         if cc>0:\n#             plt.legend('')\n#         cc+=1\n        if c == 1:\n            plt.setp(ax.get_legend().get_texts(), fontsize=10) # for legend text\n            plt.setp(ax.get_legend().get_title(), fontsize=10) # for legend title    \n            plt.legend( bbox_to_anchor = (1.15, 1.), loc='upper right')\n        else:\n            plt.legend('')\n\n\nplt.show() \n","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:14:22.211708Z","iopub.execute_input":"2023-02-20T17:14:22.212042Z","iopub.status.idle":"2023-02-20T17:14:43.722496Z","shell.execute_reply.started":"2023-02-20T17:14:22.211994Z","shell.execute_reply":"2023-02-20T17:14:43.721259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# NCVis - similar to TSNE, bust faster (from Skolkovo team)¶","metadata":{}},{"cell_type":"code","source":"!pip install ncvis\nimport ncvis\n\nreducer = ncvis.NCVis()\nr = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1 , alpha = alpha1  )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title(' ncvis ' + str_data_inf + ' n_samples='+str(len(r)) , fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:14:43.724227Z","iopub.execute_input":"2023-02-20T17:14:43.724635Z","iopub.status.idle":"2023-02-20T17:15:20.651619Z","shell.execute_reply.started":"2023-02-20T17:14:43.724599Z","shell.execute_reply":"2023-02-20T17:15:20.650350Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Trimap","metadata":{}},{"cell_type":"code","source":"%%time\n!pip install trimap \nimport trimap\n\nreducer = trimap.TRIMAP()\nr = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1 , alpha = alpha1 )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title(' trimap ' + str_data_inf  + ' n_samples='+str(len(r)), fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:15:20.654163Z","iopub.execute_input":"2023-02-20T17:15:20.654683Z","iopub.status.idle":"2023-02-20T17:15:49.423157Z","shell.execute_reply.started":"2023-02-20T17:15:20.654628Z","shell.execute_reply":"2023-02-20T17:15:49.421976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MulticoreTSNE","metadata":{}},{"cell_type":"code","source":"%%time\n!pip install MulticoreTSNE\n\nfrom MulticoreTSNE import MulticoreTSNE as TSNE\n\nreducer = TSNE(n_jobs=4)\nr = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1 , marker =  marker1 , alpha = alpha1 )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    \n    plt.title(' MulticoreTSNE ' + str_data_inf +  ' n_samples='+str(len(r)), fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:15:49.425279Z","iopub.execute_input":"2023-02-20T17:15:49.425748Z","iopub.status.idle":"2023-02-20T17:16:10.882662Z","shell.execute_reply.started":"2023-02-20T17:15:49.425700Z","shell.execute_reply":"2023-02-20T17:16:10.881406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom MulticoreTSNE import MulticoreTSNE as TSNE\nstr_model = 'MulticoreTSNE'\n\nn_x_subplots = 2; c = 0\n\nfor perplexity in [5,10,30,100]:\n    reducer = TSNE(n_jobs=4, perplexity=perplexity)\n    r = reducer.fit_transform(X)\n    \n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        plt.suptitle(str_model + ' '+ str_data_inf + ' n_samples='+str(len(r)), fontsize = 20)\n    c += 1; fig.add_subplot(1,n_x_subplots ,c)\n    i,j = 0,1\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1 , alpha = alpha1  )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title        \n       \n    plt.title('perplexity '+str(perplexity), fontsize = 20)\n    \n    plt.xlabel(str_model+str(i+1), fontsize = 20)\n    plt.ylabel(str_model+str(j+1), fontsize = 20)  \n    if c == 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title    \n        plt.legend( bbox_to_anchor = (1.15, 1.), loc='upper right')\n    else:\n        plt.legend('')\n        \n    \n        \n    \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:16:10.885083Z","iopub.execute_input":"2023-02-20T17:16:10.885610Z","iopub.status.idle":"2023-02-20T17:16:18.001695Z","shell.execute_reply.started":"2023-02-20T17:16:10.885552Z","shell.execute_reply":"2023-02-20T17:16:18.000618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# pacmap","metadata":{}},{"cell_type":"code","source":"!pip install pacmap","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:16:18.003406Z","iopub.execute_input":"2023-02-20T17:16:18.003727Z","iopub.status.idle":"2023-02-20T17:16:29.940036Z","shell.execute_reply.started":"2023-02-20T17:16:18.003698Z","shell.execute_reply":"2023-02-20T17:16:29.938597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n#!pip install trimap \nimport pacmap\n\nreducer = pacmap.PaCMAP()\nr = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:#,(0,2),(1,2),(2,3),(2,4),(3,4)]:\n    fig = plt.figure(figsize = (20,12) ); c = 0\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1 , alpha = alpha1 )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=20) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=20) # for legend title        \n    plt.title(' pacmap ' + str_data_inf  + ' n_samples='+str(len(r)), fontsize = 20)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:16:29.942561Z","iopub.execute_input":"2023-02-20T17:16:29.943080Z","iopub.status.idle":"2023-02-20T17:16:50.483924Z","shell.execute_reply.started":"2023-02-20T17:16:29.943023Z","shell.execute_reply":"2023-02-20T17:16:50.482772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport pacmap\nstr_model = 'PaCMAP'\n\n\nn_x_subplots = 2; c = 0\n\nfor param1 in [3,5,10,50]:\n    reducer = pacmap.PaCMAP( n_neighbors=param1)\n    r = reducer.fit_transform(X)\n    \n    if c % n_x_subplots == 0:\n        if c > 0:\n            plt.show()\n        fig = plt.figure(figsize = (20,5) ); c = 0\n        plt.suptitle(str_model + ' '+ str_data_inf + ' n_samples='+str(len(r)), fontsize = 20)\n    c += 1; fig.add_subplot(1,n_x_subplots ,c)\n    i,j = 0,1\n    ax = sns.scatterplot(x = r[:,i],y=r[:,j], hue = v4color,  palette = palette1, marker =  marker1 , alpha = alpha1  )\n    if 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title        \n       \n    plt.title('n_neighbors '+str(param1), fontsize = 20)\n    \n    plt.xlabel(str_model+str(i+1), fontsize = 20)\n    plt.ylabel(str_model+str(j+1), fontsize = 20)  \n    if c == 1:\n        plt.setp(ax.get_legend().get_texts(), fontsize=12) # for legend text\n        plt.setp(ax.get_legend().get_title(), fontsize=12) # for legend title    \n        plt.legend( bbox_to_anchor = (1.15, 1.), loc='upper right')\n    else:\n        plt.legend('')\n        \n    \n        \n    \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:16:50.485702Z","iopub.execute_input":"2023-02-20T17:16:50.486055Z","iopub.status.idle":"2023-02-20T17:17:11.165278Z","shell.execute_reply.started":"2023-02-20T17:16:50.486023Z","shell.execute_reply":"2023-02-20T17:17:11.164085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Direct analysis of correlations  and clustering proteins based on correlations ","metadata":{}},{"cell_type":"code","source":"%%time\ndf_corr = df.corr(method='pearson')\ndf_corr.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:17:11.167113Z","iopub.execute_input":"2023-02-20T17:17:11.168106Z","iopub.status.idle":"2023-02-20T17:17:15.090003Z","shell.execute_reply.started":"2023-02-20T17:17:11.168058Z","shell.execute_reply":"2023-02-20T17:17:15.088953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.Series( df_corr.abs().values.ravel() ).describe()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:17:15.091482Z","iopub.execute_input":"2023-02-20T17:17:15.091823Z","iopub.status.idle":"2023-02-20T17:17:15.104204Z","shell.execute_reply.started":"2023-02-20T17:17:15.091795Z","shell.execute_reply":"2023-02-20T17:17:15.103065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\nfrom scipy.spatial import distance\nfrom scipy.cluster import hierarchy\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:17:15.105759Z","iopub.execute_input":"2023-02-20T17:17:15.106153Z","iopub.status.idle":"2023-02-20T17:17:15.111549Z","shell.execute_reply.started":"2023-02-20T17:17:15.106119Z","shell.execute_reply":"2023-02-20T17:17:15.110460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"correlations = np.abs(df_corr)\ncorrelations_array = np.abs(df_corr)\n\nm_linkage = hierarchy.linkage(distance.pdist(correlations_array), method='average')  \n\nclusters = hierarchy.fcluster(m_linkage, 0.5*distance.pdist(correlations_array).max(), 'distance')\n\n\nclusters_list = []\n\nindex_clusters = []\n\nval, cnt = np.unique(clusters, return_counts=True)\n\nfor i,cluster in enumerate(clusters):        \n    clusters_list.append([df_corr.index[i], cluster])\n\n    index_cluster = f\"{df_corr.index[i]}_cl{cluster}_sz{int(cnt[val == cluster])}\"\n    index_clusters.append(index_cluster)\n\n\nindex_dict = dict(zip(df_corr.index, index_clusters))\n\n\ndf_corr.index = df_corr.index.map(index_dict.get) \n\ncorrelations_1 = np.abs(df_corr)\n\ng =  sns.clustermap(correlations_1, row_linkage=m_linkage, col_linkage=m_linkage, method=\"average\",cmap='vlag', figsize=(25, 25)) \nplt.setp(g.ax_heatmap.yaxis.get_majorticklabels(), rotation=0)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:17:15.112994Z","iopub.execute_input":"2023-02-20T17:17:15.113676Z","iopub.status.idle":"2023-02-20T17:17:18.652805Z","shell.execute_reply.started":"2023-02-20T17:17:15.113642Z","shell.execute_reply":"2023-02-20T17:17:18.651438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_clusters = pd.DataFrame(clusters_list, columns=['cd', 'cluster'])\ndf_corr_1 = df.corr(method='pearson')\n\npearson_v_counts = df_clusters['cluster'].value_counts()#.head(10)\npearson_df = df_clusters.merge(pearson_v_counts.to_frame(),\n                                left_on='cluster',\n                                right_index=True)\npearson_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    mask = df_corr_1.loc[x, x].abs()\n    result = mask.values.mean()\n    print('cluster size: ', len(x), ', mean abs correlation per cluster: ', np.round(result,3) , ', genes: ', list(x)) \n    print()\n","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:17:18.654627Z","iopub.execute_input":"2023-02-20T17:17:18.655068Z","iopub.status.idle":"2023-02-20T17:17:22.478509Z","shell.execute_reply.started":"2023-02-20T17:17:18.655027Z","shell.execute_reply":"2023-02-20T17:17:22.477245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result_list = []\nfor n in pearson:\n    mask = df_corr_1.loc[n, n].abs()\n    result = mask.values.mean()\n    result_list.append([len(n), result])\n    \ndf_cluster_mean_value = pd.DataFrame(result_list, columns=['cluster size', 'mean_correlation'])\ndf_cluster_mean_value.index.name = 'cluster ID'\ndf_cluster_mean_value = df_cluster_mean_value.sort_values(by=['mean_correlation'], ascending=False)\nprint('number of clusters:', df_cluster_mean_value.shape[0])\nprint('cluster sizes:',df_cluster_mean_value['cluster size'].tolist() )\nprint('mean correlations:',np.round(df_cluster_mean_value['mean_correlation'].tolist(),2) )\ndf_cluster_mean_value","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:17:22.479817Z","iopub.execute_input":"2023-02-20T17:17:22.480684Z","iopub.status.idle":"2023-02-20T17:17:22.504805Z","shell.execute_reply.started":"2023-02-20T17:17:22.480639Z","shell.execute_reply":"2023-02-20T17:17:22.503392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Compare with clusters obtained by other method (K-means)\n\nThe clusters from here: \nhttps://www.kaggle.com/code/alekse1pakhomov/multimodal-singlecell-correlation?scriptVersionId=119682605&cellId=1\n\nIt appears to be that the biggest (61 protein) cluster with low correlated proteins has been mainly splitted into 2:\nwith about 20 , 30 elements in each. \n\nThe cluster with 40 elements has been mainly splitted into two with 23, 22 elements.\n\nThe cluster with 24 elemtns fits cluster with 17 elements.\n\nOther clusters are small so we will not discuss them. \n\n","metadata":{}},{"cell_type":"code","source":"dict_clusters_1 = {}\nfor i,x in enumerate(pearson):\n    mask = df_corr_1.loc[x, x].abs()\n    result = mask.values.mean()\n    #print('cluster size: ', len(x), ', mean abs correlation per cluster: ', np.round(result,3) , ', genes: ', list(x)) \n    #print()\n    k = '1_cl_'+str(len(x)) + '_'+str(i)\n    dict_clusters_1[k] = list(x)\ndict_clusters_2 = {}\ndict_clusters_2['2_cl_22_4']  = ['CD3', 'CD7', 'CD194', 'CD16', 'CD25', 'Mouse-IgG1', 'Mouse-IgG2a', 'Mouse-IgG2b', 'Rat-IgG2b', 'CD196', 'CD69', 'CD27', 'CD64', 'CD35', 'CD39', 'CD42b', 'CD62P', 'Rat-IgG1', 'Rat-IgG2a', 'IgD', 'CD22', 'CD49a']\ndict_clusters_2['2_cl_30_3']  = ['CD40', 'CD8', 'CD56', 'CD19', 'CD4', 'CD14', 'CD20', 'CD146', 'IgM', 'CD195', 'CD161', 'CD1c', 'CD11b', 'CD21', 'CD79b', 'CD169', 'integrinB7', 'CD268', 'CD122', 'CD2', 'CD38', 'CD127', 'CD172a', 'CD93', 'CD73', 'TCRVa7.2', 'TCRVd2', 'CD158e1', 'CD352', 'CD94']\ndict_clusters_2['2_cl_23_5']  = ['CD154', 'CD105', 'CD45RO', 'CD279', 'TIGIT', 'CD185', 'CD223', 'KLRG1', 'CD134', 'CD141', 'CD314', 'CX3CR1', 'CD24', 'CD119', 'CD192', 'CD163', 'CD83', 'CD303', 'CD304', 'CD142', 'CD319', 'CD23', 'HLA-E']\ndict_clusters_2['2_cl_8_7']  = ['CD33', 'CD49f', 'CD58', 'CD29', 'CD63', 'CD49d', 'CD162', 'CD224']\ndict_clusters_2['2_cl_17_2']  = ['CD47', 'CD48', 'HLA-A-B-C', 'CD45RA', 'CD123', 'CD44', 'CD31', 'CD62L', 'CD95', 'HLA-DR', 'CD11a', 'CD244', 'CD13', 'CD49b', 'CD81', 'CD18', 'CD45']\ndict_clusters_2['2_cl_6_6']  = ['CD11c', 'CD54', 'CD72', 'CD9', 'CD328', 'CD101']\ndict_clusters_2['2_cl_22_1']  = ['Podoplanin', 'CD5', 'CD103', 'CD152', 'CD107a', 'CD1d', 'CD57', 'CD272', 'CD278', 'TCR', 'FceRIa', 'CD137', 'CD124', 'CD226', 'CD28', 'CD26', 'CD115', 'CD158', 'LOX-1', 'CD158b', 'CD85j', 'CD82']\ndict_clusters_2['2_cl_4_0']  = ['CD41', 'CD71', 'CD36', 'CD88']\n\nfor k2 in dict_clusters_2:\n    print(len(dict_clusters_2[k2] ))\nd = pd.DataFrame()\nfor k1 in dict_clusters_1:\n    for k2 in dict_clusters_2:\n        d.loc[k1,k2] = len( set(dict_clusters_1[k1]) &   set(dict_clusters_2[k2]) )\nsns.heatmap(d,annot=True,)        \nd        ","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:57:38.476041Z","iopub.execute_input":"2023-02-20T17:57:38.477148Z","iopub.status.idle":"2023-02-20T17:57:39.084272Z","shell.execute_reply.started":"2023-02-20T17:57:38.477106Z","shell.execute_reply":"2023-02-20T17:57:39.082859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes Sets Enrichment with GSEAPY","metadata":{}},{"cell_type":"code","source":"!pip install gseapy\nimport gseapy as gp","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:28:30.260407Z","iopub.execute_input":"2023-02-20T17:28:30.261265Z","iopub.status.idle":"2023-02-20T17:28:42.758755Z","shell.execute_reply.started":"2023-02-20T17:28:30.261200Z","shell.execute_reply":"2023-02-20T17:28:42.757546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for x in pearson:\n    mask = df_corr_1.loc[x, x].abs()\n    result = mask.values.mean()\n    if result > 0.2:\n        print('cluster size: ', len(x), ', mean correlation: ', np.round(result,3) , ', genes: ', x)\n        print()\n        list_genes = list(x)\n        #print(len(list_genes), list_genes[:50])\n        enr = gp.enrichr(\n            gene_list=list_genes, \n            gene_sets= ['KEGG_2021_Human', 'Reactome_2022','GO_Molecular_Function_2021'  ], # , 'GO_Biological_Process_2021', 'MSigDB_Oncogenic_Signatures',\n            organism='human',\n            outdir=None,\n        )\n        col = 'Adjusted P-value' #P-value # Adjusted P-value - with Bonferoni type correction for checking multiple lists, while just \"P-value\" - is without\n        display( enr.results.sort_values(col).head(5) )    \n        \n        print()\n        print()\n        ","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:28:42.760719Z","iopub.execute_input":"2023-02-20T17:28:42.761123Z","iopub.status.idle":"2023-02-20T17:29:01.088931Z","shell.execute_reply.started":"2023-02-20T17:28:42.761085Z","shell.execute_reply":"2023-02-20T17:29:01.087704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_gene_sets = ['KEGG_2021_Human', 'Reactome_2022','GO_Molecular_Function_2021'  ]\nlist_gene_sets = ['ChEA_2022',     'ENCODE_TF_ChIP-seq_2015', 'ENCODE_and_ChEA_Consensus_TFs_from_ChIP-X', \n    'TF-LOF_Expression_from_GEO', 'TF_Perturbations_Followed_by_Expression',     'TRANSFAC_and_JASPAR_PWMs',\n    'TRRUST_Transcription_Factors_2019',   ]\nfor x in pearson:\n    mask = df_corr_1.loc[x, x].abs()\n    result = mask.values.mean()\n    if result > 0.2:\n        print('cluster size: ', len(x), ', mean correlation: ', np.round(result,3) , ', genes: ', x)\n        print()\n        list_genes = list(x)\n        #print(len(list_genes), list_genes[:50])\n        enr = gp.enrichr(\n            gene_list=list_genes, \n            gene_sets= list_gene_sets, # , 'GO_Biological_Process_2021', 'MSigDB_Oncogenic_Signatures',\n            organism='human',\n            outdir=None,\n        )\n        col = 'Adjusted P-value' #P-value # Adjusted P-value - with Bonferoni type correction for checking multiple lists, while just \"P-value\" - is without\n        display( enr.results.sort_values(col).head(5) )    \n        \n        print()\n        print()\n        \n","metadata":{"execution":{"iopub.status.busy":"2023-02-20T17:29:04.903835Z","iopub.execute_input":"2023-02-20T17:29:04.904345Z","iopub.status.idle":"2023-02-20T17:29:44.182085Z","shell.execute_reply.started":"2023-02-20T17:29:04.904297Z","shell.execute_reply":"2023-02-20T17:29:44.180848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}