scanpy_pbmc3k.ipynb

{"cells":[{"cell_type":"markdown","metadata":{},"source":"# Preprocessing and clustering 3k PBMCs"},{"cell_type":"markdown","metadata":{},"source":"In May 2017, this started out as a demonstration that Scanpy would allow to reproduce most of Seurat's [guided clustering tutorial](http://satijalab.org/seurat/pbmc3k_tutorial.html) ([Satija et al., 2015](https://doi.org/10.1038/nbt.3192)).\n\nWe gratefully acknowledge Seurat's authors for the tutorial! In the meanwhile, we have added and removed a few pieces.\n\nThe data consist of *3k PBMCs from a Healthy Donor* and are freely available from 10x Genomics ([here](http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz) from this [webpage](https://support.10xgenomics.com/single-cell-gene-expression/datasets/1.1.0/pbmc3k))."},{"cell_type":"code","execution_count":null,"metadata":{"jupyter":{"source_hidden":false},"tags":[]},"outputs":[],"source":"import pandas as pd\nimport scanpy as sc\nimport warnings\n\nfrom pathlib import Path\n\nwarnings.filterwarnings(\"ignore\")"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.settings.verbosity = 3  # verbosity: errors (0), warnings (1), info (2), hints (3)\nsc.logging.print_header()\nsc.settings.set_figure_params(dpi=80, facecolor=\"white\")"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"# configuration for capsule\n\nsc._settings.ScanpyConfig.cachedir = Path(\"/data/cache/\")\nsc._settings.ScanpyConfig.datasetdir = Path(\"/data/\")\nPath(\"/results/figures\").mkdir(parents=True, exist_ok=True)\nsc._settings.ScanpyConfig.figdir = Path(\"/results/figures\")\nPath(\"/results/write\").mkdir(parents=True, exist_ok=True)\nsc._settings.ScanpyConfig.writedir = Path(\"/results/write/\")"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"results_dir = \"/results\"\nresults_file = (\n    f\"{results_dir}/write/pbmc3k.h5ad\"  # the file that will store the analysis results\n)"},{"cell_type":"markdown","metadata":{},"source":"Read in the count matrix into an [AnnData](https://anndata.readthedocs.io/en/latest/\n) object, which holds many slots for annotations and different representations of the data. It also comes with its own HDF5-based file format: `.h5ad`."},{"cell_type":"code","execution_count":null,"metadata":{"scrolled":true,"tags":[]},"outputs":[],"source":"adata = sc.read_10x_mtx(\n    \"/data/filtered_gene_bc_matrices/hg19/\",  # the directory with the `.mtx` file\n    var_names=\"gene_symbols\",  # use gene symbols for the variable names (variables-axis index)\n    cache=False,\n)  # write a cache file for faster subsequent reading"},{"cell_type":"markdown","metadata":{},"source":"<div class=\"alert alert-info\">\n\n**Note**\n    \nSee [anndata-tutorials/getting-started](https://anndata-tutorials.readthedocs.io/en/latest/getting-started.html) for a more comprehensive introduction to `AnnData`.\n\n</div>"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.var_names_make_unique()  # this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata"},{"cell_type":"markdown","metadata":{},"source":"## Preprocessing"},{"cell_type":"markdown","metadata":{},"source":"Show those genes that yield the highest fraction of counts in each single cell, across all cells."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.highest_expr_genes(\n    adata,\n    n_top=20,\n)"},{"cell_type":"markdown","metadata":{},"source":"Basic filtering:"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.filter_cells(adata, min_genes=200)\nsc.pp.filter_genes(adata, min_cells=3)"},{"cell_type":"markdown","metadata":{},"source":"Let's assemble some information about mitochondrial genes, which are important for quality control.\n\nCiting from \"Simple Single Cell\" workflows [(Lun, McCarthy & Marioni, 2017)](https://master.bioconductor.org/packages/release/workflows/html/simpleSingleCell.html#examining-gene-level-metrics):\n\n> High proportions are indicative of poor-quality cells (Islam et al. 2014; Ilicic et al. 2016), possibly because of loss of cytoplasmic RNA from perforated cells. The reasoning is that mitochondria are larger than individual transcript molecules and less likely to escape through tears in the cell membrane."},{"cell_type":"markdown","metadata":{},"source":"With `pp.calculate_qc_metrics`, we can compute many metrics very efficiently."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.var[\"mt\"] = adata.var_names.str.startswith(\n    \"MT-\"\n)  # annotate the group of mitochondrial genes as 'mt'\nsc.pp.calculate_qc_metrics(\n    adata, qc_vars=[\"mt\"], percent_top=None, log1p=False, inplace=True\n)"},{"cell_type":"markdown","metadata":{},"source":"A violin plot of some of the computed quality measures:\n\n* the number of genes expressed in the count matrix\n* the total counts per cell\n* the percentage of counts in mitochondrial genes"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.violin(\n    adata,\n    [\"n_genes_by_counts\", \"total_counts\", \"pct_counts_mt\"],\n    jitter=0.4,\n    multi_panel=True,\n)"},{"cell_type":"markdown","metadata":{},"source":"Remove cells that have too many mitochondrial genes expressed or too many total counts:"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.scatter(adata, x=\"total_counts\", y=\"pct_counts_mt\")\nsc.pl.scatter(adata, x=\"total_counts\", y=\"n_genes_by_counts\")"},{"cell_type":"markdown","metadata":{},"source":"Actually do the filtering by slicing the `AnnData` object."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata = adata[adata.obs.n_genes_by_counts < 2500, :]\nadata = adata[adata.obs.pct_counts_mt < 5, :]"},{"cell_type":"markdown","metadata":{},"source":"Total-count normalize (library-size correct) the data matrix $\\mathbf{X}$ to 10,000 reads per cell, so that counts become comparable among cells."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.normalize_total(adata, target_sum=1e4)"},{"cell_type":"markdown","metadata":{},"source":"Logarithmize the data:"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.log1p(adata)"},{"cell_type":"markdown","metadata":{},"source":"Identify highly-variable genes."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.highly_variable_genes(adata)"},{"cell_type":"markdown","metadata":{},"source":"Set the `.raw` attribute of the AnnData object to the normalized and logarithmized raw gene expression for later use in differential testing and visualizations of gene expression. This simply freezes the state of the AnnData object."},{"cell_type":"markdown","metadata":{},"source":"<div class=\"alert alert-info\">\n\n**Note**\n\nYou can get back an `AnnData` of the object in `.raw` by calling `.raw.to_adata()`.\n\n</div>    "},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.raw = adata"},{"cell_type":"markdown","metadata":{},"source":"<div class=\"alert alert-info\">\n\n**Note**\n    \nIf you don't proceed below with correcting the data with `sc.pp.regress_out` and scaling it via `sc.pp.scale`, you can also get away without using `.raw` at all. Sometimes sc.pp.regress_out can overcorrect and it can be preferable to not do that.\n\nThe result of the previous highly-variable-genes detection is stored as an annotation in `.var.highly_variable` and auto-detected by PCA and hence, `sc.pp.neighbors` and subsequent manifold/graph tools. In that case, the step *actually do the filtering* below is unnecessary, too.\n\n</div>"},{"cell_type":"markdown","metadata":{},"source":"Actually do the filtering"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata = adata[:, adata.var.highly_variable]"},{"cell_type":"markdown","metadata":{},"source":"Regress out effects of total counts per cell and the percentage of mitochondrial genes expressed. Scale the data to unit variance."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.regress_out(adata, [\"total_counts\", \"pct_counts_mt\"])"},{"cell_type":"markdown","metadata":{},"source":"Scale each gene to unit variance. Clip values exceeding standard deviation 10. "},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.scale(adata, max_value=10)"},{"cell_type":"markdown","metadata":{},"source":"## Principal component analysis"},{"cell_type":"markdown","metadata":{},"source":"Reduce the dimensionality of the data by running principal component analysis (PCA), which reveals the main axes of variation and denoises the data."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.pca(adata, svd_solver=\"arpack\")"},{"cell_type":"markdown","metadata":{},"source":"We can make a scatter plot in the PCA coordinates, but we will not use that later on."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.pca(adata, color=\"CST3\")"},{"cell_type":"markdown","metadata":{},"source":"Let us inspect the contribution of single PCs to the total variance in the data. This gives us information about how many PCs we should consider in order to compute the neighborhood relations of cells, e.g. used in the clustering function  `sc.tl.louvain()` or tSNE `sc.tl.tsne()`. In our experience, often a rough estimate of the number of PCs does fine."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.pca_variance_ratio(adata, log=True)"},{"cell_type":"markdown","metadata":{},"source":"Save the result."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.write(results_file)"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata"},{"cell_type":"markdown","metadata":{},"source":"## Computing the neighborhood graph"},{"cell_type":"markdown","metadata":{},"source":"Let us compute the neighborhood graph of cells using the PCA representation of the data matrix. You might simply use default values here. For the sake of reproducing Seurat's results, let's take the following values."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)"},{"cell_type":"markdown","metadata":{},"source":"## Embedding the neighborhood graph"},{"cell_type":"markdown","metadata":{},"source":"We suggest embedding the graph in two dimensions using UMAP ([McInnes et al., 2018](https://arxiv.org/abs/1802.03426)), see below. It is potentially more faithful to the global connectivity of the manifold than tSNE, i.e., it better preserves trajectories. In some ocassions, you might still observe disconnected clusters and similar connectivity violations. They can usually be remedied by running:\n\n```Python\nsc.tl.paga(adata)\nsc.pl.paga(adata, plot=False)  # remove `plot=False` if you want to see the coarse-grained graph\nsc.tl.umap(adata, init_pos=\"paga\")\n```"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.umap(adata)"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.var[\"gene_ids\"]"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.umap(adata, color=[\"CST3\", \"NKG7\", \"PPBP\"])"},{"cell_type":"markdown","metadata":{},"source":"As we set the `.raw` attribute of `adata`, the previous plots showed the \"raw\" (normalized, logarithmized, but uncorrected) gene expression. You can also plot the scaled and corrected gene expression by explicitly stating that you don't want to use `.raw`."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.umap(adata, color=[\"CST3\", \"NKG7\", \"PPBP\"], use_raw=False)"},{"cell_type":"markdown","metadata":{},"source":"## Clustering the neighborhood graph"},{"cell_type":"markdown","metadata":{},"source":"As with Seurat and many other frameworks, we recommend the Leiden graph-clustering method (community detection based on optimizing modularity) by [Traag *et al.* (2019)](https://www.nature.com/articles/s41598-019-41695-z). Note that Leiden clustering directly clusters the neighborhood graph of cells, which we already computed in the previous section."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.leiden(adata)"},{"cell_type":"markdown","metadata":{},"source":"Plot the clusters, which agree quite well with the result of Seurat."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.umap(adata, color=[\"leiden\", \"CST3\", \"NKG7\"])"},{"cell_type":"markdown","metadata":{},"source":"Save the result."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.write(results_file)"},{"cell_type":"markdown","metadata":{},"source":"## Finding marker genes"},{"cell_type":"markdown","metadata":{},"source":"Let us compute a ranking for the highly differential genes in each cluster. For this, by default, the `.raw` attribute of AnnData is used in case it has been initialized before. The simplest and fastest method to do so is the t-test."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.rank_genes_groups(adata, \"leiden\", method=\"t-test\")\nsc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.settings.verbosity = 2  # reduce the verbosity"},{"cell_type":"markdown","metadata":{},"source":"The result of a [Wilcoxon rank-sum (Mann-Whitney-U)](https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test) test is very similar. We recommend using the latter in publications, see e.g., [Sonison & Robinson (2018)](https://doi.org/10.1038/nmeth.4612). You might also consider much more powerful differential testing packages like MAST, limma, DESeq2 and, for python, the recent diffxpy."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.rank_genes_groups(adata, \"leiden\", method=\"wilcoxon\")\nsc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)"},{"cell_type":"markdown","metadata":{},"source":"Save the result."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.write(results_file)"},{"cell_type":"markdown","metadata":{},"source":"As an alternative, let us rank genes using logistic regression. For instance, this has been suggested by [Ntranos et al. (2018)](https://doi.org/10.1101/258566). The essential difference is that here, we use a multi-variate approach whereas conventional differential tests are uni-variate. [Clark et al. (2014)](https://doi.org/10.1186/1471-2105-15-79) has more details."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.rank_genes_groups(adata, \"leiden\", method=\"logreg\")\nsc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)"},{"cell_type":"markdown","metadata":{},"source":"With the exceptions of *IL7R*, which is only found by the t-test and *FCER1A*, which is only found by the other two approaches, all marker genes are recovered in all approaches.\n\nLouvain Group | Markers | Cell Type\n---|---|---\n0 | IL7R | CD4 T cells\n1 | CD14, LYZ | CD14+ Monocytes\n2 | MS4A1 |\tB cells\n3 | CD8A |\tCD8 T cells\n4 | GNLY, NKG7 | \tNK cells\n5 | FCGR3A, MS4A7 |\tFCGR3A+ Monocytes\n6 | FCER1A, CST3 |\tDendritic Cells\n7 | PPBP |\tMegakaryocytes"},{"cell_type":"markdown","metadata":{},"source":"Let us also define a list of marker genes for later reference."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"marker_genes = [\n    \"IL7R\",\n    \"CD79A\",\n    \"MS4A1\",\n    \"CD8A\",\n    \"CD8B\",\n    \"LYZ\",\n    \"CD14\",\n    \"LGALS3\",\n    \"S100A8\",\n    \"GNLY\",\n    \"NKG7\",\n    \"KLRB1\",\n    \"FCGR3A\",\n    \"MS4A7\",\n    \"FCER1A\",\n    \"CST3\",\n    \"PPBP\",\n]"},{"cell_type":"markdown","metadata":{},"source":"Reload the object that has been save with the Wilcoxon Rank-Sum test result."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata = sc.read(results_file)"},{"cell_type":"markdown","metadata":{},"source":"Show the 10 top ranked genes per cluster 0, 1, ..., 7 in a dataframe."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"pd.DataFrame(adata.uns[\"rank_genes_groups\"][\"names\"]).head(5)"},{"cell_type":"markdown","metadata":{},"source":"Get a table with the scores and groups."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"result = adata.uns[\"rank_genes_groups\"]\ngroups = result[\"names\"].dtype.names\npd.DataFrame(\n    {\n        group + \"_\" + key[:1]: result[key][group]\n        for group in groups\n        for key in [\"names\", \"pvals\"]\n    }\n).head(5)"},{"cell_type":"markdown","metadata":{},"source":"Compare to a single cluster: "},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.tl.rank_genes_groups(adata, \"leiden\", groups=[\"0\"], reference=\"1\", method=\"wilcoxon\")\nsc.pl.rank_genes_groups(adata, groups=[\"0\"], n_genes=20)"},{"cell_type":"markdown","metadata":{},"source":"If we want a more detailed view for a certain group, use `sc.pl.rank_genes_groups_violin`."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.rank_genes_groups_violin(adata, groups=\"0\", n_genes=8)"},{"cell_type":"markdown","metadata":{},"source":"Reload the object with the computed differential expression (i.e. DE via a comparison with the rest of the groups):"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata = sc.read(results_file)"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.rank_genes_groups_violin(adata, groups=\"0\", n_genes=8)"},{"cell_type":"markdown","metadata":{},"source":"If you want to compare a certain gene across groups, use the following."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.violin(adata, [\"CST3\", \"NKG7\", \"PPBP\"], groupby=\"leiden\")"},{"cell_type":"markdown","metadata":{},"source":"Actually mark the cell types."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":""},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.umap(\n    adata, color=\"leiden\", legend_loc=\"on data\", title=\"\", frameon=False, save=\".pdf\"\n)"},{"cell_type":"markdown","metadata":{},"source":"Now that we annotated the cell types, let us visualize the marker genes."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.dotplot(adata, marker_genes, groupby=\"leiden\")"},{"cell_type":"markdown","metadata":{},"source":"There is also a very compact violin plot."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"sc.pl.stacked_violin(adata, marker_genes, groupby=\"leiden\", rotation=90)"},{"cell_type":"markdown","metadata":{},"source":"During the course of this analysis, the AnnData accumlated the following annotations."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata"},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.write(\n    results_file, compression=\"gzip\"\n)  # `compression=\"gzip\"` saves disk space, but slows down writing and subsequent reading"},{"cell_type":"markdown","metadata":{},"source":"Get a rough overview of the file using `h5ls`, which has many options - for more details see [here](https://github.com/theislab/scanpy_usage/blob/master/170505_seurat/info_h5ad.md). The file format might still be subject to further optimization in the future. All reading functions will remain backwards-compatible, though."},{"cell_type":"markdown","metadata":{},"source":"If you want to share this file with people who merely want to use it for visualization, a simple way to reduce the file size is by removing the dense scaled and corrected data matrix. The file still contains the raw data used in the visualizations in `adata.raw`."},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":"adata.raw.to_adata().write(f\"{results_dir}/write/pbmc3k_withoutX.h5ad\")"}],"metadata":{"anaconda-cloud":{},"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"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"}},"nbformat":4,"nbformat_minor":4}