{ "cells": [ { "cell_type": "raw", "metadata": {}, "source": [ "---\n", "description: Quality control of single cell RNA-Seq data. Inspection of\n", " QC metrics including number of UMIs, number of genes expressed,\n", " mitochondrial and ribosomal expression, sex and cell cycle state.\n", "subtitle: Scanpy Toolkit\n", "title: Quality Control\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "> **Note**\n", ">\n", "> Code chunks run Python commands unless it starts with `%%bash`, in\n", "> which case, those chunks run shell commands.\n", "\n", "
\n", "\n", "## Get data\n", "\n", "In this tutorial, we will run all tutorials with a set of 8 PBMC 10x\n", "datasets from 4 covid-19 patients and 4 healthy controls, the samples\n", "have been subsampled to 1500 cells per sample. We can start by defining\n", "our paths." ] }, { "cell_type": "code", "metadata": {}, "source": [ "import os\n", "\n", "# download pre-computed annotation\n", "fetch_annotation = False\n", "\n", "path_data = \"https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq\"\n", "\n", "path_covid = \"./data/covid\"\n", "if not os.path.exists(path_covid):\n", " os.makedirs(path_covid, exist_ok=True)\n", "\n", "path_results = \"data/covid/results\"\n", "if not os.path.exists(path_results):\n", " os.makedirs(path_results, exist_ok=True)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": {}, "source": [ "import urllib.request\n", "\n", "file_list = [\n", " \"normal_pbmc_13.h5\", \"normal_pbmc_14.h5\", \"normal_pbmc_19.h5\", \"normal_pbmc_5.h5\",\n", " \"ncov_pbmc_15.h5\", \"ncov_pbmc_16.h5\", \"ncov_pbmc_17.h5\", \"ncov_pbmc_1.h5\"\n", "]\n", "\n", "for i in file_list:\n", " path_file = os.path.join(path_covid, i)\n", " if not os.path.exists(path_file):\n", " file_url = os.path.join(path_data, \"covid\", i)\n", " urllib.request.urlretrieve(file_url, path_file)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With data in place, now we can start loading libraries we will use in\n", "this tutorial." ] }, { "cell_type": "code", "metadata": {}, "source": [ "import numpy as np\n", "import pandas as pd\n", "import scanpy as sc\n", "import warnings\n", "\n", "warnings.simplefilter(action='ignore', category=Warning)\n", "\n", "# verbosity: errors (0), warnings (1), info (2), hints (3)\n", "sc.settings.verbosity = 3\n", "sc.settings.set_figure_params(dpi=80)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can first load the data individually by reading directly from HDF5\n", "file format (.h5).\n", "\n", "In Scanpy we read them into an Anndata object with the the function\n", "`read_10x_h5`" ] }, { "cell_type": "code", "metadata": {}, "source": [ "data_cov1 = sc.read_10x_h5(os.path.join(path_covid,'ncov_pbmc_1.h5'))\n", "data_cov1.var_names_make_unique()\n", "data_cov15 = sc.read_10x_h5(os.path.join(path_covid,'ncov_pbmc_15.h5'))\n", "data_cov15.var_names_make_unique()\n", "data_cov16 = sc.read_10x_h5(os.path.join(path_covid,'ncov_pbmc_16.h5'))\n", "data_cov16.var_names_make_unique()\n", "data_cov17 = sc.read_10x_h5(os.path.join(path_covid,'ncov_pbmc_17.h5'))\n", "data_cov17.var_names_make_unique()\n", "data_ctrl5 = sc.read_10x_h5(os.path.join(path_covid,'normal_pbmc_5.h5'))\n", "data_ctrl5.var_names_make_unique()\n", "data_ctrl13 = sc.read_10x_h5(os.path.join(path_covid,'normal_pbmc_13.h5'))\n", "data_ctrl13.var_names_make_unique()\n", "data_ctrl14 = sc.read_10x_h5(os.path.join(path_covid,'normal_pbmc_14.h5'))\n", "data_ctrl14.var_names_make_unique()\n", "data_ctrl19 = sc.read_10x_h5(os.path.join(path_covid,'normal_pbmc_19.h5'))\n", "data_ctrl19.var_names_make_unique()" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Collate\n", "\n", "{{< qc_collate_1 >}}" ] }, { "cell_type": "code", "metadata": {}, "source": [ "# add some metadata\n", "data_cov1.obs['type']=\"Covid\"\n", "data_cov1.obs['sample']=\"covid_1\"\n", "data_cov15.obs['type']=\"Covid\"\n", "data_cov15.obs['sample']=\"covid_15\"\n", "data_cov16.obs['type']=\"Covid\"\n", "data_cov16.obs['sample']=\"covid_16\"\n", "data_cov17.obs['type']=\"Covid\"\n", "data_cov17.obs['sample']=\"covid_17\"\n", "data_ctrl5.obs['type']=\"Ctrl\"\n", "data_ctrl5.obs['sample']=\"ctrl_5\"\n", "data_ctrl13.obs['type']=\"Ctrl\"\n", "data_ctrl13.obs['sample']=\"ctrl_13\"\n", "data_ctrl14.obs['type']=\"Ctrl\"\n", "data_ctrl14.obs['sample']=\"ctrl_14\"\n", "data_ctrl19.obs['type']=\"Ctrl\"\n", "data_ctrl19.obs['sample']=\"ctrl_19\"\n", "\n", "# merge into one object.\n", "adata = data_cov1.concatenate(data_cov15, data_cov16, data_cov17, data_ctrl5, data_ctrl13, data_ctrl14, data_ctrl19)\n", "\n", "# and delete individual datasets to save space\n", "del(data_cov1, data_cov15, data_cov16, data_cov17)\n", "del(data_ctrl5, data_ctrl13, data_ctrl14, data_ctrl19)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can print a summary of the datasets in the Scanpy object, or a\n", "summary of the whole object." ] }, { "cell_type": "code", "metadata": {}, "source": [ "print(adata.obs['sample'].value_counts())\n", "adata" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate QC\n", "\n", "Having the data in a suitable format, we can start calculating some\n", "quality metrics. We can for example calculate the percentage of\n", "mitochondrial and ribosomal genes per cell and add to the metadata. The\n", "proportion of hemoglobin genes can give an indication of red blood cell\n", "contamination. This will be helpful to visualize them across different\n", "metadata parameters (i.e. datasetID and chemistry version). There are\n", "several ways of doing this. The QC metrics are finally added to the\n", "metadata table.\n", "\n", "Citing from Simple Single Cell workflows (Lun, McCarthy & Marioni,\n", "2017): High proportions are indicative of poor-quality cells (Islam et\n", "al. 2014; Ilicic et al. 2016), possibly because of loss of cytoplasmic\n", "RNA from perforated cells. The reasoning is that mitochondria are larger\n", "than individual transcript molecules and less likely to escape through\n", "tears in the cell membrane.\n", "\n", "First, let Scanpy calculate some general qc-stats for genes and cells\n", "with the function `sc.pp.calculate_qc_metrics`, similar to\n", "`calculateQCmetrics()` in Scater. It can also calculate proportion of\n", "counts for specific gene populations, so first we need to define which\n", "genes are mitochondrial, ribosomal and hemoglobin." ] }, { "cell_type": "code", "metadata": {}, "source": [ "# mitochondrial genes\n", "adata.var['mt'] = adata.var_names.str.startswith('MT-') \n", "# ribosomal genes\n", "adata.var['ribo'] = adata.var_names.str.startswith((\"RPS\",\"RPL\"))\n", "# hemoglobin genes.\n", "adata.var['hb'] = adata.var_names.str.contains((\"^HB[^(P|E|S)]\"))\n", "\n", "adata.var" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": {}, "source": [ "sc.pp.calculate_qc_metrics(adata, qc_vars=['mt','ribo','hb'], percent_top=None, log1p=False, inplace=True)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now you can see that we have additional data in the metadata slot.\n", "\n", "Another opition to using the `calculate_qc_metrics` function is to\n", "calculate the values on your own and add to a metadata slot. An example\n", "for mito genes can be found below:" ] }, { "cell_type": "code", "metadata": {}, "source": [ "mito_genes = adata.var_names.str.startswith('MT-')\n", "# for each cell compute fraction of counts in mito genes vs. all genes\n", "# the `.A1` is only necessary as X is sparse (to transform to a dense array after summing)\n", "adata.obs['percent_mt2'] = np.sum(\n", " adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1\n", "# add the total counts per cell as observations-annotation to adata\n", "adata.obs['n_counts'] = adata.X.sum(axis=1).A1\n", "\n", "adata" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot QC\n", "\n", "Now we can plot some of the QC variables as violin plots." ] }, { "cell_type": "code", "metadata": {}, "source": [ "sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt', 'pct_counts_ribo', 'pct_counts_hb'], jitter=0.4, groupby = 'sample', rotation= 45)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, there is quite some difference in quality for these\n", "samples, with for instance the covid_15 and covid_16 samples having\n", "cells with fewer detected genes and more mitochondrial content. As the\n", "ribosomal proteins are highly expressed they will make up a larger\n", "proportion of the transcriptional landscape when fewer of the lowly\n", "expressed genes are detected. We can also plot the different QC-measures\n", "as scatter plots." ] }, { "cell_type": "code", "metadata": { "fig-height": 5, "fig-width": 5 }, "source": [ "sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', color=\"sample\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "> **Discuss**\n", ">\n", "> Plot additional QC stats that we have calculated as scatter plots. How\n", "> are the different measures correlated? Can you explain why?\n", "\n", "
\n", "\n", "## Filtering\n", "\n", "### Detection-based filtering\n", "\n", "A standard approach is to filter cells with low number of reads as well\n", "as genes that are present in at least a given number of cells. Here we\n", "will only consider cells with at least 200 detected genes and genes need\n", "to be expressed in at least 3 cells. Please note that those values are\n", "highly dependent on the library preparation method used." ] }, { "cell_type": "code", "metadata": {}, "source": [ "sc.pp.filter_cells(adata, min_genes=200)\n", "sc.pp.filter_genes(adata, min_cells=3)\n", "\n", "print(adata.n_obs, adata.n_vars)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extremely high number of detected genes could indicate doublets.\n", "However, depending on the cell type composition in your sample, you may\n", "have cells with higher number of genes (and also higher counts) from one\n", "cell type. In this case, we will run doublet prediction further down, so\n", "we will skip this step now, but the code below is an example of how it\n", "can be run:" ] }, { "cell_type": "code", "metadata": {}, "source": [ "# skip for now as we are doing doublet prediction\n", "#keep_v2 = (adata.obs['n_genes_by_counts'] < 2000) & (adata.obs['n_genes_by_counts'] > 500) & (adata.obs['lib_prep'] == 'v2')\n", "#print(sum(keep_v2))\n", "\n", "# filter for gene detection for v3\n", "#keep_v3 = (adata.obs['n_genes_by_counts'] < 4100) & (adata.obs['n_genes_by_counts'] > 1000) & (adata.obs['lib_prep'] != 'v2')\n", "#print(sum(keep_v3))\n", "\n", "# keep both sets of cells\n", "#keep = (keep_v2) | (keep_v3)\n", "#print(sum(keep))\n", "#adata = adata[keep, :]\n", "\n", "#print(\"Remaining cells %d\"%adata.n_obs)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally, we can also see which genes contribute the most to such\n", "reads. We can for instance plot the percentage of counts per gene." ] }, { "cell_type": "code", "metadata": { "fig-height": 6, "fig-width": 6 }, "source": [ "sc.pl.highest_expr_genes(adata, n_top=20)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, MALAT1 constitutes up to 30% of the UMIs from a single\n", "cell and the other top genes are mitochondrial and ribosomal genes. It\n", "is quite common that nuclear lincRNAs have correlation with quality and\n", "mitochondrial reads, so high detection of MALAT1 may be a technical\n", "issue. Let us assemble some information about such genes, which are\n", "important for quality control and downstream filtering.\n", "\n", "### Mito/Ribo filtering\n", "\n", "We also have quite a lot of cells with high proportion of mitochondrial\n", "and low proportion of ribosomal reads. It would be wise to remove those\n", "cells, if we have enough cells left after filtering. Another option\n", "would be to either remove all mitochondrial reads from the dataset and\n", "hope that the remaining genes still have enough biological signal. A\n", "third option would be to just regress out the `percent_mito` variable\n", "during scaling. In this case we had as much as 99.7% mitochondrial reads\n", "in some of the cells, so it is quite unlikely that there is much cell\n", "type signature left in those. Looking at the plots, make reasonable\n", "decisions on where to draw the cutoff. In this case, the bulk of the\n", "cells are below 20% mitochondrial reads and that will be used as a\n", "cutoff. We will also remove cells with less than 5% ribosomal reads." ] }, { "cell_type": "code", "metadata": {}, "source": [ "# filter for percent mito\n", "adata = adata[adata.obs['pct_counts_mt'] < 20, :]\n", "\n", "# filter for percent ribo > 0.05\n", "adata = adata[adata.obs['pct_counts_ribo'] > 5, :]\n", "\n", "print(\"Remaining cells %d\"%adata.n_obs)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, a large proportion of sample covid_15 is filtered out.\n", "Also, there is still quite a lot of variation in `percent_mito`, so it\n", "will have to be dealt with in the data analysis step. We can also notice\n", "that the `percent_ribo` are also highly variable, but that is expected\n", "since different cell types have different proportions of ribosomal\n", "content, according to their function.\n", "\n", "### Plot filtered QC\n", "\n", "Lets plot the same QC-stats once more." ] }, { "cell_type": "code", "metadata": {}, "source": [ "sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt','pct_counts_ribo', 'pct_counts_hb'], jitter=0.4, groupby = 'sample', rotation = 45)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Filter genes\n", "\n", "As the level of expression of mitochondrial and MALAT1 genes are judged\n", "as mainly technical, it can be wise to remove them from the dataset\n", "before any further analysis. In this case we will also remove the HB\n", "genes." ] }, { "cell_type": "code", "metadata": {}, "source": [ "malat1 = adata.var_names.str.startswith('MALAT1')\n", "# we need to redefine the mito_genes since they were first \n", "# calculated on the full object before removing low expressed genes.\n", "mito_genes = adata.var_names.str.startswith('MT-')\n", "hb_genes = adata.var_names.str.contains('^HB[^(P|E|S)]')\n", "\n", "remove = np.add(mito_genes, malat1)\n", "remove = np.add(remove, hb_genes)\n", "keep = np.invert(remove)\n", "\n", "adata = adata[:,keep]\n", "\n", "print(adata.n_obs, adata.n_vars)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sample sex\n", "\n", "When working with human or animal samples, you should ideally constrain\n", "your experiments to a single sex to avoid including sex bias in the\n", "conclusions. However this may not always be possible. By looking at\n", "reads from chromosomeY (males) and XIST (X-inactive specific transcript)\n", "expression (mainly female) it is quite easy to determine per sample\n", "which sex it is. It can also be a good way to detect if there has been\n", "any mislabelling in which case, the sample metadata sex does not agree\n", "with the computational predictions.\n", "\n", "To get choromosome information for all genes, you should ideally parse\n", "the information from the gtf file that you used in the mapping pipeline\n", "as it has the exact same annotation version/gene naming. However, it may\n", "not always be available, as in this case where we have downloaded public\n", "data. Hence, we will use biomart to fetch chromosome information." ] }, { "cell_type": "code", "metadata": {}, "source": [ "# requires pybiomart\n", "if not fetch_annotation:\n", " annot = sc.queries.biomart_annotations(\"hsapiens\", [\"ensembl_gene_id\", \"external_gene_name\", \"start_position\", \"end_position\", \"chromosome_name\"], ).set_index(\"external_gene_name\")\n", " # adata.var[annot.columns] = annot" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have the chromosome information, we can calculate the\n", "proportion of reads that comes from chromosome Y per cell." ] }, { "cell_type": "code", "metadata": {}, "source": [ "chrY_genes = adata.var_names.intersection(annot.index[annot.chromosome_name == \"Y\"])\n", "chrY_genes\n", "\n", "adata.obs['percent_chrY'] = np.sum(\n", " adata[:, chrY_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1 * 100" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then plot XIST expression vs chrY proportion. As you can see, the\n", "samples are clearly on either side, even if some cells do not have\n", "detection of either." ] }, { "cell_type": "code", "metadata": { "fig-height": 5, "fig-width": 5 }, "source": [ "# color inputs must be from either .obs or .var, so add in XIST expression to obs.\n", "adata.obs[\"XIST-counts\"] = adata.X[:,adata.var_names.str.match('XIST')].toarray()\n", "\n", "sc.pl.scatter(adata, x='XIST-counts', y='percent_chrY', color=\"sample\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot as violins." ] }, { "cell_type": "code", "metadata": { "fig-height": 5, "fig-width": 10 }, "source": [ "sc.pl.violin(adata, [\"XIST-counts\", \"percent_chrY\"], jitter=0.4, groupby = 'sample', rotation= 45)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we can see clearly that we have three males and five females, can\n", "you see which samples they are? Do you think this will cause any\n", "problems for downstream analysis? Discuss with your group: what would be\n", "the best way to deal with this type of sex bias?\n", "\n", "## Cell cycle state\n", "\n", "We here perform cell cycle scoring. To score a gene list, the algorithm\n", "calculates the difference of mean expression of the given list and the\n", "mean expression of reference genes. To build the reference, the function\n", "randomly chooses a bunch of genes matching the distribution of the\n", "expression of the given list. Cell cycle scoring adds three slots in the\n", "metadata, a score for S phase, a score for G2M phase and the predicted\n", "cell cycle phase.\n", "\n", "First read the file with cell cycle genes, from Regev lab and split into\n", "S and G2M phase genes. We first download the file." ] }, { "cell_type": "code", "metadata": {}, "source": [ "path_file = os.path.join(path_results, 'regev_lab_cell_cycle_genes.txt')\n", "if not os.path.exists(path_file):\n", " urllib.request.urlretrieve(os.path.join(path_data, 'regev_lab_cell_cycle_genes.txt'), path_file)" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": {}, "source": [ "cell_cycle_genes = [x.strip() for x in open('./data/covid/results/regev_lab_cell_cycle_genes.txt')]\n", "print(len(cell_cycle_genes))\n", "\n", "# Split into 2 lists\n", "s_genes = cell_cycle_genes[:43]\n", "g2m_genes = cell_cycle_genes[43:]\n", "\n", "cell_cycle_genes = [x for x in cell_cycle_genes if x in adata.var_names]\n", "print(len(cell_cycle_genes))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before running cell cycle we have to normalize the data. In the scanpy\n", "object, the data slot will be overwritten with the normalized data. So\n", "first, save the raw data into the slot `raw`. Then run normalization,\n", "log transformation and scale the data." ] }, { "cell_type": "code", "metadata": {}, "source": [ "# save normalized counts in raw slot.\n", "adata.raw = adata\n", "\n", "# normalize to depth 10 000\n", "sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)\n", "\n", "# logaritmize\n", "sc.pp.log1p(adata)\n", "\n", "# scale\n", "sc.pp.scale(adata)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We here perform cell cycle scoring. The function is actually a wrapper\n", "to sc.tl.score_gene_list, which is launched twice, to score separately S\n", "and G2M phases. Both sc.tl.score_gene_list and\n", "sc.tl.score_cell_cycle_genes are a port from Seurat and are supposed to\n", "work in a very similar way. To score a gene list, the algorithm\n", "calculates the difference of mean expression of the given list and the\n", "mean expression of reference genes. To build the reference, the function\n", "randomly chooses a bunch of genes matching the distribution of the\n", "expression of the given list. Cell cycle scoring adds three slots in\n", "data, a score for S phase, a score for G2M phase and the predicted cell\n", "cycle phase." ] }, { "cell_type": "code", "metadata": {}, "source": [ "sc.tl.score_genes_cell_cycle(adata, s_genes=s_genes, g2m_genes=g2m_genes)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now create a violin plot for the cell cycle scores as well." ] }, { "cell_type": "code", "metadata": { "fig-height": 5, "fig-width": 10 }, "source": [ "sc.pl.violin(adata, ['S_score', 'G2M_score'], jitter=0.4, groupby = 'sample', rotation=45)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case it looks like we only have a few cycling cells in these\n", "datasets.\n", "\n", "Scanpy does an automatic prediction of cell cycle phase with a default\n", "cutoff of the scores at zero. As you can see this does not fit this data\n", "very well, so be cautios with using these predictions. Instead we\n", "suggest that you look at the scores." ] }, { "cell_type": "code", "metadata": { "fig-height": 8, "fig-width": 10 }, "source": [ "sc.pl.scatter(adata, x='S_score', y='G2M_score', color=\"phase\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predict doublets\n", "\n", "Doublets/Multiples of cells in the same well/droplet is a common issue\n", "in scRNAseq protocols. Especially in droplet-based methods with\n", "overloading of cells. In a typical 10x experiment the proportion of\n", "doublets is linearly dependent on the amount of loaded cells. As\n", "indicated from the Chromium user guide, doublet rates are about as\n", "follows:\\\n", "![](../figs/10x_doublet_rate.png)\\\n", "Most doublet detectors simulates doublets by merging cell counts and\n", "predicts doublets as cells that have similar embeddings as the simulated\n", "doublets. Most such packages need an assumption about the\n", "number/proportion of expected doublets in the dataset. The data you are\n", "using is subsampled, but the original datasets contained about 5 000\n", "cells per sample, hence we can assume that they loaded about 9 000 cells\n", "and should have a doublet rate at about 4%.\n", "\n", "For doublet detection, we will use the package `Scrublet`, so first we\n", "need to get the raw counts from `adata.raw.X` and run scrublet with that\n", "matrix. Then we add in the doublet prediction info into our anndata\n", "object.\n", "\n", "Doublet prediction should be run for each dataset separately, so first\n", "we need to split the adata object into 6 separate objects, one per\n", "sample and then run scrublet on each of them." ] }, { "cell_type": "code", "metadata": {}, "source": [ "import scrublet as scr\n", "\n", "# split per batch into new objects.\n", "batches = adata.obs['sample'].cat.categories.tolist()\n", "alldata = {}\n", "for batch in batches:\n", " tmp = adata[adata.obs['sample'] == batch,]\n", " print(batch, \":\", tmp.shape[0], \" cells\")\n", " scrub = scr.Scrublet(tmp.raw.X)\n", " out = scrub.scrub_doublets(verbose=False, n_prin_comps = 20)\n", " alldata[batch] = pd.DataFrame({'doublet_score':out[0],'predicted_doublets':out[1]},index = tmp.obs.index)\n", " print(alldata[batch].predicted_doublets.sum(), \" predicted_doublets\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": {}, "source": [ "# add predictions to the adata object.\n", "scrub_pred = pd.concat(alldata.values())\n", "adata.obs['doublet_scores'] = scrub_pred['doublet_score'] \n", "adata.obs['predicted_doublets'] = scrub_pred['predicted_doublets'] \n", "\n", "sum(adata.obs['predicted_doublets'])" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We should expect that two cells have more detected genes than a single\n", "cell, lets check if our predicted doublets also have more detected genes\n", "in general." ] }, { "cell_type": "code", "metadata": { "fig-height": 5, "fig-width": 5 }, "source": [ "# add in column with singlet/doublet instead of True/Fals\n", "%matplotlib inline\n", "\n", "adata.obs['doublet_info'] = adata.obs[\"predicted_doublets\"].astype(str)\n", "sc.pl.violin(adata, 'n_genes_by_counts', jitter=0.4, groupby = 'doublet_info', rotation=45)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, lets run PCA and UMAP and plot doublet scores onto UMAP to check\n", "the doublet predictions." ] }, { "cell_type": "code", "metadata": { "fig-height": 4, "fig-width": 12 }, "source": [ "sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)\n", "adata = adata[:, adata.var.highly_variable]\n", "sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])\n", "sc.pp.scale(adata, max_value=10)\n", "sc.tl.pca(adata, svd_solver='arpack')\n", "sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)\n", "sc.tl.umap(adata)\n", "sc.pl.umap(adata, color=['doublet_scores','doublet_info','sample'])" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, lets remove all predicted doublets from our data." ] }, { "cell_type": "code", "metadata": {}, "source": [ "# also revert back to the raw counts as the main matrix in adata\n", "adata = adata.raw.to_adata() \n", "\n", "adata = adata[adata.obs['doublet_info'] == 'False',:]\n", "print(adata.shape)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Save data\n", "\n", "Finally, lets save the QC-filtered data for further analysis. Create\n", "output directory `data/covid/results` and save data to that folder. This\n", "will be used in downstream labs." ] }, { "cell_type": "code", "metadata": {}, "source": [ "adata.write_h5ad('data/covid/results/scanpy_covid_qc.h5ad')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Session info\n", "\n", "```{=html}\n", "
\n", "```\n", "```{=html}\n", "\n", "```\n", "Click here\n", "```{=html}\n", "\n", "```" ] }, { "cell_type": "code", "metadata": {}, "source": [ "sc.logging.print_versions()" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{=html}\n", "
\n", "```" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 4 }