{"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":"# Exploratory Data Analysis of 10X Multiome Using Scanpy\n\nThis notebook demonstrates the use of `scanpy` for doing single cell analysis. Examples include looking at gene expression, broken down by cell population and donor, as well as dimensionality reduction and clustering for identifying cell phenotypes and/or sub-phenotypes (further subcategories within the `cell_type` category provided in the dataset). This is compared to the groun truth cell phenotype data at the end.\n\nThe code in this notebook is the work of Vendekagon Labs and licensed under the [Apache 2.0 license](https://www.apache.org/licenses/LICENSE-2.0), like all public Kaggle notebooks.","metadata":{}},{"cell_type":"code","source":"%pip install scanpy\n%pip install leidenalg","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:36:23.808915Z","iopub.execute_input":"2023-09-21T02:36:23.810160Z","iopub.status.idle":"2023-09-21T02:36:55.555595Z","shell.execute_reply.started":"2023-09-21T02:36:23.810114Z","shell.execute_reply":"2023-09-21T02:36:55.554357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport scanpy as sc\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-21T02:36:55.558894Z","iopub.execute_input":"2023-09-21T02:36:55.559305Z","iopub.status.idle":"2023-09-21T02:37:04.011677Z","shell.execute_reply.started":"2023-09-21T02:36:55.559267Z","shell.execute_reply":"2023-09-21T02:37:04.010375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Above we install the [scanpy](https://scanpy.readthedocs.io/en/stable/) package for working with single cell data as well as `leidenlag` for clustering.","metadata":{}},{"cell_type":"markdown","source":"## Data Ingest\n\nWe read in the `multiome_train` dataset. This is a very, very large dataset in sparse matrix form (i.e. what you would get from the `melt` method from a `pandas` dataframe. We take a subset of 10% and then use `pivot` to create a gene x observation matrix, then put that in an `anndata` data structure, which is required to use `scanpy` functions.","metadata":{}},{"cell_type":"code","source":"m_train = pd.read_parquet(\"/kaggle/input/open-problems-single-cell-perturbations/multiome_train.parquet\")\nm_train","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:37:04.012822Z","iopub.execute_input":"2023-09-21T02:37:04.013572Z","iopub.status.idle":"2023-09-21T02:38:24.903813Z","shell.execute_reply.started":"2023-09-21T02:37:04.013525Z","shell.execute_reply":"2023-09-21T02:38:24.902244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"obs_meta_f = \"/kaggle/input/open-problems-single-cell-perturbations/multiome_obs_meta.csv\"\nobs_meta_df = pd.read_csv(obs_meta_f)\nobs_meta_df","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:38:24.905450Z","iopub.execute_input":"2023-09-21T02:38:24.905853Z","iopub.status.idle":"2023-09-21T02:38:24.986929Z","shell.execute_reply.started":"2023-09-21T02:38:24.905818Z","shell.execute_reply":"2023-09-21T02:38:24.985236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Subsetting Data\n\nThis is too large to process the entire dataset in a Kaggle notebook. We subset to 20% of the 'observations' in the dataset (using `observations` means we're not splitting apart genes, but just donor/cell/cell population, but keeping all transcript counts for each of those observations retained after sampling).","metadata":{}},{"cell_type":"code","source":"sub = obs_meta_df['obs_id'].sample(frac=0.20, random_state=10)\nlen(sub)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:38:24.991345Z","iopub.execute_input":"2023-09-21T02:38:24.991907Z","iopub.status.idle":"2023-09-21T02:38:25.005694Z","shell.execute_reply.started":"2023-09-21T02:38:24.991859Z","shell.execute_reply":"2023-09-21T02:38:25.004307Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"m_train_sub = m_train[m_train['obs_id'].isin(sub)]\nm_train_sub","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:38:25.007788Z","iopub.execute_input":"2023-09-21T02:38:25.008165Z","iopub.status.idle":"2023-09-21T02:38:56.155315Z","shell.execute_reply.started":"2023-09-21T02:38:25.008132Z","shell.execute_reply":"2023-09-21T02:38:56.153543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Scanpy Settings\n\nWe set Scanpy logging/message verbosity and plot settings here. ","metadata":{}},{"cell_type":"code","source":"sc.settings.verbosity = 3\nsc.logging.print_header()\nsc.settings.set_figure_params(dpi=80, facecolor='white')","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:38:56.157985Z","iopub.execute_input":"2023-09-21T02:38:56.159273Z","iopub.status.idle":"2023-09-21T02:39:15.626102Z","shell.execute_reply.started":"2023-09-21T02:38:56.159206Z","shell.execute_reply":"2023-09-21T02:39:15.624502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"unmelt = m_train_sub.pivot_table(index=\"obs_id\", columns=\"location\", values=\"normalized_count\").fillna(0)\nunmelt","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:39:15.628111Z","iopub.execute_input":"2023-09-21T02:39:15.629276Z","iopub.status.idle":"2023-09-21T02:40:45.722701Z","shell.execute_reply.started":"2023-09-21T02:39:15.629226Z","shell.execute_reply":"2023-09-21T02:40:45.721095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Conversion to AnnData\n\nNow that we have observation X gene ('location' in the source data), we can put this in the appropriate `AnnData` format.","metadata":{}},{"cell_type":"code","source":"adata = sc.AnnData(\n    X=unmelt,\n    obs=pd.DataFrame(unmelt.index, index=unmelt.index),\n    var=pd.DataFrame(unmelt.columns, index=unmelt.columns)\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:40:45.724451Z","iopub.execute_input":"2023-09-21T02:40:45.724831Z","iopub.status.idle":"2023-09-21T02:40:45.801996Z","shell.execute_reply.started":"2023-09-21T02:40:45.724799Z","shell.execute_reply":"2023-09-21T02:40:45.800459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Plotting Gene Expression\n\nWe can plot the genes expressed most highly at baseline from our donors using `scanpy`.","metadata":{}},{"cell_type":"code","source":"sc.pl.highest_expr_genes(adata, n_top=20, )","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:40:45.804171Z","iopub.execute_input":"2023-09-21T02:40:45.804615Z","iopub.status.idle":"2023-09-21T02:41:00.784847Z","shell.execute_reply.started":"2023-09-21T02:40:45.804577Z","shell.execute_reply":"2023-09-21T02:41:00.783502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"obs_info = adata.obs.reset_index(drop=True).merge(obs_meta_df, on='obs_id')\nobs_info","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:00.786371Z","iopub.execute_input":"2023-09-21T02:41:00.786790Z","iopub.status.idle":"2023-09-21T02:41:00.828812Z","shell.execute_reply.started":"2023-09-21T02:41:00.786755Z","shell.execute_reply":"2023-09-21T02:41:00.827395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plots by Donor and Cell Type\n\nWe should expect to see differences for donors and cell populations within donors, w/r/t gene expression. In the following cells we look at cell-population proportions by donor, then plot gene expression by donor and cell type.","metadata":{}},{"cell_type":"code","source":"from IPython.display import display, Markdown","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:00.830662Z","iopub.execute_input":"2023-09-21T02:41:00.831162Z","iopub.status.idle":"2023-09-21T02:41:00.837519Z","shell.execute_reply.started":"2023-09-21T02:41:00.831125Z","shell.execute_reply":"2023-09-21T02:41:00.835939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"donors = list(obs_info['donor_id'].unique())\ndonors","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:00.838954Z","iopub.execute_input":"2023-09-21T02:41:00.839850Z","iopub.status.idle":"2023-09-21T02:41:00.852870Z","shell.execute_reply.started":"2023-09-21T02:41:00.839810Z","shell.execute_reply":"2023-09-21T02:41:00.851351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for donor in donors:\n    print(f\"Cell type proportions for {donor}\")\n    print(f\"=================================\")\n    display(obs_info[obs_info['donor_id'] == donor]['cell_type'].value_counts(normalize=True))\n    print(f\"=================================\")","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:00.859876Z","iopub.execute_input":"2023-09-21T02:41:00.860311Z","iopub.status.idle":"2023-09-21T02:41:00.894410Z","shell.execute_reply.started":"2023-09-21T02:41:00.860278Z","shell.execute_reply":"2023-09-21T02:41:00.893443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cell_types = list(obs_info['cell_type'].unique())\ncell_types","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:00.895952Z","iopub.execute_input":"2023-09-21T02:41:00.897229Z","iopub.status.idle":"2023-09-21T02:41:00.904769Z","shell.execute_reply.started":"2023-09-21T02:41:00.897190Z","shell.execute_reply":"2023-09-21T02:41:00.903834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for donor in donors:\n    for cell_type in cell_types:\n        display(Markdown(f'## Highest expressed genes for {donor} and {cell_type} population.'))\n        sub_df = obs_info[(obs_info['donor_id'] == donor) & \\\n                          (obs_info['cell_type'] == cell_type)]\n        adata_sub = adata[list(sub_df['obs_id']), :]\n        sc.pl.highest_expr_genes(adata_sub, n_top=20, )","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:00.905953Z","iopub.execute_input":"2023-09-21T02:41:00.906530Z","iopub.status.idle":"2023-09-21T02:41:31.514519Z","shell.execute_reply.started":"2023-09-21T02:41:00.906497Z","shell.execute_reply":"2023-09-21T02:41:31.513062Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cleaning Data\n\nSingle cell data is typically quite noisy! Here we filter down cells that don't express many genes, as well as genes that are not expressed across many cells.\n\nThis next section follows examples from the [Scanpy tutorial](https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html), which also processes pbmc data.","metadata":{}},{"cell_type":"code","source":"sc.pp.filter_cells(adata, min_genes=200)\nsc.pp.filter_genes(adata, min_cells=3)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:41:31.516236Z","iopub.execute_input":"2023-09-21T02:41:31.516734Z","iopub.status.idle":"2023-09-21T02:42:08.966915Z","shell.execute_reply.started":"2023-09-21T02:41:31.516688Z","shell.execute_reply":"2023-09-21T02:42:08.965596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Quality Checks\n\nWe can use mitochondrial DNA to check the quality of our single cell data. ","metadata":{}},{"cell_type":"code","source":"adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'\nsc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:42:08.968685Z","iopub.execute_input":"2023-09-21T02:42:08.969323Z","iopub.status.idle":"2023-09-21T02:42:12.753804Z","shell.execute_reply.started":"2023-09-21T02:42:08.969280Z","shell.execute_reply":"2023-09-21T02:42:12.752538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],\n             jitter=0.4, multi_panel=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:42:12.756604Z","iopub.execute_input":"2023-09-21T02:42:12.757954Z","iopub.status.idle":"2023-09-21T02:42:14.621090Z","shell.execute_reply.started":"2023-09-21T02:42:12.757915Z","shell.execute_reply":"2023-09-21T02:42:14.620040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')\nsc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:42:14.622512Z","iopub.execute_input":"2023-09-21T02:42:14.623059Z","iopub.status.idle":"2023-09-21T02:42:15.890957Z","shell.execute_reply.started":"2023-09-21T02:42:14.623026Z","shell.execute_reply":"2023-09-21T02:42:15.889973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Remove Bad Total Count or Ratio Cells\n\nLooking at the above plots, we can see where we have calls with low total counts where a high number of counts are just mitochondrial DNA. Mitochondria are large and will retain their transcripts when the rest of the cytoplasmic RNA has been lost or drained away during the assay. So if we trim some of these, this will likely improve our remaining cell quality.","metadata":{}},{"cell_type":"code","source":"# it looks like a ceiling of 20k gets us to a better place on both plots\nadata = adata[adata.obs.n_genes_by_counts < 20000, :]\n# we can also trim out a little bit of the higher mitochondrial expression\nadata = adata[adata.obs.pct_counts_mt < 0.9, :]","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:42:15.892164Z","iopub.execute_input":"2023-09-21T02:42:15.893281Z","iopub.status.idle":"2023-09-21T02:42:15.904066Z","shell.execute_reply.started":"2023-09-21T02:42:15.893245Z","shell.execute_reply":"2023-09-21T02:42:15.902511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Renormalize\n\nRenormalize to 10,000, smaller values, log transform, and add annotations for highly variable genes.","metadata":{}},{"cell_type":"code","source":"sc.pp.normalize_total(adata, target_sum=1e4)\nsc.pp.log1p(adata)\nsc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)\nsc.pl.highly_variable_genes(adata)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:42:15.905702Z","iopub.execute_input":"2023-09-21T02:42:15.906201Z","iopub.status.idle":"2023-09-21T02:42:39.832961Z","shell.execute_reply.started":"2023-09-21T02:42:15.906167Z","shell.execute_reply":"2023-09-21T02:42:39.831398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we set the raw data basis to the cleaned, and normalized data for follow on processing. We also regress out previously calculated variables `total_counts` and `pct_counts_mt` (during the QC process), and scale the remaining data to unit variance, clipping out values more than 10 standard deviations from the mean center.","metadata":{}},{"cell_type":"code","source":"adata.raw = adata\nadata = adata[:, adata.var.highly_variable]\nsc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])\nsc.pp.scale(adata, max_value=10)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:42:39.834698Z","iopub.execute_input":"2023-09-21T02:42:39.836458Z","iopub.status.idle":"2023-09-21T02:48:48.392310Z","shell.execute_reply.started":"2023-09-21T02:42:39.836403Z","shell.execute_reply":"2023-09-21T02:48:48.390877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dimensionality Reduction\n\nBecause the space of genes expressed is quite large, and because many genes participate in the same pathways and networks (eg pro-inflammatory genes, immune suppressing genes, proliferatory genes, etc.), the intrinsic dimensionality is quite lower. Both PCA and comparable methods, and more sophisticated methods like UMAP, are used in different aspects of single cell analysis.","metadata":{}},{"cell_type":"code","source":"sc.tl.pca(adata, svd_solver='arpack')\nsc.pl.pca(adata, color='CST3')","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:48:48.394103Z","iopub.execute_input":"2023-09-21T02:48:48.394466Z","iopub.status.idle":"2023-09-21T02:49:20.535843Z","shell.execute_reply.started":"2023-09-21T02:48:48.394435Z","shell.execute_reply":"2023-09-21T02:49:20.534766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.pca_variance_ratio(adata, log=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:49:20.537622Z","iopub.execute_input":"2023-09-21T02:49:20.538501Z","iopub.status.idle":"2023-09-21T02:49:20.929994Z","shell.execute_reply.started":"2023-09-21T02:49:20.538444Z","shell.execute_reply":"2023-09-21T02:49:20.928572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"save_file = '/kaggle/working/processed_multiome_traing.h5ad'\nadata.write(save_file)\nadata","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:49:20.931874Z","iopub.execute_input":"2023-09-21T02:49:20.932250Z","iopub.status.idle":"2023-09-21T02:49:40.486128Z","shell.execute_reply.started":"2023-09-21T02:49:20.932217Z","shell.execute_reply":"2023-09-21T02:49:40.484723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### UMAP and Neighborhood Graph\n\nUMAP visualization is often a good first step toward identifying clusters by eyeball that might correspond to interesting cell phenotypes.","metadata":{}},{"cell_type":"code","source":"sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)\nsc.tl.leiden(adata)\nsc.tl.paga(adata)\nsc.pl.paga(adata, plot=True)\nsc.tl.umap(adata, init_pos='paga')\nsc.tl.umap(adata)\nsc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:49:40.488965Z","iopub.execute_input":"2023-09-21T02:49:40.490390Z","iopub.status.idle":"2023-09-21T02:50:08.240786Z","shell.execute_reply.started":"2023-09-21T02:49:40.490330Z","shell.execute_reply":"2023-09-21T02:50:08.239369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# antigen responsive CD8 signature from: https://www.frontiersin.org/articles/10.3389/fimmu.2019.02568/full\n# note: genes have many different names in different assemblies, and their names change, or genes are combined\n#       or split, etc. -- it's not always trivial to resolve a set of genes found in the literature to genes\n#       present in a dataset of interest.\ngenes = ['TNFRSF9', 'REL', 'EGR2', 'SRP14', 'FASLG', 'GZMB', 'IFNG', 'CD69', 'HMGB1', 'NAB2']\n[gene in adata.var_names for gene in genes]","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:50:08.243074Z","iopub.execute_input":"2023-09-21T02:50:08.243564Z","iopub.status.idle":"2023-09-21T02:50:08.254224Z","shell.execute_reply.started":"2023-09-21T02:50:08.243526Z","shell.execute_reply":"2023-09-21T02:50:08.252332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sc.pl.umap(adata, color=['TNFRSF9', 'EGR2', 'FASLG', 'NAB2'])","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:50:08.255970Z","iopub.execute_input":"2023-09-21T02:50:08.256441Z","iopub.status.idle":"2023-09-21T02:50:09.925673Z","shell.execute_reply.started":"2023-09-21T02:50:08.256398Z","shell.execute_reply":"2023-09-21T02:50:09.924574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Clustering\n\nAs in the Scanpy tutorial, we use the Leiden graph-clustering method by\n[Traag et al. (2018)](https://arxiv.org/abs/1810.08473). ","metadata":{}},{"cell_type":"code","source":"sc.pl.umap(adata, color=['leiden'])","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:50:09.927168Z","iopub.execute_input":"2023-09-21T02:50:09.927595Z","iopub.status.idle":"2023-09-21T02:50:10.722029Z","shell.execute_reply.started":"2023-09-21T02:50:09.927558Z","shell.execute_reply":"2023-09-21T02:50:10.720552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# persist all of our work thus far to our saved file\nadata.write(save_file)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:50:10.723841Z","iopub.execute_input":"2023-09-21T02:50:10.724713Z","iopub.status.idle":"2023-09-21T02:50:23.514442Z","shell.execute_reply.started":"2023-09-21T02:50:10.724666Z","shell.execute_reply":"2023-09-21T02:50:23.512644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Identifying Clusters with Gene Markers\n\nWhen clusters of interest are identified, it usually makes sense to next explore which genes are expressed or downregulated most strongly within each clusters. This will give us a clue as to what cell types or phenotypes are present within the cluster.","metadata":{}},{"cell_type":"code","source":"sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')\nsc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:50:23.516124Z","iopub.execute_input":"2023-09-21T02:50:23.516876Z","iopub.status.idle":"2023-09-21T02:51:57.724108Z","shell.execute_reply.started":"2023-09-21T02:50:23.516835Z","shell.execute_reply":"2023-09-21T02:51:57.722792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see several genes of interest in the groups above. Several clusters are all high for [IL7R](https://medlineplus.gov/genetics/gene/il7r/). IL7R encodes the IL-7 receptor, which is present on T and B lymphocytes and their progenitors. We also see several other highly expressed genes, like CAMK4 and RORA, though these are more broadly expressed across types. We also see some `MT-` genes, so on first viewing we cannot identify every cluster confidently with a broader lineage.\n\nNote that different methods of comparing gene expression between groups yield different results. Those found in `DESeq1` and similar packages usually yield more reliable results.","metadata":{}},{"cell_type":"code","source":"sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')\nsc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T02:51:57.726192Z","iopub.execute_input":"2023-09-21T02:51:57.726602Z","iopub.status.idle":"2023-09-21T03:01:02.213145Z","shell.execute_reply.started":"2023-09-21T02:51:57.726567Z","shell.execute_reply":"2023-09-21T03:01:02.210770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Drilling down further into phenotype\n\nWe can also look at how genes of interest are distributed across our clusters.","metadata":{}},{"cell_type":"code","source":"sc.pl.violin(adata, ['CD8A', 'FOXP3', 'EGR2', 'IL7R'], groupby='leiden')","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:01:02.215280Z","iopub.execute_input":"2023-09-21T03:01:02.216202Z","iopub.status.idle":"2023-09-21T03:01:05.226180Z","shell.execute_reply.started":"2023-09-21T03:01:02.216150Z","shell.execute_reply":"2023-09-21T03:01:05.224713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Seeing T-cell and other lymphocyte associated genes across clusters shouldn't be particularly surprising, given that this dataset has been pre-filtered to include those cells in larger numbers, especially in training data. But it is interesting to see for instance that `FOXP3` only appears to be present in one clusters.\n\nNow, if we know a priori what some marker genes are, we can use these to try and identify lineages by clusters.","metadata":{}},{"cell_type":"code","source":"marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',\n                'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',\n                'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']\nsc.pl.dotplot(adata, marker_genes, groupby='leiden');","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:01:05.228326Z","iopub.execute_input":"2023-09-21T03:01:05.228813Z","iopub.status.idle":"2023-09-21T03:01:06.132067Z","shell.execute_reply.started":"2023-09-21T03:01:05.228768Z","shell.execute_reply":"2023-09-21T03:01:06.130566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now that we've wrapped our exploratory analysis, save our file one more time.\nadata.write(save_file, compression='gzip')","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:01:06.133806Z","iopub.execute_input":"2023-09-21T03:01:06.134256Z","iopub.status.idle":"2023-09-21T03:02:30.062854Z","shell.execute_reply.started":"2023-09-21T03:01:06.134217Z","shell.execute_reply":"2023-09-21T03:02:30.061245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Comparison to Known Lineages\n\nMany of the outcomes of our analysis above, like leiden clusters, have been added to our `obs` data frame. We can extract this from our `AnnData` variable, `adata`.","metadata":{}},{"cell_type":"code","source":"adata.obs","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:02:30.065590Z","iopub.execute_input":"2023-09-21T03:02:30.066733Z","iopub.status.idle":"2023-09-21T03:02:30.091856Z","shell.execute_reply.started":"2023-09-21T03:02:30.066679Z","shell.execute_reply":"2023-09-21T03:02:30.090547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we look at the distribution of known cell types across clusters, as well as donors, as a final sanity check.","metadata":{}},{"cell_type":"code","source":"merged = obs_meta_df.merge(adata.obs.reset_index(drop=True), on='obs_id')\nmerged","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:02:30.093516Z","iopub.execute_input":"2023-09-21T03:02:30.093967Z","iopub.status.idle":"2023-09-21T03:02:30.153006Z","shell.execute_reply.started":"2023-09-21T03:02:30.093932Z","shell.execute_reply":"2023-09-21T03:02:30.151664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters = merged.groupby('leiden')['cell_type'].value_counts(normalize=True)\npd.set_option('display.max_rows', len(clusters))\nclusters","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:02:30.154985Z","iopub.execute_input":"2023-09-21T03:02:30.155751Z","iopub.status.idle":"2023-09-21T03:02:30.190127Z","shell.execute_reply.started":"2023-09-21T03:02:30.155704Z","shell.execute_reply":"2023-09-21T03:02:30.188757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The results above show that we get very good separation across phenotypes by clusters. We also see indications that we might be looking at interesting sub-populations, e.g. Tregs (regulatory CD4+ T Cells in one of the clusters) -- see the `FOXP3` violin plot several cells above.\n\nWhen we look at donors split across clusters, we can see that there does seem to be some bias of particular donors' cells being more concentrated in certain clusters instead of others. But we do have some mixture across donors present throughout the clusters.","metadata":{}},{"cell_type":"code","source":"merged.groupby('leiden')['donor_id'].value_counts(normalize=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T03:02:30.191904Z","iopub.execute_input":"2023-09-21T03:02:30.192285Z","iopub.status.idle":"2023-09-21T03:02:30.217130Z","shell.execute_reply.started":"2023-09-21T03:02:30.192253Z","shell.execute_reply":"2023-09-21T03:02:30.215810Z"},"trusted":true},"execution_count":null,"outputs":[]}]}