{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.12"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59094,"databundleVersionId":7010844,"sourceType":"competition"},{"sourceId":6636461,"sourceType":"datasetVersion","datasetId":3830903}],"dockerImageVersionId":30626,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# <a id='toc1_'></a>[Open problems - Single cell perturbations](#toc0_)\n## <a id='toc1_1_'></a>[Post-competition EDA: significance of errorneous data](#toc0_)\n### <a id='toc1_1_1_'></a>[by Jalil Nourisa](#toc0_)\nIn our previous work [1], as well as in the findings by Amrbosm [2], several issues in the training data provided for the competition were highlighted. These problems include misclassification of cell types, inadequate cell counts, and invalid differential expression (DE) values due to missing RNA counts. This study delves deeper into these errors and quantifies their impact. We find that nearly 35 percent of the data is potentially erroneous. Notably, these error-prone data display more significance compared to the remaining data. Our findings highlight the limitations of the current training dataset and emphasize the necessity of recalculating DE values from the raw expression data for future studies.  ","metadata":{}},{"cell_type":"markdown","source":"Refs: \\\n[1][https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/461159]  \\\n[2][https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/458661]","metadata":{}},{"cell_type":"code","source":"%pip install scanpy","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:09.666961Z","iopub.execute_input":"2024-01-07T13:38:09.667462Z","iopub.status.idle":"2024-01-07T13:38:20.276334Z","shell.execute_reply.started":"2024-01-07T13:38:09.667426Z","shell.execute_reply":"2024-01-07T13:38:20.274467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Imports","metadata":{}},{"cell_type":"code","source":"# The code is directly taken from Ambrosm [1]\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport scipy.stats\nimport scanpy as sc\ndef 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\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 = 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\n\nnp.set_printoptions(edgeitems=3)\npd.set_option(\"min_rows\", 10)\n\nfn = '/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet'\nde_train = pd.read_parquet(fn)# , index_col = 0)\ndata = de_train.iloc[:,5:].values\n# 18211 genes\ngene_names = de_train.columns[5:] \nt_train = de_to_t_score(de_train.set_index(['cell_type', 'sm_name'])[gene_names])\ndata_df = de_train[gene_names]\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\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\ntrain_cell_types = de_train.query(\"sm_name == 'Vorinostat'\").cell_type.sort_values().values\n\n# obs\nadata_obs = pd.read_csv('/kaggle/input/open-problems-single-cell-perturbations/adata_obs_meta.csv')\ncell_count = adata_obs.groupby(['cell_type', 'sm_name']).size().reindex_like(t_train)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:20.283574Z","iopub.execute_input":"2024-01-07T13:38:20.283883Z","iopub.status.idle":"2024-01-07T13:38:41.852492Z","shell.execute_reply.started":"2024-01-07T13:38:20.283850Z","shell.execute_reply":"2024-01-07T13:38:41.851224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## RNA count obtained from raw expression data\nWe get the data obtained by HOHO [https://www.kaggle.com/datasets/awater1223/open-problems-single-cell-perturbations-optional]\n","metadata":{"execution":{"iopub.status.busy":"2024-01-06T15:55:24.079574Z","iopub.execute_input":"2024-01-06T15:55:24.080046Z","iopub.status.idle":"2024-01-06T15:55:24.088862Z","shell.execute_reply.started":"2024-01-06T15:55:24.080009Z","shell.execute_reply":"2024-01-06T15:55:24.087295Z"}}},{"cell_type":"code","source":"# The code is directly taken from Ambrosm [1]\nfixed_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":"2024-01-07T13:38:41.854381Z","iopub.execute_input":"2024-01-07T13:38:41.854906Z","iopub.status.idle":"2024-01-07T13:38:50.520242Z","shell.execute_reply.started":"2024-01-07T13:38:41.854866Z","shell.execute_reply":"2024-01-07T13:38:50.518743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Identification of outlier samples and compounds","metadata":{"execution":{"iopub.status.busy":"2024-01-07T10:06:58.636213Z","iopub.execute_input":"2024-01-07T10:06:58.636680Z","iopub.status.idle":"2024-01-07T10:06:58.642270Z","shell.execute_reply.started":"2024-01-07T10:06:58.636639Z","shell.execute_reply":"2024-01-07T10:06:58.640844Z"}}},{"cell_type":"markdown","source":"## Cell type misclassification \nWe and Ambrosm previously highlighted concerns regarding the potential misclassification of cell types for certain compounds (refer to [1] and [2] for more details). Here, we identify the outlier compounds and the associated samples due to cell type misclassification. ","metadata":{}},{"cell_type":"markdown","source":"### Relative cell count ratio per sample","metadata":{}},{"cell_type":"code","source":"df_pivot = cell_count.unstack('cell_type')[train_cell_types]\nnormalized_ratios = df_pivot.div(df_pivot.sum(axis=1), axis=0)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:50.524340Z","iopub.execute_input":"2024-01-07T13:38:50.524649Z","iopub.status.idle":"2024-01-07T13:38:50.535346Z","shell.execute_reply.started":"2024-01-07T13:38:50.524626Z","shell.execute_reply":"2024-01-07T13:38:50.534266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Identification of the outlier compounds using histograms \nWe establish lower and upper thresholds for each cell type to detect outliers. For minority cell types, only the upper threshold is considered, acknowledging that their underrepresentation might be due to random sampling variations.\n","metadata":{}},{"cell_type":"code","source":"major_celltypes = ['NK cells', 'T cells CD4+']","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:50.536484Z","iopub.execute_input":"2024-01-07T13:38:50.536728Z","iopub.status.idle":"2024-01-07T13:38:50.545188Z","shell.execute_reply.started":"2024-01-07T13:38:50.536706Z","shell.execute_reply":"2024-01-07T13:38:50.543906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 6))\naxes = axes.flatten()\nthreshold = 3\noutliers_stack = []\n# Iterate through each cell type to create histograms\nfor idx, cell_type in enumerate(train_cell_types):\n    ax = axes[idx]\n    dataa = normalized_ratios[cell_type]\n    dataa = dataa[~np.isnan(dataa)]\n    if np.any(dataa==0):\n        raise ValueError('shouldnt be 0')\n    Q1 = np.percentile(dataa, 25)\n    Q3 = np.percentile(dataa, 75)\n    IQR = Q3 - Q1\n    upper_bound_log = Q1 + threshold * IQR\n    ax.axvline(x=upper_bound_log, color='red', linestyle='--', linewidth=1)\n    ax.hist(dataa, bins=40, label=cell_type, color='skyblue', edgecolor='black')\n    ax.set_title(f'{cell_type}')\n    ax.set_xlabel('Relative ratio of cell count (log10)')\n    ax.set_ylabel('Frequency')\n\n    if cell_type in major_celltypes:\n        lower_bound_log = Q1 - threshold * IQR\n        ax.axvline(x=lower_bound_log, color='green', linestyle='--', linewidth=1)\n        exceeding_drugs = [drug for drug, value in dataa.items() if (value < lower_bound_log) or (value > upper_bound_log)]\n    else:\n        exceeding_drugs = [drug for drug, value in dataa.items() if (value > upper_bound_log)]\n\n    outliers_stack+= exceeding_drugs\n    if True:\n        # annote outlier drugs\n        exceeding_drugs = [item[0:10] if len(item)>10 else item for item in exceeding_drugs ]\n        drugs_text = '\\n'.join(exceeding_drugs)\n        # Place the collected names in a box\n        if exceeding_drugs:\n            ax.text(0.82, 0.97, drugs_text, transform=ax.transAxes, fontsize=8, \n                    verticalalignment='top', bbox=dict(boxstyle=\"round\", alpha=0.5,  facecolor='white'))\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:50.546454Z","iopub.execute_input":"2024-01-07T13:38:50.546831Z","iopub.status.idle":"2024-01-07T13:38:51.560912Z","shell.execute_reply.started":"2024-01-07T13:38:50.546801Z","shell.execute_reply":"2024-01-07T13:38:51.559292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"outlier_compounds_1, counts_outliers = np.unique(outliers_stack, return_counts=True)\nprint(len(outlier_compounds_1))","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:51.562532Z","iopub.execute_input":"2024-01-07T13:38:51.562856Z","iopub.status.idle":"2024-01-07T13:38:51.569659Z","shell.execute_reply.started":"2024-01-07T13:38:51.562829Z","shell.execute_reply":"2024-01-07T13:38:51.568204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Identification of the outlier compounds using isolation forest analysis\nThis method uses isolation forest based on the cell count ratios to pinpoint outliers. See this tutorial for more information [https://scikit-learn.org/stable/auto_examples/ensemble/plot_isolation_forest.html#sphx-glr-auto-examples-ensemble-plot-isolation-forest-py].","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import IsolationForest\nnormalized_ratios = normalized_ratios.fillna(0)\nclf = IsolationForest(max_samples=100, random_state=0)\nclf.fit(normalized_ratios.values)\noutlier_compounds_2 = normalized_ratios.index[clf.predict(normalized_ratios.values)==-1]\noutlier_compounds_2","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:51.571105Z","iopub.execute_input":"2024-01-07T13:38:51.571498Z","iopub.status.idle":"2024-01-07T13:38:51.813540Z","shell.execute_reply.started":"2024-01-07T13:38:51.571467Z","shell.execute_reply":"2024-01-07T13:38:51.812092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Outlier compounds identified mutually by the two approaches\nThe combined results from the two methods are used to determine the outlier compounds. It's important to note that our analysis is confined to training cell types due to the limited number of samples available for test cell types.","metadata":{}},{"cell_type":"code","source":"toGoCompound_1 = np.intersect1d(outlier_compounds_1, outlier_compounds_2)\nprint(f'# of compounds: {len(outlier_compounds_1)}, {len(outlier_compounds_2)}, {len(toGoCompound_1)}')","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:51.815258Z","iopub.execute_input":"2024-01-07T13:38:51.815638Z","iopub.status.idle":"2024-01-07T13:38:51.822787Z","shell.execute_reply.started":"2024-01-07T13:38:51.815608Z","shell.execute_reply":"2024-01-07T13:38:51.821141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We take the intersection of both methods which results in 15 compounds.","metadata":{}},{"cell_type":"markdown","source":"### Samples associated with the outlier compounds \nAll samples associated with the outlier compounds are identified here which will be removed from the data later.","metadata":{}},{"cell_type":"code","source":"toGoSamples_1 = cell_count.index.get_level_values(1).isin(toGoCompound_1)\nprint(f'number of samples to go due to invalid cell type {toGoSamples_1.sum()}')","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:51.824051Z","iopub.execute_input":"2024-01-07T13:38:51.824406Z","iopub.status.idle":"2024-01-07T13:38:51.835176Z","shell.execute_reply.started":"2024-01-07T13:38:51.824378Z","shell.execute_reply":"2024-01-07T13:38:51.834006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inadequate cell count \nIn our earlier research (refer to [1]), we demonstrated that samples with low cell counts tend to exhibit abnormal DE values. We showed that those samples with 10 or less cells were noisy and potentially errorneous. Here, we identify those samples and associated compounds as outliers.","metadata":{}},{"cell_type":"code","source":"lower_bound = 10\nfig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4, 3))\nax.hist(np.log10(cell_count.values), bins=40, label=cell_type, color='skyblue', edgecolor='black')\nplt.xlabel('cell count (log10)')\nplt.ylabel('samples')\nplt.title('Distribution of cell count across all samples')\nax.axvline(x=np.log10(10), color='green', linestyle='--', linewidth=1)\nouts = cell_count[cell_count<lower_bound]\n","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:51.836888Z","iopub.execute_input":"2024-01-07T13:38:51.837304Z","iopub.status.idle":"2024-01-07T13:38:52.097328Z","shell.execute_reply.started":"2024-01-07T13:38:51.837266Z","shell.execute_reply":"2024-01-07T13:38:52.096240Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Samples associated with low cell count","metadata":{}},{"cell_type":"code","source":"toGoSamples_2_1 = (cell_count<lower_bound) \nprint(f'Number of samples to go: {toGoSamples_2_1.sum()}')","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.098765Z","iopub.execute_input":"2024-01-07T13:38:52.099147Z","iopub.status.idle":"2024-01-07T13:38:52.105531Z","shell.execute_reply.started":"2024-01-07T13:38:52.099113Z","shell.execute_reply":"2024-01-07T13:38:52.104567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Compounds as outliers due to low cell types left\nIf we remove those problamatic samples with insufficient cell counts, let see how many cell types are left for each drug.","metadata":{}},{"cell_type":"code","source":"problem_drugs = outs.index.get_level_values(1).unique()\nproblem_cell_types = outs.groupby(level=1).size()\nall_cell_types = cell_count.groupby(level=1).size()\nmerged_df = all_cell_types.to_frame().merge(problem_cell_types.to_frame(), left_index=True, right_index=True, how='inner')\nmerged_df = merged_df.rename(columns={'0_x':'initial', '0_y':'to_remove'})\nmerged_df['remaining'] = merged_df.initial - merged_df.to_remove\nmerged_df","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.110084Z","iopub.execute_input":"2024-01-07T13:38:52.110382Z","iopub.status.idle":"2024-01-07T13:38:52.134923Z","shell.execute_reply.started":"2024-01-07T13:38:52.110355Z","shell.execute_reply":"2024-01-07T13:38:52.133500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compounds with fewer than 3 cell types are irrelevant for this competition, which focuses on generalizing from 4 training cell types to 2 test cell types. Therefore, we will exclude these compounds and their samples.","metadata":{}},{"cell_type":"code","source":"toGoCompound_2 = merged_df[merged_df['remaining']<3].index.unique()\nprint(f'compounds to go with all their samples: {toGoCompound_2.shape}')","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.136473Z","iopub.execute_input":"2024-01-07T13:38:52.136815Z","iopub.status.idle":"2024-01-07T13:38:52.144799Z","shell.execute_reply.started":"2024-01-07T13:38:52.136790Z","shell.execute_reply":"2024-01-07T13:38:52.143585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"toGoSamples_2_2 = cell_count.index.get_level_values(level=1).isin(toGoCompound_2)\nprint(f'Number of samples to go: {toGoSamples_2_2.sum()}')","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.146244Z","iopub.execute_input":"2024-01-07T13:38:52.146821Z","iopub.status.idle":"2024-01-07T13:38:52.157910Z","shell.execute_reply.started":"2024-01-07T13:38:52.146784Z","shell.execute_reply":"2024-01-07T13:38:52.156840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### All samples to go due to low cell count ","metadata":{"execution":{"iopub.status.busy":"2024-01-07T10:30:11.784841Z","iopub.execute_input":"2024-01-07T10:30:11.785306Z","iopub.status.idle":"2024-01-07T10:30:11.790484Z","shell.execute_reply.started":"2024-01-07T10:30:11.785272Z","shell.execute_reply":"2024-01-07T10:30:11.789289Z"}}},{"cell_type":"code","source":"toGoSamples_2 = (toGoSamples_2_1 + toGoSamples_2_2)\nprint(f'Number of samples to go: {toGoSamples_2.sum()}')","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.159176Z","iopub.execute_input":"2024-01-07T13:38:52.159514Z","iopub.status.idle":"2024-01-07T13:38:52.172114Z","shell.execute_reply.started":"2024-01-07T13:38:52.159486Z","shell.execute_reply":"2024-01-07T13:38:52.170571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Shared outlier compounds from cell type misclassification and low cell count \nLets put together the identified compounds from cell type misclassification and low cell count and see how many are shared.\n","metadata":{}},{"cell_type":"code","source":"all_outlier_compounds = np.union1d(toGoCompound_1, toGoCompound_2)\nshared_outlier_compounds = np.intersect1d(toGoCompound_1, toGoCompound_2)\n\nprint(f'invalid compounds: {len(toGoCompound_1)=}, {len(toGoCompound_2)=}', '=> intersection: ',shared_outlier_compounds)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.173226Z","iopub.execute_input":"2024-01-07T13:38:52.173516Z","iopub.status.idle":"2024-01-07T13:38:52.182721Z","shell.execute_reply.started":"2024-01-07T13:38:52.173492Z","shell.execute_reply":"2024-01-07T13:38:52.181912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shared_outlier_compounds","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.183998Z","iopub.execute_input":"2024-01-07T13:38:52.184352Z","iopub.status.idle":"2024-01-07T13:38:52.199444Z","shell.execute_reply.started":"2024-01-07T13:38:52.184308Z","shell.execute_reply":"2024-01-07T13:38:52.198398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that these 5 compounds that resulted in low cell count also causes cell type misclassification.","metadata":{}},{"cell_type":"markdown","source":"## Compare the identified outlier compounds in this study to Ambrosm's\nTo have a comparitive analysis, we pull out the compounds which were identified by Ambrosm's work, available at [https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/458661]. He also provided the reason why these compounds are outliers.","metadata":{}},{"cell_type":"code","source":"ambrosm = {\n    \"AT13387\": \"only 7 T regulatory cells\",\n    \"Alvocidib\": \"≤10 for several cell types\",\n    \"BAY 61-3606\": \"mean t-score for CD8+ cells > 2\",\n    \"BMS-387032\": \"only 10 T cells CD8+, no Myeloid cells\",\n    \"Belinostat\": \"control compound with too many cells\",\n    \"CEP-18770 (Delanzomib)\": \"≤10 for several cell types\",\n    \"CGM-097\": \"too many T cells CD8+\",\n    \"CGP 60474\": \"≤10 for several cell types\",\n    \"Dabrafenib\": \"control compound with too many cells\",\n    \"Ganetespib (STA-9090)\": \"only 4 T regulatory cells, too many NK cells\",\n    \"I-BET151\": \"too many T cells CD8+\",\n    \"IN1451\": \"≤10 for several cell types\",\n    \"LY2090314\": \"only 6 T cells CD8+\",\n    \"MLN 2238\": \"≤10 for several cell types\",\n    \"Oprozomib (ONX 0912)\": \"≤10 for several cell types\",\n    \"Proscillaridin A;Proscillaridin-A\": \"≤10 for several cell types\",\n    \"Resminostat\": \"no T cells CD8+\",\n    \"Scriptaid\": \"only 2 T regulatory cells\",\n    \"UNII-BXU45ZH6LI\": \"only 6 T cells CD8+\",\n    \"Vorinostat\": \"only 1 T regulatory cell\"\n}","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.201240Z","iopub.execute_input":"2024-01-07T13:38:52.201599Z","iopub.status.idle":"2024-01-07T13:38:52.211497Z","shell.execute_reply.started":"2024-01-07T13:38:52.201571Z","shell.execute_reply":"2024-01-07T13:38:52.210372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from matplotlib.patches import Patch\n# the stacked bar code is adopted from Ambrosm\ndef plot_stacked_bar_chart(cell_types_in_drops, title, xticks=None, xticklabels=None, colors=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    We dont show the controls as they are too tall.\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=(25, 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    first_legend = plt.legend(title='Cell Types')\n    if xticks is not None:\n        plt.xticks(xticks, xticklabels, rotation=90)\n    \n    ax = plt.gca()\n    for ticklabel, color in zip(ax.get_xticklabels(), colors):\n        ticklabel.set_color(color)\n\n    color_legend_handles = [\n        Patch(facecolor='red', label='Shared'),\n        Patch(facecolor='blue', label='Only Ambrosm'),\n        Patch(facecolor='green', label='Only this study'),\n    ]\n    color_legend = ax.legend(handles=color_legend_handles, title='X-tick Label Colors', loc='upper right', bbox_to_anchor=(1, 1))\n    ax.add_artist(first_legend)\n    plt.show()\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')\n\nall_outliers = set(np.concatenate([all_outlier_compounds, list(ambrosm.keys())]))\n\nxticks = np.arange(len(cc))[sorted_compound_names.isin(all_outliers)]\nxticklabels = sorted_compound_names[sorted_compound_names.isin(all_outliers)]\nprint(len(xticklabels))\ncolors = []\nfor i, sm_name in enumerate(xticklabels):\n    if (sm_name in ambrosm) & (sm_name in all_outlier_compounds):\n        color = 'red'\n    elif (sm_name in ambrosm):\n        color = 'blue'\n    elif (sm_name in all_outlier_compounds):\n        color = 'green'\n    else:\n        raise ValueError('shouldnt be')\n    colors.append(color)\n\nplot_stacked_bar_chart(cc[cell_types].values, 'Composition of droplets after 24 hours', xticks=xticks, xticklabels=xticklabels, colors=colors)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:52.212844Z","iopub.execute_input":"2024-01-07T13:38:52.213147Z","iopub.status.idle":"2024-01-07T13:38:54.141378Z","shell.execute_reply.started":"2024-01-07T13:38:52.213123Z","shell.execute_reply":"2024-01-07T13:38:54.138900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"we can observe that most compounds are mutually identified as outliers between our method and Ambrosm's report. Certain compounds are only reported by him and not in this study.","metadata":{}},{"cell_type":"code","source":"# Outliers only in ambrosm\ninter_ = np.intersect1d(all_outlier_compounds, list(ambrosm.keys()))\nprint(f'from {len(all_outlier_compounds)}, there is {len(inter_)} in ambrosh, which has {len(ambrosm)}')\n\n# inter_dict = {key:value if (key in inter_) for key, value in ambrosm.items()}\nambrosm_remaining = set(ambrosm.keys()) - set(all_outlier_compounds)\n\nprint('\\n Those only in ambrosm:')\nfor key in ambrosm_remaining:\n    print(f'{key}: {ambrosm[key]}')\n# Outliers only in this study   \ndiff_outliers = list(set(all_outlier_compounds)-set(ambrosm))\nprint(diff_outliers)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:54.143728Z","iopub.execute_input":"2024-01-07T13:38:54.144114Z","iopub.status.idle":"2024-01-07T13:38:54.153781Z","shell.execute_reply.started":"2024-01-07T13:38:54.144079Z","shell.execute_reply":"2024-01-07T13:38:54.151551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our methodology differs in several aspects:\n\n- We do not exclude compounds due to an excess of cells or a high t-score (in case of BAY 61-3606, Belinostat, Dabrafenib)\n- We do not exclude compounds if there are more than 2 cell types left after removing those that have low cell count (in case of Scriptaid)","metadata":{}},{"cell_type":"markdown","source":"# Evaluation of genes without RNA count\nAs identified by Ambrosm, certain genes are not expressed as RNA count but are given in DE values. Here, we further evaluate them to understand the amount of information missing.","metadata":{}},{"cell_type":"markdown","source":"## Determine whether RNA count is present per gene","metadata":{}},{"cell_type":"code","source":"valid_genes_list = []\nfor row in t_train.iterrows():\n    ct = row[0][0]\n    sm = row[0][1]\n    this_t_score = row[1]\n    # Distinguish genes which are expressed / not expressed for this compound\n    compound_zero = (rna_count.loc[ct, sm] == 0).values.ravel()\n    valid_genes_list.append(~compound_zero)\nvalid_genes_df = pd.DataFrame(valid_genes_list, columns=de_train.columns[5:])\nvalid_genes_df.head()","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:54.155878Z","iopub.execute_input":"2024-01-07T13:38:54.156415Z","iopub.status.idle":"2024-01-07T13:38:57.267729Z","shell.execute_reply.started":"2024-01-07T13:38:54.156381Z","shell.execute_reply":"2024-01-07T13:38:57.266478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## How bad is the missingness?\nWe evaluate the distribution of invalid DE values. We show the number of missing DE values (genes) for samples as well as the count of samples missing for genes. ","metadata":{}},{"cell_type":"code","source":"# invalid genes for samples\nmissing_numbers_1 = (valid_genes_df==False).sum(axis=1).values\nfig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))\nax=axes[0]\nax.hist(missing_numbers_1, bins=40, label=cell_type, color='skyblue', edgecolor='black')\nax.set_xlabel('Number of invalid genes')\nax.set_ylabel('Sample')\n# ax.set_title('Distribution of number of invalid genes across different samples')\n\n# invalid genes stats across samples \nmissing_numbers_2 = (valid_genes_df==False).sum(axis=0).values\nax=axes[1]\nax.hist(missing_numbers_2, bins=40, label=cell_type, color='skyblue', edgecolor='black')\nplt.xlabel('Number of samples')\nplt.ylabel('Number of invalid genes')\n# plt.title('Distribution of number of genes across different samples')\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:57.269226Z","iopub.execute_input":"2024-01-07T13:38:57.269627Z","iopub.status.idle":"2024-01-07T13:38:57.790296Z","shell.execute_reply.started":"2024-01-07T13:38:57.269594Z","shell.execute_reply":"2024-01-07T13:38:57.788953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('number of genes available across all samples: ', valid_genes_df.all().sum())\nprint('maximum number of samples missing for a gene: ', missing_numbers_2.max())\nprint('maximum fractio of missing genes for a sample: ', missing_numbers_1.max()/len(gene_names))","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:57.791786Z","iopub.execute_input":"2024-01-07T13:38:57.792086Z","iopub.status.idle":"2024-01-07T13:38:57.804796Z","shell.execute_reply.started":"2024-01-07T13:38:57.792060Z","shell.execute_reply":"2024-01-07T13:38:57.803622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Only 184 genes are consistently present across all samples. Some genes are missing in 611 out of 612 samples, and certain samples are missing up to 92% of genes. Currently, we are not excluding any samples due to this missingness, but additional filtering will be necessary for refinement.","metadata":{}},{"cell_type":"markdown","source":"# Portion of missing data due to outlier samples and missing genes\nBy removing those samples identified as outliers as well as removing those invalid DE values, lets see what portion of the original data remains.","metadata":{}},{"cell_type":"code","source":"toGoSamples = (toGoSamples_1 + toGoSamples_2) \nfilter_samples_df = pd.DataFrame(np.full(data.shape, True), columns=gene_names)\nfilter_samples_df.loc[toGoSamples.values,:] = False\nfilter_genes_df = valid_genes_df\nfilter_df = filter_genes_df*filter_samples_df\nprint('portion of data remaining: ', filter_df.sum().sum()/(data.shape[0]*data.shape[1]))","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:57.806701Z","iopub.execute_input":"2024-01-07T13:38:57.807226Z","iopub.status.idle":"2024-01-07T13:38:57.855135Z","shell.execute_reply.started":"2024-01-07T13:38:57.807184Z","shell.execute_reply":"2024-01-07T13:38:57.853462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that 65 percent of data has to go due to potential problems in the reported DE values. ","metadata":{}},{"cell_type":"markdown","source":"## Importance of the removed data\nFollowing the approach used in [https://www.kaggle.com/code/alexandervc/op2-eda-housekeeping-genes?scriptVersionId=155511517&cellId=11], we present a visual comparison of the significance of data that needs to go versus remaining ones. We distinguish housekeeping genes versus other genes for better comparision.","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('https://www.tau.ac.il/~elieis/HKG/HK_genes.txt', sep=' ', header=None)\nhousekeeping_genes = df.loc[:, 0]  # Gene names\nprint(f'Number of housekeeping genes: {len(housekeeping_genes)}')\nhousekeeping_genes_filter = gene_names.isin(housekeeping_genes)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:38:57.857115Z","iopub.execute_input":"2024-01-07T13:38:57.857901Z","iopub.status.idle":"2024-01-07T13:38:59.325570Z","shell.execute_reply.started":"2024-01-07T13:38:57.857856Z","shell.execute_reply":"2024-01-07T13:38:59.323700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"d = pd.DataFrame(); IX = -1\nlist_thresholds = [0.1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16]\n\nfor t in list_thresholds:\n    # remaining all genes\n    v  = data_df.loc[:, ~housekeeping_genes_filter].values.ravel()\n    v = v[~np.isnan(v)]\n    v = 10**(-np.abs(v))\n    m = v < t\n    IX+=1\n    d.loc[IX, 'Threshold'] = t\n    d.loc[IX, 'Other genes'] = 100*m.sum()/len(v)\n\n    v  = data_df.loc[:, housekeeping_genes_filter].values.ravel()\n    v = v[~np.isnan(v)]\n    v = 10**(-np.abs(v))\n    m = v < t\n    d.loc[IX, 'House keeping genes'] = 100*m.sum()/len(v)\n\n    v  = data_df.loc[:, ~housekeeping_genes_filter][filter_df].values.ravel()\n    v = v[~np.isnan(v)]\n    v = 10**(-np.abs(v))\n    m = v < t\n    d.loc[IX, 'Other genes after filter'] = 100*m.sum()/len(v)\n\n    v  = data_df.loc[:, housekeeping_genes_filter][filter_df].values.ravel()\n    v = v[~np.isnan(v)]\n    v = 10**(-np.abs(v))\n    m = v < t\n    d.loc[IX, 'House keeping genes after filter'] = 100*m.sum()/len(v)","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:42:08.104788Z","iopub.execute_input":"2024-01-07T13:42:08.105233Z","iopub.status.idle":"2024-01-07T13:42:08.110237Z","shell.execute_reply.started":"2024-01-07T13:42:08.105195Z","shell.execute_reply":"2024-01-07T13:42:08.109379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\n# Define line styles and line width\nline_styles = ['-', '--', '-.', ':']\nline_width = 2  # Increase the line thickness\n\nplt.figure(figsize=(15, 4))\n\n# Iterate over columns and line styles\nfor col, style in zip(d.columns[1:], line_styles * (len(d.columns) // len(line_styles) + 1)):\n    plt.plot(d[col].values, label=col, linestyle=style, linewidth=line_width)\n\nplt.grid()\nplt.legend()\nplt.xticks(range(len(list_thresholds)), list_thresholds)\nplt.xlabel('Threshold on p-value')\nplt.ylabel('Percentage of significantly expressed genes')\nplt.savefig('sig.png', bbox_inches='tight')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-01-07T13:44:58.385708Z","iopub.execute_input":"2024-01-07T13:44:58.386042Z","iopub.status.idle":"2024-01-07T13:44:58.941042Z","shell.execute_reply.started":"2024-01-07T13:44:58.386019Z","shell.execute_reply":"2024-01-07T13:44:58.939633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We observe that by removing the problamatic measuerments, the overal significane of the remaining samples sharply decreases, even below housekeeping genes. ","metadata":{}}]}