{"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\nSome basic look on data.\nCorrelations, dimensional reductions, etc. \n\n## **Findings:**\n\n 1) Targets are quite noisy with potentially many outliers, which probably should be filtered out e.g. CD11c (max outlier), CD36 (min outliers)\n Other examples of targets with outliers : \n Max: CD11c, CD328, CD57, CD328, CD2\n Min: CD36, CD88, CD163, CD23\n It might be almost all targets contain outliers and that should be correct before training \n \n 2) Targets depends on days and donors - so it might be one needs to build separate models by donors\n \n 3) Some targets are NOT CD**: 'HLA-A-B-C', 'TIGIT', 'Mouse-IgG1', 'Mouse-IgG2a', 'Mouse-IgG2b',\n 'Rat-IgG2b', 'Podoplanin', 'IgM', 'KLRG1', 'HLA-DR', 'CX3CR1',\n 'integrinB7', 'TCR', 'Rat-IgG1', 'Rat-IgG2a', 'FceRIa', 'IgD',\n 'TCRVa7.2', 'TCRVd2', 'LOX-1', 'HLA-E'\n \n Some of them are probably sort of control - not so clear, what one should do with them. \n Especially: 'Rat-IgG1', 'Rat-IgG2a', 'Mouse-IgG1', 'Mouse-IgG2a', 'Mouse-IgG2b',\n We can see that they are not much correlated with other targets.\n \n 4) Analsysis of top correlated and top - see below. In particular pay attention that: CD62L is top anticorrelated with many targets \n \n 5) Dimensional reductions (PCA,UMAP,etc) do NOT show clear clusters - thus cell types (and donors, days) are not well separated from each other. Thought coloring by cell type is not random and correspond to different parts of UMAP-like figures. In ATAC-seq data cell types are seen better. \n \n 6) Some pairs shows some \"mutually exclusive\" behaviour - e.g. when one grows, another is stable \n (in spirit of anticorrlation, but a bit different and stronger). \n That might have some biological interpretation - may be these are markers of some cell types \n Some Mutually Exclusive:\n CD88 CD328\n CD36 CD328\n See: https://www.kaggle.com/code/alexandervc/mmscel-eda-targets-citeseq-02?scriptVersionId=106222383&cellId=21\n \n \nRemark: To item 2 - control antibodies - which probably : \nhttps://www.thermofisher.com/fr/fr/home/life-science/antibodies/primary-antibodies/control-antibodies/isotype-control-antibodies.html\n \n \nRemark: Discussion here questions how reasonable is the normalization used by the organizers - since it shows negative values, while concentration (of protein) cannot be negative.\nhttps://www.kaggle.com/competitions/open-problems-multimodal/discussion/346210\n \n \nPrevious version of analysis of the same part of data:\nhttps://www.kaggle.com/alexandervc/mmscel-eda-bioinfo-targets-citeseq-01\n","metadata":{"execution":{"iopub.status.busy":"2022-09-20T09:35:25.343879Z","iopub.execute_input":"2022-09-20T09:35:25.344340Z","iopub.status.idle":"2022-09-20T09:35:25.352714Z","shell.execute_reply.started":"2022-09-20T09:35:25.344300Z","shell.execute_reply":"2022-09-20T09:35:25.350815Z"}}},{"cell_type":"markdown","source":"## Findings on correlations:\n\nTop correlated pairs:\n\n CD71\tCD88\t0.817769\n CD115\tCD88\t0.812867\n CD29\tCD49d\t0.714796\n CD155\tCD29\t0.695378\n CD29\tCD81\t0.67392\n HLA-A-B-C\tCD81\t0.666818\n CD47\tHLA-A-B-C\t0.664981\n CD112\tCD31\t0.652837\n HLA-A-B-C\tCD31\t0.644859\n \nTop anti-correlated pairs:\n\n CD62L\tCD41\t-0.420337\n CD32\tCD62L\t-0.414331\n CD62L\tCD88\t-0.411008\n CD62L\tCD36\t-0.409057\n CD62L\tCD71\t-0.350056\n CD62L\tCD115\t-0.338883\n CD62L\tFceRIa\t-0.320145\n CD244\tCD41\t-0.312103\n CD244\tCD36\t-0.291549\n CD44\tCD36\t-0.286202\n CD11a\tCD41\t-0.27769\n CD44\tCD41\t-0.276898 \n\nSome groups of correlated targets: \n\n 0.8 correlation_threshold \n Top 5 cluster sizes: [3 1 1 1 1] seconds passed: 0.0\n Not more than 50 Ids in largest correlated components:\n 1-th component: ['CD71' 'CD115' 'CD88']\n\n 0.7 correlation_threshold \n Top 5 cluster sizes: [3 2 1 1 1] seconds passed: 0.0\n Not more than 50 Ids in largest correlated components:\n 1-th component: ['CD71' 'CD115' 'CD88']\n 2-th component: ['CD29' 'CD49d']\n\n 0.6 correlation_threshold \n Top 5 cluster sizes: [14 5 2 1 1] seconds passed: 0.0\n Not more than 50 Ids in largest correlated components:\n 1-th component: ['CD155' 'CD112' 'CD47' 'HLA-A-B-C' 'CD45RA' 'CD31' 'CD11a' 'CD13' 'CD29'\n 'CD81' 'CD18' 'CD45' 'CD49d' 'CD162']\n 2-th component: ['CD272' 'FceRIa' 'CD71' 'CD115' 'CD88']\n 3-th component: ['CD41' 'CD36']\n\n 0.5 correlation_threshold \n Top 5 cluster sizes: [39 4 2 1 1] seconds passed: 0.0\n Not more than 50 Ids in largest correlated components:\n 1-th component: ['CD155' 'CD112' 'CD47' 'CD48' 'CD33' 'HLA-A-B-C' 'CD45RA' 'CD123' 'CD49f'\n 'CD44' 'CD31' 'Podoplanin' 'CD32' 'CD62L' 'CD107a' 'CD95' 'HLA-DR' 'CD1d'\n 'CD272' 'CD58' 'CD11a' 'CD244' 'FceRIa' 'CD137' 'CD13' 'CD29' 'CD49b'\n 'CD81' 'CD18' 'CD45' 'CD71' 'CD26' 'CD115' 'CD63' 'CD49d' 'CD162' 'CD85j'\n 'CD88' 'CD224']\n Largest clique: ['CD31' 'HLA-A-B-C' 'CD45' 'CD13' 'CD81' 'CD29' 'HLA-DR' 'CD112'] \n Non intersecting Maximal Clique: ['CD58' 'CD107a' 'CD88' 'CD115']\n 2-th component: ['CD141' 'CX3CR1' 'CD24' 'CD83']\n 3-th component: ['CD41' 'CD36']","metadata":{}},{"cell_type":"markdown","source":"# Install/Import ","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":"2022-09-24T16:05:34.086970Z","iopub.execute_input":"2022-09-24T16:05:34.087877Z","iopub.status.idle":"2022-09-24T16:05:34.135378Z","shell.execute_reply.started":"2022-09-24T16:05:34.087762Z","shell.execute_reply":"2022-09-24T16:05:34.134065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install scanpy\nimport scanpy as sc\nimport anndata\n\nimport time\nt0start = time.time()\n\nimport pandas as pd\nimport numpy as np\nimport os\nimport sys\n\nimport matplotlib.pyplot as plt\n#plt.style.use('dark_background')\nimport seaborn as sns\n\n#If you see a urllib warning running this cell, go to \"Settings\" on the right hand side, \n#and turn on internet. Note, you need to be phone verified.\n!pip install --quiet tables\n\nimport h5py\n!pip install hdf5plugin~=2.0 # https://forum.hdfgroup.org/t/cant-open-directory-usr-local-hdf5-lib-plugin/9738/4\nimport hdf5plugin\n\nDATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")\n\ndf_cell = pd.read_csv(FP_CELL_METADATA)\ndf_cell","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:05:34.137469Z","iopub.execute_input":"2022-09-24T16:05:34.137799Z","iopub.status.idle":"2022-09-24T16:06:10.440395Z","shell.execute_reply.started":"2022-09-24T16:05:34.137773Z","shell.execute_reply":"2022-09-24T16:06:10.438447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"%%time\ndf = pd.read_hdf(FP_CITE_TRAIN_TARGETS)\nstr_data_inf = 'CITEseq Targets '\nprint(df.shape)\ndisplay(df.head() )","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:10.442270Z","iopub.execute_input":"2022-09-24T16:06:10.442675Z","iopub.status.idle":"2022-09-24T16:06:11.297194Z","shell.execute_reply.started":"2022-09-24T16:06:10.442636Z","shell.execute_reply":"2022-09-24T16:06:11.296320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('1.2G RAM consumed')","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:11.299484Z","iopub.execute_input":"2022-09-24T16:06:11.299746Z","iopub.status.idle":"2022-09-24T16:06:11.304517Z","shell.execute_reply.started":"2022-09-24T16:06:11.299723Z","shell.execute_reply":"2022-09-24T16:06:11.303612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look column names","metadata":{}},{"cell_type":"code","source":"m = np.array([t.startswith('CD') for t in df.columns])\nprint(m.sum(), (~m).sum())\nprint( df.columns[~m] )\nv = [(t[2:]).isnumeric() for t in df.columns ]\nv2 = [int(t[2:]) for t in df.columns[v] ]\n#print(np.sum(v)) ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:11.306385Z","iopub.execute_input":"2022-09-24T16:06:11.306748Z","iopub.status.idle":"2022-09-24T16:06:11.316587Z","shell.execute_reply.started":"2022-09-24T16:06:11.306718Z","shell.execute_reply":"2022-09-24T16:06:11.315374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"w = [ not((t[2:]).isnumeric()) for t in df.columns[m] ]\ndf.columns[m][w]","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:11.317732Z","iopub.execute_input":"2022-09-24T16:06:11.317995Z","iopub.status.idle":"2022-09-24T16:06:11.332073Z","shell.execute_reply.started":"2022-09-24T16:06:11.317970Z","shell.execute_reply":"2022-09-24T16:06:11.330935Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(v2), ': 99 out of 352 - roughly speaking 1/3.5 ') \nprint( list(np.sort(v2) ) )","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:11.333624Z","iopub.execute_input":"2022-09-24T16:06:11.333926Z","iopub.status.idle":"2022-09-24T16:06:11.345668Z","shell.execute_reply.started":"2022-09-24T16:06:11.333895Z","shell.execute_reply":"2022-09-24T16:06:11.344777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['CD36'].sort_values()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:11.346952Z","iopub.execute_input":"2022-09-24T16:06:11.347179Z","iopub.status.idle":"2022-09-24T16:06:11.390444Z","shell.execute_reply.started":"2022-09-24T16:06:11.347156Z","shell.execute_reply":"2022-09-24T16:06:11.389822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Statitics overall ","metadata":{}},{"cell_type":"code","source":"v = df.values.ravel()\nfig = plt.figure(figsize = (20,5)); \nfig.add_subplot(1,2,1)\nplt.hist(v , bins = 100 )\nplt.title('All elements ', fontsize = 20)\nfig.add_subplot(1,2,2)\nplt.hist(v[ (v>np.percentile(v,1)) & (v 0 '] = ( v > 0 ).sum() \nt[' < 0 '] = ( v < 0 ).sum() \nt[' == 0 '] = ( v == 0 ).sum() \nt[' NaN '] = np.isnan(v).sum()\npd.Series(t).plot(kind='barh', figsize=(15,5)); plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:15.798536Z","iopub.execute_input":"2022-09-24T16:06:15.798894Z","iopub.status.idle":"2022-09-24T16:06:15.982446Z","shell.execute_reply.started":"2022-09-24T16:06:15.798861Z","shell.execute_reply":"2022-09-24T16:06:15.981592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Basics statistics - calculate statistics ('Mean','AbsMean', 'Max', 'Min', 'Median', 'Std') for each column, then look on obtained results overall","metadata":{"execution":{"iopub.status.busy":"2022-09-20T15:48:19.988997Z","iopub.execute_input":"2022-09-20T15:48:19.989765Z","iopub.status.idle":"2022-09-20T15:48:20.017704Z","shell.execute_reply.started":"2022-09-20T15:48:19.989708Z","shell.execute_reply":"2022-09-20T15:48:20.010457Z"}}},{"cell_type":"code","source":"d = pd.DataFrame()\n\nfor str_inf in ['Mean','AbsMean', 'Max', 'Min', 'Median', 'Std']:\n if str_inf == 'Mean': v = df.mean(axis = 0)\n elif str_inf == 'AbsMean': v = df.abs().mean(axis = 0)\n elif str_inf == 'Max': v = df.max(axis = 0)\n elif str_inf == 'Min': v = df.min(axis = 0)\n elif str_inf == 'Median': v = df.median(axis = 0)\n elif str_inf == 'Std': v = df.std(axis = 0)\n \n t = np.round( np.array( [v.mean(), v.max(), v.min(), v.std()] ) , 3 )\n s =( 'Mean Max Min Std: ' + str(t) ).replace('[','').replace(']','')\n _, axs = plt.subplots(1, 2, figsize=(20, 4))\n plt.suptitle ('Statistics: '+str_inf + '. Its stat: ' + s, fontsize = 20 )\n axs[0].plot(v.values)\n axs[0].set_xlabel('Data Column Index', fontsize=14)# 'medium') \n axs[0].set_facecolor('lightblue')\n #axs[0].set_title(str_inf, fontsize = 20)\n axs[1].hist(v.values, bins = 100)\n plt.gca().set_facecolor('lightblue')\n #axs[1].set_title('Histogram for '+str_inf , fontsize = 20 )\n plt.show()\n \n v.name = str_inf\n d = d.join(v, how = 'outer')\nd['Max2AbsMean'] = d['Max']/d['AbsMean']\nd['Min2AbsMean'] = d['Min'].abs()/d['AbsMean']\nd['Max2Min'] = d['Max'].abs()/d['Min'].abs()\nd['Min2Max'] = d['Min'].abs()/d['Max'].abs()\n \ns = set(); \nk_top = 10\nfor i,str_inf in enumerate( d.columns):\n s1 = d[str_inf].abs().sort_values(ascending = False).index[:k_top]\n print(str_inf, 'Top ',k_top, ':', list(s1) )\n if i == 0: s = set(s1)\n else: s = s & set(s1)\n if i == 3: print('Intersection of tops:', s)\nprint('Intersection of tops:', s)\nd.sort_values('AbsMean', inplace = True) \ndisplay(d.head(20)) \nprint()\ndisplay(d.tail(5)) \n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:15.983578Z","iopub.execute_input":"2022-09-24T16:06:15.984114Z","iopub.status.idle":"2022-09-24T16:06:18.595001Z","shell.execute_reply.started":"2022-09-24T16:06:15.984084Z","shell.execute_reply":"2022-09-24T16:06:18.593704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d.corr()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:18.596597Z","iopub.execute_input":"2022-09-24T16:06:18.596998Z","iopub.status.idle":"2022-09-24T16:06:18.621547Z","shell.execute_reply.started":"2022-09-24T16:06:18.596964Z","shell.execute_reply":"2022-09-24T16:06:18.619798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look for columns with potential outliers","metadata":{}},{"cell_type":"code","source":"s = set(); \nk_top = 6\nn_y_subplots = 1\nn_x_subplots = 3 # int( k_top/ n_y_subplots)\nc = 0\nfor i,str_inf in enumerate( ['Max','Min','Max2AbsMean', 'Min2AbsMean', 'Std']):\n s1 = d[str_inf].abs().sort_values(ascending = False).index[:k_top]\n \n str_suptitle = str_inf + ' Top ' + str(k_top) +':'+ str( list(s1) ).replace('[','').replace(']','')\n s = s | set(s1)\n for icol, col in enumerate(s1):\n #fig = plt.figure(figsize = (20,3)); \n if c % (n_x_subplots*n_y_subplots) == 0:\n if c > 0:\n plt.show()\n fig = plt.figure(figsize = (20,3) ); c = 0\n #plt.suptitle(str_suptitle)\n #print(str_inf, 'Top ',k_top, ':', list(s1) )\n\n c += 1; fig.add_subplot(n_y_subplots,n_x_subplots ,c) \n sns.boxplot(x = df[col])\n plt.title(str(col) + ' top ' +str(icol+1) + ' for ' + str_inf, fontsize = 15)\nplt.show()\nprint('All:')\nprint(s)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:18.622817Z","iopub.execute_input":"2022-09-24T16:06:18.623135Z","iopub.status.idle":"2022-09-24T16:06:21.390040Z","shell.execute_reply.started":"2022-09-24T16:06:18.623110Z","shell.execute_reply":"2022-09-24T16:06:21.388450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Violins plots - one can see difference for some columns with outliers and more \"normal\" ones ","metadata":{"execution":{"iopub.status.busy":"2022-09-22T08:57:55.962930Z","iopub.execute_input":"2022-09-22T08:57:55.963404Z","iopub.status.idle":"2022-09-22T08:57:55.972734Z","shell.execute_reply.started":"2022-09-22T08:57:55.963365Z","shell.execute_reply":"2022-09-22T08:57:55.971279Z"}}},{"cell_type":"code","source":"s = set(); \nk_top = 4\nn_y_subplots = 1\nn_x_subplots = 4 # int( k_top/ n_y_subplots)\nc = 0\nfor i,str_inf in enumerate( ['Max2AbsMean', 'AbsMean']):\n s1 = d[str_inf].abs().sort_values(ascending = False).index[:k_top]\n \n str_suptitle = str_inf + ' Top ' + str(k_top) +':'+ str( list(s1) ).replace('[','').replace(']','')\n s = s | set(s1)\n for icol, col in enumerate(s1):\n #fig = plt.figure(figsize = (20,3)); \n if c % (n_x_subplots*n_y_subplots) == 0:\n if c > 0:\n plt.show()\n fig = plt.figure(figsize = (20,3) ); c = 0\n #plt.suptitle(str_suptitle)\n print(); print(str_inf, 'Top ',k_top, ':', list(s1) ); print()\n \n\n c += 1; fig.add_subplot(n_y_subplots,n_x_subplots ,c) \n \n sns.violinplot(x = df[col] )\n plt.title(str(col) + ' top ' +str(icol+1) + ' for ' + str_inf, fontsize = 12)\n \nplt.show()\nprint('All:')\nprint(s)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:21.391713Z","iopub.execute_input":"2022-09-24T16:06:21.392130Z","iopub.status.idle":"2022-09-24T16:06:23.557784Z","shell.execute_reply.started":"2022-09-24T16:06:21.392103Z","shell.execute_reply":"2022-09-24T16:06:23.556370Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Histograms for selected columns - similar to violin plots - one can see difference of columns with outliers comparing to more \"normal\" ones\n","metadata":{}},{"cell_type":"code","source":"s = set(); \nk_top = 3\nn_y_subplots = 1\nn_x_subplots = 3 # int( k_top/ n_y_subplots)\nc = 0\nfor i,str_inf in enumerate( ['Max2AbsMean', 'AbsMean'] ):\n s1 = d[str_inf].abs().sort_values(ascending = False).index[:k_top]\n \n str_suptitle = str_inf + ' Top ' + str(k_top) +':'+ str( list(s1) ).replace('[','').replace(']','')\n s = s | set(s1)\n for icol, col in enumerate(s1):\n #fig = plt.figure(figsize = (20,3)); \n if c % (n_x_subplots*n_y_subplots) == 0:\n if c > 0:\n plt.show()\n fig = plt.figure(figsize = (20,3) ); c = 0\n #plt.suptitle(str_suptitle)\n #print(str_inf, 'Top ',k_top, ':', list(s1) )\n\n c += 1; fig.add_subplot(n_y_subplots,n_x_subplots ,c) \n \n #sns.distplot(x = df[col] , color='g', bins=100, hist_kws={'alpha': 0.4});\n sns.histplot(x = df[col] , color='g', bins=100, kde = True)# , hist_kws={'alpha': 0.4});\n # distplot - depricacted: how to substitute: https://gist.github.com/mwaskom/de44147ed2974457ad6372750bbe5751\n plt.title(str(col) + ' top ' +str(icol+1) + ' for ' + str_inf, fontsize = 15)\n \nplt.show()\nprint('All:')\nprint(s)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:23.558657Z","iopub.execute_input":"2022-09-24T16:06:23.558879Z","iopub.status.idle":"2022-09-24T16:06:27.153930Z","shell.execute_reply.started":"2022-09-24T16:06:23.558836Z","shell.execute_reply":"2022-09-24T16:06:27.152267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pair plots for selected columns\n\nObserve some \"mutually exclusive\" patterns for some\nSome pairs shows some \"mutually exclusive\" behaviour - e.g. when one grows, another is stable \n(in spirit of anticorrlation, but a bit different and stronger). \nThat might have some biological interpretation - may be these are markers of some cell types \n\n Some Mutually Exclusive:\n CD88 CD328\n CD36 CD328\n \n \n","metadata":{}},{"cell_type":"code","source":"s = set(); \nk_top = 2\nfor i,str_inf in enumerate( ['Max','Min','AbsMean', 'Std']):\n s1 = d[str_inf].abs().sort_values(ascending = False).index[:k_top]\n print(str_inf, 'Top ',k_top, ':', list(s1) )\n s = s | set(s1)\nprint(s)\n\nsns.set()\ncolumns = list(s)\nsns.pairplot(df[columns],height = 2 ,kind ='scatter',diag_kind='kde')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:27.155535Z","iopub.execute_input":"2022-09-24T16:06:27.155875Z","iopub.status.idle":"2022-09-24T16:06:38.875130Z","shell.execute_reply.started":"2022-09-24T16:06:27.155828Z","shell.execute_reply":"2022-09-24T16:06:38.873985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mean vs Std and related plots \n\nThese figure may allow to find columns with potential outliers","metadata":{}},{"cell_type":"code","source":"sns.set()\ncolumns = list(s)\nsns.pairplot(d[['AbsMean','Max','Min','Std']],height = 2 ,kind ='scatter',diag_kind='kde')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:38.876311Z","iopub.execute_input":"2022-09-24T16:06:38.876543Z","iopub.status.idle":"2022-09-24T16:06:40.993462Z","shell.execute_reply.started":"2022-09-24T16:06:38.876520Z","shell.execute_reply":"2022-09-24T16:06:40.992429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_x_subplots = 4\nfig = plt.figure(figsize = (20,10) ); c = 0\nfor (col1, col2, col3) in [('Mean', 'Std', 'Min'), ('AbsMean', 'Max', 'Min'), ('AbsMean', 'Min', 'Max'), ('Min', 'Max', 'Std'),]:\n c += 1; fig.add_subplot(1,n_x_subplots ,c) \n ax = sns.scatterplot(x = d[col1], y=d[col2], hue = d[col3], palette = 'rainbow')\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 ax.set_xlabel(col1, fontsize=20)\n ax.set_ylabel(col2, fontsize=20)# 'medium') \n\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:40.994576Z","iopub.execute_input":"2022-09-24T16:06:40.994812Z","iopub.status.idle":"2022-09-24T16:06:42.067753Z","shell.execute_reply.started":"2022-09-24T16:06:40.994789Z","shell.execute_reply":"2022-09-24T16:06:42.066432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Simple plots and histograms by columns\n\nDay/donor dependence is seen","metadata":{}},{"cell_type":"code","source":"n_x_subplots = 5\nc = 0\nfor col in list(d.sort_values('Min').index[:140]): # + list(d.sort_values('Min').index[-10:]) :\n \n if c % n_x_subplots == 0:\n if c > 0:\n plt.show()\n fig = plt.figure(figsize = (20,3) ); c = 0\n #plt.suptitle(str(str_data_inf) + ' n_samples='+str(len(r)),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED \n c += 1; fig.add_subplot(1,n_x_subplots ,c) \n plt.plot(df[col].values, label = col)\n plt.title(col,fontsize = 14)\n plt.gca().set_facecolor('lightblue')\n #plt.legend(fontsize = 20)\nplt.show()\n \n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:06:42.068981Z","iopub.execute_input":"2022-09-24T16:06:42.069270Z","iopub.status.idle":"2022-09-24T16:07:01.636120Z","shell.execute_reply.started":"2022-09-24T16:06:42.069245Z","shell.execute_reply":"2022-09-24T16:07:01.635183Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf.iloc[:,:70].hist(figsize=(16, 20), bins=50, xlabelsize=8, ylabelsize=8); # ; avoid having the matplotlib verbose informations\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:01.637382Z","iopub.execute_input":"2022-09-24T16:07:01.638060Z","iopub.status.idle":"2022-09-24T16:07:13.450203Z","shell.execute_reply.started":"2022-09-24T16:07:01.638033Z","shell.execute_reply":"2022-09-24T16:07:13.448533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndf.iloc[:,70:140].hist(figsize=(16, 20), bins=50, xlabelsize=8, ylabelsize=8); # ; avoid having the matplotlib verbose informations\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:13.451807Z","iopub.execute_input":"2022-09-24T16:07:13.452146Z","iopub.status.idle":"2022-09-24T16:07:25.785388Z","shell.execute_reply.started":"2022-09-24T16:07:13.452121Z","shell.execute_reply":"2022-09-24T16:07:25.783821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look on missing values","metadata":{}},{"cell_type":"code","source":"v = df.isna().sum(axis = 0)\nif np.sum(v) > 0:\n #v.plot()\n _, axs = plt.subplots(1, 2, figsize=(20, 4))\n axs[0].plot(v.values)\n axs[0].set_title('Counts NA')\n axs[1].hist(v.values, bins = 100)\n axs[1].set_title('Histogram for NA counts')\n plt.show()\n display(v.describe())\n for i in range(5):\n print( i, 'missing values: ', (v==i).sum() , ' %%: ', '%.1f'%((v==i).sum()/len(v)*100) )\n print(v.sort_values(ascending = False).head(10))\nelse:\n print('No missing values ')","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:25.787149Z","iopub.execute_input":"2022-09-24T16:07:25.787461Z","iopub.status.idle":"2022-09-24T16:07:25.808013Z","shell.execute_reply.started":"2022-09-24T16:07:25.787435Z","shell.execute_reply":"2022-09-24T16:07:25.806916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlation analysis ","metadata":{}},{"cell_type":"code","source":"list_columns_ids = list( df.columns )\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:25.809767Z","iopub.execute_input":"2022-09-24T16:07:25.810234Z","iopub.status.idle":"2022-09-24T16:07:25.814897Z","shell.execute_reply.started":"2022-09-24T16:07:25.810207Z","shell.execute_reply":"2022-09-24T16:07:25.813731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif df.isna().sum(axis = 0).sum() > 0:\n cm = df.corr()\nelse: # Fast way: \n cm = pd.DataFrame(np.corrcoef(df.T) , index = list_columns_ids, columns = list_columns_ids) \ncm","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:25.816548Z","iopub.execute_input":"2022-09-24T16:07:25.816911Z","iopub.status.idle":"2022-09-24T16:07:25.953073Z","shell.execute_reply.started":"2022-09-24T16:07:25.816883Z","shell.execute_reply":"2022-09-24T16:07:25.952314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nv = cm.values[np.triu_indices(len(cm),k=1)]\nprint(len(v), len(cm)*(len(cm)-1)/2 )\nplt.hist(v, bins = 1000)# , color='k')\nplt.show()\nprint(pd.Series(v).describe() )","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:25.954441Z","iopub.execute_input":"2022-09-24T16:07:25.954988Z","iopub.status.idle":"2022-09-24T16:07:27.437449Z","shell.execute_reply.started":"2022-09-24T16:07:25.954957Z","shell.execute_reply":"2022-09-24T16:07:27.436121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Top correlated and anticorrelated ","metadata":{}},{"cell_type":"code","source":"t = np.min([30, len(v)])\nthreshold4top_corr = np.sort(v[~np.isnan(v)])[-t ]\nthreshold4bottom_corr = np.sort(v[~np.isnan(v)])[ t ]\n \nthreshold4top_corr, threshold4bottom_corr\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:27.438764Z","iopub.execute_input":"2022-09-24T16:07:27.439591Z","iopub.status.idle":"2022-09-24T16:07:27.449056Z","shell.execute_reply.started":"2022-09-24T16:07:27.439534Z","shell.execute_reply":"2022-09-24T16:07:27.448017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n\nz = np.where( np.triu(cm.values,1) >= threshold4top_corr )\ndf_corrs = pd.DataFrame()\nIX = 0\nfor (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_corrs.loc[IX, 'I1'] = i\n df_corrs.loc[IX, 'I2'] = j\n df_corrs.loc[IX, 'Id1'] = list_columns_ids[i]\n df_corrs.loc[IX, 'Id2'] = list_columns_ids[j]\n df_corrs.loc[IX, 'Correlation'] = cm.iloc[i,j]\n \ndf_corrs['I1'] = df_corrs['I1'].astype(int) \ndf_corrs['I2'] = df_corrs['I2'].astype(int) \n\nprint(df_corrs.shape)\n\ndisplay ( df_corrs.sort_values('Correlation' , ascending = False) )\n ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:27.450547Z","iopub.execute_input":"2022-09-24T16:07:27.450917Z","iopub.status.idle":"2022-09-24T16:07:27.510375Z","shell.execute_reply.started":"2022-09-24T16:07:27.450887Z","shell.execute_reply":"2022-09-24T16:07:27.508986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n\nz = np.where( np.triu(cm.values,1) <= threshold4bottom_corr )\ndf_corrs = pd.DataFrame()\nIX = 0\nfor (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_corrs.loc[IX, 'I1'] = i\n df_corrs.loc[IX, 'I2'] = j\n df_corrs.loc[IX, 'Id1'] = list_columns_ids[i]\n df_corrs.loc[IX, 'Id2'] = list_columns_ids[j]\n df_corrs.loc[IX, 'Correlation'] = cm.iloc[i,j]\n \ndf_corrs['I1'] = df_corrs['I1'].astype(int) \ndf_corrs['I2'] = df_corrs['I2'].astype(int) \n\nprint(df_corrs.shape)\n\ndisplay ( df_corrs.sort_values('Correlation' , ascending = True) )\n ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:27.516482Z","iopub.execute_input":"2022-09-24T16:07:27.516829Z","iopub.status.idle":"2022-09-24T16:07:27.576464Z","shell.execute_reply.started":"2022-09-24T16:07:27.516803Z","shell.execute_reply":"2022-09-24T16:07:27.575003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v = np.abs(cm).sum(axis = 1)\nv = v.sort_values(ascending = False)\nprint('Top and anti-top correlated:')\ndisplay(v.head(10))\nprint()\ndisplay(v.tail(10))\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:27.578649Z","iopub.execute_input":"2022-09-24T16:07:27.579008Z","iopub.status.idle":"2022-09-24T16:07:27.596038Z","shell.execute_reply.started":"2022-09-24T16:07:27.578974Z","shell.execute_reply":"2022-09-24T16:07:27.594521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlation heatmaps","metadata":{}},{"cell_type":"code","source":"N= 12; print(N, 'top correlated:')\nlist_top_cor = list( v[:N].index)\nplt.figure(figsize=(20, 12))\ncm2 = cm.loc[list_top_cor,list_top_cor]\nax = sns.heatmap(cm2[(cm2 >= 0.5) | (cm2 <= -0.1)], \n cmap='viridis', vmax=1.0, vmin=-1.0, linewidths=0.1,\n annot=True, annot_kws={\"size\": 24}, square=True);\nax.xaxis.set_tick_params(labelsize=14)\nax.yaxis.set_tick_params(labelsize=14)","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:27.598092Z","iopub.execute_input":"2022-09-24T16:07:27.598469Z","iopub.status.idle":"2022-09-24T16:07:28.041794Z","shell.execute_reply.started":"2022-09-24T16:07:27.598435Z","shell.execute_reply":"2022-09-24T16:07:28.040694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20, 16))\n\nsns.heatmap(cm[(cm >= 0.1) | (cm <= -0.1)], \n cmap='viridis', vmax=1.0, vmin=-1.0, linewidths=0.1,\n annot=False, annot_kws={\"size\": 8}, square=True);","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:28.043394Z","iopub.execute_input":"2022-09-24T16:07:28.044453Z","iopub.status.idle":"2022-09-24T16:07:30.542605Z","shell.execute_reply.started":"2022-09-24T16:07:28.044398Z","shell.execute_reply":"2022-09-24T16:07:30.541547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plt.figure(figsize = (20,14) )\n# sns.heatmap(cm)","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:30.544299Z","iopub.execute_input":"2022-09-24T16:07:30.545101Z","iopub.status.idle":"2022-09-24T16:07:30.549063Z","shell.execute_reply.started":"2022-09-24T16:07:30.545074Z","shell.execute_reply":"2022-09-24T16:07:30.547753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Avoid nan columns in cm, otherwise clustermap will crash\ns = cm.isna().sum(axis = 1)\nm = (s < cm.shape[0]) \nl = list(cm.columns[m])\n\nplt.figure(figsize=(20,12) )\nsns.clustermap(cm.loc[l,l])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:30.550224Z","iopub.execute_input":"2022-09-24T16:07:30.551388Z","iopub.status.idle":"2022-09-24T16:07:31.418465Z","shell.execute_reply.started":"2022-09-24T16:07:30.551352Z","shell.execute_reply":"2022-09-24T16:07:31.417509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Correlation analysis with \"igraph\" - show components connected by correlation greater than threshold. ","metadata":{}},{"cell_type":"code","source":"import igraph\nlist_X_column_names = list_columns_ids # = list( cm.columns )\n\nn_components2show_members = 7\n\nverbose = 0\ndf_stat = pd.DataFrame() # dict_save_largest_component_size = {} \ni = 0\nfor correlation_threshold in [0.8, 0.7, 0.6, 0.5, 0.4] :\n t0 = time.time()\n print()\n print(correlation_threshold , 'correlation_threshold ')\n corr_matr_abs_bool = np.abs(cm) > correlation_threshold\n corr_matr_abs_bool = corr_matr_abs_bool # Restrict to genes part \n corr_matr_abs_bool = np.triu(corr_matr_abs_bool,1) # Take upper triangular part \n g = igraph.Graph().Adjacency(corr_matr_abs_bool.tolist())\n g.to_undirected(mode = 'collapse')\n if verbose >= 10:\n print( corr_matr_abs_bool.astype(int) )\n print('Number of nodes ', g.vcount())\n print('Number of edges ', g.ecount() )\n print('Number of weakly connected compoenents', len( g.connected_components(mode='WEAK')))\n\n\n list_clusters_nodes_lists = list( g.connected_components(mode='WEAK') )\n list_clusers_size = [len(t) for t in list_clusters_nodes_lists ]\n list_clusers_size = np.sort(list_clusers_size)[::-1]\n print('Top 5 cluster sizes:', list_clusers_size[:5] , 'seconds passed:', np.round(time.time()-t0 , 2))\n #dict_save_largest_component_size[correlation_threshold ] = list_clusers_size[0]\n list_clusers_size_unique = np.sort(np.unique(list_clusers_size))[::-1]\n cc = 0\n print('Not more than 50 Ids in largest correlated components:')\n for i_cluser in range(np.min( [n_components2show_members, len(list_clusers_size_unique )] )): # list_clusers_size:\n for t in list_clusters_nodes_lists:\n if len(t) == list_clusers_size_unique[i_cluser]:\n if len(t) == 1: continue \n cc += 1\n if cc >= n_components2show_members: continue # Show not more than n_components2show_members components \n print(str(cc)+'-th','component:', np.array(list_X_column_names)[t[:50]])\n i += 1\n df_stat.loc[i,'correlation threshold'] = correlation_threshold\n for t in range(n_components2show_members ):\n if t == 0:\n df_stat.loc[i,'The first Largest Component Size'] = list_clusers_size[t]\n elif len(list_clusers_size) > t:\n df_stat.loc[i,str(t)+'-th Component Size'] = list_clusers_size[t]\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:31.419594Z","iopub.execute_input":"2022-09-24T16:07:31.419974Z","iopub.status.idle":"2022-09-24T16:07:31.622419Z","shell.execute_reply.started":"2022-09-24T16:07:31.419942Z","shell.execute_reply":"2022-09-24T16:07:31.620892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualization of correlated components of correlation graphs","metadata":{}},{"cell_type":"code","source":"import igraph\nlist_X_column_names = list_columns_ids # = list( cm.columns )\n\nn_components2show_members = 7\n\nverbose = 0\ndf_stat = pd.DataFrame() # dict_save_largest_component_size = {} \ni = 0\nfor correlation_threshold in [0.7,0.6, 0.5]:# 0.8, 0.7, 0.6, 0.5, 0.4] :\n t0 = time.time()\n print()\n print(correlation_threshold , 'correlation_threshold ')\n corr_matr_abs_bool = np.abs(cm) > correlation_threshold\n corr_matr_abs_bool = corr_matr_abs_bool # Restrict to genes part \n corr_matr_abs_bool = np.triu(corr_matr_abs_bool,1) # Take upper triangular part \n g = igraph.Graph().Adjacency(corr_matr_abs_bool.tolist())\n g.vs['label'] = list_columns_ids\n g.to_undirected(mode = 'collapse')\n if verbose >= 10:\n print( corr_matr_abs_bool.astype(int) )\n print('Number of nodes ', g.vcount())\n print('Number of edges ', g.ecount() )\n print('Number of weakly connected compoenents', len( g.connected_components(mode='WEAK')))\n \n connected_comps = g.connected_components(mode='WEAK')\n icc = 0\n for cc in connected_comps :\n if len(cc) <=2: continue\n subgraph = g.subgraph(cc)\n icc += 1\n cn = subgraph.clique_number() \n print('Connected component',icc,'Size', len(cc), 'Size of largest clique:', cn )\n \n ss = set()\n if 1: # https://igraph.org/python/api/latest/igraph.GraphBase.html#largest_cliques\n list_l = subgraph.largest_cliques()\n for ilc, lc in enumerate( list_l):\n print(str(ilc)+'-th', 'Largest Clique', np.array(subgraph.vs['label'])[list(lc)] )\n ss = ss | set( np.array(subgraph.vs['label'])[list(lc)] )\n if 1:# Find and output maximal clique , non intersecting with largest found above \n if cn>=4:\n list_c = subgraph.maximal_cliques(4, cn)\n for ilc, lc in enumerate( list_c):\n ttt = np.array(subgraph.vs['label'])[list(lc)]\n if len(set(ttt) & ss ) == 0 :\n print( 'Non intersecting Maximal Clique:', ttt )\n ss = ss | set(ttt)\n \n tdf = pd.DataFrame(); ii = 0 \n for (v1,v2) in subgraph.get_edgelist():\n tdf.loc[ii, 'Node1'] = subgraph.vs['label'][v1]\n tdf.loc[ii, 'Node2'] = subgraph.vs['label'][v2]\n tdf.loc[ii, 'Weight'] = np.round( cm.loc[ tdf.loc[ii, 'Node1'], tdf.loc[ii, 'Node2']],2)\n ii += 1\n #tdf \n subgraph_weights_saved = igraph.Graph.TupleList(tdf.itertuples(index=False), directed=False, weights=False, edge_attrs=[\"Weight\"])\n visual_style = {}\n visual_style[\"vertex_color\"] = ['white' for t in range(len(cc))] # ,'gray','green']\n #visual_style[\"vertex_label\"] = np.array(list_columns_ids)[ cc ]\n visual_style[\"edge_label\"] = list(subgraph_weights_saved.es['Weight'])\n #visual_style[\"edge_color\"] = list(g3.es['Color'])\n visual_style[\"vertex_size\"] = 0.3\n visual_style[\"edge_size\"] = 0.2\n \n fig, ax = plt.subplots(1,1 ,figsize = (20,10))\n igraph.plot(\n subgraph,\n mark_groups=True, palette=igraph.RainbowPalette(),\n edge_width=0.5,\n target=ax, **visual_style, \n )\n plt.axis('off')\n plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:31.623957Z","iopub.execute_input":"2022-09-24T16:07:31.624454Z","iopub.status.idle":"2022-09-24T16:07:33.524812Z","shell.execute_reply.started":"2022-09-24T16:07:31.624425Z","shell.execute_reply":"2022-09-24T16:07:33.524052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look on missing values, df2 - copy with filled missings ","metadata":{}},{"cell_type":"code","source":"df.isnull().sum().sum()/ (df.shape[0]*df.shape[1])","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:33.525910Z","iopub.execute_input":"2022-09-24T16:07:33.526967Z","iopub.status.idle":"2022-09-24T16:07:33.549746Z","shell.execute_reply.started":"2022-09-24T16:07:33.526934Z","shell.execute_reply":"2022-09-24T16:07:33.549033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df2 = df.copy()\nfor c in df2.columns:\n df2[c] =df[c].fillna(df[c].mean())\n ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:33.550769Z","iopub.execute_input":"2022-09-24T16:07:33.551785Z","iopub.status.idle":"2022-09-24T16:07:33.646462Z","shell.execute_reply.started":"2022-09-24T16:07:33.551754Z","shell.execute_reply":"2022-09-24T16:07:33.645024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA dimensional reduction and visualizations\n\nFeatures have largest std typicall the most contributing to PCA compenents ","metadata":{"execution":{"iopub.status.busy":"2022-09-20T09:16:44.970255Z","iopub.execute_input":"2022-09-20T09:16:44.971571Z","iopub.status.idle":"2022-09-20T09:16:44.977180Z","shell.execute_reply.started":"2022-09-20T09:16:44.971509Z","shell.execute_reply":"2022-09-20T09:16:44.975744Z"}}},{"cell_type":"code","source":"%%time\n\nimport numpy as np\nfrom sklearn.decomposition import PCA\n\nX = df2\npca = PCA( n_components=100)\nr = pca.fit_transform(X)\nprint(r.shape)\n\nfig = plt.figure(figsize = (15,7) )\nfig.add_subplot(1, 2, 1) \nplt.plot(pca.singular_values_,'o-')\nplt.title('Singular values')\nfig.add_subplot(1, 2, 2) \nplt.plot(pca.explained_variance_ratio_,'o-')\nplt.title('explained variance')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:33.648239Z","iopub.execute_input":"2022-09-24T16:07:33.648562Z","iopub.status.idle":"2022-09-24T16:07:35.643031Z","shell.execute_reply.started":"2022-09-24T16:07:33.648537Z","shell.execute_reply":"2022-09-24T16:07:35.642042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list(df.columns) == list_columns_ids","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:35.644417Z","iopub.execute_input":"2022-09-24T16:07:35.645704Z","iopub.status.idle":"2022-09-24T16:07:35.652370Z","shell.execute_reply.started":"2022-09-24T16:07:35.645653Z","shell.execute_reply":"2022-09-24T16:07:35.651488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(pca.components_.shape)\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = pca.components_[i,:].ravel()\n plt.plot(v,label = 'PCA'+str(i+1))\n plt.xlabel('column index', fontsize = 20)\n plt.ylabel('PCA component value', fontsize = 20)\n plt.title(str_data_inf + '. PCA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = pca.components_[i,:].ravel()\n sns.kdeplot(v,label = 'PCA'+str(i+1))\n #plt.xlabel('column index', fontsize = 20)\n #plt.ylabel('PCA component value', fontsize = 20)\n plt.title(str_data_inf + '. PCA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\ndpca = {}\nfor i in range(5):\n v = pca.components_[i,:].ravel()\n t = np.percentile(np.abs(v),90)\n IX = np.where( np.abs(v) > t)[0]\n print(len(IX))\n l = np.array(list_columns_ids)[IX]\n print(l, 'Top columns contributing to PCA'+str(i) )\n dpca[i] = l\n if i == 0:\n s = set(dpca[0])\n else:\n s = s & set(dpca[i])\n print(s, 'Intersection',i) \nplt.show()\n \nprint()\ns = set(dpca[0])\nfor i in dpca:\n s = s & set(dpca[i])\n print(s, 'Intersection',i)","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:35.653489Z","iopub.execute_input":"2022-09-24T16:07:35.654694Z","iopub.status.idle":"2022-09-24T16:07:36.128983Z","shell.execute_reply.started":"2022-09-24T16:07:35.654632Z","shell.execute_reply":"2022-09-24T16:07:36.127616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v4color = pd.Series(r[:,0], name = 'PCA1' )\nv4color = pd.Series( df.sum(axis=1).values, name = 'Sum' ) # pd.Series(r[:,0], name = 'PCA1' )\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:36.130129Z","iopub.execute_input":"2022-09-24T16:07:36.130426Z","iopub.status.idle":"2022-09-24T16:07:36.151133Z","shell.execute_reply.started":"2022-09-24T16:07:36.130402Z","shell.execute_reply":"2022-09-24T16:07:36.150056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nn_x_subplots = 3\nc = 0\nfor (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4) , (0,5), (1,5)]:\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(str_data_inf) + ' n_samples='+str(len(r)),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED \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 = 'rainbow' )\n if 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.title( ' n_cells='+str(len(r)), fontsize = 20)\n plt.xlabel('PCA'+str(i+1), fontsize = 20)\n plt.ylabel('PCA'+str(j+1), fontsize = 20) \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:07:36.152223Z","iopub.execute_input":"2022-09-24T16:07:36.152454Z","iopub.status.idle":"2022-09-24T16:08:07.780285Z","shell.execute_reply.started":"2022-09-24T16:07:36.152431Z","shell.execute_reply":"2022-09-24T16:08:07.779025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA with preprocessing and scaling ","metadata":{}},{"cell_type":"code","source":"%%time\nX2 = df2.copy()\nv1 = X2.quantile(0.95, axis=0)\nv2 = X2.quantile(0.05, axis=0)\n\n# cut off outliers:\nfor col in X2.columns:\n X2[col] = np.clip(X2[col], a_min = v2[col], a_max = v1[col])\n\n# Take top most variable \nv = X2.std(axis = 0) / X2.abs().mean(axis=0) \nIX = v.sort_values(ascending = False).iloc[:100]\nX2 = X2[IX.index]\n\n\nfrom sklearn.preprocessing import StandardScaler\ndata = [[0, 0], [0, 0], [1, 1], [1, 1]]\nscaler = StandardScaler()\nX2 = scaler.fit_transform(X2)\n\nlist_columns_ids_preprocessed = list(IX.index)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:08:07.782056Z","iopub.execute_input":"2022-09-24T16:08:07.782437Z","iopub.status.idle":"2022-09-24T16:08:08.409396Z","shell.execute_reply.started":"2022-09-24T16:08:07.782380Z","shell.execute_reply":"2022-09-24T16:08:08.407884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nimport numpy as np\nfrom sklearn.decomposition import PCA\n\npca = PCA( n_components=100)\nr = pca.fit_transform(X2)\nprint(r.shape)\n\nfig = plt.figure(figsize = (15,7) )\nfig.add_subplot(1, 2, 1) \nplt.plot(pca.singular_values_,'o-')\nplt.title('Singular values')\nfig.add_subplot(1, 2, 2) \nplt.plot(pca.explained_variance_ratio_,'o-')\nplt.title('explained variance')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:08:08.410975Z","iopub.execute_input":"2022-09-24T16:08:08.411450Z","iopub.status.idle":"2022-09-24T16:08:08.935887Z","shell.execute_reply.started":"2022-09-24T16:08:08.411415Z","shell.execute_reply":"2022-09-24T16:08:08.934010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_columns_ids_loc = list_columns_ids_preprocessed.copy()\n\nprint(pca.components_.shape)\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = pca.components_[i,:].ravel()\n plt.plot(v,label = 'PCA'+str(i+1))\n plt.xlabel('column index', fontsize = 20)\n plt.ylabel('PCA component value', fontsize = 20)\n plt.title(str_data_inf + '. PCA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = pca.components_[i,:].ravel()\n sns.kdeplot(v,label = 'PCA'+str(i+1))\n #plt.xlabel('column index', fontsize = 20)\n #plt.ylabel('PCA component value', fontsize = 20)\n plt.title(str_data_inf + '. PCA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\ndpca = {}\nfor i in range(5):\n v = pca.components_[i,:].ravel()\n t = np.percentile(np.abs(v),90)\n IX = np.where( np.abs(v) > t)[0]\n print(len(IX))\n l = np.array(list_columns_ids_loc) [IX]\n print(l, 'Top columns contributing to PCA'+str(i) )\n dpca[i] = l\n if i == 0:\n s = set(dpca[0])\n else:\n s = s & set(dpca[i])\n print(s, 'Intersection',i) \nplt.show()\n\nprint()\ns = set(dpca[0])\nfor i in dpca:\n s = s & set(dpca[i])\n print(s, 'Intersection',i)","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:08:08.938757Z","iopub.execute_input":"2022-09-24T16:08:08.939024Z","iopub.status.idle":"2022-09-24T16:08:09.421047Z","shell.execute_reply.started":"2022-09-24T16:08:08.939000Z","shell.execute_reply":"2022-09-24T16:08:09.419866Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"v4color = pd.Series(r[:,0], name = 'PCA1' )\nv4color = pd.Series( df.sum(axis=1).values, name = 'Sum' ) # pd.Series(r[:,0], name = 'PCA1' )\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:08:09.422387Z","iopub.execute_input":"2022-09-24T16:08:09.422939Z","iopub.status.idle":"2022-09-24T16:08:09.444100Z","shell.execute_reply.started":"2022-09-24T16:08:09.422911Z","shell.execute_reply":"2022-09-24T16:08:09.442824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nn_x_subplots = 3\nc = 0\nfor (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4) , (0,5), (1,5)]:\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('PCA on preprocessed data '+ str(str_data_inf) + ' n_samples='+str(len(r)),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED \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 = 'rainbow' )\n if 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.title( ' n_cells='+str(len(r)), fontsize = 20)\n plt.xlabel('PCA'+str(i+1), fontsize = 20)\n plt.ylabel('PCA'+str(j+1), fontsize = 20) \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:08:09.445401Z","iopub.execute_input":"2022-09-24T16:08:09.445725Z","iopub.status.idle":"2022-09-24T16:08:46.048464Z","shell.execute_reply.started":"2022-09-24T16:08:09.445699Z","shell.execute_reply":"2022-09-24T16:08:46.047229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# UMAP ","metadata":{}},{"cell_type":"code","source":"%%time\nimport umap\nreducer = umap.UMAP()\nr2 = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y = r2[:,j], hue = v4color, palette = 'rainbow')\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('UMAP ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:08:46.050165Z","iopub.execute_input":"2022-09-24T16:08:46.050444Z","iopub.status.idle":"2022-09-24T16:10:14.796656Z","shell.execute_reply.started":"2022-09-24T16:08:46.050419Z","shell.execute_reply":"2022-09-24T16:10:14.795589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport umap\nreducer = umap.UMAP(densmap=True, random_state=42)\nr2 = reducer.fit_transform(X)\nprint(X.shape)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y = r2[:,j], hue = v4color, palette = 'rainbow')\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('UMAP densmap On ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:10:14.798622Z","iopub.execute_input":"2022-09-24T16:10:14.799339Z","iopub.status.idle":"2022-09-24T16:12:49.151718Z","shell.execute_reply.started":"2022-09-24T16:10:14.799306Z","shell.execute_reply":"2022-09-24T16:12:49.150198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reducer = umap.UMAP()\nreducer.min_dist, reducer.n_neighbors\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:12:49.152994Z","iopub.execute_input":"2022-09-24T16:12:49.153785Z","iopub.status.idle":"2022-09-24T16:12:49.364079Z","shell.execute_reply.started":"2022-09-24T16:12:49.153757Z","shell.execute_reply":"2022-09-24T16:12:49.361814Z"},"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\nplt.suptitle('UMAP ' + str_data_inf+ ' n_samples='+str(len(r2)) , 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 = 'rainbow')\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)\nplt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:12:49.366070Z","iopub.execute_input":"2022-09-24T16:12:49.366494Z","iopub.status.idle":"2022-09-24T16:19:03.841709Z","shell.execute_reply.started":"2022-09-24T16:12:49.366459Z","shell.execute_reply":"2022-09-24T16:19:03.839932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport umap\nreducer = umap.UMAP(n_components = 4)\nr2 = reducer.fit_transform(X)\n\nfig = plt.figure(figsize = (20,16)); c = 0\nplt.suptitle('UMAP ' + str_data_inf+ ' n_samples='+str(len(r2)) , fontsize = 20 )\n\nfor (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3)]: # , (0,4),(1,4),(2,4),(3,4)]:\n c += 1; fig.add_subplot(2,3,c)\n \n ax = sns.scatterplot(x = r2[:,i],y = r2[:,j], hue = v4color, palette = 'rainbow')\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('UMAP ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n plt.ylabel('UMAP'+str(j+1), fontsize = 20)\nplt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:19:03.843691Z","iopub.execute_input":"2022-09-24T16:19:03.844080Z","iopub.status.idle":"2022-09-24T16:20:04.235583Z","shell.execute_reply.started":"2022-09-24T16:19:03.844044Z","shell.execute_reply":"2022-09-24T16:20:04.234557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport umap\nreducer = umap.UMAP()\nr2 = reducer.fit_transform(X2)\nprint(X2.shape)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y = r2[:,j], hue = v4color, palette = 'rainbow')\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('UMAP on preprocessed data ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:20:04.237088Z","iopub.execute_input":"2022-09-24T16:20:04.237358Z","iopub.status.idle":"2022-09-24T16:21:00.565724Z","shell.execute_reply.started":"2022-09-24T16:20:04.237333Z","shell.execute_reply":"2022-09-24T16:21:00.564684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport umap\nreducer = umap.UMAP(densmap=True, random_state=42)\nr2 = reducer.fit_transform(X2)\nprint(X2.shape)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y = r2[:,j], hue = v4color, palette = 'rainbow')\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('UMAP densmap On preprocessed ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('UMAP'+str(i+1), fontsize = 20)\n plt.ylabel('UMAP'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:21:00.567238Z","iopub.execute_input":"2022-09-24T16:21:00.567491Z","iopub.status.idle":"2022-09-24T16:23:47.513608Z","shell.execute_reply.started":"2022-09-24T16:21:00.567465Z","shell.execute_reply":"2022-09-24T16:23:47.512829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# NCVis - similar to UMAP, bust faster","metadata":{}},{"cell_type":"code","source":"!pip install ncvis","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:23:47.514743Z","iopub.execute_input":"2022-09-24T16:23:47.515959Z","iopub.status.idle":"2022-09-24T16:24:17.363206Z","shell.execute_reply.started":"2022-09-24T16:23:47.515930Z","shell.execute_reply":"2022-09-24T16:24:17.361985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(X.shape)\nimport ncvis\nreducer = ncvis.NCVis() # umap.UMAP()\nr2 = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y=r2[:,j], hue = v4color, palette = 'rainbow')\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('NCIVIS ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('NCIVIS'+str(i+1), fontsize = 20)\n plt.ylabel('NCIVIS'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:24:17.364666Z","iopub.execute_input":"2022-09-24T16:24:17.365107Z","iopub.status.idle":"2022-09-24T16:24:54.956169Z","shell.execute_reply.started":"2022-09-24T16:24:17.365063Z","shell.execute_reply":"2022-09-24T16:24:54.954835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport ncvis\nreducer = ncvis.NCVis() # umap.UMAP()\nr2 = reducer.fit_transform(X2)\nprint(X2.shape)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y=r2[:,j], hue = v4color, palette = 'rainbow')\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('NCIVIS on preprocessed data ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('NCIVIS'+str(i+1), fontsize = 20)\n plt.ylabel('NCIVIS'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:24:54.959363Z","iopub.execute_input":"2022-09-24T16:24:54.959898Z","iopub.status.idle":"2022-09-24T16:25:34.964391Z","shell.execute_reply.started":"2022-09-24T16:24:54.959861Z","shell.execute_reply":"2022-09-24T16:25:34.963347Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Trimap","metadata":{}},{"cell_type":"code","source":"!pip install trimap \n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:25:34.966733Z","iopub.execute_input":"2022-09-24T16:25:34.967091Z","iopub.status.idle":"2022-09-24T16:25:44.710497Z","shell.execute_reply.started":"2022-09-24T16:25:34.967061Z","shell.execute_reply":"2022-09-24T16:25:44.708862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nprint(X.shape)\nimport trimap\nreducer = trimap.TRIMAP()# ncvis.NCVis() # umap.UMAP()\nif isinstance(X, pd.DataFrame):\n r2 = reducer.fit_transform(X.values)\nelse:\n r2 = reducer.fit_transform(X)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y=r2[:,j], hue = v4color, palette = 'rainbow')\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(r2)), fontsize = 20)\n plt.xlabel('Trimap'+str(i+1), fontsize = 20)\n plt.ylabel('Trimap'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:25:44.712654Z","iopub.execute_input":"2022-09-24T16:25:44.713100Z","iopub.status.idle":"2022-09-24T16:27:03.356476Z","shell.execute_reply.started":"2022-09-24T16:25:44.713014Z","shell.execute_reply":"2022-09-24T16:27:03.355391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nimport trimap\nreducer = trimap.TRIMAP()# ncvis.NCVis() # umap.UMAP()\nr2 = reducer.fit_transform(X2)\nprint(X2.shape)\n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y=r2[:,j], hue = v4color, palette = 'rainbow')\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 on preprocessed data ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('Trimap'+str(i+1), fontsize = 20)\n plt.ylabel('Trimap'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:27:03.358009Z","iopub.execute_input":"2022-09-24T16:27:03.358424Z","iopub.status.idle":"2022-09-24T16:28:13.397771Z","shell.execute_reply.started":"2022-09-24T16:27:03.358390Z","shell.execute_reply":"2022-09-24T16:28:13.396818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# TSNE for restrictred number of samples - to speed up","metadata":{}},{"cell_type":"code","source":"%%time\nN = 2000\nprint(X.shape)\nfrom sklearn import manifold\nreducer = manifold.TSNE(n_components=2, init='pca', random_state=0) # reducer = trimap.TRIMAP()# ncvis.NCVis() # umap.UMAP()\n\nif isinstance(X, pd.DataFrame):\n r2 = reducer.fit_transform(X.iloc[:N,:])\nelse:\n r2 = reducer.fit_transform(X[:N,:])\nif isinstance(v4color, pd.Series): \n v4color2 = v4color.iloc[:N]\nelse:\n v4color2 = v4color[:N]\n \n\nfor (i,j) in [(0,1)]:# ,(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4)]:\n plt.figure(figsize = (20,10))\n ax = sns.scatterplot(x = r2[:,i],y=r2[:,j], hue = v4color2, palette = 'rainbow')\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('TSNE ' + str_data_inf + ' n_samples='+str(len(r2)), fontsize = 20)\n plt.xlabel('TSNE'+str(i+1), fontsize = 20)\n plt.ylabel('TSNE'+str(j+1), fontsize = 20)\n plt.show() ","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:28:13.399016Z","iopub.execute_input":"2022-09-24T16:28:13.399943Z","iopub.status.idle":"2022-09-24T16:28:30.891499Z","shell.execute_reply.started":"2022-09-24T16:28:13.399916Z","shell.execute_reply":"2022-09-24T16:28:30.890030Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Several other (fast) dimensional reduction methods from sklearn","metadata":{}},{"cell_type":"code","source":"# Based on: \n# https://scikit-learn.org/stable/auto_examples/manifold/plot_compare_methods.html#sphx-glr-auto-examples-manifold-plot-compare-methods-py\n# See also:\n# https://scikit-learn.org/stable/auto_examples/manifold/plot_lle_digits.html#\n\n\n\n# To speed-up reduce dimensions by PCA first\n# X_save = X.copy( )\n#r = pca.fit_transform(X)\n#X = r[:1000,:20]\n\n\n\nimport umap \nfrom sklearn import manifold\nfrom sklearn.decomposition import PCA\nfrom sklearn.decomposition import FactorAnalysis\nfrom sklearn.decomposition import NMF\nfrom sklearn.decomposition import FastICA\nfrom sklearn.decomposition import FactorAnalysis\nfrom sklearn.decomposition import LatentDirichletAllocation\nfrom sklearn.ensemble import RandomTreesEmbedding\nfrom sklearn.random_projection import SparseRandomProjection\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.decomposition import TruncatedSVD\n\n\nfrom collections import OrderedDict\nfrom functools import partial\nfrom matplotlib.ticker import NullFormatter\n\n\nn_neighbors = 10\nn_components = 2\n# Set-up manifold methods\nLLE = partial(manifold.LocallyLinearEmbedding,\n n_neighbors, n_components, eigen_solver='auto')\n\nmethods = OrderedDict()\nmethods['PCA'] = PCA()\nmethods['umap'] = umap.UMAP(n_components = n_components)\nmethods['t-SNE'] = manifold.TSNE(n_components=n_components, init='pca', random_state=0)\nmethods['ICA'] = FastICA(n_components=n_components, random_state=0)\nmethods['FA'] = FactorAnalysis(n_components=n_components, random_state=0)\n#methods['LLE'] = LLE(method='standard')\n#methods['Modified LLE'] = LLE(method='modified')\n#methods['Isomap'] = manifold.Isomap(n_neighbors, n_components)\nmethods['MDS'] = manifold.MDS(n_components, max_iter=100, n_init=1)\nmethods['SE'] = manifold.SpectralEmbedding(n_components=n_components,\n n_neighbors=n_neighbors)\nmethods['NMF'] = NMF(n_components=n_components, init='random', random_state=0) \nmethods['RandProj'] = SparseRandomProjection(n_components=n_components, random_state=42)\n\nrand_trees_embed = make_pipeline(RandomTreesEmbedding(n_estimators=200, random_state=0, max_depth=5), TruncatedSVD(n_components=n_components) )\nmethods['RandTrees'] = rand_trees_embed\nmethods['LatDirAll'] = LatentDirichletAllocation(n_components=n_components, random_state=0)\n#methods['LTSA'] = LLE(method='ltsa') \n#methods['Hessian LLE'] = LLE(method='hessian') \n\nlist_fast_methods = ['PCA','umap','FA', 'NMF','RandProj','RandTrees'] # 'ICA',\nlist_slow_methods = ['t-SNE','LLE','Modified LLE','Isomap','MDS','SE','LatDirAll','LTSA','Hessian LLE']\n\n# transformer = NeighborhoodComponentsAnalysis(init='random', n_components=2, random_state=0) # Cannot be applied since supervised - requires y \n# methods['LinDisA'] = LinearDiscriminantAnalysis(n_components=n_components)# Cannot be applied since supervised - requires y \n\n\n# Create figure\nfig = plt.figure(figsize=(25, 16))\n\n# Plot results\nc = 0\nfor i, (label, method) in enumerate(methods.items()):\n if label not in list_fast_methods: # list_slow_methods :\n continue\n \n t0 = time.time()\n try:\n r = method.fit_transform(X)\n except:\n print('Got Exception', label )\n continue \n t1 = time.time()\n print(\"%s: %.2g sec\" % (label, t1 - t0))\n c+=1\n fig.add_subplot(2, 3 , c) \n sns.scatterplot(x=r[:,0], y=r[:,1] , hue = v4color)\n plt.title(label )\n plt.legend('')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:28:30.893428Z","iopub.execute_input":"2022-09-24T16:28:30.893741Z","iopub.status.idle":"2022-09-24T16:29:52.035554Z","shell.execute_reply.started":"2022-09-24T16:28:30.893714Z","shell.execute_reply":"2022-09-24T16:29:52.034429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Other (slow) methods - consider only first N samples","metadata":{}},{"cell_type":"code","source":"N = 1000 # To speed up consider just first N samples \n\n# Based on: \n# https://scikit-learn.org/stable/auto_examples/manifold/plot_compare_methods.html#sphx-glr-auto-examples-manifold-plot-compare-methods-py\n# See also:\n# https://scikit-learn.org/stable/auto_examples/manifold/plot_lle_digits.html#\n\n\n\n# To speed-up reduce dimensions by PCA first\n# X_save = X.copy( )\n#r = pca.fit_transform(X)\n#X = r[:1000,:20]\n\n\n\nimport umap \nfrom sklearn import manifold\nfrom sklearn.decomposition import PCA\nfrom sklearn.decomposition import FactorAnalysis\nfrom sklearn.decomposition import NMF\nfrom sklearn.decomposition import FastICA\nfrom sklearn.decomposition import FactorAnalysis\nfrom sklearn.decomposition import LatentDirichletAllocation\nfrom sklearn.ensemble import RandomTreesEmbedding\nfrom sklearn.random_projection import SparseRandomProjection\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.decomposition import TruncatedSVD\n\n\nfrom collections import OrderedDict\nfrom functools import partial\nfrom matplotlib.ticker import NullFormatter\n\n\nn_neighbors = 10\nn_components = 2\n# Set-up manifold methods\nLLE = partial(manifold.LocallyLinearEmbedding,\n n_neighbors, n_components, eigen_solver='auto')\n\nmethods = OrderedDict()\nmethods['PCA'] = PCA()\nmethods['umap'] = umap.UMAP(n_components = n_components)\nmethods['t-SNE'] = manifold.TSNE(n_components=n_components, init='pca', random_state=0)\nmethods['ICA'] = FastICA(n_components=n_components, random_state=0)\nmethods['FA'] = FactorAnalysis(n_components=n_components, random_state=0)\n#methods['LLE'] = LLE(method='standard')\n#methods['Modified LLE'] = LLE(method='modified')\n#methods['Isomap'] = manifold.Isomap(n_neighbors, n_components)\nmethods['MDS'] = manifold.MDS(n_components, max_iter=100, n_init=1)\nmethods['SE'] = manifold.SpectralEmbedding(n_components=n_components,\n n_neighbors=n_neighbors)\nmethods['NMF'] = NMF(n_components=n_components, init='random', random_state=0) \nmethods['RandProj'] = SparseRandomProjection(n_components=n_components, random_state=42)\n\nrand_trees_embed = make_pipeline(RandomTreesEmbedding(n_estimators=200, random_state=0, max_depth=5), TruncatedSVD(n_components=n_components) )\nmethods['RandTrees'] = rand_trees_embed\nmethods['LatDirAll'] = LatentDirichletAllocation(n_components=n_components, random_state=0)\n#methods['LTSA'] = LLE(method='ltsa') \n#methods['Hessian LLE'] = LLE(method='hessian') \n\nlist_fast_methods = ['PCA','umap','FA', 'NMF','RandProj','RandTrees'] # 'ICA',\nlist_slow_methods = ['t-SNE','LLE','Modified LLE','Isomap','MDS','SE','LatDirAll','LTSA','Hessian LLE']\n\n# transformer = NeighborhoodComponentsAnalysis(init='random', n_components=2, random_state=0) # Cannot be applied since supervised - requires y \n# methods['LinDisA'] = LinearDiscriminantAnalysis(n_components=n_components)# Cannot be applied since supervised - requires y \n\n\n# Create figure\nfig = plt.figure(figsize=(25, 16))\n\n# Plot results\nc = 0\nfor i, (label, method) in enumerate(methods.items()):\n if label not in list_slow_methods :# list_fast_methods: # \n continue\n \n t0 = time.time()\n try:\n if isinstance(X, pd.DataFrame):\n r = method.fit_transform(X.iloc[:N,:])\n else:\n r = method.fit_transform(X[:N,:])\n except:\n print('Got Exception', label )\n continue \n t1 = time.time()\n print(\"%s: %.2g sec\" % (label, t1 - t0))\n c+=1\n fig.add_subplot(2, 3 , c) \n if isinstance(v4color, pd.Series):\n sns.scatterplot(x=r[:,0], y=r[:,1] , hue = v4color.iloc[:N])\n else:\n sns.scatterplot(x=r[:,0], y=r[:,1] , hue = v4color[:N])\n plt.title(label )\n plt.legend('')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:29:52.037300Z","iopub.execute_input":"2022-09-24T16:29:52.037580Z","iopub.status.idle":"2022-09-24T16:30:06.531943Z","shell.execute_reply.started":"2022-09-24T16:29:52.037555Z","shell.execute_reply":"2022-09-24T16:30:06.530493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ICA for dimensional reduction and analysis","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import FastICA","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:30:06.533805Z","iopub.execute_input":"2022-09-24T16:30:06.534121Z","iopub.status.idle":"2022-09-24T16:30:06.539024Z","shell.execute_reply.started":"2022-09-24T16:30:06.534091Z","shell.execute_reply":"2022-09-24T16:30:06.537835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nreducer = FastICA(n_components=6, random_state=0, whiten='unit-variance')\nr2 = reducer.fit_transform(X)\nprint(X.shape, r2.shape)\nprint(reducer.components_.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:30:06.540369Z","iopub.execute_input":"2022-09-24T16:30:06.540582Z","iopub.status.idle":"2022-09-24T16:30:10.825362Z","shell.execute_reply.started":"2022-09-24T16:30:06.540558Z","shell.execute_reply":"2022-09-24T16:30:10.824028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nn_x_subplots = 3\nc = 0\nfor (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4) , (0,5), (1,5)]:\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('ICA on preprocessed data '+ str(str_data_inf) + ' n_samples='+str(len(r)),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED \n c += 1; fig.add_subplot(1,n_x_subplots ,c) \n \n ax = sns.scatterplot(x=r2[:,i], y=r2[:,j], hue = v4color, palette = 'rainbow' )\n if 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.title( ' n_cells='+str(len(r)), fontsize = 20)\n plt.xlabel('ICA'+str(i+1), fontsize = 20)\n plt.ylabel('ICA'+str(j+1), fontsize = 20) \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:30:10.828023Z","iopub.execute_input":"2022-09-24T16:30:10.828268Z","iopub.status.idle":"2022-09-24T16:30:43.783815Z","shell.execute_reply.started":"2022-09-24T16:30:10.828244Z","shell.execute_reply":"2022-09-24T16:30:43.781685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_columns_ids_loc = list_columns_ids.copy()\nprint(reducer.components_.shape)\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = reducer.components_[i,:].ravel()\n plt.plot(v,label = 'ICA'+str(i+1))\n plt.xlabel('column index', fontsize = 20)\n plt.ylabel('ICA component value', fontsize = 20)\n plt.title(str_data_inf + '. ICA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = reducer.components_[i,:].ravel()\n sns.kdeplot(v,label = 'ICA'+str(i+1))\n #plt.xlabel('column index', fontsize = 20)\n #plt.ylabel('PCA component value', fontsize = 20)\n plt.title(str_data_inf + '. ICA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\ndpca = {}\nfor i in range(5):\n v = reducer.components_[i,:].ravel()\n t = np.percentile(np.abs(v),90)\n IX = np.where( np.abs(v) > t)[0]\n print(len(IX))\n l = np.array(list_columns_ids_loc)[IX]\n print(l, 'Top columns contributing to ICA'+str(i) )\n dpca[i] = l\n if i == 0:\n s = set(dpca[0])\n else:\n s = s & set(dpca[i])\n print(s, 'Intersection',i) \nplt.show()\n\nprint()\ns = set(dpca[0])\nfor i in dpca:\n s = s & set(dpca[i])\n print(s, 'Intersection',i)","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:30:43.786186Z","iopub.execute_input":"2022-09-24T16:30:43.786567Z","iopub.status.idle":"2022-09-24T16:30:44.457557Z","shell.execute_reply.started":"2022-09-24T16:30:43.786539Z","shell.execute_reply":"2022-09-24T16:30:44.456010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ICA on preprocessed data","metadata":{}},{"cell_type":"code","source":"%%time\nreducer = FastICA(n_components=6, random_state=0, whiten='unit-variance')\nr2 = reducer.fit_transform(X2)\nprint(X2.shape, r2.shape)\nprint(reducer.components_.shape)\n","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:30:44.459900Z","iopub.execute_input":"2022-09-24T16:30:44.460275Z","iopub.status.idle":"2022-09-24T16:30:47.124408Z","shell.execute_reply.started":"2022-09-24T16:30:44.460248Z","shell.execute_reply":"2022-09-24T16:30:47.121745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nn_x_subplots = 3\nc = 0\nfor (i,j) in [(0,1),(0,2),(1,2), (0,3),(1,3),(2,3), (0,4),(1,4),(2,4),(3,4) , (0,5), (1,5)]:\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('ICA on preprocessed data '+ str(str_data_inf) + ' n_samples='+str(len(r)),fontsize = 20 )# str_data_inf + ' n_cells: ' +str(mask.sum()) + ' RED \n c += 1; fig.add_subplot(1,n_x_subplots ,c) \n \n ax = sns.scatterplot(x=r2[:,i], y=r2[:,j], hue = v4color, palette = 'rainbow' )\n if 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.title( ' n_cells='+str(len(r)), fontsize = 20)\n plt.xlabel('ICA'+str(i+1), fontsize = 20)\n plt.ylabel('ICA'+str(j+1), fontsize = 20) \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:30:47.127561Z","iopub.execute_input":"2022-09-24T16:30:47.127909Z","iopub.status.idle":"2022-09-24T16:31:21.309610Z","shell.execute_reply.started":"2022-09-24T16:30:47.127882Z","shell.execute_reply":"2022-09-24T16:31:21.308173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list_columns_ids_loc = list_columns_ids_preprocessed.copy()\n\nprint(reducer.components_.shape)\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = reducer.components_[i,:].ravel()\n plt.plot(v,label = 'ICA'+str(i+1))\n plt.xlabel('column index', fontsize = 20)\n plt.ylabel('ICA component value', fontsize = 20)\n plt.title(str_data_inf + '. ICA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\nfor i in range(3):\n v = reducer.components_[i,:].ravel()\n sns.kdeplot(v,label = 'ICA'+str(i+1))\n #plt.xlabel('column index', fontsize = 20)\n #plt.ylabel('PCA component value', fontsize = 20)\n plt.title(str_data_inf + '. ICA components. n_samples='+str(len(df)), fontsize = 20)\nplt.legend( fontsize = 20)\nplt.show()\n\nfig = plt.figure(figsize= (20,6))\ndpca = {}\nfor i in range(5):\n v = reducer.components_[i,:].ravel()\n t = np.percentile(np.abs(v),90)\n IX = np.where( np.abs(v) > t)[0]\n print(len(IX))\n l = np.array(list_columns_ids_loc)[IX]\n print(l, 'Top columns contributing to ICA'+str(i) )\n dpca[i] = l\n if i == 0:\n s = set(dpca[0])\n else:\n s = s & set(dpca[i])\n print(s, 'Intersection',i) \nplt.show()\n\nprint()\ns = set(dpca[0])\nfor i in dpca:\n s = s & set(dpca[i])\n print(s, 'Intersection',i)","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:31:21.311742Z","iopub.execute_input":"2022-09-24T16:31:21.312065Z","iopub.status.idle":"2022-09-24T16:31:21.961565Z","shell.execute_reply.started":"2022-09-24T16:31:21.312036Z","shell.execute_reply":"2022-09-24T16:31:21.959914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Heatmap/imshow plots","metadata":{}},{"cell_type":"code","source":"\nfrom sklearn.preprocessing import MinMaxScaler\nscaler = MinMaxScaler()\n\nX2 = df2.copy()\nv1 = X2.quantile(0.95, axis=0)\nv2 = X2.quantile(0.05, axis=0)\n# cut off outliers:\nfor col in X2.columns:\n X2[col] = np.clip(X2[col], a_min = v2[col], a_max = v1[col])\n\nr = scaler.fit_transform(X2)\nprint(r.shape, X2.shape, df.shape)\n\nplt.figure(figsize = (20,12))\nplt.imshow(r, cmap = 'rainbow', aspect='auto', vmin=0, vmax=1,)\nplt.show()\nplt.figure(figsize = (20,12))\n#plt.imshow(train_inputs[:105942,I].toarray(), cmap='Blues', aspect='auto', vmin=0, vmax=1)\nplt.spy(r, aspect='auto', markersize=0.01) # vmin=0, vmax=1)\nplt.title(str_data_inf, fontsize = 20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:31:21.963955Z","iopub.execute_input":"2022-09-24T16:31:21.964296Z","iopub.status.idle":"2022-09-24T16:31:34.021169Z","shell.execute_reply.started":"2022-09-24T16:31:21.964268Z","shell.execute_reply":"2022-09-24T16:31:34.019889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('%.1f seconds passed total '%(time.time()-t0start) )","metadata":{"execution":{"iopub.status.busy":"2022-09-24T16:31:34.022715Z","iopub.execute_input":"2022-09-24T16:31:34.023187Z","iopub.status.idle":"2022-09-24T16:31:34.029118Z","shell.execute_reply.started":"2022-09-24T16:31:34.023152Z","shell.execute_reply":"2022-09-24T16:31:34.027877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}