{"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: Generating the correlations between all inputs and all targets\nThis notebook aims at generating the pearsons correlations coefficients between all inputs and all targets for both Multiome and CITEseq datasets.\n\nWARNING: For the Multiome dataset, this notebook is not supposed to run on kaggle machines as they do not have enough memory. Saturn Cloud is kindly offering us access to machines with 128GB RAM, and that is what I have been using to run this code. See [this discussion](https://www.kaggle.com/competitions/open-problems-multimodal/discussion/346999). The correlations for the CITEseq data can be generated on kaggle.\n\nWARNING 2: this version is a work in progress and has not been tested on multiome yet. Use the previous version if you want to compute multiome correlations.\n\nThe correlations matrices generated from this notebook are available in [this kaggle dataset](https://www.kaggle.com/datasets/fabiencrom/msci-correlations), so you do not have to run this notebook yourself. (Note that the multiome correlations are stored in a sparse format with weak correlations set to zero for space/memory efficiency).\n\nYou can also check these two notebooks for an analysis of the correlations:\n- https://www.kaggle.com/fabiencrom/msci-correlations-eda-multiome\n- https://www.kaggle.com/fabiencrom/msci-correlations-eda-citeseq","metadata":{"execution":{"iopub.execute_input":"2022-09-11T01:20:12.824274Z","iopub.status.busy":"2022-09-11T01:20:12.823991Z","iopub.status.idle":"2022-09-11T01:20:15.356694Z","shell.execute_reply":"2022-09-11T01:20:15.356188Z","shell.execute_reply.started":"2022-09-11T01:20:12.824250Z"},"tags":[]}},{"cell_type":"code","source":"\nimport 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\nimport scipy.sparse\n\n","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:41:45.245180Z","iopub.execute_input":"2022-09-20T11:41:45.245676Z","iopub.status.idle":"2022-09-20T11:41:48.487279Z","shell.execute_reply.started":"2022-09-20T11:41:45.245570Z","shell.execute_reply":"2022-09-20T11:41:48.485864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Path and Choosing CITEseq/Multiome\n\nFor efficiency, I use a version of the competition data encoded in sparse format that can be found in [this dataset](https://www.kaggle.com/datasets/fabiencrom/multimodal-single-cell-as-sparse-matrix). If you run this on Saturn Cloud, you will have to download the dataset and set the SPARSE_DATA_PATH variable accordingly.","metadata":{"execution":{"iopub.execute_input":"2022-09-11T11:24:10.022770Z","iopub.status.busy":"2022-09-11T11:24:10.022480Z","iopub.status.idle":"2022-09-11T11:24:10.025476Z","shell.execute_reply":"2022-09-11T11:24:10.025036Z","shell.execute_reply.started":"2022-09-11T11:24:10.022745Z"}}},{"cell_type":"code","source":"SPARSE_DATA_PATH = \"../input/multimodal-single-cell-as-sparse-matrix/\"\n\n# choose which dataset to work on:\n\n# dataset_infix = \"multi\" #multiome\ndataset_infix = \"cite\" #CITEseq \n\ninputs_fn = SPARSE_DATA_PATH + \"train_%s_inputs_values.sparse.npz\"%dataset_infix\ntargets_fn = SPARSE_DATA_PATH + \"train_%s_targets_values.sparse.npz\"%dataset_infix\n\ntargets_idxcol_names_fn = SPARSE_DATA_PATH + \"train_%s_targets_idxcol.npz\"%dataset_infix\ninputs_idxcol_names_fn = SPARSE_DATA_PATH + \"train_%s_inputs_idxcol.npz\"%dataset_infix\n\nprint(inputs_fn, targets_fn)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:41:48.490070Z","iopub.execute_input":"2022-09-20T11:41:48.490564Z","iopub.status.idle":"2022-09-20T11:41:48.498342Z","shell.execute_reply.started":"2022-09-20T11:41:48.490517Z","shell.execute_reply":"2022-09-20T11:41:48.497071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading Data","metadata":{}},{"cell_type":"code","source":"np_inputs_names = np.load(inputs_idxcol_names_fn, allow_pickle=True)[\"columns\"]\nnp_targets_names = np.load(targets_idxcol_names_fn, allow_pickle=True)[\"columns\"]","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:41:48.500098Z","iopub.execute_input":"2022-09-20T11:41:48.500828Z","iopub.status.idle":"2022-09-20T11:41:48.582649Z","shell.execute_reply.started":"2022-09-20T11:41:48.500781Z","shell.execute_reply":"2022-09-20T11:41:48.581533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nnp_tgts = scipy.sparse.load_npz(targets_fn)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:41:48.584869Z","iopub.execute_input":"2022-09-20T11:41:48.585477Z","iopub.status.idle":"2022-09-20T11:41:49.597401Z","shell.execute_reply.started":"2022-09-20T11:41:48.585440Z","shell.execute_reply":"2022-09-20T11:41:49.596284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nnp_inpts = scipy.sparse.load_npz(inputs_fn)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:41:49.598659Z","iopub.execute_input":"2022-09-20T11:41:49.598994Z","iopub.status.idle":"2022-09-20T11:42:12.611493Z","shell.execute_reply.started":"2022-09-20T11:41:49.598964Z","shell.execute_reply":"2022-09-20T11:42:12.610249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Computing correlations\n\nCorrelations are computed using the following formula for covariance: $Cov(X,Y)=E(XY)-E(X)E(Y)$. This allows us to compute the correlations from the sparse matrices. (Using the \"official\" formula $Cov(X,Y)=E((X-E(X))(Y-E(Y))$ would lead to the sparse matrices being blown up into a dense matrix).","metadata":{}},{"cell_type":"code","source":"%%time\ndef compute_mean_and_var_sparse(np_inpts):\n    np_inpts_means = np.array(np_inpts.sum(axis=0)/np_inpts.shape[0])\n    temp = np_inpts.copy()\n    temp.data = temp.data**2\n    mean_squares = np.array(temp.sum(axis=0))/np_inpts.shape[0]\n    np_inpts_vars = mean_squares - np_inpts_means**2\n    return np_inpts_means, np_inpts_vars\n    \ndef compute_mean_and_var_dense(np_inpts, only_nonzeros=False):\n    np_inpts_means = np_inpts.mean(axis=0, keepdims=True)\n    \n    if only_nonzeros:\n        nnz_inpts = np.sum(np_inpts!=0, axis=0)\n        np_inpts_means_nz = np_inpts_means*(np_inpts.shape[0]/ (nnz_inpts+1e-30))\n        \n        np_inpts_vars = np.sum(np_inpts**2, axis=0)/(nnz_inpts+1e-30) - np_inpts_means_nz**2\n#         np_inpts_vars_nz = (np_inpts_vars+np_inpts_means**2)*(np_inpts.shape[0]/ (nnz_inpts+1e-30)) - np_inpts_means_nz**2\n        np_inpts_means = np_inpts_means_nz\n#     np_inpts_vars = , np_inpts_vars_nz\n        \n    else:\n        np_inpts_vars = np_inpts.var(axis=0, keepdims=True)\n        \n    return np_inpts_means, np_inpts_vars\n\ndef compute_cor_dense(np_inpts, np_tgts, only_nonzero_inpts=False):\n    if only_nonzero_inpts:\n        np_inpts_means, np_inpts_vars = compute_mean_and_var(np_inpts, only_nonzeros=True)\n        denom = np.sum(np_inpts!=0, axis=0, keepdims=True) + 1e-30\n        np_tgts_means_nz_inpt = np_tgts.T @ (np_inpts > 0).astype(np.float32)/denom\n        np_tgts_vars_nz_inpt = (np_tgts**2).T @ (np_inpts > 0).astype(np.float32)/denom- np_tgts_means_nz_inpt**2\n        cov = (np_tgts.T @ np_inpts)/denom- (np_tgts_means_nz_inpt*np_inpts_means)\n        cor = cov/(np.sqrt(np_tgts_vars_nz_inpt+1e-30)*np.sqrt(np_inpts_vars+1e-30))\n        cor[:, np_inpts_vars[0]==0]=0\n    else:\n        denom = np_inpts.shape[0] \n        np_inpts_means, np_inpts_vars = compute_mean_and_var(np_inpts, only_nonzeros=False)\n        np_tgts_means, np_tgts_vars = compute_mean_and_var(np_tgts, only_nonzeros=False)\n        cov = (np_tgts.T @ np_inpts)/denom- (np_tgts_means.T*np_inpts_means)\n        cor = cov/(np.sqrt(np_tgts_vars.T+1e-30)*np.sqrt(np_inpts_vars+1e-30))\n        cor[:, np_inpts_vars[0]==0]=0\n    return cor\n\ndef compute_cor_sparse_slice(start, end):\n    cov_first20 = (np_tgts[:,start:end].todense().T @ np_inpts)/np_inpts.shape[0] - (np_tgts_means.T[start:end]*np_inpts_means)\n    cor_first20 = cov_first20/(np.sqrt(np_tgts_vars.T[start:end]+1e-30)*np.sqrt(np_inpts_vars+1e-30))\n    return cor_first20\n\ndef compute_cor_sparse(mb_size=1000):\n    nb_iterations = np_tgts.shape[1]//mb_size + (np_tgts.shape[1]%mb_size!=0)\n    \n    cor_list = []\n    for i in tqdm(range(nb_iterations)):\n        cor = compute_cor_sparse_slice(i*mb_size, (i+1)*mb_size).astype(np.float32)\n        cor_list.append(np.array(cor))\n    return np.vstack(cor_list)\n","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:42:12.613149Z","iopub.execute_input":"2022-09-20T11:42:12.613522Z","iopub.status.idle":"2022-09-20T11:42:12.628808Z","shell.execute_reply.started":"2022-09-20T11:42:12.613490Z","shell.execute_reply":"2022-09-20T11:42:12.627509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if dataset_infix == \"cite\":\n    np_inpts = np.array(np_inpts.todense())\n    np_tgts = np.array(np_tgts.todense())  \n    compute_mean_and_var = compute_mean_and_var_dense\n    compute_cor = compute_cor_dense\nelse:\n    np_inpts = np_inpts.tocsc()\n    np_tgts = np_tgts.tocsc()\n    compute_mean_and_var = compute_mean_and_var_sparse\n    compute_cor = compute_cor_sparse\n    ","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:42:12.630356Z","iopub.execute_input":"2022-09-20T11:42:12.630740Z","iopub.status.idle":"2022-09-20T11:42:29.814962Z","shell.execute_reply.started":"2022-09-20T11:42:12.630709Z","shell.execute_reply":"2022-09-20T11:42:29.813745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Computing correlation of targets with on/off inputs feature","metadata":{}},{"cell_type":"code","source":"\nnon_zero_binary_feature = (np_inpts!=0).astype(np.int8)\ndel np_inpts\ngc.collect()\n# posinpts = posinpts.astype(np.float32)\ncor_binary = compute_cor(non_zero_binary_feature, np_tgts, only_nonzero_inpts=False)\n# nb_nonzeros_per_inpts = np.sum(non_zero_binary_feature, axis=0)\n# del non_zero_binary_feature\ngc.collect()\n","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:42:29.816322Z","iopub.execute_input":"2022-09-20T11:42:29.816661Z","iopub.status.idle":"2022-09-20T11:42:47.532157Z","shell.execute_reply.started":"2022-09-20T11:42:29.816631Z","shell.execute_reply":"2022-09-20T11:42:47.531391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np_inpts = scipy.sparse.load_npz(inputs_fn)\nnp_inpts = np.asarray(np_inpts.todense())","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:42:47.533395Z","iopub.execute_input":"2022-09-20T11:42:47.533869Z","iopub.status.idle":"2022-09-20T11:43:18.734956Z","shell.execute_reply.started":"2022-09-20T11:42:47.533838Z","shell.execute_reply":"2022-09-20T11:43:18.733741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Computing correlation of targets with non-zero inputs ","metadata":{}},{"cell_type":"code","source":"%%time\nnp_inpts_means_nz, np_inpts_vars_nz = compute_mean_and_var(np_inpts, only_nonzeros=True)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:43:18.738082Z","iopub.execute_input":"2022-09-20T11:43:18.738465Z","iopub.status.idle":"2022-09-20T11:43:24.050134Z","shell.execute_reply.started":"2022-09-20T11:43:18.738431Z","shell.execute_reply":"2022-09-20T11:43:24.049048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nnp_tgts_means, np_tgts_vars = compute_mean_and_var(np_tgts, only_nonzeros=False)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:43:24.051414Z","iopub.execute_input":"2022-09-20T11:43:24.051979Z","iopub.status.idle":"2022-09-20T11:43:24.096912Z","shell.execute_reply.started":"2022-09-20T11:43:24.051938Z","shell.execute_reply":"2022-09-20T11:43:24.095442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncor_matrix = compute_cor(np_inpts, np_tgts, only_nonzero_inpts=True)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:43:24.098233Z","iopub.execute_input":"2022-09-20T11:43:24.098728Z","iopub.status.idle":"2022-09-20T11:43:52.386338Z","shell.execute_reply.started":"2022-09-20T11:43:24.098694Z","shell.execute_reply":"2022-09-20T11:43:52.384570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nb_nonzeros_per_inpts = np.sum(non_zero_binary_feature, axis=0)\n\nstderr_inpts = (1-cor_matrix**2)/(np.sqrt(np.maximum(nb_nonzeros_per_inpts[None]-2, 0))+1e-30)\nstderr_inpts[stderr_inpts>1e20] = 0","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:43:54.041463Z","iopub.execute_input":"2022-09-20T11:43:54.041903Z","iopub.status.idle":"2022-09-20T11:43:55.391478Z","shell.execute_reply.started":"2022-09-20T11:43:54.041866Z","shell.execute_reply":"2022-09-20T11:43:55.390061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def check_corr(np_tgts, np_inpts, i, j, only_non_zeros=False):\n    \n    if only_non_zeros:\n        X = np_tgts[np_inpts[:,i]!=0,j:j+1].astype(np.float32)\n        Y = np_inpts[np_inpts[:,i]!=0,i:i+1].astype(np.float32)\n    else:\n        X = np_tgts[:,j:j+1].astype(np.float32)\n        Y = np_inpts[:,i:i+1].astype(np.float32)\n        \n#     np.mean(X*Y)-np.mean(X)*np.mean(Y) #np_tgts_means.astype(np.float64)[0,0]*np_inpts_means.astype(np.float64)[0,0]\n    \n    return np.corrcoef(X,Y, rowvar=False)[0,1]\n    \n    cor_matrix[2,1]\n    \n    np_tgts_means[0,0], np.mean(X)\n    \n    np_inpts_means[0,0], np.mean(Y)\n    \n    np.mean(X*Y)-np_tgts_means.astype(np.float64)[0,0]*np_inpts_means.astype(np.float64)[0,0]\n    np.mean((np_tgts[np_inpts[:,0]!=0,:1] - np_tgts_means.astype(np.float64)[0,0]) * (np_inpts[np_inpts[:,0]!=0,:1]-np_inpts_means.astype(np.float64)[0,0]))\n    np_tgts[:,:1].T @ np_inpts[:,:1]\n    (np_tgts_means[:,:1].T*np_inpts_means[:,:1])\n    np_inpts[:,:1].var(axis=0)\n    \n    np.corrcoef(np_tgts[:,0], np_inpts[:,0])\n    \n    np.corrcoef(np_tgts[np_inpts[:,0]!=0,0], np_inpts[np_inpts[:,0]!=0,0])\n    \n    np_inpts_means_nz[:,0], np_inpts_vars_nz[:,0]\n    ","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:50:21.863316Z","iopub.execute_input":"2022-09-20T11:50:21.863807Z","iopub.status.idle":"2022-09-20T11:50:21.880435Z","shell.execute_reply.started":"2022-09-20T11:50:21.863766Z","shell.execute_reply":"2022-09-20T11:50:21.879108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.abs(check_corr(np_tgts, np_inpts, 2, 5, only_non_zeros=True) - cor_matrix[5,2])<1e-4","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:51:33.533151Z","iopub.execute_input":"2022-09-20T11:51:33.533737Z","iopub.status.idle":"2022-09-20T11:51:33.549288Z","shell.execute_reply.started":"2022-09-20T11:51:33.533695Z","shell.execute_reply":"2022-09-20T11:51:33.547555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.any(np.isnan(cor_matrix)), np.any(np.isnan(cor_binary)), np.any(np.isnan(stderr_inpts))","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:45:08.953659Z","iopub.execute_input":"2022-09-20T11:45:08.954070Z","iopub.status.idle":"2022-09-20T11:45:08.970593Z","shell.execute_reply.started":"2022-09-20T11:45:08.954038Z","shell.execute_reply":"2022-09-20T11:45:08.969269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Saving the correlations\nFor CITEseq, we store all correlations as well as means and variances of inputs and targets in a npz archive.\n\nFor Multiome, We store the means and variances of inputs and targets in a npz archive. But since there are more than 5 billions correlations, we do the following to keep the size of data manageable in a Kaggle notebook:\n- Threshold all correlations with absolute value smaller than 0.01 to 0. (see the EDA notebooks for a discussion on why this is a reasonable threshold)\n- Store the correlation matrix as a CSR sparse matrix","metadata":{}},{"cell_type":"code","source":"if dataset_infix == \"cite\":\n    np.savez(\"correlations_zero_nonzero_citeseq.npz\", \n         correlations_nonzero=cor_matrix.astype(np.float16),\n         correlations_nonzero_std_err = stderr_inpts.astype(np.float16),\n         correlations_zero_binary = cor_binary.astype(np.float16),\n         non_zero_binary_feature = non_zero_binary_feature.astype(np.int8),\n         inputs_means=np_inpts_means_nz, inputs_vars=np_inpts_vars_nz,\n         targets_means=np_tgts_means, targets_vars=np_tgts_vars\n        )\nelse:\n    cor_matrix = cor_matrix.astype(np.float16)\n    cor_matrix[np.abs(cor_matrix)<0.01] = 0\n    cor_matrix_sparse = scipy.sparse.csr_matrix(cor_matrix_f16)\n    scipy.sparse.save_npz(\"correlations_multiome_th001.sparse.npz\", cor_matrix_sparse)\n    np.savez(\"means_vars_multiome.npz\", \n         inputs_means=np_inpts_means, inputs_vars=np_inpts_vars,\n         targets_means=np_tgts_means, targets_vars=np_tgts_vars\n        )\n    ","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-20T11:45:27.364758Z","iopub.execute_input":"2022-09-20T11:45:27.365243Z","iopub.status.idle":"2022-09-20T11:45:31.642840Z","shell.execute_reply.started":"2022-09-20T11:45:27.365179Z","shell.execute_reply":"2022-09-20T11:45:31.641479Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!cp ../input/msci-correlations/* .","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:46:04.582119Z","iopub.execute_input":"2022-09-20T11:46:04.582922Z","iopub.status.idle":"2022-09-20T11:46:15.429723Z","shell.execute_reply.started":"2022-09-20T11:46:04.582874Z","shell.execute_reply":"2022-09-20T11:46:15.428264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!ls -lh","metadata":{"execution":{"iopub.status.busy":"2022-09-20T11:46:23.338247Z","iopub.execute_input":"2022-09-20T11:46:23.338702Z","iopub.status.idle":"2022-09-20T11:46:24.578663Z","shell.execute_reply.started":"2022-09-20T11:46:23.338661Z","shell.execute_reply":"2022-09-20T11:46:24.577339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}