<jupyter_start><jupyter_text>Welcome to Week 3 (what you've all been waiting for), Single Cell RNA! In this week, we're going to focus on the basics of single cell RNA analysis, including dimensionality reduction and opitional doublet detection and cell type identification/cell surface protein integration.First, we'll have an introduction to the main Python package used for single cell analysis - **Scanpy**, developed by <NAME> lab at the Helmholtz Zentrum München. The main publication for Scanpy is here, if you're curious:* https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1382-0?optIn=trueWe're going to start with the output of the 10X Genomics cellranger pipeline, which is used for processing of single cell RNA data generated using droplet-based capture. You can read more about how single cell RNA sequencing and the cellranger pipeline work here (and we may do a cellranger tutorial later if there is interest:  * https://www.10xgenomics.com/products/single-cell/?src=search&lss=google&cnm=sem-goog-2020-website-page-ra_g-p_brand-amr-retarget&cid=7011P000000oWys (videos at the end)* https://support.10xgenomics.com/single-cell-gene-expression/software/overview/welcome  Note though that you can also start with many other formats: the simplest being just a cells x genes counts table, but many other formats as well.  One thing that has become increasingly common is that there's a few labs that have really pushed the bounds of single-cell analysis with many, many cells - hundred thousands to millions of cells.  In many cases, these groups (e.g. Sten Linnarson, Allen Brain Atlas) have introduced their own file formats for dealing with such large datasets efficiently (.loom, .tome, respectively), and there's specific functions for importing other dataset formats:* https://icb-scanpy.readthedocs-hosted.com/en/stable/api/index.htmlreadingThe main idea is that single cells are captured in oil-water emulsion droplets along with barcoded oligos, which tag RNA transcripts with a unique cellular barcode during reverse transcription. After sequencing, the cellranger pipeline aligns the transcripts and matches them back to the cellular barcode to obtain counts for each transcript that was expressed for each cell. The main output that we will use from the cellranger pipeline is the **feature/cell matrix**, which is a matrix of cells by genes, where the values are the number of transcripts that were detected for that cell.  Single cell methods are not perfect, and one of the main things we have to deal with is **dropout**, which occurs when a transcript is expressed by a cell but is not recovered by single cell sequencing. Lowly expressed genes are more likely to be affected by dropout. For 10X datasets, typically we can expect to detect 1000-2000 genes per cell. The main approach that has been developed to deal with dropout during single cell analysis is dimensionality reduction. The main steps we will take in Seurat are:1. **Filter** cells based on QC metrics (number of genes detected, percent mitochondrial genes, etc.)2. **Normalize** data3. Identify **variable** genes to use for principal component analysis (PCA)4. **Scale** data and regress out unwanted sources of variation (i.e. cell cycle)5. **PCA** dimensionality reduction6. **Cluster** cells to find nearest neighbors in PCA space (typically corresponding to different cell types/states)7. Visualize cells using **UMAP/TSNE** non-linear dimensionality reduction8. Find **differentially expressed** genes between clusters9. **Cluster annotation**  There is obviously a lot we can do after this but these steps are common to any single cell project. To get started we will be using this publically available human PBMC dataset from 10X:* https://support.10xgenomics.com/single-cell-vdj/datasets/3.0.0/vdj_v1_hs_pbmc2_5gex_protein  Go ahead and download the **Feature / cell matrix (filtered)**. Unpack the tar.gz file either using Finder or with this command in terminal:  `tar -xzf vdj_v1_hs_pbmc2_5gex_protein_filtered_feature_bc_matrix.tar.gz`  You can also take a look at the **Summary HTML** which has the QC metrics from the cellranger pipeline. We will now be following the Scanpy Guided Clustering Tutorial found here (but using the larger PBMC dataset we just downloaded):* https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html  Installing PackagesWe're going to be using **scanpy** to do the scRNA analysis.  You'll need to install scanpy (https://icb-scanpy.readthedocs-hosted.com/en/stable/installation.html) by running (in terminal):       conda install -c bioconda scanpyYou'll also want to make sure that AnnData is installed, which is the package that contains the AnnData datatype (https://icb-anndata.readthedocs-hosted.com/en/stable/anndata.AnnData.html).  This is the datatype that we use to store single-cell data. In addition, you'll also want to install the Leiden clustering package as well:        conda install -c conda-forge leidenalg anndata        If you want to try Scrublet (doublet detection, https://github.com/AllonKleinLab/scrublet), then run:        conda install -c bioconda scrublet  What if you get stuck?If you get stuck - you don't know what the function you should use is, or what arguments/paramters it has, the first places you should look are the tutorials and documentation/API for scanpy:* https://scanpy.readthedocs.io/en/latest/api/index.html* https://scanpy.readthedocs.io/en/latest/tutorials.htmlIn the API, you can find all of the information on what tools are available, how to use them, and   A note on Scanpy vs. SeuratScanpy is probably less common than Seurat at the moment, but I think this is mostly because Seurat came first.  Scanpy has been gaining more and more usage, and I've started to see it used a little more frequently (you can sometimes tell who is using what based on how their figures look, just with default plotting formats), in addition to noticing a lot of web-based analysis tools starting to use Scanpy.  I think it's better at dealing with large datasets, and as people are starting to develop more and more tools in Python, I think it's only going to become more commonly used.**I (Kevin) have mostly used Seurat, but the more that I've used Scanpy the more that I've been impressed.**  Pros:- Scanpy runs significanlty faster, especially with larger datasets- Scanpy has implemented a lot of other useful tools (like pseudotime trajectories) and integrates with a lot of other tools (like RNA velocity)- The built-in plotting functions (https://scanpy-tutorials.readthedocs.io/en/latest/visualizing-marker-genes.html) are very impressive and I think blow Seurat's out of the water in terms of options and ease of use.- I've also really been impressed with the documentation for all of the functions and arguments, which I think is one thing that Seurat is often lacking - their vignettes are nice, but their general documentation is often pretty poor. So while I still think Seurat is the primary tool for scRNA analysis just because it's been around longer, I think Scanpy is gaining a lot of momentum and I feel really good about recommending it, and am probably going to start using it more myself.  Import dataNote that this data contains both *RNA* and *CITE-seq* information.  We're going to start with just the RNA data.<jupyter_code>import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
path = '/Users/kevin/changlab/covid19/3_scRNA/data/filtered_feature_bc_matrix/'
figpath = '/Users/kevin/changlab/covid19/3_scRNA/figures/'

# lets set the default figure settings

sc.settings.set_figure_params(dpi_save=300, dpi=300)
sc.settings.figdir = figpath
# If we wanted to import all of the data, then we would run
orig = sc.read_10x_mtx(path, gex_only=False)
print(orig)

gex_genes = orig.var.index[orig.var['feature_types'] == 'Gene Expression']
cite_genes = orig.var.index[orig.var['feature_types'] == 'Antibody Capture']

print('gex', len(gex_genes))
print('cite', len(cite_genes))
# Since we are only going to deal with the gene expression data for now, we'll say

gex = sc.read_10x_mtx(path, gex_only=True)
print(gex)<jupyter_output>View of AnnData object with n_obs × n_vars = 8258 × 33538 
    var: 'gene_ids', 'feature_types'<jupyter_text>Just to match the Seurat filtering (only including genes that are in at least 3 cells; and only including cells that have at least 200 genes), we'll say:<jupyter_code>sc.pp.filter_cells(gex, min_genes=200)
sc.pp.filter_genes(gex, min_cells=3)

print(gex)<jupyter_output>Trying to set attribute `.obs` of view, making a copy.<jupyter_text>The rest of the tutorial I will leave for you to fill in!  There are a couple of important pivot points that you can modify for your data which I will break up by processing step:1. **Filter**     a. `nFeature_RNA`: How many genes are detected per cell. I typically use 200 as a minimum but check the distribution in your data. Cells with fewer than 200 genes detected typically don't have enough information to cluster on.     b. `nCount_RNA`: How many UMIs or transcripts detected per cell. I typically don't filter using this metric as it is linearly correlated with nFeature_RNA. Some people will use an upper threshold here to try and filter out doublets but typically is not robust and we describe alternative approaches below.     c. `percent.mt`: Percent mitochondrial transcripts. Dying cells will have more mitochondrial gene content and so this can be used to filter low quality cells. Check the distribution in your data but typically upper thresholds of 5-10% mitochondrial content make sense.2. **Normalize** - I use the defaults pretty much always.3. **Variable genes** - Not a lot to modify here. If you have residual batch/unwanted variation (like cell cycle) effects after regression you can remove these genes from the variable gene set. Sometimes I also observed clustering based on variable immune genes (like TCR/Ig transcripts) so these can be removed as well.4. **Scale**     a. `vars.to.regress`: If you have batch effects/sources of unwanted variation like cell cycle effects you can regress them out here. Some other common variables to regress are nCount_RNA, percent.mt (if you observe clustering based on these features & not due to remaining low quality cells).5. **PCA**     a. `dims`: Pick how many dimensions to use for PCA. I typically use ElbowPlot to decide. Including more PCs is usually better than too few, so I would err past the elbow. Between 10-40 PCs is typically appropriate.6. **Cluster**     a. `resolution`: Argument to the FindClusters function. This is a big one you can change. Higher resolution values will give more clusters, lower resolution values give fewer clusters. How many clusters you want depends on the biological variation in your data as well as how you want to analyze it (are broad cell types enough or do you want specific cell states?). Values between 0.4-1.2 are typically appropriate. 7. **UMAP** - UMAP is superior to TSNE in many ways so I won't talk about TSNE.     a. `min.dist`: This will change how spread out the UMAP is with lower values pulling cells apart more. This is useful for closely related cells that may not seperate well in the UMAP using the default parameters.  Make QC violin plots and filter cells  Normalize the data  Find variable genes  Optional: regress out any confounding variables  Scale your data  Perform PCA  Cluster your data, get a UMAP embedding, and plot the results  Plot some common marker genes<jupyter_code>genes_to_plot = ['CD3G','CD8A','CD4',
                 'NKG7','MS4A1','CD14',
                 'FCGR3A','FCER1A','PPBP']<jupyter_output><empty_output><jupyter_text>Additional steps: Finding marker genes, annotating clusters, and making plots.Congrats on finishing the Scanpy tutorial! There are a few things you can do now to explore the data:1. Check `n_genes`, `n_counts` amd `percent_mito` for each cluster using the `pl.violin` function. Do any clusters look like low quality cells?2. Check for cell cycle effects by checking expression of cell cycle associated genes using the `pl.violin` function or looking for cell cycle genes in the cluster markers called by `tl.rank_genes_groups`. If cell cycle effects are present they can be regressed out using this vignette:  https://nbviewer.jupyter.org/github/theislab/scanpy_usage/blob/master/180209_cell_cycle/cell_cycle.ipynb3. Annotate cell clusters based on expression of marker genes. You can check expression using either the `pl.rank_genes_groups`, `pl.rank_genes_groups_stacked_violin`, `pl.stacked_violin`, `pl.umap`, and other functions. All are useful, just depends on your preference! Here are couple of good markers for PBMCs:     a. T cells: CD3G, CD3E, CD3D     b. Cytotoxic CD8+ T cells: CD8A     c. Helper CD4+ T cells: CD4     d. Naive T cells: CCR7     e. Memory CD8+ cells: EOMES     f. Effector memory CD8+ T cells: KLRD1     g. B cells: CD19, MS4A1     h. Monocytes: CD14, FCGR3A     i. NK cells: GNLY, NKG7 (CD3-, NK cells and cytotoxic T cells share many markers and may cluster together)     j. Platelets: PPBP4. Subset a broad cell type and recluster to find more accurate subtypes. For example, subset all T cell clusters and recluster starting from data scaling and variable gene detection. Do you find more clusters? How do they compare to the original clusters you found using all PBMCs?5. Make sure you save your Scanpy object (`gex`) so it can be loaded again using the .write() function built in to the anndata (gex) object. Also try saving some plots as PDFs or PNGs and tables (like the cluster markers) as text files.  The tutorial has some steps on this (https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html).6. Play with different visualizations! See here (https://scanpy-tutorials.readthedocs.io/en/latest/visualizing-marker-genes.html) for some ideas.**Extra challenge**: This dataset also has CITE-seq surface protein data! Try adding this data to your Seurat object (Hint look back at the `sc.read_10x_mtx` we read in at the beginning). Which proteins were analyzed? How does protein expression compare to RNA expression? See more details in the vignette here:  https://icb-scanpy-tutorials.readthedocs-hosted.com/en/multiomics/cite-seq/pbmc5k.html<jupyter_code># since we know markers, let's make panels for some cell types that we expect to see

t_cell_panel = ['CD3E','CD8A','CD4','CCR7','EOMES','KLRD1']
b_cell_panel = ['CD79A','CD19','MS4A1']
monocyte_panel = ['CD14','LYZ','FCGR3A','MS4A7']
nk_panel = ['GNLY','NKG7']
dc_panel = ['FCER1A','CST3']
platelet_panel = ['PPBP']<jupyter_output><empty_output><jupyter_text>Now let's deal with doublets, and also try regressing out some confounding variables.See https://github.com/AllonKleinLab/scrublet, https://github.com/AllonKleinLab/scrublet/blob/master/examples/scrublet_basics.ipynb for examples.<jupyter_code>import scrublet as scr
import scipy.io

counts_matrix = scipy.io.mmread(path + '/matrix.mtx.gz').T.tocsc()
cells = pd.read_csv(path + '/barcodes.tsv.gz', sep='\t', header=None, names=['barcode'])
# cells['barcode'] = [i.split('-')[0] for i in cells['barcode']]
cells = cells.set_index('barcode', drop=False)
cells.head()
scrub = scr.Scrublet(counts_matrix, expected_doublet_rate=0.08)
doublet_scores, predicted_doublets = scrub.scrub_doublets(min_counts=2, 
                                                          min_cells=3, 
                                                          min_gene_variability_pctl=85, 
                                                          n_prin_comps=30)
fig, _ = scrub.plot_histogram()

fig.savefig(figpath + 'scrublet_histogram.png', dpi=300, bbox_inches='tight')
predicted_doublets = scrub.call_doublets(threshold=0.25)
scrub.set_embedding('UMAP', scr.get_umap(scrub.manifold_obs_, 10, min_dist=0.3))
fig, _ = scrub.plot_embedding('UMAP', order_points=True)

fig.savefig(figpath + 'scrublet_umap.png', dpi=300, bbox_inches='tight')
cells['doublet_score'] = doublet_scores
cells['predicted_doublet'] = predicted_doublets
cells.to_csv(figpath + 'scrublet_results.txt')
cells.head()
# lets restart and add this to our data.

# Since we are only going to deal with the gene expression data for now, we'll say

gex = sc.read_10x_mtx(path, gex_only=True)
print(gex)

gex.obs['doublet_score'] = cells.loc[gex.obs.index, 'doublet_score']
gex.obs['predicted_doublet'] = cells.loc[gex.obs.index, 'predicted_doublet']

gex.obs.head()<jupyter_output>Trying to set attribute `.obs` of view, making a copy.
{'fraction_non_alphanumeric': 0.07411204245311463, 'fraction_numerical': 0.013266598323584393, 'mean_word_length': 4.55950385517935, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 19, 'https://': 18, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '196092', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4960, 'n_tokens_neox': 4542, 'n_words': 2132, 'path': 'wk4_scRNA/py_wk4_scRNA_scanpy.ipynb', 'repo_name': 'krparker/Bioinformatics-Tutorials', 'size': '16591'}