{"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 CITEseq","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 CITEseq training set. There is a [corresponding notebook for the Multiome data](https://www.kaggle.com/fabiencrom/msci-correlations-eda-multiome).\n\nThe Pearsons correlation coefficients were computed with [this notebook](https://www.kaggle.com/fabiencrom/msci-generating-all-correlations-inputs-targets/edit) 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. Some findings:\n- All targets have important correlations with at least some inputs (contrary at what happens in the Multiome data)\n- Some inputs have an important global influence on almost all targets. They should probably not be discarded when doing dimensionality reduction.\n- Contrary to what was expected, inputs whose name match a target are not always well correlated (although this can also be the case)\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 have less then 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-13T10:57:07.473084Z","iopub.execute_input":"2022-09-13T10:57:07.473597Z","iopub.status.idle":"2022-09-13T10:57:11.236526Z","shell.execute_reply.started":"2022-09-13T10:57:07.473499Z","shell.execute_reply":"2022-09-13T10:57:11.234955Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading the CITEseq inputs/targets correlations\nCorrelations are stored in a matrix with targets as rows and inputs as columns.\nThus, `correlations[i,j]` is the sample pearson correlation between CITEseq target number `i` and CITEseq input number `j`.\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 CITEseq case, n = 70988; 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.","metadata":{}},{"cell_type":"code","source":"correlations = np.load(\"../input/msci-correlations/correlations_citeseq.npz\")[\"correlations\"]","metadata":{"execution":{"iopub.status.busy":"2022-09-13T10:57:11.238859Z","iopub.execute_input":"2022-09-13T10:57:11.239248Z","iopub.status.idle":"2022-09-13T10:57:11.404040Z","shell.execute_reply.started":"2022-09-13T10:57:11.239207Z","shell.execute_reply":"2022-09-13T10:57:11.402907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np_inputs_names = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/train_cite_inputs_idxcol.npz\", allow_pickle=True)[\"columns\"]\nnp_targets_names = np.load(\"../input/multimodal-single-cell-as-sparse-matrix/train_cite_targets_idxcol.npz\", allow_pickle=True)[\"columns\"]","metadata":{"execution":{"iopub.status.busy":"2022-09-13T10:57:11.407200Z","iopub.execute_input":"2022-09-13T10:57:11.408138Z","iopub.status.idle":"2022-09-13T10:57:11.495924Z","shell.execute_reply.started":"2022-09-13T10:57:11.408084Z","shell.execute_reply":"2022-09-13T10:57:11.494823Z"},"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-13T10:57:11.499552Z","iopub.execute_input":"2022-09-13T10:57:11.500555Z","iopub.status.idle":"2022-09-13T10:57:11.696767Z","shell.execute_reply.started":"2022-09-13T10:57:11.500501Z","shell.execute_reply":"2022-09-13T10:57:11.695351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can already observe that the range of correlations are larger than in the Multiome case. This seems to indicate that at least some CITEseq targets have a strong dependence on some inputs.\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.3 with their most correlated input, which is significant (and again, much higher than what is observed in with the Multiome case).","metadata":{}},{"cell_type":"code","source":"max_abs_correl_per_target = np.max(np.abs(correlations), axis=1)\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-13T10:57:11.698429Z","iopub.execute_input":"2022-09-13T10:57:11.698831Z","iopub.status.idle":"2022-09-13T10:57:11.750062Z","shell.execute_reply.started":"2022-09-13T10:57:11.698794Z","shell.execute_reply":"2022-09-13T10:57:11.748989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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, 3e-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-13T10:57:11.751597Z","iopub.execute_input":"2022-09-13T10:57:11.752594Z","iopub.status.idle":"2022-09-13T10:57:11.979699Z","shell.execute_reply.started":"2022-09-13T10:57:11.752555Z","shell.execute_reply":"2022-09-13T10:57:11.978214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The conclusion is that few of the correlations have a correlation above 0.1 (which is in general not even considered a very strong correlation). Given that there are 22050 inputs, we can expect that, on average, each target has about 390 inputs (22050*0.0177) with which it has a correlation above 0.1 and 9 inputs with a correlation above 0.3.\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":"most_importants_inputs = np.argsort(np.squeeze(np.asarray(np.mean(np.abs(correlations), axis=0))))[::-1][:15]","metadata":{"execution":{"iopub.status.busy":"2022-09-13T10:57:11.982039Z","iopub.execute_input":"2022-09-13T10:57:11.982569Z","iopub.status.idle":"2022-09-13T10:57:12.004338Z","shell.execute_reply.started":"2022-09-13T10:57:11.982520Z","shell.execute_reply":"2022-09-13T10:57:12.002673Z"},"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-13T10:57:12.005946Z","iopub.execute_input":"2022-09-13T10:57:12.006364Z","iopub.status.idle":"2022-09-13T10:57:12.048496Z","shell.execute_reply.started":"2022-09-13T10:57:12.006327Z","shell.execute_reply":"2022-09-13T10:57:12.046965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What we can see is that some inputs, like `ENSG00000129824_RPS4Y1` seem to act like a global inhibitor, having a strong negative correlation with almost all targets. At the opposite, `ENSG00000229807_XIST` is a global \"enhancer\", with a strong positive correlation with most targets. `ENSG00000102145_GATA1` has inhibitor effect on 38% of targets but enhancing effect on 55% of them.\n\nIf someone is using dimensionality reduction, like in many notebooks so far, it seems that these inputs with strong global impact should be kept as is.","metadata":{}},{"cell_type":"markdown","source":"### Most correlated input for each target\nLet us check, for each target, the most correlated inputs. \n\nIn addition, in [this discussion](https://www.kaggle.com/competitions/open-problems-multimodal/discussion/349242), it was suggested that some of the inputs have names related to some targets, and that these inputs should be used as predictors for these targets. We are going to check the actual correlations between the targets (proteins) and its associated genes.\n\nIn addition, we indicate by a True/False flag if the input is among the global influencers we have found earlier.","metadata":{}},{"cell_type":"code","source":"# We consider that a gene is associated with a protein if its name contains the protein name. \n# This is a bit noisy, as e.g. it will associate gene ENSG00000105383_CD33 to protein CD3\n\nassociated_genes = defaultdict(list)\nfor protein_name in np_targets_names:\n    for num_gene, gene_name in enumerate(np_inputs_names):\n        if protein_name in gene_name:\n            associated_genes[protein_name].append((num_gene, gene_name))\n            \n","metadata":{"execution":{"iopub.status.busy":"2022-09-13T10:57:12.050999Z","iopub.execute_input":"2022-09-13T10:57:12.051400Z","iopub.status.idle":"2022-09-13T10:57:12.582630Z","shell.execute_reply.started":"2022-09-13T10:57:12.051363Z","shell.execute_reply":"2022-09-13T10:57:12.581355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Th following code will show the top 5 most influent gene for each target protein.","metadata":{}},{"cell_type":"code","source":"for num_tgt in range(correlations.shape[0]):\n    print(np_targets_names[num_tgt])\n    argsorted_cor = np.argsort(np.abs(correlations[num_tgt]))[::-1]\n    inv_argsorted = np.argsort(argsorted_cor)\n    for k in argsorted_cor[:5]:\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        \n    if np_targets_names[num_tgt] in associated_genes:\n        print (\"\\t******* associated genes: *******\")\n        for num_gene, gene_name in associated_genes[np_targets_names[num_tgt]]:\n            print(\"\\t%i\\t%s:   \\t%f\"%(inv_argsorted[num_gene], gene_name, correlations[num_tgt, num_gene]))\n    print()","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-09-13T10:19:58.494757Z","iopub.execute_input":"2022-09-13T10:19:58.495289Z","iopub.status.idle":"2022-09-13T10:19:59.179821Z","shell.execute_reply.started":"2022-09-13T10:19:58.495232Z","shell.execute_reply":"2022-09-13T10:19:59.178596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see again that all targets have relatively strong correlations with at least some inputs. More surprisingly, we found that the associated genes do not always have a strong correlation. `ENSG00000114013_CD86` is well correlated with `CD86`. But `ENSG00000120217_CD274` has a relatively weak correlation with `CD274`.\n\n\n","metadata":{}},{"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- All targets have important correlations with at least some inputs (contrary at what happens in the Multiome data)\n- Some inputs have an important global influence on almost all targets. They should probably not be discarded when doing dimensionality reduction.\n- Contrary to what was expected, inputs whose name match a target are not always well correlated (although this can also be the case)\n- Models could use the correlation informations for regularization (e.g. enforcing some sparsity in the matrix of a linear regression by setting to zero coefficients associated with input/target pairs with low correlations)\n","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}