{"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":"# Differential expression between annotated cell types and between leiden clusters in CITE-seq input data (train and test)","metadata":{}},{"cell_type":"markdown","source":"# Imports and set up","metadata":{}},{"cell_type":"code","source":"%%capture\n!pip3 install scanpy[leiden]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-12-14T21:10:50.345118Z","iopub.execute_input":"2022-12-14T21:10:50.346314Z","iopub.status.idle":"2022-12-14T21:11:03.179040Z","shell.execute_reply.started":"2022-12-14T21:10:50.346207Z","shell.execute_reply":"2022-12-14T21:11:03.177362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n!pip3 install apybiomart","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:03.185964Z","iopub.execute_input":"2022-12-14T21:11:03.186496Z","iopub.status.idle":"2022-12-14T21:11:14.681452Z","shell.execute_reply.started":"2022-12-14T21:11:03.186427Z","shell.execute_reply":"2022-12-14T21:11:14.679835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport sys\nimport logging","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:14.683243Z","iopub.execute_input":"2022-12-14T21:11:14.684006Z","iopub.status.idle":"2022-12-14T21:11:14.689144Z","shell.execute_reply.started":"2022-12-14T21:11:14.683965Z","shell.execute_reply":"2022-12-14T21:11:14.687996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR_RAW = \"/kaggle/input/msci-raw-counts-anndata\"\n# Cite\nFP_CITE_RAW_TRAIN_INPUTS = os.path.join(DATA_DIR_RAW,\"train_cite_inputs_raw.ann.h5\")\nFP_CITE_RAW_TRAIN_TARGETS = os.path.join(DATA_DIR_RAW,\"train_cite_targets_raw.ann.h5\")\nFP_CITE_RAW_TEST_INPUTS = os.path.join(DATA_DIR_RAW,\"test_cite_inputs_raw.ann.h5\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:14.692414Z","iopub.execute_input":"2022-12-14T21:11:14.692891Z","iopub.status.idle":"2022-12-14T21:11:14.702735Z","shell.execute_reply.started":"2022-12-14T21:11:14.692848Z","shell.execute_reply":"2022-12-14T21:11:14.701751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DONORS = [32606]\nDAYS = [2, 3, 4, 7]\n\nDAYS, DONORS = tuple(map(str, DAYS)), tuple(map(str, DONORS))","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:14.704006Z","iopub.execute_input":"2022-12-14T21:11:14.704361Z","iopub.status.idle":"2022-12-14T21:11:14.715156Z","shell.execute_reply.started":"2022-12-14T21:11:14.704331Z","shell.execute_reply":"2022-12-14T21:11:14.713764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\n\nimport anndata as ad\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport scanpy as sc\n\nsc.settings.verbosity = 3\nsc.set_figure_params(dpi=150)\nsns.set_style(\"ticks\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:14.716680Z","iopub.execute_input":"2022-12-14T21:11:14.717114Z","iopub.status.idle":"2022-12-14T21:11:17.573438Z","shell.execute_reply.started":"2022-12-14T21:11:14.717078Z","shell.execute_reply":"2022-12-14T21:11:17.572364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import apybiomart","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:17.574911Z","iopub.execute_input":"2022-12-14T21:11:17.575322Z","iopub.status.idle":"2022-12-14T21:11:17.688238Z","shell.execute_reply.started":"2022-12-14T21:11:17.575286Z","shell.execute_reply":"2022-12-14T21:11:17.687081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"metadata_fp = \"/kaggle/input/open-problems-multimodal/metadata.csv\"\nadd_metadata_fp = \"/kaggle/input/open-problems-multimodal/metadata_cite_day_2_donor_27678.csv\"\nmetadata = pd.read_csv(metadata_fp)\nadd_metadata = pd.read_csv(add_metadata_fp)\nmetadata = pd.concat([metadata, add_metadata])\nmetadata.set_index(\"cell_id\", inplace=True)\nmetadata.head()","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:17.689660Z","iopub.execute_input":"2022-12-14T21:11:17.690403Z","iopub.status.idle":"2022-12-14T21:11:18.177603Z","shell.execute_reply.started":"2022-12-14T21:11:17.690362Z","shell.execute_reply":"2022-12-14T21:11:18.176450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def annotate(adata: ad.AnnData, metadata: pd.DataFrame, f_donors: [str]):\n    adata.obs[\"cell_type\"] = metadata[\"cell_type\"][adata.obs_names]\n    adata.obs[\"day\"] = metadata[\"day\"][adata.obs_names].astype(str)\n    adata.obs[\"donor\"] = metadata[\"donor\"][adata.obs_names].astype(str)\n    adata.obs[\"gender\"] = \"m\"\n    adata.obs.loc[adata.obs[\"donor\"].isin(f_donors), \"gender\"] = \"f\"\n    adata.obs[\"n_genes_counts\"] = (adata.X > 0).sum(axis=1).A1\n    adata.obs[\"read_depth\"] = adata.X.sum(axis=1).A1","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:18.178932Z","iopub.execute_input":"2022-12-14T21:11:18.179326Z","iopub.status.idle":"2022-12-14T21:11:18.337059Z","shell.execute_reply.started":"2022-12-14T21:11:18.179294Z","shell.execute_reply":"2022-12-14T21:11:18.335605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Input raw data","metadata":{}},{"cell_type":"code","source":"cite_raw_train = sc.read_h5ad(FP_CITE_RAW_TRAIN_INPUTS)\ncite_raw_test = sc.read_h5ad(FP_CITE_RAW_TEST_INPUTS)\nadata = ad.concat([cite_raw_train, cite_raw_test], merge=\"same\")\ndel cite_raw_train\ndel cite_raw_test\ngc.collect()\nannotate(adata, metadata, f_donors = [\"13176\", ])\nadata.layers[\"counts\"] = adata.X.copy()","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:11:18.338863Z","iopub.execute_input":"2022-12-14T21:11:18.339756Z","iopub.status.idle":"2022-12-14T21:12:13.454511Z","shell.execute_reply.started":"2022-12-14T21:11:18.339704Z","shell.execute_reply":"2022-12-14T21:12:13.453321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Explore the data","metadata":{}},{"cell_type":"code","source":"adata.obs","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:12:13.456104Z","iopub.execute_input":"2022-12-14T21:12:13.456597Z","iopub.status.idle":"2022-12-14T21:12:13.480956Z","shell.execute_reply.started":"2022-12-14T21:12:13.456549Z","shell.execute_reply":"2022-12-14T21:12:13.479839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.var","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:12:13.483033Z","iopub.execute_input":"2022-12-14T21:12:13.484333Z","iopub.status.idle":"2022-12-14T21:12:13.495389Z","shell.execute_reply.started":"2022-12-14T21:12:13.484292Z","shell.execute_reply":"2022-12-14T21:12:13.494170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Annotate the genes","metadata":{}},{"cell_type":"markdown","source":"### Explore biomart attributes","metadata":{}},{"cell_type":"code","source":"biomart_datasets = apybiomart.find_datasets()\nbiomart_datasets","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:12:13.500129Z","iopub.execute_input":"2022-12-14T21:12:13.500514Z","iopub.status.idle":"2022-12-14T21:12:15.496697Z","shell.execute_reply.started":"2022-12-14T21:12:13.500485Z","shell.execute_reply":"2022-12-14T21:12:15.495897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_datasets[biomart_datasets['Dataset_name'].str.lower().str.contains('human')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:12:15.497763Z","iopub.execute_input":"2022-12-14T21:12:15.498282Z","iopub.status.idle":"2022-12-14T21:12:15.511388Z","shell.execute_reply.started":"2022-12-14T21:12:15.498247Z","shell.execute_reply":"2022-12-14T21:12:15.510274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs = apybiomart.find_attributes('hsapiens_gene_ensembl')\nbiomart_attrs","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:12:15.512953Z","iopub.execute_input":"2022-12-14T21:12:15.513296Z","iopub.status.idle":"2022-12-14T21:13:22.291765Z","shell.execute_reply.started":"2022-12-14T21:12:15.513266Z","shell.execute_reply":"2022-12-14T21:13:22.290402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs[biomart_attrs['Attribute_name'].str.lower().str.contains('gene name')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.293610Z","iopub.execute_input":"2022-12-14T21:13:22.294087Z","iopub.status.idle":"2022-12-14T21:13:22.316785Z","shell.execute_reply.started":"2022-12-14T21:13:22.294038Z","shell.execute_reply":"2022-12-14T21:13:22.315476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs[biomart_attrs['Attribute_name'].str.lower().str.contains('uniprotkb')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.318208Z","iopub.execute_input":"2022-12-14T21:13:22.318629Z","iopub.status.idle":"2022-12-14T21:13:22.338795Z","shell.execute_reply.started":"2022-12-14T21:13:22.318594Z","shell.execute_reply":"2022-12-14T21:13:22.337519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs[biomart_attrs['Attribute_name'].str.lower().str.contains('hgnc')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.340211Z","iopub.execute_input":"2022-12-14T21:13:22.340552Z","iopub.status.idle":"2022-12-14T21:13:22.361096Z","shell.execute_reply.started":"2022-12-14T21:13:22.340522Z","shell.execute_reply":"2022-12-14T21:13:22.359577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs[biomart_attrs['Attribute_name'].str.lower().str.contains('atlas')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.365316Z","iopub.execute_input":"2022-12-14T21:13:22.365708Z","iopub.status.idle":"2022-12-14T21:13:22.383678Z","shell.execute_reply.started":"2022-12-14T21:13:22.365674Z","shell.execute_reply":"2022-12-14T21:13:22.382297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs[biomart_attrs['Attribute_name'].str.lower().str.contains('description')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.385299Z","iopub.execute_input":"2022-12-14T21:13:22.385703Z","iopub.status.idle":"2022-12-14T21:13:22.405025Z","shell.execute_reply.started":"2022-12-14T21:13:22.385661Z","shell.execute_reply":"2022-12-14T21:13:22.403858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"biomart_attrs[biomart_attrs['Attribute_name'].str.lower().str.contains('gene description')]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.406803Z","iopub.execute_input":"2022-12-14T21:13:22.407247Z","iopub.status.idle":"2022-12-14T21:13:22.429346Z","shell.execute_reply.started":"2022-12-14T21:13:22.407200Z","shell.execute_reply":"2022-12-14T21:13:22.428075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Query biomart","metadata":{}},{"cell_type":"code","source":"attrs_to_get = [\n    'ensembl_gene_id',\n    'ensembl_gene_id_version',\n    'external_gene_name',\n    'description',\n    'phenotype_description',\n    'uniprot_gn_symbol',\n    'uniprot_gn_id',\n    'hgnc_symbol',\n    'hgnc_id',\n    'arrayexpress',\n    'hpa_accession',\n    'hpa_id'\n]\n\n# Batches are needed due to \"too many attributes selected\" biomart exception\nattrs_to_get_batches = [attrs_to_get[i:i + 3] for i in range(0, len(attrs_to_get), 3)]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.430967Z","iopub.execute_input":"2022-12-14T21:13:22.431815Z","iopub.status.idle":"2022-12-14T21:13:22.439257Z","shell.execute_reply.started":"2022-12-14T21:13:22.431775Z","shell.execute_reply":"2022-12-14T21:13:22.438105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"logger = logging.getLogger('biomart')\nlogger.setLevel(logging.INFO)\n\n# to log into console\nch = logging.StreamHandler()\nlogger.addHandler(ch)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.440988Z","iopub.execute_input":"2022-12-14T21:13:22.441472Z","iopub.status.idle":"2022-12-14T21:13:22.450348Z","shell.execute_reply.started":"2022-12-14T21:13:22.441425Z","shell.execute_reply":"2022-12-14T21:13:22.449208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df = None\n\nfor attrs_batch in attrs_to_get_batches:\n    \n    batch_df = apybiomart.query(\n        attributes=attrs_batch,\n        dataset='hsapiens_gene_ensembl',\n        filters=None\n    )\n    \n    if gene_ids_df is None:\n        gene_ids_df = batch_df\n    else:\n        gene_ids_df = pd.concat([gene_ids_df, batch_df], axis=1)\n        \n    logger.info(f'Queried:{*attrs_batch,}')","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:13:22.452414Z","iopub.execute_input":"2022-12-14T21:13:22.452789Z","iopub.status.idle":"2022-12-14T21:16:38.711500Z","shell.execute_reply.started":"2022-12-14T21:13:22.452752Z","shell.execute_reply":"2022-12-14T21:16:38.709963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.713923Z","iopub.execute_input":"2022-12-14T21:16:38.714446Z","iopub.status.idle":"2022-12-14T21:16:38.742909Z","shell.execute_reply.started":"2022-12-14T21:16:38.714395Z","shell.execute_reply":"2022-12-14T21:16:38.741655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df['Gene stable ID'].isnull().value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.744761Z","iopub.execute_input":"2022-12-14T21:16:38.745167Z","iopub.status.idle":"2022-12-14T21:16:38.764993Z","shell.execute_reply.started":"2022-12-14T21:16:38.745130Z","shell.execute_reply":"2022-12-14T21:16:38.763524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df.dropna(subset=['Gene stable ID'], inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.766953Z","iopub.execute_input":"2022-12-14T21:16:38.768114Z","iopub.status.idle":"2022-12-14T21:16:38.871043Z","shell.execute_reply.started":"2022-12-14T21:16:38.768062Z","shell.execute_reply":"2022-12-14T21:16:38.869916Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df['Gene name'].isnull().value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.872516Z","iopub.execute_input":"2022-12-14T21:16:38.872856Z","iopub.status.idle":"2022-12-14T21:16:38.887447Z","shell.execute_reply.started":"2022-12-14T21:16:38.872827Z","shell.execute_reply":"2022-12-14T21:16:38.886178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df['Gene stable ID'] + '_' + gene_ids_df['Gene name']","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.888773Z","iopub.execute_input":"2022-12-14T21:16:38.889163Z","iopub.status.idle":"2022-12-14T21:16:38.912420Z","shell.execute_reply.started":"2022-12-14T21:16:38.889129Z","shell.execute_reply":"2022-12-14T21:16:38.911363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df['id_name'] = gene_ids_df['Gene stable ID'] + '_' + gene_ids_df['Gene name']","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.913700Z","iopub.execute_input":"2022-12-14T21:16:38.914072Z","iopub.status.idle":"2022-12-14T21:16:38.936031Z","shell.execute_reply.started":"2022-12-14T21:16:38.914040Z","shell.execute_reply":"2022-12-14T21:16:38.934636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df.set_index('id_name', inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.937518Z","iopub.execute_input":"2022-12-14T21:16:38.937900Z","iopub.status.idle":"2022-12-14T21:16:38.944924Z","shell.execute_reply.started":"2022-12-14T21:16:38.937848Z","shell.execute_reply":"2022-12-14T21:16:38.943604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_ids_df","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.946229Z","iopub.execute_input":"2022-12-14T21:16:38.946555Z","iopub.status.idle":"2022-12-14T21:16:38.971569Z","shell.execute_reply.started":"2022-12-14T21:16:38.946527Z","shell.execute_reply":"2022-12-14T21:16:38.970355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Add information to the adata","metadata":{}},{"cell_type":"code","source":"cols_to_join = ['Gene name', 'Gene description', 'Phenotype description']","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.973091Z","iopub.execute_input":"2022-12-14T21:16:38.973465Z","iopub.status.idle":"2022-12-14T21:16:38.979247Z","shell.execute_reply.started":"2022-12-14T21:16:38.973433Z","shell.execute_reply":"2022-12-14T21:16:38.977675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.var = adata.var.join(gene_ids_df[cols_to_join])","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:38.980315Z","iopub.execute_input":"2022-12-14T21:16:38.980670Z","iopub.status.idle":"2022-12-14T21:16:39.027258Z","shell.execute_reply.started":"2022-12-14T21:16:38.980639Z","shell.execute_reply":"2022-12-14T21:16:39.025859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Basic vizualizations","metadata":{}},{"cell_type":"code","source":"ax = sc.pl.highest_expr_genes(adata, n_top=10, gene_symbols='Gene name', show=False)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:39.029068Z","iopub.execute_input":"2022-12-14T21:16:39.029574Z","iopub.status.idle":"2022-12-14T21:16:51.686036Z","shell.execute_reply.started":"2022-12-14T21:16:39.029520Z","shell.execute_reply":"2022-12-14T21:16:51.685045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_10_names = [ticklabel._text for ticklabel in ax.get_ymajorticklabels()]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:51.687204Z","iopub.execute_input":"2022-12-14T21:16:51.687714Z","iopub.status.idle":"2022-12-14T21:16:51.693474Z","shell.execute_reply.started":"2022-12-14T21:16:51.687682Z","shell.execute_reply":"2022-12-14T21:16:51.692310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_10_genes = gene_ids_df[gene_ids_df['Gene name'].isin(top_10_names)]\ntop_10_genes","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:51.695912Z","iopub.execute_input":"2022-12-14T21:16:51.696681Z","iopub.status.idle":"2022-12-14T21:16:51.729493Z","shell.execute_reply.started":"2022-12-14T21:16:51.696632Z","shell.execute_reply":"2022-12-14T21:16:51.728456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, gene in top_10_genes.iterrows():\n    print(gene['Gene name'], gene['Gene description'], gene['Phenotype description'], '---', sep='\\n\\n')","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:51.730991Z","iopub.execute_input":"2022-12-14T21:16:51.731371Z","iopub.status.idle":"2022-12-14T21:16:51.739426Z","shell.execute_reply.started":"2022-12-14T21:16:51.731338Z","shell.execute_reply":"2022-12-14T21:16:51.738147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Surprise: these are the houskeeping genes :D**","metadata":{}},{"cell_type":"markdown","source":"# Filter by donors and days","metadata":{}},{"cell_type":"code","source":"adata = adata[adata.obs[adata.obs.donor.isin(DONORS) & adata.obs.day.isin(DAYS)].index, :]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:51.740995Z","iopub.execute_input":"2022-12-14T21:16:51.741330Z","iopub.status.idle":"2022-12-14T21:16:51.786331Z","shell.execute_reply.started":"2022-12-14T21:16:51.741300Z","shell.execute_reply":"2022-12-14T21:16:51.785078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Get HVG and normalise variance with analytical pearson residuals, count PCA","metadata":{}},{"cell_type":"code","source":"sc.experimental.pp.highly_variable_genes(adata, n_top_genes=3000, batch_key=\"donor\")\nsc.experimental.pp.normalize_pearson_residuals_pca(adata)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:16:51.794501Z","iopub.execute_input":"2022-12-14T21:16:51.794912Z","iopub.status.idle":"2022-12-14T21:17:39.512978Z","shell.execute_reply.started":"2022-12-14T21:16:51.794874Z","shell.execute_reply":"2022-12-14T21:17:39.511742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot PCA on analytical pearson residuals","metadata":{}},{"cell_type":"code","source":"sc.pl.pca(adata, color=[\"cell_type\", \"donor\", \"day\", \"n_genes_counts\"])","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:17:39.514595Z","iopub.execute_input":"2022-12-14T21:17:39.515354Z","iopub.status.idle":"2022-12-14T21:17:41.826576Z","shell.execute_reply.started":"2022-12-14T21:17:39.515312Z","shell.execute_reply":"2022-12-14T21:17:41.825141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot UMAP","metadata":{}},{"cell_type":"code","source":"sc.pp.neighbors(adata, n_pcs=15, n_neighbors=20)\nsc.tl.umap(adata, min_dist=0.2)\nsc.pl.umap(adata, color=[\"cell_type\", \"donor\", \"day\"], frameon=False)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:17:41.828261Z","iopub.execute_input":"2022-12-14T21:17:41.828730Z","iopub.status.idle":"2022-12-14T21:18:57.827308Z","shell.execute_reply.started":"2022-12-14T21:17:41.828685Z","shell.execute_reply":"2022-12-14T21:18:57.825894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.obs[\"EryP\"] = (adata.obs.cell_type == \"EryP\").astype(\"str\")\nadata.obs[\"MkP\"] = (adata.obs.cell_type == \"MkP\").astype(\"str\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:18:57.829450Z","iopub.execute_input":"2022-12-14T21:18:57.830262Z","iopub.status.idle":"2022-12-14T21:18:57.869600Z","shell.execute_reply.started":"2022-12-14T21:18:57.830218Z","shell.execute_reply":"2022-12-14T21:18:57.868387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"G1S_genes_Tirosh = ['ENSG00000100297', 'ENSG00000132646', 'ENSG00000176890', 'ENSG00000168496', 'ENSG00000073111', 'ENSG00000104738', 'ENSG00000167325', 'ENSG00000076248', 'ENSG00000131153', 'ENSG00000076003', 'ENSG00000144354', 'ENSG00000143476', 'ENSG00000198056', 'ENSG00000276043', 'ENSG00000151725', 'ENSG00000119969', 'ENSG00000049541', 'ENSG00000117748', 'ENSG00000132780', 'ENSG00000111247', 'ENSG00000112312', 'ENSG00000092470', 'ENSG00000163950', 'ENSG00000175305', 'ENSG00000012963', 'ENSG00000077514', 'ENSG00000095002', 'ENSG00000156802', 'ENSG00000051180', 'ENSG00000171848', 'ENSG00000093009', 'ENSG00000094804', 'ENSG00000174371', 'ENSG00000075131', 'ENSG00000136982', 'ENSG00000197299', 'ENSG00000118412', 'ENSG00000162607', 'ENSG00000092853', 'ENSG00000101868', 'ENSG00000159259', 'ENSG00000136492', 'ENSG00000129173']\nG2M_genes_Tirosh = ['ENSG00000164104', 'ENSG00000170312', 'ENSG00000137804', 'ENSG00000175063', 'ENSG00000089685', 'ENSG00000088325', 'ENSG00000131747', 'ENSG00000080986', 'ENSG00000123975', 'ENSG00000143228', 'ENSG00000173207', 'ENSG00000148773', 'ENSG00000120802', 'ENSG00000117724', 'ENSG00000013810', 'ENSG00000129195', 'ENSG00000113810', 'ENSG00000157456', 'ENSG00000169607', 'ENSG00000136108', 'ENSG00000178999', 'ENSG00000169679', 'ENSG00000138160', 'ENSG00000143401', 'ENSG00000188229', 'ENSG00000075218', 'ENSG00000138182', 'ENSG00000123485', 'ENSG00000111665', 'ENSG00000189159', 'ENSG00000117399', 'ENSG00000112742', 'ENSG00000158402', 'ENSG00000142945', 'ENSG00000100401', 'ENSG00000010292', 'ENSG00000126787', 'ENSG00000184661', 'ENSG00000134690', 'ENSG00000114346', 'ENSG00000137807', 'ENSG00000072571', 'ENSG00000087586', 'ENSG00000134222', 'ENSG00000011426', 'ENSG00000143815', 'ENSG00000175216', 'ENSG00000138778', 'ENSG00000102974', 'ENSG00000117650', 'ENSG00000092140', 'ENSG00000139354', 'ENSG00000094916', 'ENSG00000115163']","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:18:57.871394Z","iopub.execute_input":"2022-12-14T21:18:57.872258Z","iopub.status.idle":"2022-12-14T21:18:57.884631Z","shell.execute_reply.started":"2022-12-14T21:18:57.872207Z","shell.execute_reply":"2022-12-14T21:18:57.883118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.var[\"ensembl_ids\"] = [ t.split('_')[0] for t in adata.var_names]\nadata.var[\"gene_names\"] = [ t.split('_')[1] for t in adata.var_names]","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:18:57.886467Z","iopub.execute_input":"2022-12-14T21:18:57.886987Z","iopub.status.idle":"2022-12-14T21:18:57.920817Z","shell.execute_reply.started":"2022-12-14T21:18:57.886920Z","shell.execute_reply":"2022-12-14T21:18:57.919674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adata.obs[\"G1S_genes_Tirosh_mean\"] = adata[:,adata.var[\"ensembl_ids\"].isin(G1S_genes_Tirosh)].X.sum(axis=1).A1 / len(G1S_genes_Tirosh)\nadata.obs[\"G2M_genes_Tirosh_mean\"] = adata[:,adata.var[\"ensembl_ids\"].isin(G2M_genes_Tirosh)].X.sum(axis=1).A1 / len(G2M_genes_Tirosh)\nadata.obs[\"G1S+2M_genes_Tirosh_mean\"] = adata[:,adata.var[\"ensembl_ids\"].isin(G1S_genes_Tirosh + G2M_genes_Tirosh)].X.sum(axis=1).A1 / len(G1S_genes_Tirosh + G2M_genes_Tirosh)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:18:57.922769Z","iopub.execute_input":"2022-12-14T21:18:57.923289Z","iopub.status.idle":"2022-12-14T21:18:59.102836Z","shell.execute_reply.started":"2022-12-14T21:18:57.923242Z","shell.execute_reply":"2022-12-14T21:18:59.101612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.umap(adata, color=[\"cell_type\", \"EryP\", \"MkP\", \"G1S+2M_genes_Tirosh_mean\", \"n_genes_counts\", \"read_depth\", \"day\", \"donor\"], frameon=False, ncols=3)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:18:59.104226Z","iopub.execute_input":"2022-12-14T21:18:59.104578Z","iopub.status.idle":"2022-12-14T21:19:03.129651Z","shell.execute_reply.started":"2022-12-14T21:18:59.104547Z","shell.execute_reply":"2022-12-14T21:19:03.128216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.umap(adata, color=[\"cell_type\", \"G1S+2M_genes_Tirosh_mean\", \"n_genes_counts\", \"read_depth\", \"day\", \"donor\"], frameon=False, ncols=3)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:19:03.131433Z","iopub.execute_input":"2022-12-14T21:19:03.131833Z","iopub.status.idle":"2022-12-14T21:19:06.472510Z","shell.execute_reply.started":"2022-12-14T21:19:03.131794Z","shell.execute_reply":"2022-12-14T21:19:06.470699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Normalise data for differential expression","metadata":{}},{"cell_type":"code","source":"# Для отрисовки и дифференциальной экспрессии\nsc.pp.normalize_total(adata, target_sum=1e4)\nsc.pp.log1p(adata)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:19:06.474513Z","iopub.execute_input":"2022-12-14T21:19:06.475448Z","iopub.status.idle":"2022-12-14T21:19:10.352657Z","shell.execute_reply.started":"2022-12-14T21:19:06.475408Z","shell.execute_reply":"2022-12-14T21:19:10.351378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Differential expression between annotated cell types","metadata":{}},{"cell_type":"markdown","source":"## Wilcoxon","metadata":{}},{"cell_type":"code","source":"sc.tl.rank_genes_groups(adata, groupby=\"cell_type\", method=\"wilcoxon\", key_added=\"rank_genes_groups_cell_type\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:19:10.354110Z","iopub.execute_input":"2022-12-14T21:19:10.354469Z","iopub.status.idle":"2022-12-14T21:20:49.423840Z","shell.execute_reply.started":"2022-12-14T21:19:10.354439Z","shell.execute_reply":"2022-12-14T21:20:49.422580Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.rank_genes_groups(adata, key=\"rank_genes_groups_cell_type\", gene_symbols=\"gene_names\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:20:49.425722Z","iopub.execute_input":"2022-12-14T21:20:49.426183Z","iopub.status.idle":"2022-12-14T21:20:51.450285Z","shell.execute_reply.started":"2022-12-14T21:20:49.426150Z","shell.execute_reply":"2022-12-14T21:20:51.449134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.rank_genes_groups_heatmap(adata, key=\"rank_genes_groups_cell_type\", gene_symbols=\"gene_names\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:20:51.451859Z","iopub.execute_input":"2022-12-14T21:20:51.452383Z","iopub.status.idle":"2022-12-14T21:20:56.253871Z","shell.execute_reply.started":"2022-12-14T21:20:51.452344Z","shell.execute_reply":"2022-12-14T21:20:56.252105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.rank_genes_groups_matrixplot(\n    adata,\n    n_genes=5, key=\"rank_genes_groups_cell_type\",\n    values_to_plot=\"logfoldchanges\",\n    cmap='bwr',\n    vmin=-4,\n    vmax=4,\n    min_logfoldchange=2,\n    colorbar_title='log fold change',\n)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:20:56.255600Z","iopub.execute_input":"2022-12-14T21:20:56.256036Z","iopub.status.idle":"2022-12-14T21:21:04.127262Z","shell.execute_reply.started":"2022-12-14T21:20:56.255985Z","shell.execute_reply":"2022-12-14T21:21:04.126344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for cell_type in adata.obs.cell_type.unique():\n    de_df = sc.get.rank_genes_groups_df(adata, group=cell_type, key=\"rank_genes_groups_cell_type\")\n    de_df.to_csv(f\"wilcoxon_{cell_type}_vs_rest.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:21:04.128674Z","iopub.execute_input":"2022-12-14T21:21:04.129205Z","iopub.status.idle":"2022-12-14T21:21:08.012691Z","shell.execute_reply.started":"2022-12-14T21:21:04.129172Z","shell.execute_reply":"2022-12-14T21:21:08.011416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## logreg","metadata":{}},{"cell_type":"code","source":"sc.tl.rank_genes_groups(adata, groupby=\"cell_type\", method=\"logreg\", key_added=\"rank_genes_groups_logreg_cell_type\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:21:08.014599Z","iopub.execute_input":"2022-12-14T21:21:08.016385Z","iopub.status.idle":"2022-12-14T21:25:27.835094Z","shell.execute_reply.started":"2022-12-14T21:21:08.016328Z","shell.execute_reply":"2022-12-14T21:25:27.833862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.rank_genes_groups(adata, key=\"rank_genes_groups_logreg_cell_type\", gene_symbols=\"gene_names\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:25:27.836778Z","iopub.execute_input":"2022-12-14T21:25:27.837158Z","iopub.status.idle":"2022-12-14T21:25:29.945060Z","shell.execute_reply.started":"2022-12-14T21:25:27.837124Z","shell.execute_reply":"2022-12-14T21:25:29.944111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Leiden clustering","metadata":{}},{"cell_type":"code","source":"for resolution in [0.2, 0.5, 1, 2]:\n    sc.tl.leiden(adata, resolution=resolution, key_added=f\"leiden_{resolution}\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:25:29.946751Z","iopub.execute_input":"2022-12-14T21:25:29.947127Z","iopub.status.idle":"2022-12-14T21:26:19.233685Z","shell.execute_reply.started":"2022-12-14T21:25:29.947092Z","shell.execute_reply":"2022-12-14T21:26:19.232153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cluster_keys = [\"cell_type\", \"G1S+2M_genes_Tirosh_mean\", \"read_depth\", \"day\", \"donor\"] + [f\"leiden_{resolution}\" for resolution in [0.2, 0.5]]\nsc.pl.umap(\n    adata,\n    color=cluster_keys,\n    frameon=False,\n#     legend_loc=\"on data\",\n    legend_fontoutline=2,\n    ncols=3\n)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:26:19.235534Z","iopub.execute_input":"2022-12-14T21:26:19.235894Z","iopub.status.idle":"2022-12-14T21:26:23.186974Z","shell.execute_reply.started":"2022-12-14T21:26:19.235853Z","shell.execute_reply":"2022-12-14T21:26:23.185777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.tl.dendrogram(adata, groupby=\"leiden_0.2\", n_pcs=15, use_rep=\"X_pca\")\nsc.pl.dendrogram(adata, groupby=\"leiden_0.2\", orientation=\"left\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:26:23.188582Z","iopub.execute_input":"2022-12-14T21:26:23.189592Z","iopub.status.idle":"2022-12-14T21:26:23.450881Z","shell.execute_reply.started":"2022-12-14T21:26:23.189547Z","shell.execute_reply":"2022-12-14T21:26:23.449757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.tl.dendrogram(adata, groupby=\"cell_type\", n_pcs=15, use_rep=\"X_pca\")\nsc.pl.dendrogram(adata, groupby=\"cell_type\", orientation=\"left\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:26:23.452535Z","iopub.execute_input":"2022-12-14T21:26:23.452887Z","iopub.status.idle":"2022-12-14T21:26:23.712761Z","shell.execute_reply.started":"2022-12-14T21:26:23.452854Z","shell.execute_reply":"2022-12-14T21:26:23.711585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cluster_keys = [\"cell_type\", ] + [f\"leiden_{resolution}\" for resolution in [0.2,]]\nsc.pl.umap(\n    adata,\n    color=cluster_keys,\n    frameon=False,\n    legend_loc=\"on data\",\n    legend_fontoutline=2,\n    ncols=3\n)","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:26:23.714334Z","iopub.execute_input":"2022-12-14T21:26:23.714667Z","iopub.status.idle":"2022-12-14T21:26:24.676558Z","shell.execute_reply.started":"2022-12-14T21:26:23.714638Z","shell.execute_reply":"2022-12-14T21:26:24.675258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Differential expression between leiden clusters","metadata":{}},{"cell_type":"code","source":"sc.tl.rank_genes_groups(adata, groupby=\"leiden_0.2\", method=\"wilcoxon\", key_added=\"rank_genes_groups_leiden_0.2\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:26:24.678354Z","iopub.execute_input":"2022-12-14T21:26:24.679650Z","iopub.status.idle":"2022-12-14T21:28:03.787144Z","shell.execute_reply.started":"2022-12-14T21:26:24.679603Z","shell.execute_reply":"2022-12-14T21:28:03.785836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.rank_genes_groups(adata, key=\"rank_genes_groups_leiden_0.2\", gene_symbols=\"gene_names\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:28:03.788825Z","iopub.execute_input":"2022-12-14T21:28:03.789610Z","iopub.status.idle":"2022-12-14T21:28:06.585651Z","shell.execute_reply.started":"2022-12-14T21:28:03.789570Z","shell.execute_reply":"2022-12-14T21:28:06.584435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for cluster in adata.obs[\"leiden_0.2\"].unique():\n    de_df = sc.get.rank_genes_groups_df(adata, group=cluster, key=\"rank_genes_groups_leiden_0.2\")\n    de_df.to_csv(f\"wilcoxon_leiden_0.2_{cluster}_vs_rest.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:28:06.587549Z","iopub.execute_input":"2022-12-14T21:28:06.588486Z","iopub.status.idle":"2022-12-14T21:28:11.669735Z","shell.execute_reply.started":"2022-12-14T21:28:06.588425Z","shell.execute_reply":"2022-12-14T21:28:11.668696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.tl.rank_genes_groups(adata, groupby=\"leiden_0.2\", method=\"wilcoxon\", key_added=\"rank_genes_groups_leiden_0.2_5_6vs1\", groups=[\"5\", \"6\"], reference=\"1\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:28:11.671475Z","iopub.execute_input":"2022-12-14T21:28:11.671797Z","iopub.status.idle":"2022-12-14T21:28:46.791827Z","shell.execute_reply.started":"2022-12-14T21:28:11.671769Z","shell.execute_reply":"2022-12-14T21:28:46.790362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.rank_genes_groups(adata, key=\"rank_genes_groups_leiden_0.2_5_6vs1\", gene_symbols=\"gene_names\")","metadata":{"execution":{"iopub.status.busy":"2022-12-14T21:28:46.793604Z","iopub.execute_input":"2022-12-14T21:28:46.794059Z","iopub.status.idle":"2022-12-14T21:28:47.390307Z","shell.execute_reply.started":"2022-12-14T21:28:46.794007Z","shell.execute_reply":"2022-12-14T21:28:47.389080Z"},"trusted":true},"execution_count":null,"outputs":[]}]}