diff --git a/docs/06-data-cleaning.md b/docs/06-data-cleaning.md
index e4af639d..812539ca 100644
--- a/docs/06-data-cleaning.md
+++ b/docs/06-data-cleaning.md
@@ -198,8 +198,8 @@ dec.zeisel.qc[c(1,500,1000,2000),1:6]
 ##       <numeric> <numeric> <numeric> <numeric>   <numeric>   <numeric>
 ## Plp1    3.97379  16.80426  0.275699 16.528564 0.00000e+00 0.00000e+00
 ## Shfm1   2.64467   1.55363  0.534672  1.018959 5.54040e-47 2.04392e-45
-## Cdc42   3.34431   1.04001  0.385442  0.654563 9.97095e-38 2.96240e-36
-## Vdac2   1.81697   1.08138  0.720031  0.361348 7.90102e-05 2.55550e-04
+## Cdc42   3.34431   1.04001  0.385442  0.654563 9.97095e-38 2.96239e-36
+## Vdac2   1.81697   1.08138  0.720031  0.361348 7.90101e-05 2.55550e-04
 ```
 
 In this dataframe, the `bio` column represents the excess variation in gene expression (that is, the difference between the observed and expected expression). We are looking at four different genes - the gene with the greatest excess variation, as well as the last genes included in the top 500, top 1000, and top 2000 genes.
diff --git a/docs/07-clustering_files/figure-html/unnamed-chunk-3-1.png b/docs/07-clustering_files/figure-html/unnamed-chunk-3-1.png
index 0ddcd750..19fbe668 100644
Binary files a/docs/07-clustering_files/figure-html/unnamed-chunk-3-1.png and b/docs/07-clustering_files/figure-html/unnamed-chunk-3-1.png differ
diff --git a/docs/07-clustering_files/figure-html/unnamed-chunk-3-2.png b/docs/07-clustering_files/figure-html/unnamed-chunk-3-2.png
index 52fa8a34..bb2cd29d 100644
Binary files a/docs/07-clustering_files/figure-html/unnamed-chunk-3-2.png and b/docs/07-clustering_files/figure-html/unnamed-chunk-3-2.png differ
diff --git a/docs/07-clustering_files/figure-html/unnamed-chunk-3-3.png b/docs/07-clustering_files/figure-html/unnamed-chunk-3-3.png
index 8048872a..34a7f5de 100644
Binary files a/docs/07-clustering_files/figure-html/unnamed-chunk-3-3.png and b/docs/07-clustering_files/figure-html/unnamed-chunk-3-3.png differ
diff --git a/docs/07-clustering_files/figure-html/unnamed-chunk-5-1.png b/docs/07-clustering_files/figure-html/unnamed-chunk-5-1.png
index 52fa8a34..bb2cd29d 100644
Binary files a/docs/07-clustering_files/figure-html/unnamed-chunk-5-1.png and b/docs/07-clustering_files/figure-html/unnamed-chunk-5-1.png differ
diff --git a/docs/07-clustering_files/figure-html/unnamed-chunk-7-1.png b/docs/07-clustering_files/figure-html/unnamed-chunk-7-1.png
index e8c311a9..ba1f5dcc 100644
Binary files a/docs/07-clustering_files/figure-html/unnamed-chunk-7-1.png and b/docs/07-clustering_files/figure-html/unnamed-chunk-7-1.png differ
diff --git a/docs/404.html b/docs/404.html
index ca709a4c..30a0c8dd 100644
--- a/docs/404.html
+++ b/docs/404.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/About.md b/docs/About.md
index 6002daaa..8cdb30bf 100644
--- a/docs/About.md
+++ b/docs/About.md
@@ -42,7 +42,7 @@ These credits are based on our [course contributors table guidelines](https://gi
 ##  collate  en_US.UTF-8
 ##  ctype    en_US.UTF-8
 ##  tz       Etc/UTC
-##  date     2023-03-22
+##  date     2023-05-25
 ##  pandoc   2.17.1.1 @ /usr/local/bin/ (via rmarkdown)
 ## 
 ## ─ Packages ───────────────────────────────────────────────────────────────────
diff --git a/docs/about-the-authors.html b/docs/about-the-authors.html
index e1567074..ba67e8dd 100644
--- a/docs/about-the-authors.html
+++ b/docs/about-the-authors.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
@@ -311,7 +321,7 @@ <h1>About the Authors</h1>
 ##  collate  en_US.UTF-8
 ##  ctype    en_US.UTF-8
 ##  tz       Etc/UTC
-##  date     2023-03-22
+##  date     2023-05-25
 ##  pandoc   2.17.1.1 @ /usr/local/bin/ (via rmarkdown)
 ## 
 ## ─ Packages ───────────────────────────────────────────────────────────────────
diff --git a/docs/cell-type-annotation.html b/docs/cell-type-annotation.html
index f5c04880..b5a64954 100644
--- a/docs/cell-type-annotation.html
+++ b/docs/cell-type-annotation.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/choosing-features-for-the-analysis.html b/docs/choosing-features-for-the-analysis.html
index 79bec49b..cf4e46e3 100644
--- a/docs/choosing-features-for-the-analysis.html
+++ b/docs/choosing-features-for-the-analysis.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
@@ -258,8 +268,8 @@ <h2><span class="header-section-number">4.2</span> Selecting highly variable gen
 ##       &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt;   &lt;numeric&gt;   &lt;numeric&gt;
 ## Plp1    3.97379  16.80426  0.275699 16.528564 0.00000e+00 0.00000e+00
 ## Shfm1   2.64467   1.55363  0.534672  1.018959 5.54040e-47 2.04392e-45
-## Cdc42   3.34431   1.04001  0.385442  0.654563 9.97095e-38 2.96240e-36
-## Vdac2   1.81697   1.08138  0.720031  0.361348 7.90102e-05 2.55550e-04</code></pre>
+## Cdc42   3.34431   1.04001  0.385442  0.654563 9.97095e-38 2.96239e-36
+## Vdac2   1.81697   1.08138  0.720031  0.361348 7.90101e-05 2.55550e-04</code></pre>
 <p>In this dataframe, the <code>bio</code> column represents the excess variation in gene expression (that is, the difference between the observed and expected expression). We are looking at four different genes - the gene with the greatest excess variation, as well as the last genes included in the top 500, top 1000, and top 2000 genes.</p>
 <div class="reflection">
 <p>QUESTIONS</p>
diff --git a/docs/clustering.html b/docs/clustering.html
index bebd3bed..d81741af 100644
--- a/docs/clustering.html
+++ b/docs/clustering.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/dimensionality-reduction.html b/docs/dimensionality-reduction.html
index d32abc06..5001329a 100644
--- a/docs/dimensionality-reduction.html
+++ b/docs/dimensionality-reduction.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/index.html b/docs/index.html
index 7e19078d..d449db3a 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
@@ -223,7 +233,7 @@ <h1>
 </div>
 <div id="header">
 <h1 class="title">Statistics for Genomics: scRNAseq</h1>
-<p class="date"><em>March 22, 2023</em></p>
+<p class="date"><em>May 25, 2023</em></p>
 </div>
 <div id="about-this-book" class="section level1 unnumbered">
 <h1>About this Book</h1>
diff --git a/docs/index.md b/docs/index.md
index 1f057770..f48b8770 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,6 +1,6 @@
 ---
 title: "Statistics for Genomics: scRNAseq"
-date: "March 22, 2023"
+date: "May 25, 2023"
 site: bookdown::bookdown_site
 documentclass: book
 bibliography: book.bib
@@ -29,3 +29,4 @@ _Programming skills_
 This exercise is part of a collection of teaching resources developed through the *Genomic Data Science Community Network* (GDSCN). GDSCN works towards a vision where researchers, educators, and students from diverse backgrounds are able to fully participate in genomic data science research.  Learn more about GDSCN by visiting https://www.gdscn.org/home or reading the [article in Genome Research](https://doi.org/10.1101/gr.276496.121).
 
 Please check out our full collection of AnVIL and related resources: https://hutchdatascience.org/AnVIL_Collection/
+
diff --git a/docs/marker-gene-detection.html b/docs/marker-gene-detection.html
index 64fbb982..b54f209e 100644
--- a/docs/marker-gene-detection.html
+++ b/docs/marker-gene-detection.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/normalization.html b/docs/normalization.html
index bd862064..cf22d58c 100644
--- a/docs/normalization.html
+++ b/docs/normalization.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/overview.html b/docs/overview.html
index f530cbf4..dbd110a2 100644
--- a/docs/overview.html
+++ b/docs/overview.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/quality-control-of-the-library.html b/docs/quality-control-of-the-library.html
index 6930c357..d6039622 100644
--- a/docs/quality-control-of-the-library.html
+++ b/docs/quality-control-of-the-library.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/references.html b/docs/references.html
index 06da2403..3f60a6e5 100644
--- a/docs/references.html
+++ b/docs/references.html
@@ -23,7 +23,7 @@
 
 
 
-<meta name="date" content="2023-03-22" />
+<meta name="date" content="2023-05-25" />
 
   <meta name="viewport" content="width=device-width, initial-scale=1" />
   <meta name="apple-mobile-web-app-capable" content="yes" />
@@ -52,6 +52,16 @@
 
 <link href="libs/anchor-sections-1.0.1/anchor-sections.css" rel="stylesheet" />
 <script src="libs/anchor-sections-1.0.1/anchor-sections.js"></script>
+  
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=G-TWH8T48JSX"></script>
+<script>
+  window.dataLayer = window.dataLayer || [];
+  function gtag(){dataLayer.push(arguments);}
+  gtag('js', new Date());
+
+  gtag('config', 'G-TWH8T48JSX');
+</script>
 
 
 <style type="text/css">
diff --git a/docs/search_index.json b/docs/search_index.json
index b004e7c4..baeabeb7 100644
--- a/docs/search_index.json
+++ b/docs/search_index.json
@@ -1 +1 @@
-[["index.html", "Statistics for Genomics: scRNAseq About this Book Skills Level GDSCN Collection", " Statistics for Genomics: scRNAseq March 22, 2023 About this Book This book is intended to introduce concepts in single cell RNA-seq analysis (scRNA-seq). Skills Level Genetics Programming skills GDSCN Collection This exercise is part of a collection of teaching resources developed through the Genomic Data Science Community Network (GDSCN). GDSCN works towards a vision where researchers, educators, and students from diverse backgrounds are able to fully participate in genomic data science research. Learn more about GDSCN by visiting https://www.gdscn.org/home or reading the article in Genome Research. Please check out our full collection of AnVIL and related resources: https://hutchdatascience.org/AnVIL_Collection/ "],["overview.html", "Chapter 1 OVERVIEW 1.1 Scope 1.2 Bioconductor reference book 1.3 SingleCellExperiment 1.4 Obtaining Data", " Chapter 1 OVERVIEW 1.1 Scope There are many excellent programs available for scRNA-seq analysis. In this lesson, we will be using the scran and scater packages available from Bioconductor. We do not cover .fastq processing, nor do we explore the Seurat and scanpy pipelines. 1.2 Bioconductor reference book The lessons and activities are adapted from the Orchestrating Single-Cell Analysis with Bioconductor book by Robert Amezquita, Aaron Lun, Stephanie Hicks, and Raphael Gottardo. 1.3 SingleCellExperiment Bioconductor stores data from single-cell experiments in a special class (the SingleCellExperiment class). The rows in this data class represent features of interest (such as genes, transcripts, or genomic regions), while the columns represent individual cells. This data class is the main data structure used in the scran and scater packages. In addition to storing count data, it can also store the results of dimensionality reduction methods as well as data for additional, alternative features (like spike-in transcripts). SingleCellExperiment objects are created using a command of the same name. # AnVIL::install(c(&quot;SingleCellExperiment&quot;)) library(SingleCellExperiment) #create an empty data frame to use as an example counts &lt;- matrix(rpois(100, lambda = 10), ncol = 10, nrow = 10) pretend.cell.labels &lt;- sample(letters, ncol(counts), replace = TRUE) counts &lt;- matrix(rpois(100, lambda = 10), ncol = 10, nrow = 10) pretend.cell.labels &lt;- sample(letters, ncol(counts), replace = TRUE) pretend.gene.lengths &lt;- sample(10000, nrow(counts)) #create SingleCellExperiment object sce &lt;- SingleCellExperiment(list(counts = counts), colData=DataFrame(label = pretend.cell.labels), rowData=DataFrame(length = pretend.gene.lengths), metadata=list(study = &quot;GSE111111&quot;)) 1.4 Obtaining Data For these exercises, we are using a dataset of cells from a mouse brain from Zeisel et al. 2015. This dataset contains about 3000 brain cells of a wide variety of types (such as neurons, glia, and oligodendrocytes). The library was prepared using a UMI-based (unique molecular identifiers) protocol, and the expression data you are working with is UMI counts that have been mapped to each gene. You can learn more about UMI protocols and why they are useful here. This dataset has already had low-quality cells removed prior to publication, but we will still run QC steps to verify that the remaining cells are of good quality. The dataset is stored as a SingleCellExperiment object in the scRNAseq package. Here we retrieve the dataset, as well as do a little bit of extra formatting to merge together redundant rows (which are the result of alternative genomic coordinates for the same gene). # AnVIL::install(c(&quot;scRNAseq&quot;, &quot;scater&quot;)) library(scRNAseq) sce.zeisel &lt;- ZeiselBrainData() #library(scater) #sce.zeisel &lt;- aggregateAcrossFeatures(sce.zeisel, id=sub(&quot;_loc[0-9]+$&quot;, &quot;&quot;, rownames(sce.zeisel))) "],["quality-control-of-the-library.html", "Chapter 2 Quality control of the library 2.1 Visualizing distributions of diagnostic metrics and identifying low-quality cells 2.2 Identifying and removing low-quality cells", " Chapter 2 Quality control of the library Laboratory and technical difficulties during the scRNA-seq data generation can result in poor quality libraries. Before any analysis can be run, the data itself must be checked and cleaned. If these low quality samples are not removed from the data, they can cause issues in the downstream analysis. 2.1 Visualizing distributions of diagnostic metrics and identifying low-quality cells There are several common measures researchers use to determine sample quality in scRNA-seq: library size, the total sum of counts across all relevant features (in this case, genes) for each cell. A small library size indicates the RNA has been lost at some point during the library preparation. number of expressed features per cell, the number of genes with non-zero counts. Few expressed genes suggests the transcript population was not successfully captured during the library prep. proportion of reads mapped to mitochondrial genes relative to total count across all features. A high proportion suggests cytoplasmic RNA was lost from cells due to perforation of the cell membrane. These holes are large enough to let individual transcripts to escape but too small for mitochondria to escape, leading to the library being enriched for mitochondrial genes. Some experiments will use other metrics for quality control, such as RNA spike-ins (RNA transcripts of known sequence and quantity). In these experiments, researchers may choose to use proportion of spike-in reads relative to total counts as a QC (quality control) metric. This proportion should be similar across all samples, so any enrichment of the spike-in indicates a poor-quality sample. The Zeisel brain dataset does include spike-ins (s standardized control mix, labeled ‘ERCC’ [for the External RNA Controls Consortium, which sets standards]), so we will use this metric. It is always a good idea to plot any metrics you’re using as indicators of sample quality as the first step. The easiest way to identify low-quality cells is to set a quality threshold and remove any samples that fail to reach it. Some researchers use fixed thresholds, although this approach requires prior knowledge and experience to know what appropriate thresholds might be. We instead can define thresholds based on deviation from the median value of each metric. This approach assumes most of the samples in a dataset are high-quality. Deviation in this case is defined as the median absolute deviation (MAD). Samples are considered outliers if any of the following are true: library size &lt; 3 MADs from the median number of expressed reads &lt; 3 MADs from the median proportion of mitochondrial reads &gt; 3 MADs from the median (If your dataset uses spike-ins, treat the proportion of spike-ins similar to the proportion of mitochondrial reads.) library(scater) # identifying low-quality cells stats &lt;- perCellQCMetrics(sce.zeisel, subsets = list(Mt = rowData(sce.zeisel)$featureType == &quot;mito&quot;)) qc &lt;- quickPerCellQC(stats, percent_subsets = c(&quot;altexps_ERCC_percent&quot;, &quot;subsets_Mt_percent&quot;)) # final sample size of samples failing QC colSums(as.matrix(qc)) ## low_lib_size low_n_features high_altexps_ERCC_percent ## 0 3 65 ## high_subsets_Mt_percent discard ## 128 189 It would also be informative to see where the samples that failed QC (quality control) fall in the distributions of all of our quality markers. # visualizing low-quality cells in distributions unfiltered &lt;- sce.zeisel # we&#39;ll keep a copy of the unfiltered data, so we can see where on the distributions the eliminated samples lie colData(unfiltered) &lt;- cbind(colData(unfiltered), stats) unfiltered$discard &lt;- qc$discard # this gives us a list of the samples that have failed QC # plotting the distributions, coloring the samples that will be removed in orange gridExtra::grid.arrange( plotColData(unfiltered, y = &quot;sum&quot;, colour_by = &quot;discard&quot;) + scale_y_log10() + ggtitle(&quot;Total count&quot;), plotColData(unfiltered, y = &quot;detected&quot;, colour_by = &quot;discard&quot;) + scale_y_log10() + ggtitle(&quot;Detected features&quot;), plotColData(unfiltered, y = &quot;altexps_ERCC_percent&quot;, colour_by = &quot;discard&quot;) + ggtitle(&quot;ERCC percent&quot;), plotColData(unfiltered, y = &quot;subsets_Mt_percent&quot;, colour_by = &quot;discard&quot;) + ggtitle(&quot;Mito percent&quot;), ncol=2 ) 2.2 Identifying and removing low-quality cells After you have identified the low-quality cells, you can easily remove them in R. # creating a dataset with only the samples that pass QC sce.zeisel.qc &lt;- sce.zeisel[,!qc$discard] QUESTIONS How many samples did you remove as suspected low-quality cells? Do all the quality metrics (total count, number of expressed features, percentage of mitochondrial reads, and percentage of ERCC reads) agree when it comes to identifying low-quality samples? Which metric resulted in the removal of the greatest number of samples? Some researchers mark the suspected low-quality cells instead of removing them from the analysis. Ideally, the low-quality cells form their own cluster that can be ignored when interpreting the scRNA-seq results. This approach prevents the researcher from discarding cells or cell types that represent a true biological state and happen to have poor QC metrics. "],["normalization.html", "Chapter 3 Normalization 3.1 Calculating scaling factors 3.2 Applying scaling factors and computing normalized expression values", " Chapter 3 Normalization Single-cell RNA-seq libraries often vary in sequence coverage due to technical differences in cDNA capture or PCR amplification. We use normalization to correct for these differences. 3.1 Calculating scaling factors The simplest approach is to use library size as a scaling factor for each cell. One assumption underlying the library size normalization method is that library size is an unbiased estimate of capture or amplification issues. While this assumption is not necessarily biologically realistic, in practice the accuracy of the normalization step is unlikely to cause issues in exploratory scRNA-seq analyses. Sometimes unbalanced differential expression exists between samples, which results in composition biases in the dataset. Library size normalization can result in the non-differentially-expressed genes falsely appearing to be downregulated in one cell compared to another. If this sounds familiar, it’s because this topic has been well-studied for bulk RNA-seq analysis. Several bulk normalization methods have been created to deal with this problem in bulk RNA-seq analyses. In order to translate these algorithms to scRNA-seq and the large number of zero and low counts, we pool counts from many cells for accurate sizing factor estimation. These estimated sizing factors are then “deconvolved” into cell-based factors. This process is called normalization by deconvolution. # AnVIL::install(&quot;scran&quot;) library(scran) # calculating scaling factors set.seed(1000) clusters &lt;- quickCluster(sce.zeisel.qc) sce.zeisel.qc &lt;- computeSumFactors(sce.zeisel.qc, cluster = clusters) # comparing scaling factors based on library size (librarySizeFactors) and deconvolution (sizeFactors) plot(librarySizeFactors(sce.zeisel.qc), sizeFactors(sce.zeisel.qc), pch = 16, xlab = &quot;Library size factors&quot;, ylab = &quot;Deconvolution factors&quot;, log = &quot;xy&quot;) QUESTIONS How do the scaling factors using library size normalization compare to the scaling factors using deconvolution normalization? 3.2 Applying scaling factors and computing normalized expression values Once we have calculated scaling factors, we can compute normalized expression values for each sample by dividing the count for each gene by the size factor for the cell. The function we are using also log-transforms these normalized values, which will allow us to calculate log-fold changes in expression (a much more intuitive measure of differential expression!) This lets us focus on contributions from genes with strong relative differences. # computing normalized expression values sce.zeisel.qc &lt;- logNormCounts(sce.zeisel.qc) "],["choosing-features-for-the-analysis.html", "Chapter 4 Choosing features for the analysis 4.1 Examining per-gene variation 4.2 Selecting highly variable genes", " Chapter 4 Choosing features for the analysis The goal with many scRNA-seq analyses is to identify differences in gene expression across multiple cell types; in order to do this, cells are clustered together into groups based on a single similarity (or dissimilarity) value between a pair of cells. Which genes you choose to include in the calculation of the similarity metric can greatly influence your downstream analyses. In general, you will choose to include those genes with the greatest variability in gene expression across all the samples. The underlying assumption is this increased variability is “true” biological variation as opposed to technical noise. 4.1 Examining per-gene variation We can calculate the per-gene variation across all samples using the variance of the log-normalized expression values you computed in the Normalization section. Because we are using a log-transformation, the total variance of a gene is partly driven by its abundance (not underlying biological heterogeneity). To correct for this relationship, we use a model that calculates the expected amount of variation for each abundance value. We then use the difference between the observed variation and the expected variation for each gene to identify highly-variable genes while controlling for cell abundance. The Zeisel data was generated across multiple plates, so we need to be wary of any possible batch effects that could be driving highly-variable genes. Normally we would address this by blocking. However, each plate only contained 20-40 cells, and the cell population as a whole is highly heterogeneous, making it unlikely the sampled cell types on each plate is the same (one of the assumptions of blocking on plate). Thus, we will ignore blocking for this analysis. If you work with a dataset that does not contain spike-ins, you would use the modelGeneVar() command instead of modelGeneVarWithSpikes(). # calculating per-gene variation for each samples using the log-normalized expression value; the &quot;ERCC&quot; option tells the algorithm that information about technical (that is, non-biological) variation can be calculated using the data column labeled &quot;ERCC&quot; dec.zeisel.qc &lt;- modelGeneVarWithSpikes(sce.zeisel.qc, &quot;ERCC&quot;) # plotting the expected variance (blue curve and orange points), as well as the observed variance (black points) plot(dec.zeisel.qc$mean, dec.zeisel.qc$total, pch = 16, cex = 0.5, xlab = &quot;Mean of log-expression&quot;, ylab = &quot;Variance of log-expression&quot;) curfit &lt;- metadata(dec.zeisel.qc) points(curfit$mean, curfit$var, col = &quot;orange&quot;, pch = 16) curve(curfit$trend(x), col = &#39;dodgerblue&#39;, add = TRUE, lwd = 2) The expected variance for each log-expression value is plotted in orange. Observed values (the black points) much higher than this trend line will be used for downstream analyses. 4.2 Selecting highly variable genes When choosing the genes for downstream analysis, you will need to balance choosing as many genes as possible (so as not to exclude important biological variation) with limiting the amount of random noise. The most straightforward approach is to simply take the top n genes after ranking based on the biological component of variance. We are creating objects with the top 500, 1000, and 2000 genes, so that we can examine how our choice of n impacts downstream calculations and steps. # getting the top 500, 1000, and 2000 most-variable genes top500.hvgs &lt;- getTopHVGs(dec.zeisel.qc, n = 500) top1000.hvgs &lt;- getTopHVGs(dec.zeisel.qc, n = 1000) top2000.hvgs &lt;- getTopHVGs(dec.zeisel.qc, n = 2000) # look at the cut-offs dec.zeisel.qc &lt;- dec.zeisel.qc[order(dec.zeisel.qc$bio, decreasing = TRUE),] dec.zeisel.qc[c(1,500,1000,2000),1:6] ## DataFrame with 4 rows and 6 columns ## mean total tech bio p.value FDR ## &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Plp1 3.97379 16.80426 0.275699 16.528564 0.00000e+00 0.00000e+00 ## Shfm1 2.64467 1.55363 0.534672 1.018959 5.54040e-47 2.04392e-45 ## Cdc42 3.34431 1.04001 0.385442 0.654563 9.97095e-38 2.96240e-36 ## Vdac2 1.81697 1.08138 0.720031 0.361348 7.90102e-05 2.55550e-04 In this dataframe, the bio column represents the excess variation in gene expression (that is, the difference between the observed and expected expression). We are looking at four different genes - the gene with the greatest excess variation, as well as the last genes included in the top 500, top 1000, and top 2000 genes. QUESTIONS What is the range of log-fold expression changes (the excess variation) when you choose the top 500 genes? What about when you choose the top 1000 genes? The top 2000 genes? sessionInfo() ## R version 4.1.3 (2022-03-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.5 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] scran_1.22.1 scater_1.22.0 ## [3] ggplot2_3.3.5 scuttle_1.4.0 ## [5] scRNAseq_2.8.0 SingleCellExperiment_1.16.0 ## [7] SummarizedExperiment_1.24.0 Biobase_2.54.0 ## [9] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 ## [11] IRanges_2.28.0 S4Vectors_0.32.4 ## [13] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 ## [15] matrixStats_0.61.0 ## ## loaded via a namespace (and not attached): ## [1] AnnotationHub_3.2.2 BiocFileCache_2.2.1 ## [3] igraph_1.3.1 lazyeval_0.2.2 ## [5] BiocParallel_1.28.3 digest_0.6.29 ## [7] ensembldb_2.18.4 htmltools_0.5.2 ## [9] viridis_0.6.2 fansi_1.0.3 ## [11] magrittr_2.0.3 memoise_2.0.1 ## [13] ScaledMatrix_1.2.0 cluster_2.1.2 ## [15] limma_3.50.3 Biostrings_2.62.0 ## [17] prettyunits_1.1.1 colorspace_2.0-3 ## [19] blob_1.2.3 rappdirs_0.3.3 ## [21] ggrepel_0.9.1 xfun_0.26 ## [23] dplyr_1.0.8 crayon_1.5.1 ## [25] RCurl_1.98-1.6 jsonlite_1.8.0 ## [27] glue_1.6.2 gtable_0.3.0 ## [29] zlibbioc_1.40.0 XVector_0.34.0 ## [31] DelayedArray_0.20.0 BiocSingular_1.10.0 ## [33] scales_1.2.1 edgeR_3.36.0 ## [35] DBI_1.1.2 Rcpp_1.0.8.3 ## [37] viridisLite_0.4.0 xtable_1.8-4 ## [39] progress_1.2.2 dqrng_0.3.0 ## [41] bit_4.0.4 rsvd_1.0.5 ## [43] metapod_1.2.0 httr_1.4.2 ## [45] ellipsis_0.3.2 pkgconfig_2.0.3 ## [47] XML_3.99-0.9 farver_2.1.0 ## [49] sass_0.4.1 dbplyr_2.1.1 ## [51] locfit_1.5-9.5 utf8_1.2.2 ## [53] tidyselect_1.1.2 labeling_0.4.2 ## [55] rlang_1.0.2 later_1.3.0 ## [57] AnnotationDbi_1.56.2 munsell_0.5.0 ## [59] BiocVersion_3.14.0 tools_4.1.3 ## [61] cachem_1.0.6 cli_3.2.0 ## [63] generics_0.1.2 RSQLite_2.2.12 ## [65] ExperimentHub_2.2.1 evaluate_0.15 ## [67] stringr_1.4.0 fastmap_1.1.0 ## [69] yaml_2.3.5 knitr_1.33 ## [71] bit64_4.0.5 purrr_0.3.4 ## [73] KEGGREST_1.34.0 AnnotationFilter_1.18.0 ## [75] sparseMatrixStats_1.6.0 mime_0.12 ## [77] xml2_1.3.3 biomaRt_2.50.3 ## [79] compiler_4.1.3 beeswarm_0.4.0 ## [81] filelock_1.0.2 curl_4.3.2 ## [83] png_0.1-7 interactiveDisplayBase_1.32.0 ## [85] statmod_1.4.36 tibble_3.1.6 ## [87] bslib_0.3.1 stringi_1.7.6 ## [89] highr_0.9 GenomicFeatures_1.46.5 ## [91] lattice_0.20-45 bluster_1.4.0 ## [93] ProtGenerics_1.26.0 Matrix_1.4-0 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.16 ## [99] jquerylib_0.1.4 BiocNeighbors_1.12.0 ## [101] bitops_1.0-7 irlba_2.3.5 ## [103] httpuv_1.6.5 rtracklayer_1.54.0 ## [105] R6_2.5.1 BiocIO_1.4.0 ## [107] bookdown_0.24 promises_1.2.0.1 ## [109] gridExtra_2.3 vipor_0.4.5 ## [111] assertthat_0.2.1 rjson_0.2.21 ## [113] withr_2.5.0 GenomicAlignments_1.30.0 ## [115] Rsamtools_2.10.0 GenomeInfoDbData_1.2.7 ## [117] parallel_4.1.3 hms_1.1.1 ## [119] grid_4.1.3 beachmat_2.10.0 ## [121] rmarkdown_2.10 DelayedMatrixStats_1.16.0 ## [123] shiny_1.7.1 ggbeeswarm_0.6.0 ## [125] restfulr_0.0.13 "],["dimensionality-reduction.html", "Chapter 5 Dimensionality reduction 5.1 Calculating and choosing PCs 5.2 Applying non-linear visualization methods to PCs", " Chapter 5 Dimensionality reduction Each gene in the data represents a different dimension of the data. Reducing the number of dimensions in our data has multiple benefits, including reducing the computational work needed for downstream analyses. It also reduces the noise in the data through averaging the signal across multiple genes. Dimensionality reduction is a very common technique used in data science in general, not just in scRNA-seq analysis. You will find yourself using it over and over whenever you work with high-dimensional data. Dimensionality reduction is possible for genomic expression methods because so many genes have correlated expression. This is a consequence of different genes being involved in the same biological processes. 5.1 Calculating and choosing PCs We first use principal component analysis (PCA), which is a dimensionality reduction method that maximizes the amount of variation captured by each component, or PC. It’s up to the researcher to choose how many PCs to use for downstream analyses. More PCs mean that more biological signal is retained in the data, but it also increases the potential for noise. In our analyses, we will use 50 (the default). We will also use the top 1000 genes to calculate the PCs. library(scran) # calculating PCA # the denoisePCA command calculates PCs and removes those that primarily capture technical noise set.seed(101011001) sce.zeisel.PCA.1000 &lt;- denoisePCA(sce.zeisel.qc, technical = dec.zeisel.qc, subset.row = top1000.hvgs) #the technical option tells R where to find information about how much of the variation is attributed to &quot;technical&quot;, or non-biological, sources # visualizing the percentage of variation explained by each PC percent.var &lt;- attr(reducedDim(sce.zeisel.PCA.1000), &quot;percentVar&quot;) plot(percent.var, log = &quot;y&quot;, xlab = &quot;PC&quot;, ylab = &quot;Variance explained (%)&quot;) In PCA, the total amount of variation captured decreases for each subsequent PC. By the 10th PC, each additional PC is contributing only a small fraction to the total amount of variation explained in the dataset. Excluding them from downstream analyses has no major effect, and researchers will typically choose to include somewhere between 10 to 50 PCs. Including more PCs in the downstream analyses could theoretically cause the calculations to take longer, but in reality most calculations are fast enough that any slowdown isn’t really noticeable. More detailed information on calculating and choosing PCs for genomic analyses can be found in the Statistics for Genomics: PCA book. 5.2 Applying non-linear visualization methods to PCs In scRNA-seq analysis, plotting PCs generally does not offer enough resolution to visualize cell clusters. Instead, we rely on additional dimensionality reduction methods that can use non-linear data transformation algorithms. The most common approach is the t-stochastic neighbor embedding (t-SNE) method (Van der Maaten and Hinton 2008). t-SNE maps high-dimensional data in a low-dimensional space by first calculating the Euclidean distance between each set of points, then converting those distances into the probability that given pair of points are neighbors. On t-SNE plots, points that are members of the same cluster have a high probability of being neighbors. However, you can’t judge the similarity of different clusters based on their position on the final plot, because the t-SNE algorithm does not retain that information. The t-SNE approach is computationally complex. In practice, we reduce the computational complexity in scRNA-seq analysis by performing t-SNE calculations on the top PCs in a dataset (this both decreases the amount of computational power and time needed for analysis). We also need to set an initial starting seed and the perplexity parameter (the number of effective neighbors for each point). This parameter will determine the resolution of the plot. Lower perplexity values allow for finer resolution of population structure but can also be noisy. It’s a good idea to test multiple perplexity values when running your t-SNE analysis. We’re going to use the PCs calculated using the top 1000 highly-variable genes (a common threshold) for the rest of the analyses. library(BiocSingular) # this code first calculates the t-SNE values using PCs, and then creates a plot of the first two t-SNE dimensions set.seed(100) sce.zeisel.tsne5 &lt;- runTSNE(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;, perplexity = 5) out5 &lt;- plotReducedDim(sce.zeisel.tsne5, dimred = &quot;TSNE&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;perplexity = 5&quot;) set.seed(100) sce.zeisel.tsne20 &lt;- runTSNE(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;, perplexity = 20) out20 &lt;- plotReducedDim(sce.zeisel.tsne20, dimred = &quot;TSNE&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;perplexity = 20&quot;) set.seed(100) sce.zeisel.tsne80 &lt;- runTSNE(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;, perplexity = 80) out80 &lt;- plotReducedDim(sce.zeisel.tsne80, dimred = &quot;TSNE&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;perplexity = 80&quot;) out5 out20 out80 Some researchers will use a different method for the non-linear visualization step in their analysis. This algorithm, called uniform manifold approximation and projection (UMAP), is faster and preserves more of the global data structure when reducing dimensions compared to t-SNE (that is, you get more information about the similarity of clusters, not just of points). As a result, though, the resolution within each cluster is reduced. UMAP is becoming the method of choice as scRNA-seq datasets become larger and larger. # calculate the UMAP values from the PCs, then plot the first two UMAP dimensions set.seed(1100101001) sce.zeisel.umap &lt;- runUMAP(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;) out.umap &lt;- plotReducedDim(sce.zeisel.umap, dimred = &quot;UMAP&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;UMAP&quot;) out20 out.umap QUESTIONS 1. How does changing the perplexity parameter affect the t-SNE plot? How does the t-SNE plot compare to the UMAP plot? "],["clustering.html", "Chapter 6 Clustering 6.1 Clustering using graph-based methods", " Chapter 6 Clustering At its core, clustering is a tool that allows us to examine structure and patterns in our data. There isn’t really a “true” answer to how the data should be clustered. Instead, we can change the algorithms and parameters to explore a variety of possibilities that work best for each dataset and question the researcher is trying to answer. 6.1 Clustering using graph-based methods Graph-based clustering is based on identifying the nearest neighbors of each cell in high-dimensional space. The connections between a cell and its neighbors (called edges) are weighted based on the similarity of the two cells connected. An edge is assigned a higher weight if the two cells it connects are more closely related. After all cells have been connected to their neighbors, we apply an algorithm to identify clusters, or communities, of related cells. Each cell within a community will be more closely related to any cell within the same community than to cells outside the community. Graph-based clustering scales easily, because it only used a k-nearest neighbor search. These searches run more quickly than other methods (like hierarchical clustering). Unfortunately, no information is retained about relationships beyond the neighboring cells. This effect also means that clustering resolution depends on cell density. # let the R algorithm define and label our clusters nn.clusters &lt;- clusterCells(sce.zeisel.tsne20, use.dimred=&quot;PCA&quot;) # this command tells us how many clusters were identified (the top row) and how many cells belong to each cluster (the bottom row) table(nn.clusters) ## nn.clusters ## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ## 284 147 107 195 604 166 475 62 119 256 221 67 53 28 32 We assigned the cluster assignments back into our SingleCellExperiment object as a factor in the column metadata, which allows us to visualize the cluster assignment on a t-SNE plot. # create a t-SNE plot showing the identified clusters colLabels(sce.zeisel.tsne20) &lt;- nn.clusters plotReducedDim(sce.zeisel.tsne20, &quot;TSNE&quot;, colour_by=&quot;label&quot;) QUESTIONS 1. How many clusters were identified using graph-based clustering? Which cluster contained the most cells, and how many cells did it have? sessionInfo() ## R version 4.1.3 (2022-03-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.5 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] uwot_0.1.11 Matrix_1.4-0 ## [3] BiocSingular_1.10.0 scran_1.22.1 ## [5] scater_1.22.0 ggplot2_3.3.5 ## [7] scuttle_1.4.0 scRNAseq_2.8.0 ## [9] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0 ## [11] Biobase_2.54.0 GenomicRanges_1.46.1 ## [13] GenomeInfoDb_1.30.1 IRanges_2.28.0 ## [15] S4Vectors_0.32.4 BiocGenerics_0.40.0 ## [17] MatrixGenerics_1.6.0 matrixStats_0.61.0 ## ## loaded via a namespace (and not attached): ## [1] AnnotationHub_3.2.2 BiocFileCache_2.2.1 ## [3] igraph_1.3.1 lazyeval_0.2.2 ## [5] BiocParallel_1.28.3 digest_0.6.29 ## [7] ensembldb_2.18.4 htmltools_0.5.2 ## [9] viridis_0.6.2 fansi_1.0.3 ## [11] magrittr_2.0.3 memoise_2.0.1 ## [13] ScaledMatrix_1.2.0 cluster_2.1.2 ## [15] limma_3.50.3 Biostrings_2.62.0 ## [17] prettyunits_1.1.1 colorspace_2.0-3 ## [19] blob_1.2.3 rappdirs_0.3.3 ## [21] ggrepel_0.9.1 xfun_0.26 ## [23] dplyr_1.0.8 crayon_1.5.1 ## [25] RCurl_1.98-1.6 jsonlite_1.8.0 ## [27] glue_1.6.2 gtable_0.3.0 ## [29] zlibbioc_1.40.0 XVector_0.34.0 ## [31] DelayedArray_0.20.0 scales_1.2.1 ## [33] DBI_1.1.2 edgeR_3.36.0 ## [35] Rcpp_1.0.8.3 viridisLite_0.4.0 ## [37] xtable_1.8-4 progress_1.2.2 ## [39] dqrng_0.3.0 bit_4.0.4 ## [41] rsvd_1.0.5 metapod_1.2.0 ## [43] httr_1.4.2 FNN_1.1.3 ## [45] ellipsis_0.3.2 farver_2.1.0 ## [47] pkgconfig_2.0.3 XML_3.99-0.9 ## [49] sass_0.4.1 dbplyr_2.1.1 ## [51] locfit_1.5-9.5 utf8_1.2.2 ## [53] labeling_0.4.2 tidyselect_1.1.2 ## [55] rlang_1.0.2 later_1.3.0 ## [57] AnnotationDbi_1.56.2 munsell_0.5.0 ## [59] BiocVersion_3.14.0 tools_4.1.3 ## [61] cachem_1.0.6 cli_3.2.0 ## [63] generics_0.1.2 RSQLite_2.2.12 ## [65] ExperimentHub_2.2.1 evaluate_0.15 ## [67] stringr_1.4.0 fastmap_1.1.0 ## [69] yaml_2.3.5 knitr_1.33 ## [71] bit64_4.0.5 purrr_0.3.4 ## [73] KEGGREST_1.34.0 AnnotationFilter_1.18.0 ## [75] sparseMatrixStats_1.6.0 mime_0.12 ## [77] xml2_1.3.3 biomaRt_2.50.3 ## [79] compiler_4.1.3 beeswarm_0.4.0 ## [81] filelock_1.0.2 curl_4.3.2 ## [83] png_0.1-7 interactiveDisplayBase_1.32.0 ## [85] statmod_1.4.36 tibble_3.1.6 ## [87] bslib_0.3.1 stringi_1.7.6 ## [89] highr_0.9 RSpectra_0.16-0 ## [91] GenomicFeatures_1.46.5 lattice_0.20-45 ## [93] bluster_1.4.0 ProtGenerics_1.26.0 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.16 ## [99] jquerylib_0.1.4 BiocNeighbors_1.12.0 ## [101] bitops_1.0-7 irlba_2.3.5 ## [103] httpuv_1.6.5 rtracklayer_1.54.0 ## [105] R6_2.5.1 BiocIO_1.4.0 ## [107] bookdown_0.24 promises_1.2.0.1 ## [109] gridExtra_2.3 vipor_0.4.5 ## [111] assertthat_0.2.1 rjson_0.2.21 ## [113] withr_2.5.0 GenomicAlignments_1.30.0 ## [115] Rsamtools_2.10.0 GenomeInfoDbData_1.2.7 ## [117] parallel_4.1.3 hms_1.1.1 ## [119] grid_4.1.3 beachmat_2.10.0 ## [121] rmarkdown_2.10 DelayedMatrixStats_1.16.0 ## [123] Rtsne_0.16 shiny_1.7.1 ## [125] ggbeeswarm_0.6.0 restfulr_0.0.13 "],["marker-gene-detection.html", "Chapter 7 Marker gene detection 7.1 Calculating and ranking effect size summary statistics 7.2 Comparing gene expression levels across clusters", " Chapter 7 Marker gene detection We have managed to identify clusters of related cells in the data. However, these clusters aren’t very useful until we can identify the biology meaning of each group. This is where functional annotation comes in. The real art, and the greatest challenge, in an scRNA-seq analysis comes when interpreting the results. Up to this point (cleaning and clustering the data), the analysis and computation has been straightforward. Figuring out the biological state that each cluster represents, on the other hand, is more difficult, as it requires applying prior biological knowledge to the dataset. Thanks to previous research, we know many marker genes, or genes can be used to identify particular cell types. These genes are differentially expressed across cell types, and by examining the expression profiles of multiple marker genes across all the clusters, we can assign particular cell type identities to each cluster. 7.1 Calculating and ranking effect size summary statistics We begin by comparing each pair of clusters and calculating scores for expression differences between the two for each gene. We have multiple options for the statistics used to compare expression values. AUC (area under the curve) is the probability that a randomly chosen observation from cluster A is greater than a randomly chosen observation from cluster B. This statistic is a way to quantify how well we can distinguish between two distributions (clusters) in a pairwise comparison. An AUC of 1 means all values in cluster A are greater than any value from cluster B and suggests upregulation. An AUC of of 0.5 means the two clusters are indistinguishable from each other, while an AUC of 0 suggests the marker gene observations in cluster A are downregulated compared to those in cluster B. Cohen’s d is a standardized log-fold change, and can be thought of as the number of standard deviations that separate the two groups. Positive values of Cohen’s d suggest that our cluster of interest (cluster A) are upregulated compared to cluster B, while negative values suggest the marker gene observations in cluster A are downregulated compared to cluster B. log-fold change (logFC) is a measure of whether there is a difference in expression between clusters. Keep in mind that these values ignore the magnitude of the change. As with the others, positive values indicate upregulation in the cluster of interest (cluster A), while negative values indicate downregulation. For each of these statistics, scoreMarkers calculates mean, median, minimum value (min), maximum value (max), and minimum rank (rank; the smallest rank of each gene across all pairwise comparisons). For most of these measures, a larger number indicates upregulation. For minimum rank, however, a small value means the gene is one of the top upregulated genes. AUC or Cohen’s d are effective regardless of the magnitude of the expression values and thus are good choices for general marker detection. The log-fold change in the detected proportion is specifically useful for identifying binary changes in expression. For this exercise, we’re going to focus on upregulated markers, since those are particularly useful for identifying cell types in a heterogeneous population like the Zeisel dataset. We use the findMarkers command, which quickly identifies potential marker genes. # identify those genes which are upregulated in some clusters compared to others markers &lt;- findMarkers(sce.zeisel.tsne20, direction = &quot;up&quot;) marker.set &lt;- markers[[&quot;1&quot;]] head(marker.set, 10) ## DataFrame with 10 rows and 18 columns ## Top p.value FDR summary.logFC logFC.2 logFC.3 ## &lt;integer&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Syngr3 1 5.71433e-140 1.90535e-137 2.44603 1.156187 1.1330628 ## Mllt11 1 3.31857e-249 7.37681e-246 2.87647 0.971725 1.4504728 ## Ndrg4 1 0.00000e+00 0.00000e+00 3.83936 1.146289 0.9309904 ## Slc6a1 1 1.75709e-160 8.57373e-158 3.46313 3.063161 3.2821288 ## Gad1 1 4.70797e-234 8.56251e-231 4.54282 4.024699 4.3901265 ## Gad2 1 1.32400e-208 1.76586e-205 4.24931 3.756435 3.9505814 ## Atp1a3 1 3.54623e-279 1.41892e-275 3.45285 1.213985 0.0213418 ## Slc32a1 2 1.50335e-109 2.48562e-107 1.91750 1.791641 1.8496954 ## Rcan2 2 9.25451e-129 2.50197e-126 2.21812 1.432185 2.1127708 ## Rab3a 2 1.61475e-202 1.61523e-199 2.52266 0.884615 0.6329209 ## logFC.4 logFC.5 logFC.6 logFC.7 logFC.8 logFC.9 logFC.10 ## &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Syngr3 2.38928 2.44603 0.7085503 0.8470822 1.159959 2.525116 1.030908 ## Mllt11 2.94289 2.87647 0.4897958 0.0882716 1.030983 3.262597 0.471414 ## Ndrg4 3.67150 3.83936 0.3680955 0.7169171 1.103046 3.646328 0.791827 ## Slc6a1 2.92969 2.88778 3.6063318 3.5292982 2.218826 0.704948 3.571574 ## Gad1 4.37025 4.46268 4.6478826 4.5428243 4.477666 4.386587 4.572889 ## Gad2 4.08621 4.15971 4.2916168 4.2222656 4.229715 4.207704 4.249309 ## Atp1a3 3.18678 3.45285 0.6498389 0.1656598 0.426762 3.436859 -0.228370 ## Slc32a1 1.93798 1.91658 1.9393715 1.9092434 1.913338 1.917505 1.910826 ## Rcan2 1.81884 1.25488 1.1051970 2.1730716 2.170402 2.102633 2.218119 ## Rab3a 2.73747 2.52266 0.0267138 0.3133349 0.695216 2.921335 0.163635 ## logFC.11 logFC.12 logFC.13 logFC.14 logFC.15 ## &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Syngr3 0.684417 2.56482 2.55195 2.58830 2.42180 ## Mllt11 0.229326 3.04261 2.77725 3.08889 2.01497 ## Ndrg4 0.565654 3.82140 4.14743 3.97144 4.01655 ## Slc6a1 3.463131 3.43495 3.34170 3.05042 3.25178 ## Gad1 4.568481 4.51547 4.50591 4.39482 4.63836 ## Gad2 4.263975 4.16831 4.29808 4.09088 4.35014 ## Atp1a3 1.118434 3.31854 3.78426 3.67087 3.60817 ## Slc32a1 1.916869 1.94707 1.94707 1.94707 1.90680 ## Rcan2 1.142480 2.09515 1.26318 1.46319 2.24552 ## Rab3a 0.508261 3.16930 2.54753 3.06401 2.30063 This dataframe shows us the in log-fold expression change for each potential marker gene between cluster 1 and every other cluster. 7.2 Comparing gene expression levels across clusters Once we’ve identified potential marker genes, we can use a heatmap to compare gene expression in each cell between clusters. # pull the top-most upregulated markers from cluster 1 (compared to the rest of the clusters) and look at their expression in all clusters top.markers &lt;- rownames(marker.set)[marker.set$Top &lt;= 10] plotHeatmap(sce.zeisel.tsne20, features = top.markers, order_columns_by = &quot;label&quot;) In this heatmap, clusters are on the horizontal, while the top upregulated genes in cluster 1 are on the vertical. The magnitude of the log-fold expression change is indicated by color of each cell. We can also create a heatmap that shows the mean log-fold change of cluster 1 cells compared to the mean of each other cluster. This can simplify the heatmap and is useful when dealing with many clusters. # AnVIL::install(&quot;pheatmap&quot;) library(pheatmap) #this heatmap lets us compare the average expression of the gene within a cluster compared to the other clusters logFCs &lt;- getMarkerEffects(marker.set[1:50,]) pheatmap(logFCs, breaks = seq(-5, 5, length.out = 101)) Here we see that three genes are generally upregulated in Cluster 1 compared to the other clusters: Gad1, Gad2, and Slc6a1. This is where prior biological knowledge comes in handy, as both Gad1 and Slc6a1 are known interneuron markers (Zeng et al. 2012). QUESTION 1. Are there any groups or patterns you see in the second heatmap that look interesting? "],["cell-type-annotation.html", "Chapter 8 Cell type annotation 8.1 Annotating cell types", " Chapter 8 Cell type annotation Unless you are already an expert in neuronal cell expression, you probably didn’t know Gad1 and Slc6a1 are known interneuron markers until you were told. Luckily, we have several high-quality references databases that can be used for annotating scRNA-seq datasets. Some of these references are useful for identifying cell types (the Zeisel dataset we have been using is a well-known reference for identifying neuronal cell types). Others, such as the Gene Ontology (GO) or Kyoto Encyclopedia of Genes and Genomes (KEGG) collections, are useful for identifying biological processes associated with each cluster. 8.1 Annotating cell types There are three basic strategies for annotating datasets: match the expression profile of each individual cell to the expression profile of cells from a reference dataset (the reference dataset approach); identify sets of marker genes highly expressed in each cell and match to gene sets from known cell types (the gene set approach); or perform a gene-set enrichment analysis on the marker genes that define each cluster. The Zeisel dataset we have been working with has actually been annotated all this time. # calculate the top marker genes assigned to each cell type (level1class) in the Zeisel dataset wilcox.z &lt;- pairwiseWilcox(sce.zeisel.qc, sce.zeisel.qc$level1class, lfc = 1, direction = &quot;up&quot;) markers.z &lt;- getTopMarkers(wilcox.z$statistics, wilcox.z$pairs, pairwise = FALSE, n = 50) # look at how many cell-type categories there are, as well as how many cells assigned to each category lengths(markers.z) ## astrocytes_ependymal endothelial-mural interneurons ## 78 83 119 ## microglia oligodendrocytes pyramidal CA1 ## 69 80 124 ## pyramidal SS ## 147 This dataset has grouped the individual cells into 7 different categories of neuronal subtypes. We can use the gene sets that define these categories to annotate a second brain scRNA-seq dataset from Tasic et al. (2016). (We will not go through all the data cleaning and clustering steps for this new dataset - having made it this far, we trust you can do that!) #load the new dataset sce.tasic &lt;- TasicBrainData() We first create the gene set lists using the GSEABase package. The AUCell package identifies and ranks marker sets highly expressed in each cell using an area under the curve (AUC) approach. We can assign cell type identity in the Tasic dataset by taking the marker set with the highest AUC as the label for the cell. # AnVIL::install(&quot;GSEABase&quot;) library(GSEABase) #create a dataset that contains just the information about the Zeisel cell-type categories and the marker genes that define each cell-type all.sets &lt;- lapply(names(markers.z), function(x) { GeneSet(markers.z[[x]], setName = x) }) all.sets &lt;- GeneSetCollection(all.sets) # AnVIL::install(&quot;AUCell&quot;) library(AUCell) # rank genes by expression values within each cell rankings &lt;- AUCell_buildRankings(counts(sce.tasic), plotStats = FALSE, verbose = FALSE) # calculate AUC for each previously-defined marker set (from Zeisel) in the Tasic data cell.aucs &lt;- AUCell_calcAUC(all.sets, rankings) results &lt;- t(assay(cell.aucs)) After assigning cell type identities to each cluster, a researcher should always verify that the identities make sense. Since the Tasic dataset has also been annotated, we can compare our annotation to the researcher-provided annotation as a sort of sanity check. # assign cell type identity in the Tasic dataset by assumig the marker set with the top AUC is the proper label new.labels &lt;- colnames(results)[max.col(results)] # compare our annotations to the annotations provided by the Tasic dataset tab &lt;- table(new.labels, sce.tasic$broad_type) tab ## ## new.labels Astrocyte Endothelial Cell GABA-ergic Neuron ## astrocytes_ependymal 43 2 0 ## endothelial-mural 0 27 0 ## interneurons 0 0 759 ## microglia 0 0 0 ## oligodendrocytes 0 0 1 ## pyramidal SS 0 0 1 ## ## new.labels Glutamatergic Neuron Microglia Oligodendrocyte ## astrocytes_ependymal 0 0 0 ## endothelial-mural 0 0 0 ## interneurons 2 0 0 ## microglia 0 22 0 ## oligodendrocytes 0 0 38 ## pyramidal SS 810 0 0 ## ## new.labels Oligodendrocyte Precursor Cell Unclassified ## astrocytes_ependymal 19 4 ## endothelial-mural 0 2 ## interneurons 0 15 ## microglia 0 1 ## oligodendrocytes 3 0 ## pyramidal SS 0 60 As you can see, the labels applied by the researchers and by our annotation mostly match. (Pyramidal SS nerves are primarily glutamatergic, so although the categories are labeled differently, those two categories do indeed match!) The major exception is the oligodendrocyte precursor cells, which our annotation called astrocytes. However, this mismatch isn’t as concerning as you might think, once you know that both astrocytes and oligodendrocytes come from the same precursor cell lineage. QUESTION 1. Should we be concerned when 1 or 2 cells are assigned a different cell type identity by our annotation than by the researchers? Why or why not? sessionInfo() ## R version 4.1.3 (2022-03-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.5 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] AUCell_1.16.0 GSEABase_1.56.0 ## [3] graph_1.72.0 annotate_1.72.0 ## [5] XML_3.99-0.9 AnnotationDbi_1.56.2 ## [7] pheatmap_1.0.12 BiocSingular_1.10.0 ## [9] scran_1.22.1 scater_1.22.0 ## [11] ggplot2_3.3.5 scuttle_1.4.0 ## [13] scRNAseq_2.8.0 SingleCellExperiment_1.16.0 ## [15] SummarizedExperiment_1.24.0 Biobase_2.54.0 ## [17] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 ## [19] IRanges_2.28.0 S4Vectors_0.32.4 ## [21] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 ## [23] matrixStats_0.61.0 ## ## loaded via a namespace (and not attached): ## [1] AnnotationHub_3.2.2 BiocFileCache_2.2.1 ## [3] igraph_1.3.1 lazyeval_0.2.2 ## [5] BiocParallel_1.28.3 digest_0.6.29 ## [7] ensembldb_2.18.4 htmltools_0.5.2 ## [9] viridis_0.6.2 fansi_1.0.3 ## [11] magrittr_2.0.3 memoise_2.0.1 ## [13] ScaledMatrix_1.2.0 cluster_2.1.2 ## [15] limma_3.50.3 Biostrings_2.62.0 ## [17] R.utils_2.12.2 prettyunits_1.1.1 ## [19] colorspace_2.0-3 blob_1.2.3 ## [21] rappdirs_0.3.3 ggrepel_0.9.1 ## [23] xfun_0.26 dplyr_1.0.8 ## [25] crayon_1.5.1 RCurl_1.98-1.6 ## [27] jsonlite_1.8.0 glue_1.6.2 ## [29] gtable_0.3.0 zlibbioc_1.40.0 ## [31] XVector_0.34.0 DelayedArray_0.20.0 ## [33] scales_1.2.1 DBI_1.1.2 ## [35] edgeR_3.36.0 Rcpp_1.0.8.3 ## [37] viridisLite_0.4.0 xtable_1.8-4 ## [39] progress_1.2.2 dqrng_0.3.0 ## [41] bit_4.0.4 rsvd_1.0.5 ## [43] metapod_1.2.0 httr_1.4.2 ## [45] RColorBrewer_1.1-3 ellipsis_0.3.2 ## [47] R.methodsS3_1.8.1 farver_2.1.0 ## [49] pkgconfig_2.0.3 sass_0.4.1 ## [51] dbplyr_2.1.1 locfit_1.5-9.5 ## [53] utf8_1.2.2 tidyselect_1.1.2 ## [55] rlang_1.0.2 later_1.3.0 ## [57] munsell_0.5.0 BiocVersion_3.14.0 ## [59] tools_4.1.3 cachem_1.0.6 ## [61] cli_3.2.0 generics_0.1.2 ## [63] RSQLite_2.2.12 ExperimentHub_2.2.1 ## [65] evaluate_0.15 stringr_1.4.0 ## [67] fastmap_1.1.0 yaml_2.3.5 ## [69] knitr_1.33 bit64_4.0.5 ## [71] purrr_0.3.4 KEGGREST_1.34.0 ## [73] AnnotationFilter_1.18.0 sparseMatrixStats_1.6.0 ## [75] mime_0.12 R.oo_1.24.0 ## [77] xml2_1.3.3 biomaRt_2.50.3 ## [79] compiler_4.1.3 beeswarm_0.4.0 ## [81] filelock_1.0.2 curl_4.3.2 ## [83] png_0.1-7 interactiveDisplayBase_1.32.0 ## [85] statmod_1.4.36 tibble_3.1.6 ## [87] bslib_0.3.1 stringi_1.7.6 ## [89] highr_0.9 GenomicFeatures_1.46.5 ## [91] lattice_0.20-45 bluster_1.4.0 ## [93] ProtGenerics_1.26.0 Matrix_1.4-0 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.16 ## [99] jquerylib_0.1.4 BiocNeighbors_1.12.0 ## [101] data.table_1.14.2 bitops_1.0-7 ## [103] irlba_2.3.5 httpuv_1.6.5 ## [105] rtracklayer_1.54.0 R6_2.5.1 ## [107] BiocIO_1.4.0 bookdown_0.24 ## [109] promises_1.2.0.1 gridExtra_2.3 ## [111] vipor_0.4.5 assertthat_0.2.1 ## [113] rjson_0.2.21 withr_2.5.0 ## [115] GenomicAlignments_1.30.0 Rsamtools_2.10.0 ## [117] GenomeInfoDbData_1.2.7 parallel_4.1.3 ## [119] hms_1.1.1 grid_4.1.3 ## [121] beachmat_2.10.0 rmarkdown_2.10 ## [123] DelayedMatrixStats_1.16.0 Rtsne_0.16 ## [125] shiny_1.7.1 ggbeeswarm_0.6.0 ## [127] restfulr_0.0.13 "],["references.html", "Chapter 9 References", " Chapter 9 References Tasic B et al. Adult mouse cortical cell taxonomy revealed by single cell transcriptomics. Nat Neurosci. 2016 Feb;19(2):335-46. doi: 10.1038/nn.4216. Epub 2016 Jan 4. PMID: 26727548; PMCID: PMC4985242. Van der Maaten L, Hinton G. Visualizing data using t-SNE. J. Mach. Learn. Res. 2008 Nov;9:2579-2605. Zeisel A et al. Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science. 2015 Mar 6;347(6226):1138-42. doi: 10.1126/science.aaa1934. Epub 2015 Feb 19. PMID: 25700174. Zeng H et al. Large-scale cellular-resolution gene profiling in human neocortex reveals species-specific molecular signatures. Cell. 2012 Apr 13;149(2):483-96. doi: 10.1016/j.cell.2012.02.052. PMID: 22500809; PMCID: PMC3328777. "],["about-the-authors.html", "About the Authors", " About the Authors These credits are based on our course contributors table guidelines.     Credits Names Pedagogy Lead Content Instructor Elizabeth Humphries Content Idea Contributor Ava Hoffman, Frederick Tan Content Directors Kasper Hansen, Jeff Leek, Frederick Tan Production Content Publisher Ira Gooding Technical Template Publishing Engineers Candace Savonen, Carrie Wright, Ava Hoffman Publishing Maintenance Engineer Candace Savonen Technical Publishing Stylists Carrie Wright, Candace Savonen Package Developers (ottrpal) John Muschelli, Candace Savonen, Carrie Wright Funding Funder National Human Genome Research Institute (NHGRI) Funding Staff Fallon Bachman, Jennifer Vessio, Emily Voeglein   ## ─ Session info ─────────────────────────────────────────────────────────────── ## setting value ## version R version 4.1.3 (2022-03-10) ## os Ubuntu 20.04.5 LTS ## system x86_64, linux-gnu ## ui X11 ## language (EN) ## collate en_US.UTF-8 ## ctype en_US.UTF-8 ## tz Etc/UTC ## date 2023-03-22 ## pandoc 2.17.1.1 @ /usr/local/bin/ (via rmarkdown) ## ## ─ Packages ─────────────────────────────────────────────────────────────────── ## package * version date (UTC) lib source ## bookdown 0.24 2023-03-22 [1] Github (rstudio/bookdown@88bc4ea) ## bslib 0.3.1 2021-10-06 [1] RSPM (R 4.1.0) ## cachem 1.0.6 2021-08-19 [1] RSPM (R 4.1.0) ## callr 3.7.0 2021-04-20 [1] RSPM (R 4.1.0) ## cli 3.2.0 2022-02-14 [1] RSPM (R 4.1.0) ## crayon 1.5.1 2022-03-26 [1] RSPM (R 4.1.0) ## desc 1.4.1 2022-03-06 [1] RSPM (R 4.1.0) ## devtools 2.4.3 2021-11-30 [1] RSPM (R 4.1.0) ## digest 0.6.29 2021-12-01 [1] RSPM (R 4.1.0) ## ellipsis 0.3.2 2021-04-29 [1] RSPM (R 4.1.0) ## evaluate 0.15 2022-02-18 [1] RSPM (R 4.1.0) ## fastmap 1.1.0 2021-01-25 [1] RSPM (R 4.1.0) ## fs 1.5.2 2021-12-08 [1] RSPM (R 4.1.0) ## glue 1.6.2 2022-02-24 [1] RSPM (R 4.1.0) ## htmltools 0.5.2 2021-08-25 [1] RSPM (R 4.1.0) ## jquerylib 0.1.4 2021-04-26 [1] RSPM (R 4.1.0) ## jsonlite 1.8.0 2022-02-22 [1] RSPM (R 4.1.0) ## knitr 1.33 2023-03-22 [1] Github (yihui/knitr@a1052d1) ## lifecycle 1.0.1 2021-09-24 [1] RSPM (R 4.1.0) ## magrittr 2.0.3 2022-03-30 [1] RSPM (R 4.1.0) ## memoise 2.0.1 2021-11-26 [1] RSPM (R 4.1.0) ## pkgbuild 1.3.1 2021-12-20 [1] RSPM (R 4.1.0) ## pkgload 1.2.4 2021-11-30 [1] RSPM (R 4.1.0) ## prettyunits 1.1.1 2020-01-24 [1] RSPM (R 4.1.0) ## processx 3.5.3 2022-03-25 [1] RSPM (R 4.1.0) ## ps 1.6.0 2021-02-28 [1] RSPM (R 4.1.0) ## purrr 0.3.4 2020-04-17 [1] RSPM (R 4.1.0) ## R6 2.5.1 2021-08-19 [1] RSPM (R 4.1.0) ## remotes 2.4.2 2021-11-30 [1] RSPM (R 4.1.0) ## rlang 1.0.2 2023-03-22 [1] Github (r-lib/rlang@34b04a8) ## rmarkdown 2.10 2023-03-22 [1] Github (rstudio/rmarkdown@02d3c25) ## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.1.3) ## sass 0.4.1 2022-03-23 [1] RSPM (R 4.1.0) ## sessioninfo 1.2.2 2021-12-06 [1] RSPM (R 4.1.0) ## stringi 1.7.6 2021-11-29 [1] RSPM (R 4.1.0) ## stringr 1.4.0 2019-02-10 [1] RSPM (R 4.1.0) ## testthat 3.0.1 2023-03-22 [1] Github (R-lib/testthat@e99155a) ## usethis 2.1.6.9000 2023-03-22 [1] Github (r-lib/usethis@b4b8623) ## withr 2.5.0 2022-03-03 [1] RSPM (R 4.1.0) ## xfun 0.26 2023-03-22 [1] Github (yihui/xfun@74c2a66) ## yaml 2.3.5 2022-02-21 [1] RSPM (R 4.1.0) ## ## [1] /usr/local/lib/R/site-library ## [2] /usr/local/lib/R/library ## ## ────────────────────────────────────────────────────────────────────────────── "],["404.html", "Page not found", " Page not found The page you requested cannot be found (perhaps it was moved or renamed). You may want to try searching to find the page's new location, or use the table of contents to find the page you are looking for. "]]
+[["index.html", "Statistics for Genomics: scRNAseq About this Book Skills Level GDSCN Collection", " Statistics for Genomics: scRNAseq May 25, 2023 About this Book This book is intended to introduce concepts in single cell RNA-seq analysis (scRNA-seq). Skills Level Genetics Programming skills GDSCN Collection This exercise is part of a collection of teaching resources developed through the Genomic Data Science Community Network (GDSCN). GDSCN works towards a vision where researchers, educators, and students from diverse backgrounds are able to fully participate in genomic data science research. Learn more about GDSCN by visiting https://www.gdscn.org/home or reading the article in Genome Research. Please check out our full collection of AnVIL and related resources: https://hutchdatascience.org/AnVIL_Collection/ "],["overview.html", "Chapter 1 OVERVIEW 1.1 Scope 1.2 Bioconductor reference book 1.3 SingleCellExperiment 1.4 Obtaining Data", " Chapter 1 OVERVIEW 1.1 Scope There are many excellent programs available for scRNA-seq analysis. In this lesson, we will be using the scran and scater packages available from Bioconductor. We do not cover .fastq processing, nor do we explore the Seurat and scanpy pipelines. 1.2 Bioconductor reference book The lessons and activities are adapted from the Orchestrating Single-Cell Analysis with Bioconductor book by Robert Amezquita, Aaron Lun, Stephanie Hicks, and Raphael Gottardo. 1.3 SingleCellExperiment Bioconductor stores data from single-cell experiments in a special class (the SingleCellExperiment class). The rows in this data class represent features of interest (such as genes, transcripts, or genomic regions), while the columns represent individual cells. This data class is the main data structure used in the scran and scater packages. In addition to storing count data, it can also store the results of dimensionality reduction methods as well as data for additional, alternative features (like spike-in transcripts). SingleCellExperiment objects are created using a command of the same name. # AnVIL::install(c(&quot;SingleCellExperiment&quot;)) library(SingleCellExperiment) #create an empty data frame to use as an example counts &lt;- matrix(rpois(100, lambda = 10), ncol = 10, nrow = 10) pretend.cell.labels &lt;- sample(letters, ncol(counts), replace = TRUE) counts &lt;- matrix(rpois(100, lambda = 10), ncol = 10, nrow = 10) pretend.cell.labels &lt;- sample(letters, ncol(counts), replace = TRUE) pretend.gene.lengths &lt;- sample(10000, nrow(counts)) #create SingleCellExperiment object sce &lt;- SingleCellExperiment(list(counts = counts), colData=DataFrame(label = pretend.cell.labels), rowData=DataFrame(length = pretend.gene.lengths), metadata=list(study = &quot;GSE111111&quot;)) 1.4 Obtaining Data For these exercises, we are using a dataset of cells from a mouse brain from Zeisel et al. 2015. This dataset contains about 3000 brain cells of a wide variety of types (such as neurons, glia, and oligodendrocytes). The library was prepared using a UMI-based (unique molecular identifiers) protocol, and the expression data you are working with is UMI counts that have been mapped to each gene. You can learn more about UMI protocols and why they are useful here. This dataset has already had low-quality cells removed prior to publication, but we will still run QC steps to verify that the remaining cells are of good quality. The dataset is stored as a SingleCellExperiment object in the scRNAseq package. Here we retrieve the dataset, as well as do a little bit of extra formatting to merge together redundant rows (which are the result of alternative genomic coordinates for the same gene). # AnVIL::install(c(&quot;scRNAseq&quot;, &quot;scater&quot;)) library(scRNAseq) sce.zeisel &lt;- ZeiselBrainData() #library(scater) #sce.zeisel &lt;- aggregateAcrossFeatures(sce.zeisel, id=sub(&quot;_loc[0-9]+$&quot;, &quot;&quot;, rownames(sce.zeisel))) "],["quality-control-of-the-library.html", "Chapter 2 Quality control of the library 2.1 Visualizing distributions of diagnostic metrics and identifying low-quality cells 2.2 Identifying and removing low-quality cells", " Chapter 2 Quality control of the library Laboratory and technical difficulties during the scRNA-seq data generation can result in poor quality libraries. Before any analysis can be run, the data itself must be checked and cleaned. If these low quality samples are not removed from the data, they can cause issues in the downstream analysis. 2.1 Visualizing distributions of diagnostic metrics and identifying low-quality cells There are several common measures researchers use to determine sample quality in scRNA-seq: library size, the total sum of counts across all relevant features (in this case, genes) for each cell. A small library size indicates the RNA has been lost at some point during the library preparation. number of expressed features per cell, the number of genes with non-zero counts. Few expressed genes suggests the transcript population was not successfully captured during the library prep. proportion of reads mapped to mitochondrial genes relative to total count across all features. A high proportion suggests cytoplasmic RNA was lost from cells due to perforation of the cell membrane. These holes are large enough to let individual transcripts to escape but too small for mitochondria to escape, leading to the library being enriched for mitochondrial genes. Some experiments will use other metrics for quality control, such as RNA spike-ins (RNA transcripts of known sequence and quantity). In these experiments, researchers may choose to use proportion of spike-in reads relative to total counts as a QC (quality control) metric. This proportion should be similar across all samples, so any enrichment of the spike-in indicates a poor-quality sample. The Zeisel brain dataset does include spike-ins (s standardized control mix, labeled ‘ERCC’ [for the External RNA Controls Consortium, which sets standards]), so we will use this metric. It is always a good idea to plot any metrics you’re using as indicators of sample quality as the first step. The easiest way to identify low-quality cells is to set a quality threshold and remove any samples that fail to reach it. Some researchers use fixed thresholds, although this approach requires prior knowledge and experience to know what appropriate thresholds might be. We instead can define thresholds based on deviation from the median value of each metric. This approach assumes most of the samples in a dataset are high-quality. Deviation in this case is defined as the median absolute deviation (MAD). Samples are considered outliers if any of the following are true: library size &lt; 3 MADs from the median number of expressed reads &lt; 3 MADs from the median proportion of mitochondrial reads &gt; 3 MADs from the median (If your dataset uses spike-ins, treat the proportion of spike-ins similar to the proportion of mitochondrial reads.) library(scater) # identifying low-quality cells stats &lt;- perCellQCMetrics(sce.zeisel, subsets = list(Mt = rowData(sce.zeisel)$featureType == &quot;mito&quot;)) qc &lt;- quickPerCellQC(stats, percent_subsets = c(&quot;altexps_ERCC_percent&quot;, &quot;subsets_Mt_percent&quot;)) # final sample size of samples failing QC colSums(as.matrix(qc)) ## low_lib_size low_n_features high_altexps_ERCC_percent ## 0 3 65 ## high_subsets_Mt_percent discard ## 128 189 It would also be informative to see where the samples that failed QC (quality control) fall in the distributions of all of our quality markers. # visualizing low-quality cells in distributions unfiltered &lt;- sce.zeisel # we&#39;ll keep a copy of the unfiltered data, so we can see where on the distributions the eliminated samples lie colData(unfiltered) &lt;- cbind(colData(unfiltered), stats) unfiltered$discard &lt;- qc$discard # this gives us a list of the samples that have failed QC # plotting the distributions, coloring the samples that will be removed in orange gridExtra::grid.arrange( plotColData(unfiltered, y = &quot;sum&quot;, colour_by = &quot;discard&quot;) + scale_y_log10() + ggtitle(&quot;Total count&quot;), plotColData(unfiltered, y = &quot;detected&quot;, colour_by = &quot;discard&quot;) + scale_y_log10() + ggtitle(&quot;Detected features&quot;), plotColData(unfiltered, y = &quot;altexps_ERCC_percent&quot;, colour_by = &quot;discard&quot;) + ggtitle(&quot;ERCC percent&quot;), plotColData(unfiltered, y = &quot;subsets_Mt_percent&quot;, colour_by = &quot;discard&quot;) + ggtitle(&quot;Mito percent&quot;), ncol=2 ) 2.2 Identifying and removing low-quality cells After you have identified the low-quality cells, you can easily remove them in R. # creating a dataset with only the samples that pass QC sce.zeisel.qc &lt;- sce.zeisel[,!qc$discard] QUESTIONS How many samples did you remove as suspected low-quality cells? Do all the quality metrics (total count, number of expressed features, percentage of mitochondrial reads, and percentage of ERCC reads) agree when it comes to identifying low-quality samples? Which metric resulted in the removal of the greatest number of samples? Some researchers mark the suspected low-quality cells instead of removing them from the analysis. Ideally, the low-quality cells form their own cluster that can be ignored when interpreting the scRNA-seq results. This approach prevents the researcher from discarding cells or cell types that represent a true biological state and happen to have poor QC metrics. "],["normalization.html", "Chapter 3 Normalization 3.1 Calculating scaling factors 3.2 Applying scaling factors and computing normalized expression values", " Chapter 3 Normalization Single-cell RNA-seq libraries often vary in sequence coverage due to technical differences in cDNA capture or PCR amplification. We use normalization to correct for these differences. 3.1 Calculating scaling factors The simplest approach is to use library size as a scaling factor for each cell. One assumption underlying the library size normalization method is that library size is an unbiased estimate of capture or amplification issues. While this assumption is not necessarily biologically realistic, in practice the accuracy of the normalization step is unlikely to cause issues in exploratory scRNA-seq analyses. Sometimes unbalanced differential expression exists between samples, which results in composition biases in the dataset. Library size normalization can result in the non-differentially-expressed genes falsely appearing to be downregulated in one cell compared to another. If this sounds familiar, it’s because this topic has been well-studied for bulk RNA-seq analysis. Several bulk normalization methods have been created to deal with this problem in bulk RNA-seq analyses. In order to translate these algorithms to scRNA-seq and the large number of zero and low counts, we pool counts from many cells for accurate sizing factor estimation. These estimated sizing factors are then “deconvolved” into cell-based factors. This process is called normalization by deconvolution. # AnVIL::install(&quot;scran&quot;) library(scran) # calculating scaling factors set.seed(1000) clusters &lt;- quickCluster(sce.zeisel.qc) sce.zeisel.qc &lt;- computeSumFactors(sce.zeisel.qc, cluster = clusters) # comparing scaling factors based on library size (librarySizeFactors) and deconvolution (sizeFactors) plot(librarySizeFactors(sce.zeisel.qc), sizeFactors(sce.zeisel.qc), pch = 16, xlab = &quot;Library size factors&quot;, ylab = &quot;Deconvolution factors&quot;, log = &quot;xy&quot;) QUESTIONS How do the scaling factors using library size normalization compare to the scaling factors using deconvolution normalization? 3.2 Applying scaling factors and computing normalized expression values Once we have calculated scaling factors, we can compute normalized expression values for each sample by dividing the count for each gene by the size factor for the cell. The function we are using also log-transforms these normalized values, which will allow us to calculate log-fold changes in expression (a much more intuitive measure of differential expression!) This lets us focus on contributions from genes with strong relative differences. # computing normalized expression values sce.zeisel.qc &lt;- logNormCounts(sce.zeisel.qc) "],["choosing-features-for-the-analysis.html", "Chapter 4 Choosing features for the analysis 4.1 Examining per-gene variation 4.2 Selecting highly variable genes", " Chapter 4 Choosing features for the analysis The goal with many scRNA-seq analyses is to identify differences in gene expression across multiple cell types; in order to do this, cells are clustered together into groups based on a single similarity (or dissimilarity) value between a pair of cells. Which genes you choose to include in the calculation of the similarity metric can greatly influence your downstream analyses. In general, you will choose to include those genes with the greatest variability in gene expression across all the samples. The underlying assumption is this increased variability is “true” biological variation as opposed to technical noise. 4.1 Examining per-gene variation We can calculate the per-gene variation across all samples using the variance of the log-normalized expression values you computed in the Normalization section. Because we are using a log-transformation, the total variance of a gene is partly driven by its abundance (not underlying biological heterogeneity). To correct for this relationship, we use a model that calculates the expected amount of variation for each abundance value. We then use the difference between the observed variation and the expected variation for each gene to identify highly-variable genes while controlling for cell abundance. The Zeisel data was generated across multiple plates, so we need to be wary of any possible batch effects that could be driving highly-variable genes. Normally we would address this by blocking. However, each plate only contained 20-40 cells, and the cell population as a whole is highly heterogeneous, making it unlikely the sampled cell types on each plate is the same (one of the assumptions of blocking on plate). Thus, we will ignore blocking for this analysis. If you work with a dataset that does not contain spike-ins, you would use the modelGeneVar() command instead of modelGeneVarWithSpikes(). # calculating per-gene variation for each samples using the log-normalized expression value; the &quot;ERCC&quot; option tells the algorithm that information about technical (that is, non-biological) variation can be calculated using the data column labeled &quot;ERCC&quot; dec.zeisel.qc &lt;- modelGeneVarWithSpikes(sce.zeisel.qc, &quot;ERCC&quot;) # plotting the expected variance (blue curve and orange points), as well as the observed variance (black points) plot(dec.zeisel.qc$mean, dec.zeisel.qc$total, pch = 16, cex = 0.5, xlab = &quot;Mean of log-expression&quot;, ylab = &quot;Variance of log-expression&quot;) curfit &lt;- metadata(dec.zeisel.qc) points(curfit$mean, curfit$var, col = &quot;orange&quot;, pch = 16) curve(curfit$trend(x), col = &#39;dodgerblue&#39;, add = TRUE, lwd = 2) The expected variance for each log-expression value is plotted in orange. Observed values (the black points) much higher than this trend line will be used for downstream analyses. 4.2 Selecting highly variable genes When choosing the genes for downstream analysis, you will need to balance choosing as many genes as possible (so as not to exclude important biological variation) with limiting the amount of random noise. The most straightforward approach is to simply take the top n genes after ranking based on the biological component of variance. We are creating objects with the top 500, 1000, and 2000 genes, so that we can examine how our choice of n impacts downstream calculations and steps. # getting the top 500, 1000, and 2000 most-variable genes top500.hvgs &lt;- getTopHVGs(dec.zeisel.qc, n = 500) top1000.hvgs &lt;- getTopHVGs(dec.zeisel.qc, n = 1000) top2000.hvgs &lt;- getTopHVGs(dec.zeisel.qc, n = 2000) # look at the cut-offs dec.zeisel.qc &lt;- dec.zeisel.qc[order(dec.zeisel.qc$bio, decreasing = TRUE),] dec.zeisel.qc[c(1,500,1000,2000),1:6] ## DataFrame with 4 rows and 6 columns ## mean total tech bio p.value FDR ## &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Plp1 3.97379 16.80426 0.275699 16.528564 0.00000e+00 0.00000e+00 ## Shfm1 2.64467 1.55363 0.534672 1.018959 5.54040e-47 2.04392e-45 ## Cdc42 3.34431 1.04001 0.385442 0.654563 9.97095e-38 2.96239e-36 ## Vdac2 1.81697 1.08138 0.720031 0.361348 7.90101e-05 2.55550e-04 In this dataframe, the bio column represents the excess variation in gene expression (that is, the difference between the observed and expected expression). We are looking at four different genes - the gene with the greatest excess variation, as well as the last genes included in the top 500, top 1000, and top 2000 genes. QUESTIONS What is the range of log-fold expression changes (the excess variation) when you choose the top 500 genes? What about when you choose the top 1000 genes? The top 2000 genes? sessionInfo() ## R version 4.1.3 (2022-03-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.5 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] scran_1.22.1 scater_1.22.0 ## [3] ggplot2_3.3.5 scuttle_1.4.0 ## [5] scRNAseq_2.8.0 SingleCellExperiment_1.16.0 ## [7] SummarizedExperiment_1.24.0 Biobase_2.54.0 ## [9] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 ## [11] IRanges_2.28.0 S4Vectors_0.32.4 ## [13] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 ## [15] matrixStats_0.61.0 ## ## loaded via a namespace (and not attached): ## [1] AnnotationHub_3.2.2 BiocFileCache_2.2.1 ## [3] igraph_1.3.1 lazyeval_0.2.2 ## [5] BiocParallel_1.28.3 digest_0.6.29 ## [7] ensembldb_2.18.4 htmltools_0.5.2 ## [9] viridis_0.6.2 fansi_1.0.3 ## [11] magrittr_2.0.3 memoise_2.0.1 ## [13] ScaledMatrix_1.2.0 cluster_2.1.2 ## [15] limma_3.50.3 Biostrings_2.62.0 ## [17] prettyunits_1.1.1 colorspace_2.0-3 ## [19] blob_1.2.3 rappdirs_0.3.3 ## [21] ggrepel_0.9.1 xfun_0.26 ## [23] dplyr_1.0.8 crayon_1.5.1 ## [25] RCurl_1.98-1.6 jsonlite_1.8.0 ## [27] glue_1.6.2 gtable_0.3.0 ## [29] zlibbioc_1.40.0 XVector_0.34.0 ## [31] DelayedArray_0.20.0 BiocSingular_1.10.0 ## [33] scales_1.2.1 edgeR_3.36.0 ## [35] DBI_1.1.2 Rcpp_1.0.8.3 ## [37] viridisLite_0.4.0 xtable_1.8-4 ## [39] progress_1.2.2 dqrng_0.3.0 ## [41] bit_4.0.4 rsvd_1.0.5 ## [43] metapod_1.2.0 httr_1.4.2 ## [45] ellipsis_0.3.2 pkgconfig_2.0.3 ## [47] XML_3.99-0.9 farver_2.1.0 ## [49] sass_0.4.1 dbplyr_2.1.1 ## [51] locfit_1.5-9.5 utf8_1.2.2 ## [53] tidyselect_1.1.2 labeling_0.4.2 ## [55] rlang_1.0.2 later_1.3.0 ## [57] AnnotationDbi_1.56.2 munsell_0.5.0 ## [59] BiocVersion_3.14.0 tools_4.1.3 ## [61] cachem_1.0.6 cli_3.2.0 ## [63] generics_0.1.2 RSQLite_2.2.12 ## [65] ExperimentHub_2.2.1 evaluate_0.15 ## [67] stringr_1.4.0 fastmap_1.1.0 ## [69] yaml_2.3.5 knitr_1.33 ## [71] bit64_4.0.5 purrr_0.3.4 ## [73] KEGGREST_1.34.0 AnnotationFilter_1.18.0 ## [75] sparseMatrixStats_1.6.0 mime_0.12 ## [77] xml2_1.3.3 biomaRt_2.50.3 ## [79] compiler_4.1.3 beeswarm_0.4.0 ## [81] filelock_1.0.2 curl_4.3.2 ## [83] png_0.1-7 interactiveDisplayBase_1.32.0 ## [85] statmod_1.4.36 tibble_3.1.6 ## [87] bslib_0.3.1 stringi_1.7.6 ## [89] highr_0.9 GenomicFeatures_1.46.5 ## [91] lattice_0.20-45 bluster_1.4.0 ## [93] ProtGenerics_1.26.0 Matrix_1.4-0 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.16 ## [99] jquerylib_0.1.4 BiocNeighbors_1.12.0 ## [101] bitops_1.0-7 irlba_2.3.5 ## [103] httpuv_1.6.5 rtracklayer_1.54.0 ## [105] R6_2.5.1 BiocIO_1.4.0 ## [107] bookdown_0.24 promises_1.2.0.1 ## [109] gridExtra_2.3 vipor_0.4.5 ## [111] assertthat_0.2.1 rjson_0.2.21 ## [113] withr_2.5.0 GenomicAlignments_1.30.0 ## [115] Rsamtools_2.10.0 GenomeInfoDbData_1.2.7 ## [117] parallel_4.1.3 hms_1.1.1 ## [119] grid_4.1.3 beachmat_2.10.0 ## [121] rmarkdown_2.10 DelayedMatrixStats_1.16.0 ## [123] shiny_1.7.1 ggbeeswarm_0.6.0 ## [125] restfulr_0.0.13 "],["dimensionality-reduction.html", "Chapter 5 Dimensionality reduction 5.1 Calculating and choosing PCs 5.2 Applying non-linear visualization methods to PCs", " Chapter 5 Dimensionality reduction Each gene in the data represents a different dimension of the data. Reducing the number of dimensions in our data has multiple benefits, including reducing the computational work needed for downstream analyses. It also reduces the noise in the data through averaging the signal across multiple genes. Dimensionality reduction is a very common technique used in data science in general, not just in scRNA-seq analysis. You will find yourself using it over and over whenever you work with high-dimensional data. Dimensionality reduction is possible for genomic expression methods because so many genes have correlated expression. This is a consequence of different genes being involved in the same biological processes. 5.1 Calculating and choosing PCs We first use principal component analysis (PCA), which is a dimensionality reduction method that maximizes the amount of variation captured by each component, or PC. It’s up to the researcher to choose how many PCs to use for downstream analyses. More PCs mean that more biological signal is retained in the data, but it also increases the potential for noise. In our analyses, we will use 50 (the default). We will also use the top 1000 genes to calculate the PCs. library(scran) # calculating PCA # the denoisePCA command calculates PCs and removes those that primarily capture technical noise set.seed(101011001) sce.zeisel.PCA.1000 &lt;- denoisePCA(sce.zeisel.qc, technical = dec.zeisel.qc, subset.row = top1000.hvgs) #the technical option tells R where to find information about how much of the variation is attributed to &quot;technical&quot;, or non-biological, sources # visualizing the percentage of variation explained by each PC percent.var &lt;- attr(reducedDim(sce.zeisel.PCA.1000), &quot;percentVar&quot;) plot(percent.var, log = &quot;y&quot;, xlab = &quot;PC&quot;, ylab = &quot;Variance explained (%)&quot;) In PCA, the total amount of variation captured decreases for each subsequent PC. By the 10th PC, each additional PC is contributing only a small fraction to the total amount of variation explained in the dataset. Excluding them from downstream analyses has no major effect, and researchers will typically choose to include somewhere between 10 to 50 PCs. Including more PCs in the downstream analyses could theoretically cause the calculations to take longer, but in reality most calculations are fast enough that any slowdown isn’t really noticeable. More detailed information on calculating and choosing PCs for genomic analyses can be found in the Statistics for Genomics: PCA book. 5.2 Applying non-linear visualization methods to PCs In scRNA-seq analysis, plotting PCs generally does not offer enough resolution to visualize cell clusters. Instead, we rely on additional dimensionality reduction methods that can use non-linear data transformation algorithms. The most common approach is the t-stochastic neighbor embedding (t-SNE) method (Van der Maaten and Hinton 2008). t-SNE maps high-dimensional data in a low-dimensional space by first calculating the Euclidean distance between each set of points, then converting those distances into the probability that given pair of points are neighbors. On t-SNE plots, points that are members of the same cluster have a high probability of being neighbors. However, you can’t judge the similarity of different clusters based on their position on the final plot, because the t-SNE algorithm does not retain that information. The t-SNE approach is computationally complex. In practice, we reduce the computational complexity in scRNA-seq analysis by performing t-SNE calculations on the top PCs in a dataset (this both decreases the amount of computational power and time needed for analysis). We also need to set an initial starting seed and the perplexity parameter (the number of effective neighbors for each point). This parameter will determine the resolution of the plot. Lower perplexity values allow for finer resolution of population structure but can also be noisy. It’s a good idea to test multiple perplexity values when running your t-SNE analysis. We’re going to use the PCs calculated using the top 1000 highly-variable genes (a common threshold) for the rest of the analyses. library(BiocSingular) # this code first calculates the t-SNE values using PCs, and then creates a plot of the first two t-SNE dimensions set.seed(100) sce.zeisel.tsne5 &lt;- runTSNE(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;, perplexity = 5) out5 &lt;- plotReducedDim(sce.zeisel.tsne5, dimred = &quot;TSNE&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;perplexity = 5&quot;) set.seed(100) sce.zeisel.tsne20 &lt;- runTSNE(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;, perplexity = 20) out20 &lt;- plotReducedDim(sce.zeisel.tsne20, dimred = &quot;TSNE&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;perplexity = 20&quot;) set.seed(100) sce.zeisel.tsne80 &lt;- runTSNE(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;, perplexity = 80) out80 &lt;- plotReducedDim(sce.zeisel.tsne80, dimred = &quot;TSNE&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;perplexity = 80&quot;) out5 out20 out80 Some researchers will use a different method for the non-linear visualization step in their analysis. This algorithm, called uniform manifold approximation and projection (UMAP), is faster and preserves more of the global data structure when reducing dimensions compared to t-SNE (that is, you get more information about the similarity of clusters, not just of points). As a result, though, the resolution within each cluster is reduced. UMAP is becoming the method of choice as scRNA-seq datasets become larger and larger. # calculate the UMAP values from the PCs, then plot the first two UMAP dimensions set.seed(1100101001) sce.zeisel.umap &lt;- runUMAP(sce.zeisel.PCA.1000, dimred = &quot;PCA&quot;) out.umap &lt;- plotReducedDim(sce.zeisel.umap, dimred = &quot;UMAP&quot;, colour_by = &quot;level1class&quot;) + ggtitle(&quot;UMAP&quot;) out20 out.umap QUESTIONS 1. How does changing the perplexity parameter affect the t-SNE plot? How does the t-SNE plot compare to the UMAP plot? "],["clustering.html", "Chapter 6 Clustering 6.1 Clustering using graph-based methods", " Chapter 6 Clustering At its core, clustering is a tool that allows us to examine structure and patterns in our data. There isn’t really a “true” answer to how the data should be clustered. Instead, we can change the algorithms and parameters to explore a variety of possibilities that work best for each dataset and question the researcher is trying to answer. 6.1 Clustering using graph-based methods Graph-based clustering is based on identifying the nearest neighbors of each cell in high-dimensional space. The connections between a cell and its neighbors (called edges) are weighted based on the similarity of the two cells connected. An edge is assigned a higher weight if the two cells it connects are more closely related. After all cells have been connected to their neighbors, we apply an algorithm to identify clusters, or communities, of related cells. Each cell within a community will be more closely related to any cell within the same community than to cells outside the community. Graph-based clustering scales easily, because it only used a k-nearest neighbor search. These searches run more quickly than other methods (like hierarchical clustering). Unfortunately, no information is retained about relationships beyond the neighboring cells. This effect also means that clustering resolution depends on cell density. # let the R algorithm define and label our clusters nn.clusters &lt;- clusterCells(sce.zeisel.tsne20, use.dimred=&quot;PCA&quot;) # this command tells us how many clusters were identified (the top row) and how many cells belong to each cluster (the bottom row) table(nn.clusters) ## nn.clusters ## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ## 284 147 107 195 604 166 475 62 119 256 221 67 53 28 32 We assigned the cluster assignments back into our SingleCellExperiment object as a factor in the column metadata, which allows us to visualize the cluster assignment on a t-SNE plot. # create a t-SNE plot showing the identified clusters colLabels(sce.zeisel.tsne20) &lt;- nn.clusters plotReducedDim(sce.zeisel.tsne20, &quot;TSNE&quot;, colour_by=&quot;label&quot;) QUESTIONS 1. How many clusters were identified using graph-based clustering? Which cluster contained the most cells, and how many cells did it have? sessionInfo() ## R version 4.1.3 (2022-03-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.5 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] uwot_0.1.11 Matrix_1.4-0 ## [3] BiocSingular_1.10.0 scran_1.22.1 ## [5] scater_1.22.0 ggplot2_3.3.5 ## [7] scuttle_1.4.0 scRNAseq_2.8.0 ## [9] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0 ## [11] Biobase_2.54.0 GenomicRanges_1.46.1 ## [13] GenomeInfoDb_1.30.1 IRanges_2.28.0 ## [15] S4Vectors_0.32.4 BiocGenerics_0.40.0 ## [17] MatrixGenerics_1.6.0 matrixStats_0.61.0 ## ## loaded via a namespace (and not attached): ## [1] AnnotationHub_3.2.2 BiocFileCache_2.2.1 ## [3] igraph_1.3.1 lazyeval_0.2.2 ## [5] BiocParallel_1.28.3 digest_0.6.29 ## [7] ensembldb_2.18.4 htmltools_0.5.2 ## [9] viridis_0.6.2 fansi_1.0.3 ## [11] magrittr_2.0.3 memoise_2.0.1 ## [13] ScaledMatrix_1.2.0 cluster_2.1.2 ## [15] limma_3.50.3 Biostrings_2.62.0 ## [17] prettyunits_1.1.1 colorspace_2.0-3 ## [19] blob_1.2.3 rappdirs_0.3.3 ## [21] ggrepel_0.9.1 xfun_0.26 ## [23] dplyr_1.0.8 crayon_1.5.1 ## [25] RCurl_1.98-1.6 jsonlite_1.8.0 ## [27] glue_1.6.2 gtable_0.3.0 ## [29] zlibbioc_1.40.0 XVector_0.34.0 ## [31] DelayedArray_0.20.0 scales_1.2.1 ## [33] DBI_1.1.2 edgeR_3.36.0 ## [35] Rcpp_1.0.8.3 viridisLite_0.4.0 ## [37] xtable_1.8-4 progress_1.2.2 ## [39] dqrng_0.3.0 bit_4.0.4 ## [41] rsvd_1.0.5 metapod_1.2.0 ## [43] httr_1.4.2 FNN_1.1.3 ## [45] ellipsis_0.3.2 farver_2.1.0 ## [47] pkgconfig_2.0.3 XML_3.99-0.9 ## [49] sass_0.4.1 dbplyr_2.1.1 ## [51] locfit_1.5-9.5 utf8_1.2.2 ## [53] labeling_0.4.2 tidyselect_1.1.2 ## [55] rlang_1.0.2 later_1.3.0 ## [57] AnnotationDbi_1.56.2 munsell_0.5.0 ## [59] BiocVersion_3.14.0 tools_4.1.3 ## [61] cachem_1.0.6 cli_3.2.0 ## [63] generics_0.1.2 RSQLite_2.2.12 ## [65] ExperimentHub_2.2.1 evaluate_0.15 ## [67] stringr_1.4.0 fastmap_1.1.0 ## [69] yaml_2.3.5 knitr_1.33 ## [71] bit64_4.0.5 purrr_0.3.4 ## [73] KEGGREST_1.34.0 AnnotationFilter_1.18.0 ## [75] sparseMatrixStats_1.6.0 mime_0.12 ## [77] xml2_1.3.3 biomaRt_2.50.3 ## [79] compiler_4.1.3 beeswarm_0.4.0 ## [81] filelock_1.0.2 curl_4.3.2 ## [83] png_0.1-7 interactiveDisplayBase_1.32.0 ## [85] statmod_1.4.36 tibble_3.1.6 ## [87] bslib_0.3.1 stringi_1.7.6 ## [89] highr_0.9 RSpectra_0.16-0 ## [91] GenomicFeatures_1.46.5 lattice_0.20-45 ## [93] bluster_1.4.0 ProtGenerics_1.26.0 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.16 ## [99] jquerylib_0.1.4 BiocNeighbors_1.12.0 ## [101] bitops_1.0-7 irlba_2.3.5 ## [103] httpuv_1.6.5 rtracklayer_1.54.0 ## [105] R6_2.5.1 BiocIO_1.4.0 ## [107] bookdown_0.24 promises_1.2.0.1 ## [109] gridExtra_2.3 vipor_0.4.5 ## [111] assertthat_0.2.1 rjson_0.2.21 ## [113] withr_2.5.0 GenomicAlignments_1.30.0 ## [115] Rsamtools_2.10.0 GenomeInfoDbData_1.2.7 ## [117] parallel_4.1.3 hms_1.1.1 ## [119] grid_4.1.3 beachmat_2.10.0 ## [121] rmarkdown_2.10 DelayedMatrixStats_1.16.0 ## [123] Rtsne_0.16 shiny_1.7.1 ## [125] ggbeeswarm_0.6.0 restfulr_0.0.13 "],["marker-gene-detection.html", "Chapter 7 Marker gene detection 7.1 Calculating and ranking effect size summary statistics 7.2 Comparing gene expression levels across clusters", " Chapter 7 Marker gene detection We have managed to identify clusters of related cells in the data. However, these clusters aren’t very useful until we can identify the biology meaning of each group. This is where functional annotation comes in. The real art, and the greatest challenge, in an scRNA-seq analysis comes when interpreting the results. Up to this point (cleaning and clustering the data), the analysis and computation has been straightforward. Figuring out the biological state that each cluster represents, on the other hand, is more difficult, as it requires applying prior biological knowledge to the dataset. Thanks to previous research, we know many marker genes, or genes can be used to identify particular cell types. These genes are differentially expressed across cell types, and by examining the expression profiles of multiple marker genes across all the clusters, we can assign particular cell type identities to each cluster. 7.1 Calculating and ranking effect size summary statistics We begin by comparing each pair of clusters and calculating scores for expression differences between the two for each gene. We have multiple options for the statistics used to compare expression values. AUC (area under the curve) is the probability that a randomly chosen observation from cluster A is greater than a randomly chosen observation from cluster B. This statistic is a way to quantify how well we can distinguish between two distributions (clusters) in a pairwise comparison. An AUC of 1 means all values in cluster A are greater than any value from cluster B and suggests upregulation. An AUC of of 0.5 means the two clusters are indistinguishable from each other, while an AUC of 0 suggests the marker gene observations in cluster A are downregulated compared to those in cluster B. Cohen’s d is a standardized log-fold change, and can be thought of as the number of standard deviations that separate the two groups. Positive values of Cohen’s d suggest that our cluster of interest (cluster A) are upregulated compared to cluster B, while negative values suggest the marker gene observations in cluster A are downregulated compared to cluster B. log-fold change (logFC) is a measure of whether there is a difference in expression between clusters. Keep in mind that these values ignore the magnitude of the change. As with the others, positive values indicate upregulation in the cluster of interest (cluster A), while negative values indicate downregulation. For each of these statistics, scoreMarkers calculates mean, median, minimum value (min), maximum value (max), and minimum rank (rank; the smallest rank of each gene across all pairwise comparisons). For most of these measures, a larger number indicates upregulation. For minimum rank, however, a small value means the gene is one of the top upregulated genes. AUC or Cohen’s d are effective regardless of the magnitude of the expression values and thus are good choices for general marker detection. The log-fold change in the detected proportion is specifically useful for identifying binary changes in expression. For this exercise, we’re going to focus on upregulated markers, since those are particularly useful for identifying cell types in a heterogeneous population like the Zeisel dataset. We use the findMarkers command, which quickly identifies potential marker genes. # identify those genes which are upregulated in some clusters compared to others markers &lt;- findMarkers(sce.zeisel.tsne20, direction = &quot;up&quot;) marker.set &lt;- markers[[&quot;1&quot;]] head(marker.set, 10) ## DataFrame with 10 rows and 18 columns ## Top p.value FDR summary.logFC logFC.2 logFC.3 ## &lt;integer&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Syngr3 1 5.71433e-140 1.90535e-137 2.44603 1.156187 1.1330628 ## Mllt11 1 3.31857e-249 7.37681e-246 2.87647 0.971725 1.4504728 ## Ndrg4 1 0.00000e+00 0.00000e+00 3.83936 1.146289 0.9309904 ## Slc6a1 1 1.75709e-160 8.57373e-158 3.46313 3.063161 3.2821288 ## Gad1 1 4.70797e-234 8.56251e-231 4.54282 4.024699 4.3901265 ## Gad2 1 1.32400e-208 1.76586e-205 4.24931 3.756435 3.9505814 ## Atp1a3 1 3.54623e-279 1.41892e-275 3.45285 1.213985 0.0213418 ## Slc32a1 2 1.50335e-109 2.48562e-107 1.91750 1.791641 1.8496954 ## Rcan2 2 9.25451e-129 2.50197e-126 2.21812 1.432185 2.1127708 ## Rab3a 2 1.61475e-202 1.61523e-199 2.52266 0.884615 0.6329209 ## logFC.4 logFC.5 logFC.6 logFC.7 logFC.8 logFC.9 logFC.10 ## &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Syngr3 2.38928 2.44603 0.7085503 0.8470822 1.159959 2.525116 1.030908 ## Mllt11 2.94289 2.87647 0.4897958 0.0882716 1.030983 3.262597 0.471414 ## Ndrg4 3.67150 3.83936 0.3680955 0.7169171 1.103046 3.646328 0.791827 ## Slc6a1 2.92969 2.88778 3.6063318 3.5292982 2.218826 0.704948 3.571574 ## Gad1 4.37025 4.46268 4.6478826 4.5428243 4.477666 4.386587 4.572889 ## Gad2 4.08621 4.15971 4.2916168 4.2222656 4.229715 4.207704 4.249309 ## Atp1a3 3.18678 3.45285 0.6498389 0.1656598 0.426762 3.436859 -0.228370 ## Slc32a1 1.93798 1.91658 1.9393715 1.9092434 1.913338 1.917505 1.910826 ## Rcan2 1.81884 1.25488 1.1051970 2.1730716 2.170402 2.102633 2.218119 ## Rab3a 2.73747 2.52266 0.0267138 0.3133349 0.695216 2.921335 0.163635 ## logFC.11 logFC.12 logFC.13 logFC.14 logFC.15 ## &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; &lt;numeric&gt; ## Syngr3 0.684417 2.56482 2.55195 2.58830 2.42180 ## Mllt11 0.229326 3.04261 2.77725 3.08889 2.01497 ## Ndrg4 0.565654 3.82140 4.14743 3.97144 4.01655 ## Slc6a1 3.463131 3.43495 3.34170 3.05042 3.25178 ## Gad1 4.568481 4.51547 4.50591 4.39482 4.63836 ## Gad2 4.263975 4.16831 4.29808 4.09088 4.35014 ## Atp1a3 1.118434 3.31854 3.78426 3.67087 3.60817 ## Slc32a1 1.916869 1.94707 1.94707 1.94707 1.90680 ## Rcan2 1.142480 2.09515 1.26318 1.46319 2.24552 ## Rab3a 0.508261 3.16930 2.54753 3.06401 2.30063 This dataframe shows us the in log-fold expression change for each potential marker gene between cluster 1 and every other cluster. 7.2 Comparing gene expression levels across clusters Once we’ve identified potential marker genes, we can use a heatmap to compare gene expression in each cell between clusters. # pull the top-most upregulated markers from cluster 1 (compared to the rest of the clusters) and look at their expression in all clusters top.markers &lt;- rownames(marker.set)[marker.set$Top &lt;= 10] plotHeatmap(sce.zeisel.tsne20, features = top.markers, order_columns_by = &quot;label&quot;) In this heatmap, clusters are on the horizontal, while the top upregulated genes in cluster 1 are on the vertical. The magnitude of the log-fold expression change is indicated by color of each cell. We can also create a heatmap that shows the mean log-fold change of cluster 1 cells compared to the mean of each other cluster. This can simplify the heatmap and is useful when dealing with many clusters. # AnVIL::install(&quot;pheatmap&quot;) library(pheatmap) #this heatmap lets us compare the average expression of the gene within a cluster compared to the other clusters logFCs &lt;- getMarkerEffects(marker.set[1:50,]) pheatmap(logFCs, breaks = seq(-5, 5, length.out = 101)) Here we see that three genes are generally upregulated in Cluster 1 compared to the other clusters: Gad1, Gad2, and Slc6a1. This is where prior biological knowledge comes in handy, as both Gad1 and Slc6a1 are known interneuron markers (Zeng et al. 2012). QUESTION 1. Are there any groups or patterns you see in the second heatmap that look interesting? "],["cell-type-annotation.html", "Chapter 8 Cell type annotation 8.1 Annotating cell types", " Chapter 8 Cell type annotation Unless you are already an expert in neuronal cell expression, you probably didn’t know Gad1 and Slc6a1 are known interneuron markers until you were told. Luckily, we have several high-quality references databases that can be used for annotating scRNA-seq datasets. Some of these references are useful for identifying cell types (the Zeisel dataset we have been using is a well-known reference for identifying neuronal cell types). Others, such as the Gene Ontology (GO) or Kyoto Encyclopedia of Genes and Genomes (KEGG) collections, are useful for identifying biological processes associated with each cluster. 8.1 Annotating cell types There are three basic strategies for annotating datasets: match the expression profile of each individual cell to the expression profile of cells from a reference dataset (the reference dataset approach); identify sets of marker genes highly expressed in each cell and match to gene sets from known cell types (the gene set approach); or perform a gene-set enrichment analysis on the marker genes that define each cluster. The Zeisel dataset we have been working with has actually been annotated all this time. # calculate the top marker genes assigned to each cell type (level1class) in the Zeisel dataset wilcox.z &lt;- pairwiseWilcox(sce.zeisel.qc, sce.zeisel.qc$level1class, lfc = 1, direction = &quot;up&quot;) markers.z &lt;- getTopMarkers(wilcox.z$statistics, wilcox.z$pairs, pairwise = FALSE, n = 50) # look at how many cell-type categories there are, as well as how many cells assigned to each category lengths(markers.z) ## astrocytes_ependymal endothelial-mural interneurons ## 78 83 119 ## microglia oligodendrocytes pyramidal CA1 ## 69 80 124 ## pyramidal SS ## 147 This dataset has grouped the individual cells into 7 different categories of neuronal subtypes. We can use the gene sets that define these categories to annotate a second brain scRNA-seq dataset from Tasic et al. (2016). (We will not go through all the data cleaning and clustering steps for this new dataset - having made it this far, we trust you can do that!) #load the new dataset sce.tasic &lt;- TasicBrainData() We first create the gene set lists using the GSEABase package. The AUCell package identifies and ranks marker sets highly expressed in each cell using an area under the curve (AUC) approach. We can assign cell type identity in the Tasic dataset by taking the marker set with the highest AUC as the label for the cell. # AnVIL::install(&quot;GSEABase&quot;) library(GSEABase) #create a dataset that contains just the information about the Zeisel cell-type categories and the marker genes that define each cell-type all.sets &lt;- lapply(names(markers.z), function(x) { GeneSet(markers.z[[x]], setName = x) }) all.sets &lt;- GeneSetCollection(all.sets) # AnVIL::install(&quot;AUCell&quot;) library(AUCell) # rank genes by expression values within each cell rankings &lt;- AUCell_buildRankings(counts(sce.tasic), plotStats = FALSE, verbose = FALSE) # calculate AUC for each previously-defined marker set (from Zeisel) in the Tasic data cell.aucs &lt;- AUCell_calcAUC(all.sets, rankings) results &lt;- t(assay(cell.aucs)) After assigning cell type identities to each cluster, a researcher should always verify that the identities make sense. Since the Tasic dataset has also been annotated, we can compare our annotation to the researcher-provided annotation as a sort of sanity check. # assign cell type identity in the Tasic dataset by assumig the marker set with the top AUC is the proper label new.labels &lt;- colnames(results)[max.col(results)] # compare our annotations to the annotations provided by the Tasic dataset tab &lt;- table(new.labels, sce.tasic$broad_type) tab ## ## new.labels Astrocyte Endothelial Cell GABA-ergic Neuron ## astrocytes_ependymal 43 2 0 ## endothelial-mural 0 27 0 ## interneurons 0 0 759 ## microglia 0 0 0 ## oligodendrocytes 0 0 1 ## pyramidal SS 0 0 1 ## ## new.labels Glutamatergic Neuron Microglia Oligodendrocyte ## astrocytes_ependymal 0 0 0 ## endothelial-mural 0 0 0 ## interneurons 2 0 0 ## microglia 0 22 0 ## oligodendrocytes 0 0 38 ## pyramidal SS 810 0 0 ## ## new.labels Oligodendrocyte Precursor Cell Unclassified ## astrocytes_ependymal 19 4 ## endothelial-mural 0 2 ## interneurons 0 15 ## microglia 0 1 ## oligodendrocytes 3 0 ## pyramidal SS 0 60 As you can see, the labels applied by the researchers and by our annotation mostly match. (Pyramidal SS nerves are primarily glutamatergic, so although the categories are labeled differently, those two categories do indeed match!) The major exception is the oligodendrocyte precursor cells, which our annotation called astrocytes. However, this mismatch isn’t as concerning as you might think, once you know that both astrocytes and oligodendrocytes come from the same precursor cell lineage. QUESTION 1. Should we be concerned when 1 or 2 cells are assigned a different cell type identity by our annotation than by the researchers? Why or why not? sessionInfo() ## R version 4.1.3 (2022-03-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.5 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] stats4 stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] AUCell_1.16.0 GSEABase_1.56.0 ## [3] graph_1.72.0 annotate_1.72.0 ## [5] XML_3.99-0.9 AnnotationDbi_1.56.2 ## [7] pheatmap_1.0.12 BiocSingular_1.10.0 ## [9] scran_1.22.1 scater_1.22.0 ## [11] ggplot2_3.3.5 scuttle_1.4.0 ## [13] scRNAseq_2.8.0 SingleCellExperiment_1.16.0 ## [15] SummarizedExperiment_1.24.0 Biobase_2.54.0 ## [17] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 ## [19] IRanges_2.28.0 S4Vectors_0.32.4 ## [21] BiocGenerics_0.40.0 MatrixGenerics_1.6.0 ## [23] matrixStats_0.61.0 ## ## loaded via a namespace (and not attached): ## [1] AnnotationHub_3.2.2 BiocFileCache_2.2.1 ## [3] igraph_1.3.1 lazyeval_0.2.2 ## [5] BiocParallel_1.28.3 digest_0.6.29 ## [7] ensembldb_2.18.4 htmltools_0.5.2 ## [9] viridis_0.6.2 fansi_1.0.3 ## [11] magrittr_2.0.3 memoise_2.0.1 ## [13] ScaledMatrix_1.2.0 cluster_2.1.2 ## [15] limma_3.50.3 Biostrings_2.62.0 ## [17] R.utils_2.12.2 prettyunits_1.1.1 ## [19] colorspace_2.0-3 blob_1.2.3 ## [21] rappdirs_0.3.3 ggrepel_0.9.1 ## [23] xfun_0.26 dplyr_1.0.8 ## [25] crayon_1.5.1 RCurl_1.98-1.6 ## [27] jsonlite_1.8.0 glue_1.6.2 ## [29] gtable_0.3.0 zlibbioc_1.40.0 ## [31] XVector_0.34.0 DelayedArray_0.20.0 ## [33] scales_1.2.1 DBI_1.1.2 ## [35] edgeR_3.36.0 Rcpp_1.0.8.3 ## [37] viridisLite_0.4.0 xtable_1.8-4 ## [39] progress_1.2.2 dqrng_0.3.0 ## [41] bit_4.0.4 rsvd_1.0.5 ## [43] metapod_1.2.0 httr_1.4.2 ## [45] RColorBrewer_1.1-3 ellipsis_0.3.2 ## [47] R.methodsS3_1.8.1 farver_2.1.0 ## [49] pkgconfig_2.0.3 sass_0.4.1 ## [51] dbplyr_2.1.1 locfit_1.5-9.5 ## [53] utf8_1.2.2 tidyselect_1.1.2 ## [55] rlang_1.0.2 later_1.3.0 ## [57] munsell_0.5.0 BiocVersion_3.14.0 ## [59] tools_4.1.3 cachem_1.0.6 ## [61] cli_3.2.0 generics_0.1.2 ## [63] RSQLite_2.2.12 ExperimentHub_2.2.1 ## [65] evaluate_0.15 stringr_1.4.0 ## [67] fastmap_1.1.0 yaml_2.3.5 ## [69] knitr_1.33 bit64_4.0.5 ## [71] purrr_0.3.4 KEGGREST_1.34.0 ## [73] AnnotationFilter_1.18.0 sparseMatrixStats_1.6.0 ## [75] mime_0.12 R.oo_1.24.0 ## [77] xml2_1.3.3 biomaRt_2.50.3 ## [79] compiler_4.1.3 beeswarm_0.4.0 ## [81] filelock_1.0.2 curl_4.3.2 ## [83] png_0.1-7 interactiveDisplayBase_1.32.0 ## [85] statmod_1.4.36 tibble_3.1.6 ## [87] bslib_0.3.1 stringi_1.7.6 ## [89] highr_0.9 GenomicFeatures_1.46.5 ## [91] lattice_0.20-45 bluster_1.4.0 ## [93] ProtGenerics_1.26.0 Matrix_1.4-0 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.16 ## [99] jquerylib_0.1.4 BiocNeighbors_1.12.0 ## [101] data.table_1.14.2 bitops_1.0-7 ## [103] irlba_2.3.5 httpuv_1.6.5 ## [105] rtracklayer_1.54.0 R6_2.5.1 ## [107] BiocIO_1.4.0 bookdown_0.24 ## [109] promises_1.2.0.1 gridExtra_2.3 ## [111] vipor_0.4.5 assertthat_0.2.1 ## [113] rjson_0.2.21 withr_2.5.0 ## [115] GenomicAlignments_1.30.0 Rsamtools_2.10.0 ## [117] GenomeInfoDbData_1.2.7 parallel_4.1.3 ## [119] hms_1.1.1 grid_4.1.3 ## [121] beachmat_2.10.0 rmarkdown_2.10 ## [123] DelayedMatrixStats_1.16.0 Rtsne_0.16 ## [125] shiny_1.7.1 ggbeeswarm_0.6.0 ## [127] restfulr_0.0.13 "],["references.html", "Chapter 9 References", " Chapter 9 References Tasic B et al. Adult mouse cortical cell taxonomy revealed by single cell transcriptomics. Nat Neurosci. 2016 Feb;19(2):335-46. doi: 10.1038/nn.4216. Epub 2016 Jan 4. PMID: 26727548; PMCID: PMC4985242. Van der Maaten L, Hinton G. Visualizing data using t-SNE. J. Mach. Learn. Res. 2008 Nov;9:2579-2605. Zeisel A et al. Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science. 2015 Mar 6;347(6226):1138-42. doi: 10.1126/science.aaa1934. Epub 2015 Feb 19. PMID: 25700174. Zeng H et al. Large-scale cellular-resolution gene profiling in human neocortex reveals species-specific molecular signatures. Cell. 2012 Apr 13;149(2):483-96. doi: 10.1016/j.cell.2012.02.052. PMID: 22500809; PMCID: PMC3328777. "],["about-the-authors.html", "About the Authors", " About the Authors These credits are based on our course contributors table guidelines.     Credits Names Pedagogy Lead Content Instructor Elizabeth Humphries Content Idea Contributor Ava Hoffman, Frederick Tan Content Directors Kasper Hansen, Jeff Leek, Frederick Tan Production Content Publisher Ira Gooding Technical Template Publishing Engineers Candace Savonen, Carrie Wright, Ava Hoffman Publishing Maintenance Engineer Candace Savonen Technical Publishing Stylists Carrie Wright, Candace Savonen Package Developers (ottrpal) John Muschelli, Candace Savonen, Carrie Wright Funding Funder National Human Genome Research Institute (NHGRI) Funding Staff Fallon Bachman, Jennifer Vessio, Emily Voeglein   ## ─ Session info ─────────────────────────────────────────────────────────────── ## setting value ## version R version 4.1.3 (2022-03-10) ## os Ubuntu 20.04.5 LTS ## system x86_64, linux-gnu ## ui X11 ## language (EN) ## collate en_US.UTF-8 ## ctype en_US.UTF-8 ## tz Etc/UTC ## date 2023-05-25 ## pandoc 2.17.1.1 @ /usr/local/bin/ (via rmarkdown) ## ## ─ Packages ─────────────────────────────────────────────────────────────────── ## package * version date (UTC) lib source ## bookdown 0.24 2023-03-22 [1] Github (rstudio/bookdown@88bc4ea) ## bslib 0.3.1 2021-10-06 [1] RSPM (R 4.1.0) ## cachem 1.0.6 2021-08-19 [1] RSPM (R 4.1.0) ## callr 3.7.0 2021-04-20 [1] RSPM (R 4.1.0) ## cli 3.2.0 2022-02-14 [1] RSPM (R 4.1.0) ## crayon 1.5.1 2022-03-26 [1] RSPM (R 4.1.0) ## desc 1.4.1 2022-03-06 [1] RSPM (R 4.1.0) ## devtools 2.4.3 2021-11-30 [1] RSPM (R 4.1.0) ## digest 0.6.29 2021-12-01 [1] RSPM (R 4.1.0) ## ellipsis 0.3.2 2021-04-29 [1] RSPM (R 4.1.0) ## evaluate 0.15 2022-02-18 [1] RSPM (R 4.1.0) ## fastmap 1.1.0 2021-01-25 [1] RSPM (R 4.1.0) ## fs 1.5.2 2021-12-08 [1] RSPM (R 4.1.0) ## glue 1.6.2 2022-02-24 [1] RSPM (R 4.1.0) ## htmltools 0.5.2 2021-08-25 [1] RSPM (R 4.1.0) ## jquerylib 0.1.4 2021-04-26 [1] RSPM (R 4.1.0) ## jsonlite 1.8.0 2022-02-22 [1] RSPM (R 4.1.0) ## knitr 1.33 2023-03-22 [1] Github (yihui/knitr@a1052d1) ## lifecycle 1.0.1 2021-09-24 [1] RSPM (R 4.1.0) ## magrittr 2.0.3 2022-03-30 [1] RSPM (R 4.1.0) ## memoise 2.0.1 2021-11-26 [1] RSPM (R 4.1.0) ## pkgbuild 1.3.1 2021-12-20 [1] RSPM (R 4.1.0) ## pkgload 1.2.4 2021-11-30 [1] RSPM (R 4.1.0) ## prettyunits 1.1.1 2020-01-24 [1] RSPM (R 4.1.0) ## processx 3.5.3 2022-03-25 [1] RSPM (R 4.1.0) ## ps 1.6.0 2021-02-28 [1] RSPM (R 4.1.0) ## purrr 0.3.4 2020-04-17 [1] RSPM (R 4.1.0) ## R6 2.5.1 2021-08-19 [1] RSPM (R 4.1.0) ## remotes 2.4.2 2021-11-30 [1] RSPM (R 4.1.0) ## rlang 1.0.2 2023-03-22 [1] Github (r-lib/rlang@34b04a8) ## rmarkdown 2.10 2023-03-22 [1] Github (rstudio/rmarkdown@02d3c25) ## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.1.3) ## sass 0.4.1 2022-03-23 [1] RSPM (R 4.1.0) ## sessioninfo 1.2.2 2021-12-06 [1] RSPM (R 4.1.0) ## stringi 1.7.6 2021-11-29 [1] RSPM (R 4.1.0) ## stringr 1.4.0 2019-02-10 [1] RSPM (R 4.1.0) ## testthat 3.0.1 2023-03-22 [1] Github (R-lib/testthat@e99155a) ## usethis 2.1.6.9000 2023-03-22 [1] Github (r-lib/usethis@b4b8623) ## withr 2.5.0 2022-03-03 [1] RSPM (R 4.1.0) ## xfun 0.26 2023-03-22 [1] Github (yihui/xfun@74c2a66) ## yaml 2.3.5 2022-02-21 [1] RSPM (R 4.1.0) ## ## [1] /usr/local/lib/R/site-library ## [2] /usr/local/lib/R/library ## ## ────────────────────────────────────────────────────────────────────────────── "],["404.html", "Page not found", " Page not found The page you requested cannot be found (perhaps it was moved or renamed). You may want to try searching to find the page's new location, or use the table of contents to find the page you are looking for. "]]