{"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":"# MSCI Correlations EDA for Multiome","metadata":{}},{"cell_type":"markdown","source":"WARNING: It seems that some of my interpretations in this notebooks were wrong. Please check [this topic first](https://www.kaggle.com/competitions/open-problems-multimodal/discussion/351725).\n\nIn this notebook, we are going to analyze the correlations between the inputs and targets of the Multiome training set. There is a [corresponding notebook for the CITEseq data](https://www.kaggle.com/fabiencrom/msci-correlations-eda-citeseq).\n\nThe Pearsons correlation coefficients were computed with [this notebook](https://www.kaggle.com/fabiencrom/msci-generating-all-correlations-inputs-targets/edit) on a Saturn Cloud machine with 128GB memory and are availables in [this dataset](https://www.kaggle.com/datasets/fabiencrom/msci-correlations).\n\nOur goal will be to find out which targets are dependent on which inputs and by how much. Notably, we observe that 560 of the targets are constants.\n\nFirst, a warning on what are and are not Pearson correlation coefficients. I think this figure and this explanation taken from [this wikipedia entry](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient) do a good job of showing what to expect from such coefficients:\n\n<img src=\"https://upload.wikimedia.org/wikipedia/commons/d/d4/Correlation_examples2.svg\">\n\n> Several sets of (x, y) points, with the correlation coefficient of x and y for each set. Note that the correlation reflects the strength and direction of a linear relationship (top row), but not the slope of that relationship (middle), nor many aspects of nonlinear relationships (bottom). N.B.: the figure in the center has a slope of 0 but in that case the correlation coefficient is undefined because the variance of Y is zero.\n\nSo, as you can see, it is perfectly possible for two variables to have a strong dependency while having a near-zero correlation (see the circle example on the bottom right). Pearson coefficients only measure *linear* correlation (both positive and negative). In practice, however, I think it is uncommon for two variables to have a strong dependency without any linear correlation.\n\nIn addition, we can only compute the *empirical* correlation (or *sample* correlation) which can be different from the *real* correlation (or *population* correlation). In short, the population correlation is the one we could compute if we had an unlimited number of observations. Here, we only have about 100k observations (or samples). This means the number we compute is only an approximation of the population correlation. We might therefore find that some input and target that are in reality independent have a non-zero correlation in our computations. An estimate of the standard error for this dataset is 0.003 (see below). Therefore I expect correlation coefficients whose absolute value is above 0.01 to have a good chance to indicate real correlations.","metadata":{}},{"cell_type":"code","source":"import os\nimport copy\nimport gc\nimport math\nimport itertools\nimport pickle\nimport glob\nimport joblib\nimport json\nimport random\nimport re\nimport operator\n\nfrom collections import defaultdict\nfrom operator import itemgetter, attrgetter\n\nfrom tqdm.notebook import tqdm\n\nimport torch\nimport torch.nn as nn\n\nimport numpy as np\nimport pandas as pd\nimport plotly.express as px\n\nimport scipy\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:17:36.732005Z","iopub.execute_input":"2022-09-11T12:17:36.733376Z","iopub.status.idle":"2022-09-11T12:17:42.054049Z","shell.execute_reply.started":"2022-09-11T12:17:36.733252Z","shell.execute_reply":"2022-09-11T12:17:42.052781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading the Multiome inputs/targets correlations\nCorrelations are stored in a sparse csr matrix, with targets as rows and inputs as columns.\nThus, `correlations[i,j]` is the sample pearson correlation between multiome target number `i` and multiome input number `j`.\n\nWe remind that in this datset, correlations with absolute value smaller than 0.01 have been thresholded to zero. This is mainly for memory optimization (otherwise the full correlation matrix would use around 10GB of memory and would be difficult to manipulate without running into Out-of-Memory errors on kaggle notebooks).\n\nNote that the standard error on sample pearson correlation can be estimated by this formula (see [this wikipedia entry](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Inference)):\n\n$\\sigma _{r}={\\frac {1-r^{2}}{\\sqrt {n-2}}}$\n\nIn the Multiome case, n = 105942; which lead to a standard error of about 0.003 for small correlations. The threshold of 0.01 is therefore about 3 standard errors away from the null correlation. This should ensure that a large numbers of the correlations above 0.01 in absolute value have some significance. But keep in mind that in the Multiome case, we are computing 23418 targets X 228942 inputs = 5.3 billions correlations. It is therefore certain that sampling noise will still give us lot of false positives (ie. relatively large correlations between a feature and a target that are actually independent).","metadata":{}},{"cell_type":"code","source":"correlations = scipy.sparse.load_npz(\"../input/msci-correlations/correlations_multiome.sparse.npz\")","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:17:42.056581Z","iopub.execute_input":"2022-09-11T12:17:42.057046Z","iopub.status.idle":"2022-09-11T12:18:00.523895Z","shell.execute_reply.started":"2022-09-11T12:17:42.056982Z","shell.execute_reply":"2022-09-11T12:18:00.522652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np_inputs_names = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_inputs_idxcol.npz\", allow_pickle=True)[\"columns\"]\nnp_targets_names = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_idxcol.npz\", allow_pickle=True)[\"columns\"]","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:00.525435Z","iopub.execute_input":"2022-09-11T12:18:00.525814Z","iopub.status.idle":"2022-09-11T12:18:00.787243Z","shell.execute_reply.started":"2022-09-11T12:18:00.525777Z","shell.execute_reply":"2022-09-11T12:18:00.785965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Analysis","metadata":{}},{"cell_type":"markdown","source":"### Range of correlations\nFirst let's get an idea of the range of correlations we have in our correlation matrix:","metadata":{}},{"cell_type":"code","source":"print(\"range:\", np.min(correlations), \"-\", np.max(correlations))\nprint(\"absolute values range:\", np.min(np.abs(correlations)), \"-\", np.max(np.abs(correlations)))","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:00.790678Z","iopub.execute_input":"2022-09-11T12:18:00.791224Z","iopub.status.idle":"2022-09-11T12:18:03.728537Z","shell.execute_reply.started":"2022-09-11T12:18:00.791171Z","shell.execute_reply":"2022-09-11T12:18:03.727014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can already observe that the range of correlations is smaller than in the CITEseq case. This seems to indicate that Multiome targets usually do not have strong dependencies with a specific input.\n\n### Maximum correlation per target\nWe can look into that a bit more precisely. For each target, we check the input with which it has the maximum correlation. We find that, on average, targets have a correlation of 0.03-0.04 with their most correlated input, which is quite low (and again, lower than with the CITEseq case). We note that some targets even have a zero correlation with all the inputs!","metadata":{}},{"cell_type":"code","source":"max_abs_correl_per_target = np.array(np.max(np.abs(correlations), axis=1).todense())[:,0]\nprint(\"min:\", np.min(max_abs_correl_per_target), \n      \"mean:\", np.mean(max_abs_correl_per_target), \n      \"median:\", np.median(max_abs_correl_per_target), \n      \"max:\", np.max(max_abs_correl_per_target))","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:03.730103Z","iopub.execute_input":"2022-09-11T12:18:03.730481Z","iopub.status.idle":"2022-09-11T12:18:04.606616Z","shell.execute_reply.started":"2022-09-11T12:18:03.730450Z","shell.execute_reply":"2022-09-11T12:18:04.605311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Inputs and targets with only zero correlations\nTherefore let us check the targets that have zero correlation with all inputs and the inputs that have zero correlation with all targets.","metadata":{}},{"cell_type":"code","source":"# Using `correlations==0' on a sparse matrix would blow up memory. Hence the  lightly convoluted way.\ninputs_with_all_zeros_correlations = np.sum(correlations!=0, axis=0)==0\nnp.sum(inputs_with_all_zeros_correlations)","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:04.608353Z","iopub.execute_input":"2022-09-11T12:18:04.608711Z","iopub.status.idle":"2022-09-11T12:18:06.817701Z","shell.execute_reply.started":"2022-09-11T12:18:04.608678Z","shell.execute_reply":"2022-09-11T12:18:06.816497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"targets_with_all_zeros_correlations = np.sum(correlations!=0, axis=1)==0\nnp.sum(targets_with_all_zeros_correlations)","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:06.819331Z","iopub.execute_input":"2022-09-11T12:18:06.819791Z","iopub.status.idle":"2022-09-11T12:18:08.170908Z","shell.execute_reply.started":"2022-09-11T12:18:06.819747Z","shell.execute_reply":"2022-09-11T12:18:08.169672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Therefore, we have 10 inputs that have no correlations with any target. And 560 targets that have no correlations with any input. Now this could be due to the thresholding we have applied to the Multiome correlations. But actually, the result is the same with the non-thresholded correlations I obtained on Saturn Cloud. Therefore the correlations are really all exactly zero.\n\nThe only realistic way this can happen with real data is that the corresponding entries are all constant. Actually, after checking, they are indeed all equal to zero. Let us prove this.","metadata":{}},{"cell_type":"code","source":"# Let us get the positions of targets with all zeros correlations\nnz_tgts = np.nonzero(np.squeeze(np.asarray(targets_with_all_zeros_correlations)))[0]\n\n# Loading the original targets values\nnp_targets = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_targets_values.sparse.npz\")\n\nprint(\"Number of non-zero elements for the selected columns:\", np.sum(np_targets[:, nz_tgts]!=0))\n\ndel np_targets","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:08.172161Z","iopub.execute_input":"2022-09-11T12:18:08.173182Z","iopub.status.idle":"2022-09-11T12:18:33.348872Z","shell.execute_reply.started":"2022-09-11T12:18:08.173143Z","shell.execute_reply":"2022-09-11T12:18:33.347687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let us get the positions of targets with all zeros correlations\nnz_inpts = np.nonzero(np.squeeze(np.asarray(inputs_with_all_zeros_correlations)))[0]\n\n# Loading the original targets values\nnp_inputs = scipy.sparse.load_npz(\"../input/multimodal-single-cell-as-sparse-matrix/train_multi_inputs_values.sparse.npz\")\n\nprint(\"Number of non-zero elements for the selected columns:\", np.sum(np_inputs[:, nz_inpts]!=0))\ndel np_inputs","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:18:33.350355Z","iopub.execute_input":"2022-09-11T12:18:33.351361Z","iopub.status.idle":"2022-09-11T12:19:31.998072Z","shell.execute_reply.started":"2022-09-11T12:18:33.351322Z","shell.execute_reply":"2022-09-11T12:19:31.996999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"CONCLUSION: There are 10 Multiome inputs constantly equal to zero. And more importantly, there are 560 targets constantly equal to zeros. This mean we should probably force these 560 targets to be equal to zero in our predictions.","metadata":{}},{"cell_type":"markdown","source":"### Distribution of the correlations\nWhich percentage of correlation can be considered significants? It depends on how strict we are for the level of significance. In any case, I expect correlations with absolute values <0.01 to be mostly sampling noise, as explained in the above sections. Let us see how many correlations are left if we suppress the one whose absolute value is below a given significance threshold.","metadata":{}},{"cell_type":"code","source":"for threshold_of_significance in [1e-3, 1e-2, 5e-2, 1e-1]:\n    print(\"threshold\", threshold_of_significance,\n      \"ratio of significant correlations:\", \n          np.sum(np.abs(correlations) > threshold_of_significance)/(correlations.shape[0]*correlations.shape[1])*100, \"%\")\n    ","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:21:19.620432Z","iopub.execute_input":"2022-09-11T12:21:19.620856Z","iopub.status.idle":"2022-09-11T12:21:27.575800Z","shell.execute_reply.started":"2022-09-11T12:21:19.620821Z","shell.execute_reply":"2022-09-11T12:21:27.574588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Of course, since we are using a correlation matrix that has been thresholded at 0.01, it is normal that results for threshold = 0.001 and threshold=0.01 are the same. FYI, here would be the output on the original non-thresholded matrix:\n```\nthreshold 0.0001 ratio of significant correlations: 95.65890430882376 %\nthreshold 0.001 ratio of significant correlations: 77.23980017893044 %\nthreshold 0.005 ratio of significant correlations: 19.388668635600034 %\nthreshold 0.01 ratio of significant correlations: 3.9919700982885518 %\nthreshold 0.05 ratio of significant correlations: 0.033951604159723425 %\nthreshold 0.1 ratio of significant correlations: 0.0014207952205211275 %'```","metadata":{}},{"cell_type":"markdown","source":"The conclusion is that very few of the >5 billion correlations have a correlation above 0.1 (which is in general not even considered a very strong correlation). Given that there are 228942 inputs, we can expect that, on average, each target has about 3 inputs (228942*0.000014) with which it has a correlation above 0.1.\n\nLet us see that in more details.","metadata":{}},{"cell_type":"markdown","source":"### Most globally correlated inputs on average.\n\nLet us check the inputs that have the highest average correlation with the target. For each of them, we display:\n- their average absolute correlation\n- their average correlation\n- What percentage of targets they have a significantly positive or negative correlation\n- Min and max correlation with all targets","metadata":{}},{"cell_type":"code","source":"correlations = correlations.tocsc()\nmost_importants_inputs = np.argsort(np.squeeze(np.asarray(np.mean(np.abs(correlations), axis=0))))[::-1][:15]","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:50:36.552921Z","iopub.execute_input":"2022-09-11T12:50:36.553470Z","iopub.status.idle":"2022-09-11T12:50:38.635458Z","shell.execute_reply.started":"2022-09-11T12:50:36.553428Z","shell.execute_reply":"2022-09-11T12:50:38.634276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print(\"Input             \\tavg_abs_cor\\tavg_cor\\tpc_pos_cor\\tpc_neg_cor\")\n\nname_col = []\nmean_abs_cor_col = []\nmean_cor_col = []\npc_positive_cor_col = []\npc_negative_cor_col = []\nmin_col = []\nmax_col = []\n\nfor num_input in most_importants_inputs:\n    mean_abs_cor = np.mean(np.abs(correlations[:,num_input]))\n    mean_cor = np.mean(correlations[:,num_input])\n    pc_positive_cor = np.mean(correlations[:,num_input]>0.01)*100\n    pc_negative_cor = np.mean(correlations[:,num_input]<-0.01)*100\n    min_val = np.min(correlations[:,num_input])\n    max_val = np.max(correlations[:,num_input])\n    \n    name_col.append(np_inputs_names[num_input])\n    mean_abs_cor_col.append(mean_abs_cor)\n    mean_cor_col.append(mean_cor)\n    pc_positive_cor_col.append(pc_positive_cor)\n    pc_negative_cor_col.append(pc_negative_cor)\n    \n    min_col.append(min_val)\n    max_col.append(max_val)\n    \n#     print(\"%s:    \\t%.3f\\t%.3f\\t%.1f\\t%.1f\"%(np_inputs_names[num_input], mean_abs_cor, mean_cor, pc_positive_cor, pc_negative_cor))\n\npd.DataFrame({\"Input\":name_col, \n              \"avg abs cor\": mean_abs_cor_col,\n              \"avg cor\": mean_cor_col,\n              \"%positive\": pc_positive_cor_col,\n              \"%negative\": pc_negative_cor_col,\n              \"min\":min_col,\n              \"max\":max_col\n             } )","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:50:38.637311Z","iopub.execute_input":"2022-09-11T12:50:38.637757Z","iopub.status.idle":"2022-09-11T12:50:38.711067Z","shell.execute_reply.started":"2022-09-11T12:50:38.637725Z","shell.execute_reply":"2022-09-11T12:50:38.709940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What we can see is that, contrary to the CITEseq dataset, there is no single input with a very high correlation with all targets. We can still see some kind of global trends. For example `chr1:47180897-47181792` seem to be a enhancer for 35% of targets, and an inhibitor for 11% of them. On the contrary, `chr2:28557270-28558187` seem to be an inhibitor for 30% of targets and an enhancer for 11% of them. \n\nThese approximate ratios of  30% / 10% of negative/positive correlations appear surprisingly often. There is probably an explanation to that, but I have not looked into it yet. Possibly there are two subgroups of highly correlated targets representing about 30% and 10% of all targets and that have a very similar response to the same inputs.","metadata":{}},{"cell_type":"markdown","source":"### Most correlated input for each target\nIn contrast with CITEseq, Multiome has way too many targets to display them exhaustively. Let us just look at the first ones to get an idea. Again, for each target, we display the inputs with the highest absolute correlation.","metadata":{}},{"cell_type":"code","source":"for num_tgt in range(20): #cor_matrix.shape[0]):\n    print(np_targets_names[num_tgt])\n    argsorted_cor = np.argsort(np.squeeze(np.asarray(np.abs(correlations[num_tgt]).todense())))[::-1]\n    inv_argsorted = np.argsort(argsorted_cor)\n    for k in argsorted_cor[:10]:\n        is_globally_influential = k in most_importants_inputs\n        print(\"\\t%i\\t%s:   \\t%f\\t%s\"%(inv_argsorted[k], np_inputs_names[k], correlations[num_tgt, k], is_globally_influential))\n        ","metadata":{"execution":{"iopub.status.busy":"2022-09-11T12:59:25.654193Z","iopub.execute_input":"2022-09-11T12:59:25.655195Z","iopub.status.idle":"2022-09-11T12:59:31.882274Z","shell.execute_reply.started":"2022-09-11T12:59:25.655150Z","shell.execute_reply":"2022-09-11T12:59:31.880923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Again, we observe rather low correlations. At best 0.03 for most of them. But we have, for example `ENSG00000115977` that has a not so small correlation of 0.1 with `chr2:28557270-28558187:`.\n\nWe can also observe that, while for a given target the best correlated input has a correlation ranging from 0.01 to  0.1, each target appear to have very similar correlation with a whole set of inputs. This seem to point to the facts that there are gtoups of very correlated inputs and that when a target is correlated with one input of the group, it is correlated similarly with the others.\n","metadata":{}},{"cell_type":"code","source":"most_importants_targets = np.argsort(np.squeeze(np.asarray(np.mean(np.abs(correlations), axis=1))))[::-1][:10]","metadata":{"execution":{"iopub.status.busy":"2022-09-11T13:14:30.194153Z","iopub.execute_input":"2022-09-11T13:14:30.194659Z","iopub.status.idle":"2022-09-11T13:14:32.809450Z","shell.execute_reply.started":"2022-09-11T13:14:30.194621Z","shell.execute_reply":"2022-09-11T13:14:32.807979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for num_tgt in most_importants_targets: #cor_matrix.shape[0]):\n    print(np_targets_names[num_tgt])\n    argsorted_cor = np.argsort(np.squeeze(np.asarray(np.abs(correlations[num_tgt]).todense())))[::-1]\n    inv_argsorted = np.argsort(argsorted_cor)\n    for k in argsorted_cor[:10]:\n        is_globally_influential = k in most_importants_inputs\n        print(\"\\t%i\\t%s:   \\t%f\\t%s\"%(inv_argsorted[k], np_inputs_names[k], correlations[num_tgt, k], is_globally_influential))\n        ","metadata":{"execution":{"iopub.status.busy":"2022-09-11T13:14:32.823775Z","iopub.execute_input":"2022-09-11T13:14:32.824629Z","iopub.status.idle":"2022-09-11T13:14:36.286662Z","shell.execute_reply.started":"2022-09-11T13:14:32.824576Z","shell.execute_reply":"2022-09-11T13:14:36.285365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here, we can see that at least some of the targets have rather important correlations with some inputs.","metadata":{}},{"cell_type":"markdown","source":"# Conclusions\nHere are some of the main conclusion that can be drawn from this analysis:\n- On average, the individual correlations between a single target and a single input are rather small, and quite smaller than for CITEseq.\n- 560 targets are constantly equal to zero in the train dataset. In addition, 10 inputs are constantly equal to zero.\n- Many hints point to the existence of highly correlated subgroups of targets as well as highly correlated subgroups of inputs. This should be investigated furthers.","metadata":{}}]}