{"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":"# EDA Multiome targets, which is zero and which is non-zero?\n\n## Purpose of this notebook \n[Fabien Crom calculated all correlations](https://www.kaggle.com/competitions/open-problems-multimodal/discussion/351725) , but I found strange features in scatter graphs as I posted. Data can be split into 2 or 3 groups, and there is gaps between them. Which implied 0 is NA, not measured 0.\nI'm afraid zero data are not deterministic.(Even if donor, day and cell type are the same, some have statistical distribution around some mean, but some are just 0). So I took a closer look at raw data.<br>\nOn the other hand, NORIFUMI IRIE and ALEXANDER CHERVOV use heatmaps to study similar things about multiome inputs and implied that chromatin accessibility depends on days or cell types. ( [Heatmap of Multiome Data](https://www.kaggle.com/code/norifumiirie/heatmap-of-multiome-data) ,  [Heatmap of ATAC-seq (Multiome) Data](https://www.kaggle.com/code/alexandervc/heatmap-of-atac-seq-multiome-data) ) So I applied similar heatmaps to multiome ouptputs(gene expressions.) \n\n## Result\nEven if donor, day and cell type are the same, if gene expression is zero or non-zero are (seemingly) random <br>\nThe question is, the correct answers are chosen from non-zero gene expressions, or contain 0 as training data contains?<br>\n<br>\nAbout heatmaps, although I only plotted donor 32606, I could recognize day 7 is thinner than other days. But I could not recognize a difference between cell types. And I could recognize irregularity within the same donor, day, and cell-types.\n\n### P.S.\nI had created metadata & sparse data set [convert-data-for-speed-up](https://www.kaggle.com/code/konomuabe/convert-data-for-speed-up) or [Open Problems - Multimodal Single-Cell Metadata](https://www.kaggle.com/datasets/konomuabe/open-problems-multimodal-singlecell-metadata). This is useful for these analysis.\n","metadata":{}},{"cell_type":"code","source":"import os, gc, pickle\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom colorama import Fore, Back, Style\nfrom matplotlib.ticker import MaxNLocator\nimport time\nimport tqdm\n\nfrom sklearn.base import BaseEstimator, TransformerMixin\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import StandardScaler, scale\nfrom sklearn.decomposition import PCA\nfrom sklearn.dummy import DummyRegressor\nfrom sklearn.pipeline import make_pipeline, Pipeline\nfrom sklearn.linear_model import Ridge, LinearRegression\nfrom sklearn.metrics import mean_squared_error\n\nfrom scipy.sparse import *\nfrom scipy.sparse import linalg\nfrom sklearn.decomposition import TruncatedSVD\nfrom scipy import interpolate\nimport random as rd\n\nimport psutil\n\ntry:\n    os.environ['SATURN_IMAGE']\n    kernel = \"saturn_cloud\"\nexcept KeyError:\n    kernel = \"kaggle\"\n\nDATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nif kernel==\"saturn_cloud\":\n    SPARSE_DIR = \"./\"\n    META_DIR = \"./\"\nelse: # kaggle\n#    SPARSE_DIR = \"../input/opmsci-sparse-data/\"\n    META_DIR = \"../input/open-problems-multimodal-singlecell-metadata/\"\n    SPARSE_DIR = \"../input/open-problems-multimodal-singlecell-metadata/\"\n\nMETA_CITE_TRAIN = os.path.join(META_DIR,\"cite_train_meta.csv\")\nMETA_CITE_TEST = os.path.join(META_DIR,\"cite_test_meta.csv\")\nMETA_MULTI_TRAIN = os.path.join(META_DIR,\"multiome_train_meta.csv\")\nMETA_MULTI_TEST = os.path.join(META_DIR,\"multiome_test_meta.csv\")\n\nSPS_CITE_TRAIN_INPUTS = os.path.join(SPARSE_DIR,\"train_cite_inputs.npz\")\nSPS_CITE_TRAIN_TARGETS = os.path.join(SPARSE_DIR,\"train_cite_targets.npz\")\nSPS_CITE_TEST_INPUTS = os.path.join(SPARSE_DIR,\"test_cite_inputs.npz\")\n\nSPS_MULTIOME_TRAIN_INPUTS = os.path.join(SPARSE_DIR,\"train_multi_inputs.npz\")\nSPS_MULTIOME_TRAIN_TARGETS = os.path.join(SPARSE_DIR,\"train_multi_targets.npz\")\nSPS_MULTIOME_TEST_INPUTS = os.path.join(SPARSE_DIR,\"test_multi_inputs.npz\")\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\nif not os.path.exists('/opt/conda/lib/python3.7/site-packages/tables'):  ## AMBROSM's trick\n    !pip install --quiet tables","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-09-17T10:59:09.108839Z","iopub.execute_input":"2022-09-17T10:59:09.109324Z","iopub.status.idle":"2022-09-17T10:59:25.043016Z","shell.execute_reply.started":"2022-09-17T10:59:09.109225Z","shell.execute_reply":"2022-09-17T10:59:25.041647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## load variables","metadata":{}},{"cell_type":"code","source":"use_rows_percent = 100  ## how much percent is used for calculation\nrd.seed(1)\n\nmulti_x = load_npz(SPS_MULTIOME_TRAIN_INPUTS)\nif use_rows_percent< 100:\n    rows_to_use = (multi_x.shape[0]*use_rows_percent)//100\n    rows_to_use = rd.sample(range(multi_x.shape[0]),rows_to_use)\nelse:\n    rows_to_use = list(range(multi_x.shape[0]))\n\nprint(\"Used %.1f Gbyte RAM\"%(psutil.virtual_memory()[3]/(1024**3)))\nmulti_x = multi_x[rows_to_use,:]\ngc.collect()\nprint(\"Used %.1f Gbyte RAM\"%(psutil.virtual_memory()[3]/(1024**3)))\n\nmulti_y = load_npz(SPS_MULTIOME_TRAIN_TARGETS)\n\nprint(\"Used %.1f Gbyte RAM\"%(psutil.virtual_memory()[3]/(1024**3)))\nmulti_y = multi_y[rows_to_use,:]\ngc.collect()\nprint(\"Used %.1f Gbyte RAM\"%(psutil.virtual_memory()[3]/(1024**3)))\n\nmulti_ym = pd.read_csv(META_MULTI_TRAIN)\nmulti_ym = multi_ym.iloc[rows_to_use,:]  ## Bug fixed at Version 2 (but conclusion was not affected)\n\nexpression_names = pd.read_hdf(FP_MULTIOME_TRAIN_TARGETS,start=0,stop=1).columns.values\naccess_names = pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS,start=0,stop=1).columns.values\n\n#multi_y = multi_y.toarray()\nmulti_y = multi_y.tocsc()\nmulti_x = multi_x.tocsc()  \nprint(\"Used %.1f Gbyte RAM\"%(psutil.virtual_memory()[3]/(1024**3)))","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-09-17T10:59:25.045593Z","iopub.execute_input":"2022-09-17T10:59:25.046004Z","iopub.status.idle":"2022-09-17T11:01:48.021773Z","shell.execute_reply.started":"2022-09-17T10:59:25.045964Z","shell.execute_reply":"2022-09-17T11:01:48.020478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## First, plot scatter of some input/target pairs to see how they are distributed.\nLet's choose random 4x4 pairs.","metadata":{}},{"cell_type":"code","source":"access_n ,express_n = 4,4\naccess_samples = rd.sample(range(multi_x.shape[1]),access_n)\nexpress_samples = rd.sample(range(multi_y.shape[1]),express_n)","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:01:48.023226Z","iopub.execute_input":"2022-09-17T11:01:48.023580Z","iopub.status.idle":"2022-09-17T11:01:48.030038Z","shell.execute_reply.started":"2022-09-17T11:01:48.023546Z","shell.execute_reply":"2022-09-17T11:01:48.028803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig,ax =plt.subplots(express_n,access_n,figsize=(15,15))\nplt.subplots_adjust(wspace=0.4, hspace=0.6)\nfor ei in range(express_n):\n    yidx = express_samples[ei]\n    for ai in range(access_n):\n        xidx = access_samples[ai]\n        ax[ei,ai].scatter(multi_x[:,xidx].toarray(), multi_y[:,yidx].toarray())\n        _ = ax[ei,ai].set_xlabel('chromatine accesibility')\n        _ = ax[ei,ai].set_ylabel('gene expression')\n        _ = ax[ei,ai].set_title('%s \\n %s'%(expression_names[yidx],access_names[xidx]))\n        corr = np.corrcoef(multi_x[:,xidx].toarray().T, multi_y[:,yidx].toarray().T)\n        print(\"%s %.2f %% non-zero , %s %.2f %% non-zero, correlation %f \"%(\n            expression_names[yidx],\n                (multi_y[:,yidx].getnnz()/multi_y[:,yidx].shape[0]) * 100,\n            access_names[xidx],\n                (multi_x[:,xidx].getnnz()/multi_x[:,xidx].shape[0]) * 100,\n            corr[0][1]\n        ) )","metadata":{"tags":[],"execution":{"iopub.status.busy":"2022-09-17T11:01:48.031518Z","iopub.execute_input":"2022-09-17T11:01:48.031884Z","iopub.status.idle":"2022-09-17T11:01:56.035169Z","shell.execute_reply.started":"2022-09-17T11:01:48.031853Z","shell.execute_reply":"2022-09-17T11:01:56.033939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Because the data is sparse, only a limited number of cells are scattered. Others are aligned on the vertical or horizontal lines, because the other parameter is 0.<br>\nFor Multiome cases, there are 228942 inputs and 23418 targets. I guess the technology can measure a small part of them. If the gene expression is not caught for the cell, the value is filled with 0 (my guess). So there are gaps between naturally distributed data and vertically or horizontally aligned data.\n## Distribution of non-zero inputs/targets in one cell\nIn submission, we are requested to submit 3512 gene expressions per cell. So I wonder how much of gene expressions in one cell are non-zero.\nLet's plot the distribution of non-zero numbers for inputs and targets","metadata":{}},{"cell_type":"code","source":"fig,axes = plt.subplots(1,2,figsize=(15,3))\n_ = axes[1].hist(multi_y.getnnz(axis=1),1000)\n_ = axes[0].hist(multi_x.getnnz(axis=1),1000)\n_ =axes[0].set_xlabel('Non-zero accesibilities in one cell')\n_ = axes[0].set_ylabel('number of cells')\n_ = axes[0].set_title('Distribution of non-zero inputs in one cell')\n_ = axes[0].set_xlim([0,228942])\n_ = axes[1].set_xlabel('Non-zero gene expressions in one cell')\n_ = axes[1].set_ylabel('number of cells')\n_ = axes[1].set_title('Distribution of non-zero targets in one cell')\n_ = axes[1].set_xlim([0,23418])\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:01:56.039037Z","iopub.execute_input":"2022-09-17T11:01:56.039945Z","iopub.status.idle":"2022-09-17T11:02:06.926530Z","shell.execute_reply.started":"2022-09-17T11:01:56.039895Z","shell.execute_reply":"2022-09-17T11:02:06.925244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, the numbers of accessibilities or gene expressions within the same cell_id are not fixed, Rather, statistically distributed. I guess this is due to the measurement principal of 10X technologies.<br> \nPeak for gene expression distribution is about 3700, this may be why the competition hosts chose 3512 for numuber of gene expressions submission. <br>\n<br>\nNext thing we are interested is, which gene expression is non-zero is random or day/donor dependent. I select donor 32606 days 2 as a sample and plot again.","metadata":{}},{"cell_type":"markdown","source":"## Non-zero gene expressions in the same donor same day cell (donor 32606 days 2)\n### Distribution","metadata":{}},{"cell_type":"code","source":"#multi_ym.donor.unique()  # [32606, 13176, 31800]\n#multi_ym.day.unique()  # [2,3,4,7]\ndo32606d2 = multi_y[(multi_ym.day==2) & (multi_ym.donor==32606)]\n\nfig,axes = plt.subplots(1,2,figsize=(15,3))\n_ = axes[0].hist(do32606d2.getnnz(axis=1),1000)\n_ = axes[0].set_xlabel('Non-zero accesibilities in one cell')\n_ = axes[0].set_ylabel('number of cells')\n_ = axes[0].set_title('Distribution of non-zero inputs in one cell')\n_ = axes[0].set_xlim([0,228942])\n_ = axes[1].hist(do32606d2.getnnz(axis=1),1000)\n_ = axes[1].set_xlabel('Non-zero gene expressions in one cell')\n_ = axes[1].set_ylabel('number of cells')\n_ = axes[1].set_title('Distribution of non-zero targets in one cell')\n_ = axes[1].set_xlim([0,23418])\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:02:06.927952Z","iopub.execute_input":"2022-09-17T11:02:06.928438Z","iopub.status.idle":"2022-09-17T11:02:12.298381Z","shell.execute_reply.started":"2022-09-17T11:02:06.928405Z","shell.execute_reply":"2022-09-17T11:02:12.297101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Distribution of donor 32606 day 2 is not so different with over all.<br> From now, I focuse on targets, because targets are what we have to predict. ","metadata":{}},{"cell_type":"markdown","source":"### Which gene expressions are non-zero\nFirst I plotted if gene expression is non-zero(1) or 0(0) against columns index. From what I can see in a small graph is non-zero covers all gene expression.(Althoug they are sparsely covered) Then I magnified the horizontal axis (columns index) to 50. Surprised or not, non zero-genes appear randomly. As measured cell number grows, measurement covers denser. ","metadata":{}},{"cell_type":"code","source":"fig,axes = plt.subplots(1,3,figsize=(15,3))\nfor i in range(3):\n    if i==0:\n        _ = axes[i].plot(do32606d2[1:20].toarray().T !=0)\n        _ = axes[i].set_title('plotted with 20 cells')\n    if i==1:\n        _ = axes[i].plot(do32606d2[1:20,0:50].toarray().T !=0)\n        _ = axes[i].set_title('plotted with 20 cells')\n    if i==2:\n        _ = axes[i].plot(do32606d2[1:200,0:50].toarray().T !=0)\n        _ = axes[i].set_title('plotted with 200 cells')\n    _ = axes[i].set_xlabel('column index')\n    _ = axes[i].set_ylabel('Gene expression non-zero')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:02:12.299596Z","iopub.execute_input":"2022-09-17T11:02:12.299969Z","iopub.status.idle":"2022-09-17T11:02:15.124413Z","shell.execute_reply.started":"2022-09-17T11:02:12.299938Z","shell.execute_reply":"2022-09-17T11:02:15.123180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next question is, is this because of cell_type? So, I narrow down the data to donor 32606 days 2, cell type NeuP\n## Non-zero gene expressions in the same donor same day cell, same cell-type <br>(donor 32606 days 2 cell-type NeuP)","metadata":{}},{"cell_type":"code","source":"do32606d2_NeuP = multi_y[(multi_ym.day==2) & (multi_ym.donor==32606) &(multi_ym.cell_type=='NeuP')]\nfig,axes = plt.subplots(1,3,figsize=(15,3))\nfor i in range(3):\n    if i==0:\n        _ = axes[i].plot(do32606d2_NeuP[1:20].toarray().T !=0)\n        _ = axes[i].set_title('20 cells')\n    if i==1:\n        _ = axes[i].plot(do32606d2_NeuP[1:20,0:50].toarray().T !=0)\n        _ = axes[i].set_title('20 cells')\n    if i==2:\n        _ = axes[i].plot(do32606d2_NeuP[1:200,0:50].toarray().T !=0)\n        _ = axes[i].set_title('200 cells')\n    _ = axes[i].set_xlabel('column index')\n    _ = axes[i].set_ylabel('Gene expression non-zero')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:02:15.126316Z","iopub.execute_input":"2022-09-17T11:02:15.126751Z","iopub.status.idle":"2022-09-17T11:02:18.882889Z","shell.execute_reply.started":"2022-09-17T11:02:15.126713Z","shell.execute_reply":"2022-09-17T11:02:18.881577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Even if donor day cell_type conditions are the same, measured gene expressions are randomly changed. Some gene expression return non-zero values for some cell_ids, but return zero for other cell_ids, even if conditions are the same.","metadata":{}},{"cell_type":"markdown","source":"Now that,we took a very clos look at non-zero data, let's look at wider area.\n# Heatmaps\nThis time I use matplotlib' spy plot, which is a visualization tool for sparse matrix. Visualize the same data as above (donor 32606 days 2 cell-type NeuP) but cells are 500.","metadata":{}},{"cell_type":"code","source":"fix,axes = plt.subplots(1,1,figsize = (30,12))\n_ = axes.spy(do32606d2_NeuP[0:500,0:1000],aspect='auto',  markersize=1.5)\n_ = axes.set_title(\"donor 32606 , days 2 , cell type NeuP\")\n_ = axes.set_xlabel('gene expression[gene_id index]')\n_ = axes.set_ylabel('cells ')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:02:18.884742Z","iopub.execute_input":"2022-09-17T11:02:18.885090Z","iopub.status.idle":"2022-09-17T11:02:19.396602Z","shell.execute_reply.started":"2022-09-17T11:02:18.885060Z","shell.execute_reply":"2022-09-17T11:02:19.395425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We could not see the density with overlaid plots. Now we can see there are dense gene expressions and sparse gene expressions.<br>\nNext, look at the difference in days.\n## Dependency on days","metadata":{}},{"cell_type":"code","source":"fix,axes = plt.subplots(4,1,figsize = (30,15))\nplt.subplots_adjust(hspace=0.6)\nax_i = 0\nfor d in multi_ym.day.unique():\n    axes[ax_i].spy(multi_y[(multi_ym.day==d) & (multi_ym.donor==32606)][0:1200,:],\n                aspect='auto',  markersize=0.01)\n    print(multi_y[(multi_ym.day==d) & (multi_ym.donor==32606)].shape)\n    axes[ax_i].set_title(\"donor 32606 , days %s\"%d)\n    axes[ax_i].set_xlabel('gene expression[gene_id index]')\n    axes[ax_i].set_ylabel('cells ')\n    ax_i += 1\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:02:19.397992Z","iopub.execute_input":"2022-09-17T11:02:19.398363Z","iopub.status.idle":"2022-09-17T11:02:42.510460Z","shell.execute_reply.started":"2022-09-17T11:02:19.398322Z","shell.execute_reply":"2022-09-17T11:02:42.508795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is not a big difference between days. But day 7 is a little sparser than other days. I am not sure if this is just a batch effect or if cells' states really changed.\nWithin the same day, we can recognize horizontal lines caused by ununiformity. There is a possibility that thick or thin lines belong to some specific cell types. Let's see the same data grouped by cell types.\n## Dependency on cell types","metadata":{}},{"cell_type":"code","source":"fix,axes = plt.subplots(7,1,figsize = (30,20))\nplt.subplots_adjust(hspace=0.8)\nax_i = 0\nfor ct in multi_ym.cell_type.unique():\n    axes[ax_i].spy(multi_y[(multi_ym.day==2) & (multi_ym.donor==32606) &(multi_ym.cell_type==ct),:][0:700,:],\n                aspect='auto',  markersize=0.01)\n    print(multi_y[(multi_ym.day==2) & (multi_ym.donor==32606) &(multi_ym.cell_type==ct),:].shape)\n    axes[ax_i].set_title(ct)\n    axes[ax_i].set_xlabel('gene expression[gene_id index]')\n    axes[ax_i].set_ylabel('cells ')\n    ax_i += 1\n","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:02:42.512132Z","iopub.execute_input":"2022-09-17T11:02:42.512466Z","iopub.status.idle":"2022-09-17T11:03:08.212821Z","shell.execute_reply.started":"2022-09-17T11:02:42.512436Z","shell.execute_reply":"2022-09-17T11:03:08.210398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"MoP and BP look very sparse, but those have one order smaller data and are plotted with lower density. Let's ignore these two. I can not recognize differences between cell types, but I still recognize horizontal lines. I am not sure if these are noise like batch effects or real data reflecting the cell's state. My guess is noise because the stripe pattern looks artificial.","metadata":{}},{"cell_type":"markdown","source":"# Targets self-correlation\nOne more thing I am curious about is the mixture of barcodes.\nIn [\"Chromium\nNext GEM\nSingle Cell Multiome\nATAC + Gene Expression\" USER GUIDE](https://assets.ctfassets.net/an68im79xiti/30evMTksKKnv78Ii69VlUo/0553cd05226dfd6853ee9f4d3980073b/CG000338_ChromiumNextGEM_Multiome_ATAC_GEX_User_Guide_RevF.pdf) P17, 10x Barcoded Gel Beads are encapsulated with a cell. Gel Beads are shown in many colors, which possibly implies Gel Beads are not all the same.<br>\nI have two assumptions. (Although I am a Laymon for biology) <br> \nA:  Barcoded Gel Beads have limited numbers of variations. Say, each variation has 3700 kinds of barcodes, and there are 6 variations of Gel Beads, 3700 X 6 = 22200 and this covers all gene expressions.<br> \nB:  Each Barcoded Gel Bead has a randomly selected unique mixture of barcodes, and many Gel Beads eventually cover all gene expressions.<br>\nThere is one thing we can do to get a hint, which is to calculate the self-correlation of gene expression.\nIf we can find a strong correlation between some cells, it is highly probable that Barcoded Gel Beards have a limited number of variations.<br>\nLet's see. I chose the number of cells to 100, considering 6 variations might cover all gene expression. If we choose this number unnecesarily big, we have trouble seeing the dots.","metadata":{}},{"cell_type":"code","source":"do32606d2_NeuP = multi_y[(multi_ym.day==2) & (multi_ym.donor==32606) &(multi_ym.cell_type=='NeuP')]\ndo32606d2_NeuP = do32606d2_NeuP[:100,:]\ncor_matrix = np.corrcoef(do32606d2_NeuP.toarray()!=0)  ### correlation of only zero or non-zero\nfig,ax = plt.subplots(figsize=(12,12))\npos = ax.imshow(cor_matrix)\nfig.colorbar(pos, ax=ax)\nax.set_title('donor 32606 days 2 cell-type NeuP : Gene expression non-zero self-correlation')\n_ = ax.set_xlabel('cell number in donor 32606 days 2 NeuP')\n_ = ax.set_ylabel('cell number in donor 32606 days 2 NeuP')            ","metadata":{"execution":{"iopub.status.busy":"2022-09-17T11:04:18.180774Z","iopub.execute_input":"2022-09-17T11:04:18.181186Z","iopub.status.idle":"2022-09-17T11:04:19.531264Z","shell.execute_reply.started":"2022-09-17T11:04:18.181155Z","shell.execute_reply":"2022-09-17T11:04:19.530123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can not find bright dots between cells, which implies Barcodes are ramdomly mixed into Gel Beads. <br>\nBut the first figure in section heatmaps shows vertical lines, which means some gene expressions are measured in most of the cells. This denies above guess.\nOne thing we should notice is there are brighter lines and darker lines, which have strong or weak correlations with the other cells. But, I have no idea what this means. Or, my assumpsion is all wrong.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}