{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59094,"databundleVersionId":7010844,"sourceType":"competition"},{"sourceId":6636461,"sourceType":"datasetVersion","datasetId":3830903}],"dockerImageVersionId":30587,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# SCP: EDA which makes sense ⭐️⭐️⭐️⭐️⭐️\n\nThis notebook points out some little-known aspects of the competition data:\n- t-scores are better to work with than log10pvalues.\n- Cell types come in certain proportions, and it looks like some cells are annotated with the wrong cell type.\n- 20 of the 145 compounds should be seen as outliers.\n- The data contains artefacts produced by Limma. The artefacts can be exploited for the competition, but they distract from the real objective, which is *helping to accurately predict chemical perturbations in new cell types*.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"!pip install scanpy\n","metadata":{"scrolled":true,"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:22:30.860297Z","iopub.execute_input":"2023-12-01T00:22:30.861013Z","iopub.status.idle":"2023-12-01T00:22:44.903897Z","shell.execute_reply.started":"2023-12-01T00:22:30.860984Z","shell.execute_reply":"2023-12-01T00:22:44.902886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport scipy.stats\nfrom scipy.stats import norm, multinomial\n\nimport scanpy as sc\n\nnp.set_printoptions(edgeitems=3)\npd.set_option(\"min_rows\", 10)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:22:44.905930Z","iopub.execute_input":"2023-12-01T00:22:44.906190Z","iopub.status.idle":"2023-12-01T00:22:50.499364Z","shell.execute_reply.started":"2023-12-01T00:22:44.906167Z","shell.execute_reply":"2023-12-01T00:22:50.498411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some functions","metadata":{}},{"cell_type":"code","source":"def mean_rowwise_rmse(y_true, y_pred):\n    \"\"\"Competition metric\n    \n    Calling convention like in sklearn.metrics\n    \"\"\"\n    mrrmse = np.sqrt(np.square(y_true - y_pred).mean(axis=1)).mean()\n    return mrrmse\n\ndef de_to_t_score(de):\n    \"\"\"Convert log10pvalues to t-scores\n    \n    Parameter:\n    de: array or DataFrame of log10pvalues\n    \n    Return value:\n    t_score: array or DataFrame of t-scores\n    \"\"\"\n    p_value = 10 ** (-np.abs(de))\n    return - scipy.stats.t.ppf(p_value / 2, df=420) * np.sign(de)\n#     return - norm.ppf(p_value / 2) * np.sign(de)\n\ndef t_score_to_de(t_score):\n    \"\"\"Convert t-scores to log10pvalues (inverse of de_to_t_score)\n    \n    Parameter:\n    t_score: array or DataFrame of t-scores\n    \n    Return value:\n    de: array or DataFrame of log10pvalues\n    \"\"\"\n    p_value = scipy.stats.t.cdf(- np.abs(t_score), df=420) * 2\n#     p_value = norm.cdf(- np.abs(t_score)) * 2\n    p_value = p_value.clip(1e-180, None)\n    return - np.log10(p_value) * np.sign(t_score)\n\ndef mode(series):\n    \"\"\"Mode of a pandas series\"\"\"\n    return series.value_counts().index[0]\n","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-12-01T00:22:50.500322Z","iopub.execute_input":"2023-12-01T00:22:50.500810Z","iopub.status.idle":"2023-12-01T00:22:50.508675Z","shell.execute_reply.started":"2023-12-01T00:22:50.500787Z","shell.execute_reply":"2023-12-01T00:22:50.507417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading the data\n\nWe read three files:\n\n- `de_train.parquet` --> transformed into `t_train`\n- `adata_obs_meta.csv` containing the metadata for every cell --> transformed into `cell_count` of shape (614,)\n- `train_or_control_bulk_by_cell_type_adata.h5ad` containing the pseudobulked RNA counts --> `fixed_bulk_adata` of shape 2558 × 18211 --> `rna_count` of shape (614, 18211) and `rna_count_dmso` of shape (6, 18211)","metadata":{}},{"cell_type":"code","source":"fn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\n# fn = '/kaggle/input/scp-3-merge/de_train.parquet'\nde_train = pd.read_parquet(fn)# , index_col = 0)\n# display(de_train)\n\n# fn = '/kaggle/input/open-problems-single-cell-perturbations/id_map.csv'\n# id_map = pd.read_csv(fn, index_col = 0)\n\n# 18211 genes\ngenes = de_train.columns[5:] \nt_train = de_to_t_score(de_train.set_index(['cell_type', 'sm_name'])[genes])\n\n# All 146 sm_names\nsm_names = sorted(de_train.sm_name.unique())\n# Determine the 17 compounds (including the two control compounds) \n# with data for almost all cell types\ntrain_sm_names = de_train.query(\"cell_type == 'B cells'\").sm_name.sort_values().values\n# The other 129 sm_names\ntest_sm_names = [sm for sm in sm_names if sm not in train_sm_names]\n# The three control sm_names\ncontrols3 = ['Dabrafenib', 'Belinostat', 'Dimethyl Sulfoxide']\n\n# All 6 cell types\n# cell_types = sorted(de_train.cell_type.unique())\ncell_types = ['NK cells', 'T cells CD4+', 'T cells CD8+', \n              'T regulatory cells', 'B cells', 'Myeloid cells']\n# Determine the 4 cell types with data for almost all compounds\n# ['NK cells', 'T cells CD4+', 'T cells CD8+', 'T regulatory cells']\ntrain_cell_types = de_train.query(\"sm_name == 'Vorinostat'\").cell_type.sort_values().values\n# The other 2 cell types: ['B cells', 'Myeloid cells']\ntest_cell_types = [ct for ct in cell_types if ct not in train_cell_types]\n","metadata":{"execution":{"iopub.status.busy":"2023-12-01T00:22:50.510937Z","iopub.execute_input":"2023-12-01T00:22:50.511811Z","iopub.status.idle":"2023-12-01T00:23:08.855674Z","shell.execute_reply.started":"2023-12-01T00:22:50.511789Z","shell.execute_reply":"2023-12-01T00:23:08.855057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata_obs = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv')\n\ncell_count = adata_obs.groupby(['cell_type', 'sm_name']).size().reindex_like(t_train)\n","metadata":{"execution":{"iopub.status.busy":"2023-12-01T00:23:08.856710Z","iopub.execute_input":"2023-12-01T00:23:08.857687Z","iopub.status.idle":"2023-12-01T00:23:09.757787Z","shell.execute_reply.started":"2023-12-01T00:23:08.857645Z","shell.execute_reply":"2023-12-01T00:23:09.756892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fixed_bulk_adata = sc.read_h5ad('/kaggle/input/open-problems-single-cell-perturbations-optional/train_or_control_bulk_by_cell_type_adata.h5ad')\nfixed_bulk_adata.X = fixed_bulk_adata.layers['counts']\n\nlincs_id_mapping_df = pd.read_parquet('/kaggle/input/open-problems-single-cell-perturbations-optional/lincs_id_compound_mapping.parquet')\nfixed_bulk_adata.obs['sm_name'] = lincs_id_mapping_df.set_index('compound_id')['sm_name'].reindex(fixed_bulk_adata.obs.compound_id).values\n\nadata = pd.DataFrame(fixed_bulk_adata.X,\n                     index=pd.MultiIndex.from_frame(fixed_bulk_adata.obs[['cell_type', 'sm_name']]),\n                     columns=fixed_bulk_adata.var_names).astype(int)\nrna_count = adata.groupby(['cell_type', 'sm_name'], observed=True).sum().rename_axis(columns='gene')\nrna_count_dmso = rna_count.query(\"sm_name == 'Dimethyl Sulfoxide'\")\nrna_count = rna_count.reindex(de_train[['cell_type', 'sm_name']])\ndisplay(rna_count_dmso)","metadata":{"execution":{"iopub.status.busy":"2023-12-01T00:23:09.759324Z","iopub.execute_input":"2023-12-01T00:23:09.759776Z","iopub.status.idle":"2023-12-01T00:23:17.566671Z","shell.execute_reply.started":"2023-12-01T00:23:09.759748Z","shell.execute_reply":"2023-12-01T00:23:17.565673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# T-scores are better than log10pvalues\n\nLimma performs t-tests. t-scores are (almost) normally distributed, which is good for machine learning inputs. For this competition, the t-scores were nonlinearly transformed to log10pvalues. The transformation squeezes the nice bell shape into a distribution with a much higher kurtosis.\n\nMy machine learning models perform better if I transform the log10pvalues into t-score in a preprocessing step, predict t-scores, and transform the predictions back afterwards. Perhaps working with log-fold changes or even RNA counts would be even better.\n","metadata":{}},{"cell_type":"code","source":"ts = np.linspace(-28, 28, 61)\nls = t_score_to_de(ts) # between -198 and +198\nspace = 1.03\nplt.figure(figsize=(10, 6))\nplt.subplot(3, 1, 1)\nplt.hist(de_train[genes].values.ravel(),\n         bins=500,\n         color='b',\n         label='Histogram of log10pvalues')\nplt.xticks(np.linspace(-180, 180, 7))\nplt.yticks([])\nplt.xlim(-198*space, 198*space)\nplt.legend()\n\nplt.subplot(3, 1, 2)\nfor t, l in zip(list(ts), list(ls / 198 * 30)):\n    plt.plot([t, l], [0, 1], c='orange')\nplt.yticks([0, 1], ['t-score', 'signed -log10(pvalue)'])\nplt.xticks([])\nplt.xlim(-30*space, 30*space)\nplt.twiny()\nplt.xticks([])\nplt.xlim(-198*space, 198*space)\n\nplt.subplot(3, 1, 3)\nplt.hist(t_train.values.ravel(),\n         bins=500,\n         color='g',\n         label='Histogram of t-scores')\nplt.yticks([])\nplt.xlim(-30*space, 30*space)\nplt.legend()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:17.567735Z","iopub.execute_input":"2023-12-01T00:23:17.568003Z","iopub.status.idle":"2023-12-01T00:23:20.195785Z","shell.execute_reply.started":"2023-12-01T00:23:17.567982Z","shell.execute_reply":"2023-12-01T00:23:20.195087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Don't trust the cell types!\n\nLet's recapitulate the course of the experiment in a simplified form. We can imagine an experimenter who is in front of a large pot of human peripheral blood mononuclear cells. The pot contains a mixture of six cell types in certain proportions. T cells CD4+ are most abundant (42 %), only 2 % are T regulatory cells:","metadata":{}},{"cell_type":"code","source":"# Cell type ratios (extrapolated from 17 train_sm_names)\ntemp = adata_obs.groupby(['cell_type', 'sm_name']).size().unstack().loc[cell_types]\ncell_type_ratio = temp[list(train_sm_names) + ['Dimethyl Sulfoxide']].sum(axis=1)\ncell_type_ratio /= cell_type_ratio.sum()\nplt.pie(cell_type_ratio, labels=cell_type_ratio.index, autopct=\"%.0f%%\")\nplt.show()\n\ndrops = temp[list(train_sm_names) + ['Dimethyl Sulfoxide']].dropna(axis=1).sum(axis=0)\ndrops = drops[['CHIR-99021', 'Crizotinib', 'Dactolisib', 'Foretinib',\n               'Idelalisib', 'LDN 193189', 'Linagliptin', 'MLN 2238', \n               'O-Demethylated Adapalene', 'Palbociclib', 'Penfluridol',\n               'Porcn Inhibitor III', 'R428']]\nprint(f\"Mean drop size: {drops.mean():.0f}, standard deviation of drop size: {drops.std():.0f}\")\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:20.196685Z","iopub.execute_input":"2023-12-01T00:23:20.197050Z","iopub.status.idle":"2023-12-01T00:23:20.364224Z","shell.execute_reply.started":"2023-12-01T00:23:20.197020Z","shell.execute_reply":"2023-12-01T00:23:20.363366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The experimenter now takes 145 droplets out of the large pot. Every droplet contains 1551 ± 244 cells. If we counted the cells of every cell type in a droplet, we'd see a multinomial distribution. The 145 droplets might be composed like in the following bar chart (fictitious data, sorted from smallest to largest droplet):","metadata":{}},{"cell_type":"code","source":"drops = norm.rvs(loc=1551, scale=244, size=145).astype(int)\ndrops.sort()\ncell_types_in_drops = np.array([multinomial.rvs(n, cell_type_ratio) for n in drops])\n\ndef plot_stacked_bar_chart(cell_types_in_drops, title, xticks=None, xticklabels=None):\n    \"\"\"Plot a stacked bar chart of cell counts\n    \n    The plot has one vertical bar per drop, and every cell type gets\n    its own color.\n    \n    Parameters:\n    cell_types_in_drops: array of shape (n_drops, n_cell_types)\n    xticks, xticklabels: parameters for plt.xticks, shape (n_drops,)\n    \"\"\"\n    # Add a column of zeros to the left and compute the cumulative sums\n    cc = np.hstack([np.zeros((len(cell_types_in_drops), 1)), cell_types_in_drops])\n    cc_cs = cc.cumsum(axis=1)\n\n    plt.figure(figsize=(15, 4))\n    for i in range(len(cell_types)):\n        plt.bar(np.arange(len(cc_cs)),\n                cc_cs[:,i+1] - cc_cs[:,i],\n                bottom=cc_cs[:,i],\n                label=cell_types[i])\n    plt.legend()\n    plt.title(title)\n    plt.xlabel('Drop')\n    plt.ylabel('Cell count')\n    plt.ylim(0, 2300)\n    if xticks is not None:\n        plt.xticks(xticks, xticklabels, rotation=90)\n    plt.show()\n\nplot_stacked_bar_chart(cell_types_in_drops, 'Composition of droplets when experiment starts')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:20.365210Z","iopub.execute_input":"2023-12-01T00:23:20.365465Z","iopub.status.idle":"2023-12-01T00:23:21.872474Z","shell.execute_reply.started":"2023-12-01T00:23:20.365443Z","shell.execute_reply":"2023-12-01T00:23:21.871564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the next step, the experimenter adds 145 substances to the 145 droplets and waits 24 hours. After 24 hours the cells are analyzed. If we count the cells again, we get the following picture, as taken from the competition's training data (cell counts for the test data are hidden):","metadata":{}},{"cell_type":"code","source":"xtick_compounds = ['Ganetespib (STA-9090)', 'LY2090314', 'CGM-097',\n                   'MLN 2238', 'CEP-18770 (Delanzomib)', \n                   'Oprozomib (ONX 0912)', 'IN1451']\n\ncc = cell_count.unstack('cell_type') # shape (n_compounds, n_cell_types)\ncc = cc.query(\"~sm_name.isin(@controls3)\") # remove the long bars of the control compounds\ncc.loc[:, train_cell_types] = cc.loc[:, train_cell_types].fillna(0)\n\n# Sort by the number of cells of the four training cell types\ncc['total'] = cc[train_cell_types].sum(axis=1)\ncc.sort_values('total', inplace=True)\n\n# Plot\nsorted_compound_names = cc.index.get_level_values('sm_name')\nxticks = np.arange(len(cc))[sorted_compound_names.isin(xtick_compounds)]\nxticklabels = sorted_compound_names[sorted_compound_names.isin(xtick_compounds)]\nplot_stacked_bar_chart(cc[cell_types].values, 'Composition of droplets after 24 hours', xticks=xticks, xticklabels=xticklabels)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:21.874920Z","iopub.execute_input":"2023-12-01T00:23:21.875264Z","iopub.status.idle":"2023-12-01T00:23:23.343570Z","shell.execute_reply.started":"2023-12-01T00:23:21.875243Z","shell.execute_reply":"2023-12-01T00:23:23.342728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this diagram we first see that some compounds are so toxic that in some droplets less than 100 cells survive. These droplets are represented by the leftmost bars in the bar chart.\n\nThe second observation is much more important: The long red part in the bars for Oprozomib and IN1451 show that these droplets contain several hundred T regulatory cells — much more than at the start of the experiment. Other compounds (e.g., CGM-079) have too many T cells CD8+. How can we interpret this observation?\n1. Does IN1451 incite the T regulatory cells to multiply so that we have five times more of them after 24 hours? No.\n1. Does IN1451 magically convert NK cells into T regulatory cells? No.\n1. Does IN1451 affect the cells in such a way that they are misclassified? Maybe.\n\nDiscussing differential gene expression for specific cell types becomes pointless if the cells change their type during the experiment. For the Kaggle competition this means that we have to deal with many outliers: Beyond the at least five toxic compounds, there are at least seven compounds which change the cells' types. Differential expression for these outliers is hard to model. They make cross-validation unreliable, and the outliers in the private leaderboard can't even be predicted by probing the public leaderboard.","metadata":{}},{"cell_type":"markdown","source":"# Cell count shouldn't affect differential gene expression\n\nDoes gene expression in a cell depend on how many cells are in the experiment? Theoretically, it doesn't. A cell behaves the same way whether there are 10 cells in the experiment or 10000. We'd expect, however, a difference in the significance of the experimental results: An experiment with 10000 cells should give more precise measurements than a 10-cell experiment: As the cell count grows, variance of the measurements should decrease, t-score should be farther away from zero, and pvalues should decrease.\n\nThe competition data don't fulfill this expectation. If we plot the mean t-scores versus the cell count for the 602 cell type–compound combinations (excluding the control compounds), we see a linear relationship: For every cell type, compounds with lower cell counts have positive t-score means, and compounds with higher cell counts have negative t-score means. This correlation between cell counts and t-scores shouldn't exist. It is an artefact of Limma rather than a biological effect.\n\nYou can plot the diagram with median or variance instead of mean — it will look similar. You can even compare the cell counts to the first principal component of the t-scores and see the same correlation. \n\n","metadata":{}},{"cell_type":"code","source":"without_controls = de_train[~de_train.sm_name.isin(controls3).values]\nt_scores_without_controls = de_to_t_score(without_controls.set_index(['cell_type', 'sm_name'])[genes])\nmeans = t_scores_without_controls.mean(axis=1)\n# means = t_scores_without_controls.median(axis=1)\n# means = t_scores_without_controls.var(axis=1)\ndf = pd.DataFrame({'cell count': adata_obs.groupby(['cell_type', 'sm_name']).size().reindex(means.index),\n                   'mean t-score': means})\ndf.reset_index(inplace=True)\n\n# sns.lmplot(\n#     data=df,\n#     x=\"cell count\", y=\"mean t-score\", hue=\"cell_type\",\n#     scatter_kws={'s': 5},\n#     fit_reg=False,\n#     height=4, aspect=2\n# )\n# plt.title('Mean t-scores depend on cell count')\n# plt.show()\n\ncolor_discrete_map = {ct: col for (ct, col) in zip(cell_types, plt.rcParams['axes.prop_cycle'].by_key()['color'])}\nfig = px.scatter(df, x=\"cell count\", y=\"mean t-score\",\n                 color=\"cell_type\",\n                 color_discrete_map=color_discrete_map,\n                 title='Mean t-scores depend on cell count',\n                 hover_data=['sm_name'])\nfig.show()\n\n# https://plotly.com/python-api-reference/generated/plotly.express.scatter.html","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:23.344623Z","iopub.execute_input":"2023-12-01T00:23:23.344891Z","iopub.status.idle":"2023-12-01T00:23:39.254061Z","shell.execute_reply.started":"2023-12-01T00:23:23.344865Z","shell.execute_reply":"2023-12-01T00:23:39.253206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can now put together a list of 20 compounds which are to be considered outliers. Notice that we don't declare single rows of the dataset to be outliers, but all 86 rows related to the 20 compounds:\n\n```\nAT13387                           only 7 T regulatory cells\nAlvocidib                         ≤10 for several cell types\nBAY 61-3606                       mean t-score for CD8+ cells > 2\nBMS-387032                        only 10 T cells CD8+, no Myeloid cells\nBelinostat                        control compound with too many cells\nCEP-18770 (Delanzomib)            ≤10 for several cell types\nCGM-097                           too many T cells CD8+\nCGP 60474                         ≤10 for several cell types\nDabrafenib                        control compound with too many cells\nGanetespib (STA-9090)             only 4 T regulatory cells, too many NK cells\nI-BET151                          too many T cells CD8+\nIN1451                            ≤10 for several cell types\nLY2090314                         only 6 T cells CD8+\nMLN 2238                          ≤10 for several cell types\nOprozomib (ONX 0912)              ≤10 for several cell types\nProscillaridin A;Proscillaridin-A ≤10 for several cell types\nResminostat                       no T cells CD8+\nScriptaid                         only 2 T regulatory cells\nUNII-BXU45ZH6LI                   only 6 T cells CD8+\nVorinostat                        only 1 T regulatory cell\n```","metadata":{}},{"cell_type":"markdown","source":"After removing the outliers, the diagram looks much cleaner. The variance of the cell counts remains. It is a source of noise which impedes the correct interpretation (and prediction) of differential expressions. Maybe we'd get cleaner data if we equalized the cell counts before library size normalization. This would amount to throwing away a part of the measurements, which isn't desirable either.","metadata":{}},{"cell_type":"code","source":"removed_compounds = ['AT13387', 'Alvocidib', 'BAY 61-3606', 'BMS-387032', \n                     'Belinostat', 'CEP-18770 (Delanzomib)', 'CGM-097', 'CGP 60474', \n                     'Dabrafenib', 'Ganetespib (STA-9090)', 'I-BET151', 'IN1451', \n                     'LY2090314', 'MLN 2238', 'Oprozomib (ONX 0912)', \n                     'Proscillaridin A;Proscillaridin-A', 'Resminostat',\n                     'Scriptaid', 'UNII-BXU45ZH6LI', 'Vorinostat']\ndf = df.query('~sm_name.isin(@removed_compounds)')\n\n# sns.lmplot(\n#     data=df,\n#     x=\"cell count\", y=\"mean t-score\", hue=\"cell_type\",\n#     scatter_kws={'s': 5},\n#     line_kws={'lw': 2},\n#     ci=None, robust=True,\n#     height=4, aspect=2\n# )\n# plt.title('After outlier removal')\n# plt.show()\n\nfig = px.scatter(df, x=\"cell count\", y=\"mean t-score\",\n                 color=\"cell_type\",\n                 color_discrete_map=color_discrete_map,\n                 trendline='ols',\n                 title='After outlier removal',\n                 hover_data=['sm_name'])\nfig.show()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:39.255319Z","iopub.execute_input":"2023-12-01T00:23:39.255568Z","iopub.status.idle":"2023-12-01T00:23:40.430494Z","shell.execute_reply.started":"2023-12-01T00:23:39.255548Z","shell.execute_reply":"2023-12-01T00:23:40.429538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# A mixture of distributions\n\nA histogram of a single row of the training data (18211 t-scores for T cells CD8+ treated with Scriptaid) shows that the distribution is multimodal.\n\nThe highest mode consists of 269 genes with a t-score of -3.769. It turns out that these are the 269 genes which are never expressed in T cells CD8+, neither with the negative control nor with any other compound. Isn't this strange? A gene which is never expressed in the whole experiment should have a log-fold change of zero and should not get a t-score at all (because t-score computation involves a division by the variance, and the variance of a never-expressed gene is zero).\n\nFor Myeloid cells treated with Foretinib, 3856 genes are not expressed (RNA count of zero), yet most of them have a positive t-score. Their highest t-score is 6.228 (resulting in a pvalue of 4e-10 and a log10pvalue of 9.33). If an RNA count is zero, the corresponding log-fold-change (and t-score) should never be positive.\n\nWe may say that the distribution of the values is a mixture of two distributions:\n1. The values for the genes which are expressed (blue) have a more or less bell-shaped distribution.\n2. The values for the genes which are not expressed (orange) have a distribution with an unusual shape, and it is strange that positive differential expressions are reported when not a single piece of RNA is counted.\n\nWhat we see here is an artefact of Limma, which affects every row of the datset. It suggests that Limma output can be biased and is not ideal for investigating cell-type translation of differential expressions.\n","metadata":{}},{"cell_type":"code","source":"def plot_histogram(ct, sm):\n    \"\"\"Plot the histogram of gene t-scores for given cell type and compound\"\"\"\n    \n    # Get the row of the matrix\n    this_t_score = t_train.loc[(ct, sm)] # 1 row of the t-score matrix comprising 18211 genes\n    \n    # Distinguish genes which are expressed / not expressed for this compound\n    compound_zero = (rna_count.loc[ct, sm] == 0).values.ravel()\n    print(f\"Genes expressed in {ct} {sm}:     {(~compound_zero).sum():5}\")\n    print(f\"Genes not expressed in {ct} {sm}: {compound_zero.sum():5}\")\n\n    # Look at the mode\n    t_mode = mode(this_t_score)\n    this_equal_mode = this_t_score == t_mode # expression is zero for all compounds\n    print(f'Mode: {t_mode:.3f} for {this_equal_mode.sum()} genes not expressed at all in {ct}')\n    \n    cc = cell_count.loc[(ct, sm)] # integer which is shown in the diagram title\n    \n    plt.figure(figsize=(12, 2))\n    plt.hist(this_t_score.values.ravel()[~compound_zero],\n             bins=np.linspace(-15, 15, 601), \n             density=False, \n             label='gene expressed for this compound')\n    plt.hist(this_t_score.values.ravel()[compound_zero], \n             bins=np.linspace(-15, 15, 601),\n             density=False,\n             alpha=0.7,\n             label='gene not expressed for this compound')\n    plt.xlim(-15, 15)\n    plt.legend()\n    plt.title(f\"{ct} {sm} t-score histogram ({cc} cells)\")\n    plt.xlabel('t-score')\n    plt.yticks([])\n    plt.show()\n    \nplot_histogram('T cells CD8+', 'Scriptaid')\nplot_histogram('Myeloid cells', 'Foretinib')\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-01T00:23:40.432068Z","iopub.execute_input":"2023-12-01T00:23:40.432418Z","iopub.status.idle":"2023-12-01T00:23:44.153825Z","shell.execute_reply.started":"2023-12-01T00:23:40.432388Z","shell.execute_reply":"2023-12-01T00:23:44.152869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}