{"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":"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.24518Z","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.24557Z","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.02277Z","iopub.status.busy":"2022-09-11T11:24:10.02248Z","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.49007Z","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.58544Z","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.61349Z","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.63074Z","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.38457Z"},"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.95407Z","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.64284Z","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":[]}]}