<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Differential abundance testing with Milo - Mouse gastrulation example • miloR</title>
<!-- jquery --><script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.4.1/jquery.min.js" integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script><!-- Bootstrap --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/css/bootstrap.min.css" integrity="sha256-bZLfwXAP04zRMK2BjiO8iu9pf4FbLqX6zitd+tIvLhE=" crossorigin="anonymous">
<script src="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/js/bootstrap.min.js" integrity="sha256-nuL8/2cJ5NDSSwnKD8VqreErSWHtnEP9E7AySL+1ev4=" crossorigin="anonymous"></script><!-- bootstrap-toc --><link rel="stylesheet" href="../bootstrap-toc.css">
<script src="../bootstrap-toc.js"></script><!-- Font Awesome icons --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous">
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous">
<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><!-- headroom.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../pkgdown.js"></script><meta property="og:title" content="Differential abundance testing with Milo - Mouse gastrulation example">
<meta property="og:description" content="miloR">
<!-- mathjax --><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body data-spy="scroll" data-target="#toc">
    <div class="container template-article">
      <header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
        <span class="sr-only">Toggle navigation</span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <span class="navbar-brand">
        <a class="navbar-link" href="../index.html">miloR</a>
        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.1.0</span>
      </span>
    </div>

    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
<li>
  <a href="../index.html">
    <span class="fas fa fas fa-home fa-lg"></span>
     
  </a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
    Articles
     
    <span class="caret"></span>
  </a>
  <ul class="dropdown-menu" role="menu">
<li>
      <a href="../articles/milo_demo.html">Differential abundance testing with Milo</a>
    </li>
    <li>
      <a href="../articles/milo_gastrulation.html">Differential abundance testing with Milo - Mouse gastrulation example</a>
    </li>
  </ul>
</li>
      </ul>
<ul class="nav navbar-nav navbar-right"></ul>
</div>
<!--/.nav-collapse -->
  </div>
<!--/.container -->
</div>
<!--/.navbar -->

      

      </header><script src="milo_gastrulation_files/header-attrs-2.5/header-attrs.js"></script><script src="milo_gastrulation_files/accessible-code-block-0.0.1/empty-anchor.js"></script><link href="milo_gastrulation_files/anchor-sections-1.0/anchor-sections.css" rel="stylesheet">
<script src="milo_gastrulation_files/anchor-sections-1.0/anchor-sections.js"></script><div class="row">
  <div class="col-md-9 contents">
    <div class="page-header toc-ignore">
      <h1 data-toc-skip>Differential abundance testing with Milo - Mouse gastrulation example</h1>
                        <h4 class="author">Emma Dann</h4>
                        <h4 class="author">Mike Morgan</h4>
            
      
      
      <div class="hidden name"><code>milo_gastrulation.Rmd</code></div>

    </div>

    
    
<div class="sourceCode" id="cb1"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">miloR</span>)
<span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">SingleCellExperiment</span>)
<span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">scater</span>)
<span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">dplyr</span>)
<span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">patchwork</span>)</pre></body></html></div>
<div id="load-data" class="section level1">
<h1 class="hasAnchor">
<a href="#load-data" class="anchor"></a>Load data</h1>
<p>For this vignette we will use the mouse gastrulation single-cell data from <a href="https://www.nature.com/articles/s41586-019-0933-9">Pijuan-Sala et al. 2019</a>. The dataset can be downloaded as a <code>SingleCellExperiment</code> object from the <a href="https://bioconductor.org/packages/3.12/data/experiment/html/MouseGastrulationData.html"><code>MouseGastrulationData</code></a> package on Bioconductor. To make computations faster, here we will download just a subset of samples, 4 samples at stage E7 and 4 samples at stage E7.5.</p>
<p>This dataset has already been pre-processed and contains a <code>pca.corrected</code> dimensionality reduction, which was built after batch correction using <a href="https://bioconductor.org/packages/release/bioc/vignettes/batchelor/inst/doc/correction.html"><code>fastMNN</code></a>.</p>
<div class="sourceCode" id="cb2"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">MouseGastrulationData</span>)
<span class="no">select_samples</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fl">2</span>,  <span class="fl">3</span>,  <span class="fl">6</span>, <span class="fl">15</span>,
                    <span class="co"># 4, 19, </span>
                    <span class="fl">10</span>, <span class="fl">14</span>, <span class="fl">20</span>, <span class="fl">30</span>
                    <span class="co">#31, 32</span>
                    )
<span class="no">embryo_data</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/pkg/MouseGastrulationData/man/EmbryoAtlasData.html">EmbryoAtlasData</a></span>(<span class="kw">samples</span> <span class="kw">=</span> <span class="no">select_samples</span>)</pre></body></html></div>
<pre><code>## Warning: `select_()` is deprecated as of dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.</code></pre>
<pre><code>## Warning: `filter_()` is deprecated as of dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.</code></pre>
<div class="sourceCode" id="cb5"><html><body><pre class="r"><span class="no">embryo_data</span></pre></body></html></div>
<pre><code>## class: SingleCellExperiment 
## dim: 29452 14679 
## metadata(0):
## assays(1): counts
## rownames(29452): ENSMUSG00000051951 ENSMUSG00000089699 ...
##   ENSMUSG00000096730 ENSMUSG00000095742
## rowData names(2): ENSEMBL SYMBOL
## colnames(14679): cell_361 cell_362 ... cell_107079 cell_107080
## colData names(16): cell barcode ... celltype colour
## reducedDimNames(2): pca.corrected umap
## spikeNames(0):
## altExpNames(0):</code></pre>
</div>
<div id="visualize-the-data" class="section level1">
<h1 class="hasAnchor">
<a href="#visualize-the-data" class="anchor"></a>Visualize the data</h1>
<p>We recompute PCA and UMAP embedding for this subset of cells to visualize the data.</p>
<div class="sourceCode" id="cb7"><html><body><pre class="r"><span class="no">embryo_data</span> <span class="kw">&lt;-</span> <span class="no">embryo_data</span>[,<span class="fu"><a href="https://rdrr.io/r/base/apply.html">apply</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/SingleCellExperiment/man/reducedDims.html">reducedDim</a></span>(<span class="no">embryo_data</span>, <span class="st">"pca.corrected"</span>), <span class="fl">1</span>, <span class="kw">function</span>(<span class="no">x</span>) !<span class="fu"><a href="https://rdrr.io/r/base/all.html">all</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/NA.html">is.na</a></span>(<span class="no">x</span>)))]
<span class="no">embryo_data</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/pkg/scater/man/runUMAP.html">runUMAP</a></span>(<span class="no">embryo_data</span>, <span class="kw">dimred</span> <span class="kw">=</span> <span class="st">"pca.corrected"</span>, <span class="kw">name</span> <span class="kw">=</span> <span class="st">'umap'</span>)

<span class="fu"><a href="https://rdrr.io/pkg/scater/man/plotReducedDim.html">plotReducedDim</a></span>(<span class="no">embryo_data</span>, <span class="kw">colour_by</span><span class="kw">=</span><span class="st">"stage"</span>, <span class="kw">dimred</span> <span class="kw">=</span> <span class="st">"umap"</span>)</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-3-1.png" width="700"></p>
<p>We will test for significant differences in abundance of cells between these stages of development, and the associated gene signatures.</p>
</div>
<div id="differential-abundance-testing" class="section level1">
<h1 class="hasAnchor">
<a href="#differential-abundance-testing" class="anchor"></a>Differential abundance testing</h1>
<div id="create-a-milo-object" class="section level2">
<h2 class="hasAnchor">
<a href="#create-a-milo-object" class="anchor"></a>Create a Milo object</h2>
<p>For differential abundance analysis on graph neighbourhoods we first construct a <code>Milo</code> object. This extends the <a href="https://bioconductor.org/packages/release/bioc/html/SingleCellExperiment.html"><code>SingleCellExperiment</code></a> class to store information about neighbourhoods on the KNN graph.</p>
<div class="sourceCode" id="cb8"><html><body><pre class="r"><span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/Milo.html">Milo</a></span>(<span class="no">embryo_data</span>)
<span class="no">embryo_milo</span></pre></body></html></div>
<pre><code>## class: Milo 
## dim: 29452 13002 
## metadata(0):
## assays(1): counts
## rownames(29452): ENSMUSG00000051951 ENSMUSG00000089699 ...
##   ENSMUSG00000096730 ENSMUSG00000095742
## rowData names(2): ENSEMBL SYMBOL
## colnames(13002): cell_361 cell_362 ... cell_107079 cell_107080
## colData names(16): cell barcode ... celltype colour
## reducedDimNames(2): pca.corrected umap
## spikeNames(0):
## altExpNames(0):
## nhoods dimensions(1): 0
## nhoodCounts dimensions(2): 1 1
## nhoodDistances dimension(1): 0
## graph names(0):
## nhoodIndex names(1): 0
## nhoodExpression dimension(2): 1 1
## nhoodReducedDim names(0):
## nhoodGraph names(0):
## nhoodAdjacency dimension(0):</code></pre>
</div>
<div id="construct-knn-graph" class="section level2">
<h2 class="hasAnchor">
<a href="#construct-knn-graph" class="anchor"></a>Construct KNN graph</h2>
<p>We need to add the KNN graph to the Milo object. This is stored in the <code>graph</code> slot, in <a href="https://igraph.org/r/"><code>igraph</code></a> format. The <code>miloR</code> package includes functionality to build and store the graph from the PCA dimensions stored in the <code>reducedDim</code> slot. In this case, we specify that we want to build the graph from the MNN corrected PCA dimensions.</p>
<div class="sourceCode" id="cb10"><html><body><pre class="r"><span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/buildGraph.html">buildGraph</a></span>(<span class="no">embryo_milo</span>, <span class="kw">k</span> <span class="kw">=</span> <span class="fl">30</span>, <span class="kw">d</span> <span class="kw">=</span> <span class="fl">30</span>, <span class="kw">reduced.dim</span> <span class="kw">=</span> <span class="st">"pca.corrected"</span>)</pre></body></html></div>
<p>Alternatively, one can add a precomputed KNN graph (for example constructed with Seurat or scanpy) to the <code>graph</code> slot using the adjacency matrix, through the helper function <code>buildFromAdjacency</code>.</p>
<!-- Alternatively, if you already have a KNN graph (for example constructed with Seurat/scanpy) you can add it from the adjacency matrix. -->
<!-- ```{r} -->
<!-- # ## Build up a mock SNN graph made with Seurat -->
<!-- # pca_df <- reducedDim(traj_milo, "PCA") -->
<!-- # rownames(pca_df) <- traj_milo$cell_id -->
<!-- # snn_graph <- FindNeighbors(pca_df)[["snn"]] -->
<!-- #  -->
<!-- # graph(traj_milo) <-  graph(buildFromAdjacency(snn_graph, k=10)) -->
<!-- ``` -->
</div>
<div id="defining-representative-neighbourhoods-on-the-knn-graph" class="section level2">
<h2 class="hasAnchor">
<a href="#defining-representative-neighbourhoods-on-the-knn-graph" class="anchor"></a>Defining representative neighbourhoods on the KNN graph</h2>
<p>We define the neighbourhood of a cell, the index, as the group of cells connected by an edge in the KNN graph to the index cell. For efficiency, we don’t test for DA in the neighbourhood of every cell, but we sample as indices a subset of representative cells, using a KNN sampling algorithm used by <a href="https://www.nature.com/articles/nmeth.3545">Gut et al. 2015</a>.</p>
<p>For sampling you need to define a few parameters:</p>
<ul>
<li>
<code>prop</code>: the proportion of cells to randomly sample to start with (usually 0.1 - 0.2 is sufficient)</li>
<li>
<code>k</code>: the k to use for KNN refinement (we recommend using the same k used for KNN graph building)</li>
<li>
<code>d</code>: the number of reduced dimensions to use for KNN refinement (we recommend using the same d used for KNN graph building)</li>
<li>
<code>refined</code>: indicates whether you want to use the sampling refinement algorith, or just pick cells at random. The default and recommended way to go is to use refinement. The only situation in which you might consider using <code>random</code> instead, is if you have batch corrected your data with a graph based correction algorithm, such as <a href="https://github.com/Teichlab/bbknn">BBKNN</a>, but the results of DA testing will be suboptimal.</li>
</ul>
<div class="sourceCode" id="cb11"><html><body><pre class="r"><span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/makeNhoods.html">makeNhoods</a></span>(<span class="no">embryo_milo</span>, <span class="kw">prop</span> <span class="kw">=</span> <span class="fl">0.1</span>, <span class="kw">k</span> <span class="kw">=</span> <span class="fl">30</span>, <span class="kw">d</span><span class="kw">=</span><span class="fl">30</span>, <span class="kw">refined</span> <span class="kw">=</span> <span class="fl">TRUE</span>, <span class="kw">reduced_dims</span> <span class="kw">=</span> <span class="st">"pca.corrected"</span>)</pre></body></html></div>
<p>Once we have defined neighbourhoods, it’s good to take a look at how big the neighbourhoods are (i.e. how many cells form each neighbourhood). This affects the power of DA testing. We can check this out using the <code>plotNhoodSizeHist</code> function. Empirically, we found it’s best to have a distribution peaking above 20. Otherwise you might consider rerunning <code>makeNhoods</code> increasing <code>k</code> and/or <code>prop</code>.</p>
<div class="sourceCode" id="cb12"><html><body><pre class="r"><span class="fu"><a href="../reference/plotNhoodSizeHist.html">plotNhoodSizeHist</a></span>(<span class="no">embryo_milo</span>)</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-7-1.png" width="700"></p>
</div>
<div id="counting-cells-in-neighbourhoods" class="section level2">
<h2 class="hasAnchor">
<a href="#counting-cells-in-neighbourhoods" class="anchor"></a>Counting cells in neighbourhoods</h2>
<p><em>Milo</em> leverages the variation in cell numbers between replicates for the same experimental condition to test for differential abundance. Therefore we have to count how many cells from each sample are in each neighbourhood. We need to use the cell metadata and specify which column contains the sample information.</p>
<div class="sourceCode" id="cb13"><html><body><pre class="r"><span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/countCells.html">countCells</a></span>(<span class="no">embryo_milo</span>, <span class="kw">meta.data</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html">data.frame</a></span>(<span class="fu">colData</span>(<span class="no">embryo_milo</span>)), <span class="kw">sample</span><span class="kw">=</span><span class="st">"sample"</span>)</pre></body></html></div>
<p>This adds to the <code>Milo</code> object a <span class="math inline">\(n \times m\)</span> matrix, where <span class="math inline">\(n\)</span> is the number of neighbourhoods and <span class="math inline">\(m\)</span> is the number of experimental samples. Values indicate the number of cells from each sample counted in a neighbourhood. This count matrix will be used for DA testing.</p>
<div class="sourceCode" id="cb14"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span>(<span class="fu"><a href="../reference/methods.html">nhoodCounts</a></span>(<span class="no">embryo_milo</span>))</pre></body></html></div>
<pre><code>## 6 x 8 sparse Matrix of class "dgCMatrix"
##   2  3  6 15 10 14 20 30
## 1 .  2 22  7 39  1 12 55
## 2 .  1 12 24 73 16  3 53
## 3 7  6 20  .  .  . 13  3
## 4 1  6 34 14  6  6 30 30
## 5 8 11 38  1  .  2 19  .
## 6 .  .  3  8  7  5  2 29</code></pre>
</div>
<div id="defining-experimental-design" class="section level2">
<h2 class="hasAnchor">
<a href="#defining-experimental-design" class="anchor"></a>Defining experimental design</h2>
<p>Now we are all set to test for differential abundance in neighbourhoods. We implement this hypothesis testing in a generalized linear model (GLM) framework, specifically using the Negative Binomial GLM implementation in <a href="https://bioconductor.org/packages/release/bioc/html/edgeR.html"><code>edgeR</code></a>.</p>
<p>We first need to think about our experimental design. The design matrix should match each sample to the experimental condition of interest for DA testing. In this case, we want to detect DA between embryonic stages, stored in the <code>stage</code> column of the dataset <code>colData</code>. We also include the <code>sequencing.batch</code> column in the design matrix. This represents a known technical covariate that we want to account for in DA testing.</p>
<div class="sourceCode" id="cb16"><html><body><pre class="r"><span class="no">embryo_design</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html">data.frame</a></span>(<span class="fu">colData</span>(<span class="no">embryo_milo</span>))[,<span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"sample"</span>, <span class="st">"stage"</span>, <span class="st">"sequencing.batch"</span>)]
<span class="co">## Convert batch info from integer to factor</span>
<span class="no">embryo_design</span>$<span class="no">sequencing.batch</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/factor.html">as.factor</a></span>(<span class="no">embryo_design</span>$<span class="no">sequencing.batch</span>)
<span class="no">embryo_design</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://dplyr.tidyverse.org/reference/distinct.html">distinct</a></span>(<span class="no">embryo_design</span>)
<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">rownames</a></span>(<span class="no">embryo_design</span>) <span class="kw">&lt;-</span> <span class="no">embryo_design</span>$<span class="no">sample</span>

<span class="no">embryo_design</span></pre></body></html></div>
<pre><code>##    sample stage sequencing.batch
## 2       2  E7.5                1
## 3       3  E7.5                1
## 6       6  E7.5                1
## 15     15  E7.0                2
## 10     10  E7.0                1
## 14     14  E7.0                2
## 20     20  E7.5                2
## 30     30  E7.0                3</code></pre>
</div>
<div id="computing-neighbourhood-connectivity" class="section level2">
<h2 class="hasAnchor">
<a href="#computing-neighbourhood-connectivity" class="anchor"></a>Computing neighbourhood connectivity</h2>
<p>Milo uses an adaptation of the Spatial FDR correction introduced by <a href="https://bioconductor.org/packages/release/bioc/html/cydar.html">cydar</a>, which accounts for the overlap between neighbourhoods. Specifically, each hypothesis test P-value is weighted by the reciprocal of the kth nearest neighbour distance. To use this statistic we first need to store the distances between nearest neighbors in the Milo object. This is done by the <code>calcNhoodDistance</code> function (N.B. this step is the most time consuming of the analysis workflow and might take a couple of minutes for large datasets).</p>
<div class="sourceCode" id="cb18"><html><body><pre class="r"><span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/calcNhoodDistance.html">calcNhoodDistance</a></span>(<span class="no">embryo_milo</span>, <span class="kw">d</span><span class="kw">=</span><span class="fl">30</span>, <span class="kw">reduced.dim</span> <span class="kw">=</span> <span class="st">"pca.corrected"</span>)</pre></body></html></div>
</div>
<div id="testing" class="section level2">
<h2 class="hasAnchor">
<a href="#testing" class="anchor"></a>Testing</h2>
<p>Now we can do the DA test, explicitly defining our experimental design. In this case, we want to dest for differences between experimental stages, while accounting for the variability between technical batches (You can find more info on how to use formulas to define a testing design in R <a href="https://r4ds.had.co.nz/model-basics.html#formulas-and-model-families">here</a>)</p>
<div class="sourceCode" id="cb19"><html><body><pre class="r"><span class="no">da_results</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/testNhoods.html">testNhoods</a></span>(<span class="no">embryo_milo</span>, <span class="kw">design</span> <span class="kw">=</span> ~ <span class="no">sequencing.batch</span> + <span class="no">stage</span>, <span class="kw">design.df</span> <span class="kw">=</span> <span class="no">embryo_design</span>)</pre></body></html></div>
<p>This calculates a Fold-change and corrected P-value for each neighbourhood, which indicates wheather there is significant differential abundance between developmental stages.</p>
<div class="sourceCode" id="cb20"><html><body><pre class="r"><span class="no">da_results</span> <span class="kw">%&gt;%</span>
  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/arrange.html">arrange</a></span>(<span class="no">SpatialFDR</span>) <span class="kw">%&gt;%</span>
  <span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span>()</pre></body></html></div>
<pre><code>##         logFC   logCPM        F       PValue         FDR Nhood  SpatialFDR
## 45  -7.481702 19.97088 13.94271 0.0001914740 0.008423267    45 0.007301757
## 78  -8.286835 20.48252 14.97000 0.0001112280 0.008423267    78 0.007301757
## 91  -7.772064 20.03656 14.33882 0.0001552522 0.008423267    91 0.007301757
## 95   7.521388 19.68571 13.67331 0.0002208687 0.008423267    95 0.007301757
## 198 -7.579592 19.99991 14.05889 0.0001800453 0.008423267   198 0.007301757
## 274 -7.623243 19.86048 13.95762 0.0001899676 0.008423267   274 0.007301757</code></pre>
</div>
</div>
<div id="inspecting-da-testing-results" class="section level1">
<h1 class="hasAnchor">
<a href="#inspecting-da-testing-results" class="anchor"></a>Inspecting DA testing results</h1>
<p>We can start inspecting the results of our DA analysis from a couple of standard diagnostic plots. We first inspect the distribution of uncorrected P values, to verify that the test was balanced.</p>
<div class="sourceCode" id="cb22"><html><body><pre class="r"><span class="fu">ggplot</span>(<span class="no">da_results</span>, <span class="fu">aes</span>(<span class="no">PValue</span>)) + <span class="fu">geom_histogram</span>(<span class="kw">bins</span><span class="kw">=</span><span class="fl">50</span>)</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-14-1.png" width="700"></p>
<p>Then we visualize the test results with a volcano plot (remember that each point here represents a neighbourhood, <em>not</em> a cell).</p>
<div class="sourceCode" id="cb23"><html><body><pre class="r"><span class="fu">ggplot</span>(<span class="no">da_results</span>, <span class="fu">aes</span>(<span class="no">logFC</span>, -<span class="fu"><a href="https://rdrr.io/r/base/Log.html">log10</a></span>(<span class="no">SpatialFDR</span>))) +
  <span class="fu">geom_point</span>() +
  <span class="fu">geom_hline</span>(<span class="kw">yintercept</span> <span class="kw">=</span> <span class="fl">1</span>) <span class="co">## Mark significance threshold (10% FDR)</span></pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-15-1.png" width="700"></p>
<p>Looks like we have detected several neighbourhoods were there is a significant difference in cell abundances between developmental stages.</p>
<p>To visualize DA results relating them to the embedding of single cells, we can build an abstracted graph of neighbourhoods that we can superimpose on the single-cell embedding. Here each node represents a neighbourhood, while edges indicate how many cells two neighbourhoods have in common. Here the layout of nodes is determined by the position of the index cell in the UMAP embedding of all single-cells. The neighbourhoods displaying singificant DA are colored by their log-Fold Change.</p>
<div class="sourceCode" id="cb24"><html><body><pre class="r"><span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/buildNhoodGraph.html">buildNhoodGraph</a></span>(<span class="no">embryo_milo</span>)

<span class="co">## Plot single-cell UMAP</span>
<span class="no">umap_pl</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/pkg/scater/man/plotReducedDim.html">plotReducedDim</a></span>(<span class="no">embryo_milo</span>, <span class="kw">dimred</span> <span class="kw">=</span> <span class="st">"umap"</span>, <span class="kw">colour_by</span><span class="kw">=</span><span class="st">"celltype"</span>, <span class="kw">text_by</span> <span class="kw">=</span> <span class="st">"celltype"</span>, <span class="kw">text_size</span> <span class="kw">=</span> <span class="fl">3</span>) +
  <span class="fu">guides</span>(<span class="kw">fill</span><span class="kw">=</span><span class="st">"none"</span>)

<span class="co">## Plot neighbourhood graph</span>
<span class="no">nh_graph_pl</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/plotNhoodGraphDA.html">plotNhoodGraphDA</a></span>(<span class="no">embryo_milo</span>, <span class="no">da_results</span>, <span class="kw">layout</span><span class="kw">=</span><span class="st">"umap"</span>,<span class="kw">alpha</span><span class="kw">=</span><span class="fl">0.05</span>)

<span class="no">umap_pl</span> + <span class="no">nh_graph_pl</span> +
  <span class="fu"><a href="https://patchwork.data-imaginist.com/reference/plot_layout.html">plot_layout</a></span>(<span class="kw">guides</span><span class="kw">=</span><span class="st">"collect"</span>)</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-16-1.png" width="1440"></p>
<p>We might also be interested in visualizing wheather DA is particularly evident in certain cell types. To do this, we assign a cell type label to each neighbourhood by finding the most abundant cell type within cells in each neighbourhood. We can label neighbourhoods in the results <code>data.frame</code> using the function <code>annotateNhoods</code>. This also saves the fraction of cells harbouring the label.</p>
<div class="sourceCode" id="cb25"><html><body><pre class="r"><span class="no">da_results</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/annotateNhoods.html">annotateNhoods</a></span>(<span class="no">embryo_milo</span>, <span class="no">da_results</span>, <span class="kw">coldata_col</span> <span class="kw">=</span> <span class="st">"celltype"</span>)
<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span>(<span class="no">da_results</span>)</pre></body></html></div>
<pre><code>##        logFC   logCPM          F      PValue        FDR Nhood SpatialFDR
## 1 -0.6106328 20.42087  0.1788746 0.672368256 0.73556076     1 0.73343784
## 2 -3.4251846 20.79253  4.5355029 0.033268444 0.07663154     2 0.07569422
## 3  6.6199781 19.28856 10.5760794 0.001156497 0.01284465     3 0.01232507
## 4  1.3986635 20.35858  0.7999356 0.371174279 0.45256305     4 0.45113575
## 5  4.5521075 19.83243  6.4076675 0.011406385 0.03999106     5 0.03912235
## 6 -2.1349494 19.30941  2.2261867 0.135778701 0.20863903     6 0.20738645
##                    celltype celltype_fraction
## 1                  Epiblast         0.9927536
## 2              ExE ectoderm         1.0000000
## 3      Rostral neurectoderm         0.8163265
## 4            Mixed mesoderm         0.8031496
## 5                Mesenchyme         0.9873418
## 6 Anterior Primitive Streak         0.9814815</code></pre>
<p>While neighbourhoods tend to be homogeneous, we can define a threshold for <code>celltype_fraction</code> to exclude neighbourhoods that are a mix of cell types.</p>
<div class="sourceCode" id="cb27"><html><body><pre class="r"><span class="fu">ggplot</span>(<span class="no">da_results</span>, <span class="fu">aes</span>(<span class="no">celltype_fraction</span>)) + <span class="fu">geom_histogram</span>(<span class="kw">bins</span><span class="kw">=</span><span class="fl">50</span>)</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-18-1.png" width="700"></p>
<div class="sourceCode" id="cb28"><html><body><pre class="r"><span class="no">da_results</span>$<span class="no">celltype</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/ifelse.html">ifelse</a></span>(<span class="no">da_results</span>$<span class="no">celltype_fraction</span> <span class="kw">&lt;</span> <span class="fl">0.7</span>, <span class="st">"Mixed"</span>, <span class="no">da_results</span>$<span class="no">celltype</span>)</pre></body></html></div>
<p>Now we can visualize the distribution of DA Fold Changes in different cell types</p>
<div class="sourceCode" id="cb29"><html><body><pre class="r"><span class="fu"><a href="../reference/plotDAbeeswarm.html">plotDAbeeswarm</a></span>(<span class="no">da_results</span>, <span class="kw">group.by</span> <span class="kw">=</span> <span class="st">"celltype"</span>)</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-20-1.png" width="672"></p>
<p>This is already quite informative: we can see that certain early development cell types, such as epiblast and primitive streak, are enriched in the earliest time stage, while others are enriched later in development, such as ectoderm cells. Interestingly, we also see plenty of DA neighbourhood with a mixed label. This could indicate that transitional states show changes in abundance in time.</p>
</div>
<div id="identifying-signatures-of-da-subpopulations" class="section level1">
<h1 class="hasAnchor">
<a href="#identifying-signatures-of-da-subpopulations" class="anchor"></a>Identifying signatures of DA subpopulations</h1>
<p>Once we have identified neighbourhoods with significant DA, we might be interested in learning more about the gene expression signatures that define these subpopulations. Here the analyst might get creative, depending on the specific characteristics of their dataset and the biological question of interest. In the <code>miloR</code> package, we provide functionality to explore neighbourhood markers through the function <code>findNhoodMarkers</code>. This groups significantly DA neighbourhoods that show concordant Fold-Change and performs a test for differential expression between cells in those neighbourhoods.</p>
<p>In practice, it might be convenient to subset a selected number of neighbourhoods of interest for gene-level downstream analysis. For the sake of demonstration, here we focus on identifying signatures of DA subpopulations in the endoderm development lineage.</p>
<div class="sourceCode" id="cb30"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/pkg/SingleCellExperiment/man/assays.html">logcounts</a></span>(<span class="no">embryo_milo</span>) <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Log.html">log1p</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/SingleCellExperiment/man/assays.html">counts</a></span>(<span class="no">embryo_milo</span>))
<span class="no">dge_smp</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/findNhoodMarkers.html">findNhoodMarkers</a></span>(<span class="no">embryo_milo</span>, <span class="no">da_results</span>,
                                     <span class="kw">assay</span> <span class="kw">=</span> <span class="st">"counts"</span>, <span class="kw">gene.offset</span> <span class="kw">=</span> <span class="fl">FALSE</span>, <span class="kw">da.fdr</span> <span class="kw">=</span> <span class="fl">0.1</span>,
                                     <span class="kw">aggregate.samples</span> <span class="kw">=</span> <span class="fl">TRUE</span>, <span class="kw">sample_col</span> <span class="kw">=</span> <span class="st">"sample"</span>,
                                     <span class="kw">subset.nhoods</span> <span class="kw">=</span> <span class="no">da_results</span>$<span class="no">celltype</span> <span class="kw">%in%</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"Anterior Primitive Streak"</span>, <span class="st">"Def. endoderm"</span>, <span class="st">"Gut"</span>, <span class="st">"Visceral endoderm"</span>)
                                     )

<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span>(<span class="no">dge_smp</span>)</pre></body></html></div>
<pre><code>##                       logFC_1 adj.P.Val_1       logFC_2 adj.P.Val_2    logFC_3
## ENSMUSG00000051951 0.00000000           1 -1.281371e-15           1  0.0000000
## ENSMUSG00000089699 0.00000000           1 -1.281371e-15           1  0.0000000
## ENSMUSG00000102343 0.00000000           1 -1.281371e-15           1  0.0000000
## ENSMUSG00000025900 0.00000000           1 -1.281371e-15           1  0.0000000
## ENSMUSG00000025902 0.01878518           1  4.214069e-01           1 -0.4762632
## ENSMUSG00000104328 0.00000000           1 -1.281371e-15           1  0.0000000
##                    adj.P.Val_3             GeneID
## ENSMUSG00000051951           1 ENSMUSG00000051951
## ENSMUSG00000089699           1 ENSMUSG00000089699
## ENSMUSG00000102343           1 ENSMUSG00000102343
## ENSMUSG00000025900           1 ENSMUSG00000025900
## ENSMUSG00000025902           1 ENSMUSG00000025902
## ENSMUSG00000104328           1 ENSMUSG00000104328</code></pre>
<p>This identifies n marker genes at FDR 10% that distinguish two main groups within the epiblast neighbourhoods, one significantly depleted in the early stage and one significantly enriched. We can visualize expression of the detected marker genes using the function <code>plotNhoodExpressionDA</code>. This shows the average expression in each neighbourhood, ranked by log-Fold Change in the DA test. Note that the gene x nhood expression matrix can be pre-computed and stored using the <code>calcNhoodExpression</code> function, to avoid repeating the computation every time you need to plot.</p>
<p>In this case we mainly identified negative markers of the epiblast neighbourhoods enriched with age.</p>
<div class="sourceCode" id="cb32"><html><body><pre class="r"><span class="no">markers</span> <span class="kw">&lt;-</span> <span class="no">dge_smp</span>[<span class="fu"><a href="https://rdrr.io/r/base/which.html">which</a></span>(<span class="no">dge_smp</span>$<span class="no">adj.P.Val_1</span> <span class="kw">&lt;</span> <span class="fl">0.1</span> ), <span class="st">"GeneID"</span>]
<span class="fu"><a href="https://rdrr.io/pkg/SingleCellExperiment/man/assays.html">logcounts</a></span>(<span class="no">embryo_milo</span>) <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Log.html">log1p</a></span>(<span class="fu"><a href="https://rdrr.io/pkg/SingleCellExperiment/man/assays.html">counts</a></span>(<span class="no">embryo_milo</span>))
<span class="no">embryo_milo</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/calcNhoodExpression.html">calcNhoodExpression</a></span>(<span class="no">embryo_milo</span>, <span class="kw">subset.row</span><span class="kw">=</span><span class="no">markers</span>)

<span class="fu"><a href="../reference/plotNhoodExpressionDA.html">plotNhoodExpressionDA</a></span>(<span class="no">embryo_milo</span>, <span class="no">da_results</span>, <span class="kw">features</span> <span class="kw">=</span> <span class="no">markers</span>,
                      <span class="kw">subset.nhoods</span> <span class="kw">=</span> <span class="no">da_results</span>$<span class="no">celltype</span> <span class="kw">%in%</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="st">"Anterior Primitive Streak"</span>, <span class="st">"Def. endoderm"</span>, <span class="st">"Gut"</span>, <span class="st">"Visceral endoderm"</span>),
                      <span class="kw">assay</span><span class="kw">=</span><span class="st">"logcounts"</span>,
                      <span class="kw">scale_to_1</span> <span class="kw">=</span> <span class="fl">TRUE</span>, <span class="kw">cluster_features</span> <span class="kw">=</span> <span class="fl">TRUE</span>
                      )</pre></body></html></div>
<p><img src="milo_gastrulation_files/figure-html/unnamed-chunk-22-1.png" width="864"></p>
<details><p><summary><strong>Session Info</strong></summary></p>
<div class="sourceCode" id="cb33"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/utils/sessionInfo.html">sessionInfo</a></span>()</pre></body></html></div>
<pre><code>## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas/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] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] MouseGastrulationData_1.0.0 patchwork_1.0.1            
##  [3] dplyr_1.0.2                 scater_1.14.6              
##  [5] ggplot2_3.3.2               SingleCellExperiment_1.8.0 
##  [7] SummarizedExperiment_1.16.1 DelayedArray_0.12.3        
##  [9] BiocParallel_1.20.1         matrixStats_0.56.0         
## [11] Biobase_2.46.0              GenomicRanges_1.38.0       
## [13] GenomeInfoDb_1.22.1         IRanges_2.20.2             
## [15] S4Vectors_0.24.4            BiocGenerics_0.32.0        
## [17] miloR_0.1.0                 edgeR_3.28.1               
## [19] limma_3.42.2                BiocStyle_2.14.4           
## 
## loaded via a namespace (and not attached):
##   [1] ggbeeswarm_0.6.0              colorspace_1.4-1             
##   [3] ellipsis_0.3.1                rprojroot_1.3-2              
##   [5] XVector_0.26.0                BiocNeighbors_1.4.2          
##   [7] fs_1.5.0                      rstudioapi_0.11              
##   [9] farver_2.0.3                  graphlayouts_0.7.0           
##  [11] ggrepel_0.8.2                 bit64_4.0.2                  
##  [13] AnnotationDbi_1.48.0          interactiveDisplayBase_1.24.0
##  [15] splines_3.6.3                 codetools_0.2-16             
##  [17] knitr_1.30                    polyclip_1.10-0              
##  [19] dbplyr_1.4.4                  uwot_0.1.8                   
##  [21] ggforce_0.3.2                 shiny_1.5.0                  
##  [23] BiocManager_1.30.10           compiler_3.6.3               
##  [25] httr_1.4.2                    backports_1.1.9              
##  [27] fastmap_1.0.1                 assertthat_0.2.1             
##  [29] Matrix_1.2-18                 later_1.1.0.1                
##  [31] tweenr_1.0.1                  BiocSingular_1.2.2           
##  [33] htmltools_0.5.0               tools_3.6.3                  
##  [35] rsvd_1.0.3                    igraph_1.2.5                 
##  [37] gtable_0.3.0                  glue_1.4.1                   
##  [39] GenomeInfoDbData_1.2.2        rappdirs_0.3.1               
##  [41] Rcpp_1.0.5                    pkgdown_1.5.1                
##  [43] vctrs_0.3.2                   ExperimentHub_1.12.0         
##  [45] DelayedMatrixStats_1.8.0      ggraph_2.0.3                 
##  [47] xfun_0.19                     stringr_1.4.0                
##  [49] mime_0.9                      lifecycle_0.2.0              
##  [51] irlba_2.3.3                   gtools_3.8.2                 
##  [53] statmod_1.4.34                AnnotationHub_2.18.0         
##  [55] zlibbioc_1.32.0               MASS_7.3-51.5                
##  [57] scales_1.1.1                  tidygraph_1.2.0              
##  [59] promises_1.1.1                yaml_2.2.1                   
##  [61] curl_4.3                      memoise_1.1.0                
##  [63] gridExtra_2.3                 stringi_1.4.6                
##  [65] RSQLite_2.2.0                 BiocVersion_3.10.1           
##  [67] desc_1.2.0                    rlang_0.4.7                  
##  [69] pkgconfig_2.0.3               bitops_1.0-6                 
##  [71] evaluate_0.14                 lattice_0.20-40              
##  [73] purrr_0.3.4                   labeling_0.3                 
##  [75] cowplot_1.0.0                 bit_4.0.4                    
##  [77] tidyselect_1.1.0              RcppAnnoy_0.0.16             
##  [79] magrittr_1.5                  bookdown_0.20                
##  [81] R6_2.4.1                      generics_0.0.2               
##  [83] DBI_1.1.0                     pillar_1.4.6                 
##  [85] withr_2.2.0                   RCurl_1.98-1.2               
##  [87] tibble_3.0.3                  crayon_1.3.4                 
##  [89] BiocFileCache_1.10.2          rmarkdown_2.5                
##  [91] viridis_0.5.1                 locfit_1.5-9.4               
##  [93] grid_3.6.3                    blob_1.2.1                   
##  [95] digest_0.6.25                 xtable_1.8-4                 
##  [97] tidyr_1.1.1                   httpuv_1.5.4                 
##  [99] munsell_0.5.0                 beeswarm_0.2.3               
## [101] viridisLite_0.3.0             vipor_0.4.5</code></pre>
</details>
</div>
  </div>

  <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">

        <nav id="toc" data-toggle="toc"><h2 data-toc-skip>Contents</h2>
    </nav>
</div>

</div>



      <footer><div class="copyright">
  <p>Developed by Mike Morgan, Emma Dann.</p>
</div>

<div class="pkgdown">
  <p>Site built with <a href="https://pkgdown.r-lib.org/">pkgdown</a> 1.5.1.</p>
</div>

      </footer>
</div>

  


  </body>
</html>
