<!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, shrink-to-fit=no">
<title>Gene Expression Analysis with `myTAI` • myTAI</title>
<!-- favicons --><link rel="icon" type="image/png" sizes="16x16" href="../favicon-16x16.png">
<link rel="icon" type="image/png" sizes="32x32" href="../favicon-32x32.png">
<link rel="apple-touch-icon" type="image/png" sizes="180x180" href="../apple-touch-icon.png">
<link rel="apple-touch-icon" type="image/png" sizes="120x120" href="../apple-touch-icon-120x120.png">
<link rel="apple-touch-icon" type="image/png" sizes="76x76" href="../apple-touch-icon-76x76.png">
<link rel="apple-touch-icon" type="image/png" sizes="60x60" href="../apple-touch-icon-60x60.png">
<script src="../deps/jquery-3.6.0/jquery-3.6.0.min.js"></script><meta name="viewport" content="width=device-width, initial-scale=1, shrink-to-fit=no">
<link href="../deps/bootstrap-5.3.1/bootstrap.min.css" rel="stylesheet">
<script src="../deps/bootstrap-5.3.1/bootstrap.bundle.min.js"></script><link href="../deps/font-awesome-6.4.2/css/all.min.css" rel="stylesheet">
<link href="../deps/font-awesome-6.4.2/css/v4-shims.min.css" rel="stylesheet">
<script src="../deps/headroom-0.11.0/headroom.min.js"></script><script src="../deps/headroom-0.11.0/jQuery.headroom.min.js"></script><script src="../deps/bootstrap-toc-1.0.1/bootstrap-toc.min.js"></script><script src="../deps/clipboard.js-2.0.11/clipboard.min.js"></script><script src="../deps/search-1.0.0/autocomplete.jquery.min.js"></script><script src="../deps/search-1.0.0/fuse.min.js"></script><script src="../deps/search-1.0.0/mark.min.js"></script><!-- pkgdown --><script src="../pkgdown.js"></script><meta property="og:title" content="Gene Expression Analysis with `myTAI`">
</head>
<body>
    <a href="#main" class="visually-hidden-focusable">Skip to contents</a>


    <nav class="navbar navbar-expand-lg fixed-top bg-light" data-bs-theme="light" aria-label="Site navigation"><div class="container">

    <a class="navbar-brand me-2" href="../index.html">myTAI</a>

    <small class="nav-text text-muted me-auto" data-bs-toggle="tooltip" data-bs-placement="bottom" title="">2.0.0</small>


    <button class="navbar-toggler" type="button" data-bs-toggle="collapse" data-bs-target="#navbar" aria-controls="navbar" aria-expanded="false" aria-label="Toggle navigation">
      <span class="navbar-toggler-icon"></span>
    </button>

    <div id="navbar" class="collapse navbar-collapse ms-3">
      <ul class="navbar-nav me-auto">
<li class="nav-item"><a class="nav-link" href="../reference/index.html">Reference</a></li>
<li class="active nav-item dropdown">
  <button class="nav-link dropdown-toggle" type="button" id="dropdown-articles" data-bs-toggle="dropdown" aria-expanded="false" aria-haspopup="true">Articles</button>
  <ul class="dropdown-menu" aria-labelledby="dropdown-articles">
<li><a class="dropdown-item" href="../articles/Advanced.html">Advanced Phylotranscriptomics Analyses</a></li>
    <li><a class="dropdown-item" href="../articles/Enrichment.html">Enrichment Analyses</a></li>
    <li><a class="dropdown-item" href="../articles/Expression.html">Gene Expression Analysis with `myTAI`</a></li>
    <li><a class="dropdown-item" href="../articles/Intermediate.html">Intermediate Phylotranscriptomics Analyses</a></li>
    <li><a class="dropdown-item" href="../articles/Introduction.html">Introduction to the myTAI package</a></li>
    <li><a class="dropdown-item" href="../articles/Phylostratigraphy.html">Discussing Genomic Phylostratigraphy</a></li>
  </ul>
</li>
<li class="nav-item"><a class="nav-link" href="../news/index.html">Changelog</a></li>
      </ul>
<ul class="navbar-nav">
<li class="nav-item"><form class="form-inline" role="search">
 <input class="form-control" type="search" name="search-input" id="search-input" autocomplete="off" aria-label="Search site" placeholder="Search for" data-search-index="../search.json">
</form></li>
<li class="nav-item"><a class="external-link nav-link" href="https://github.com/drostlab/myTAI/" aria-label="GitHub"><span class="fa fab fa-github fa-lg"></span></a></li>
      </ul>
</div>


  </div>
</nav><div class="container template-article">




<div class="row">
  <main id="main" class="col-md-9"><div class="page-header">
      <img src="../logo.png" class="logo" alt=""><h1>Gene Expression Analysis with `myTAI`</h1>
            
            <h4 data-toc-skip class="date">2024-10-03</h4>
      
      <small class="dont-index">Source: <a href="https://github.com/drostlab/myTAI/blob/HEAD/vignettes/Expression.Rmd" class="external-link"><code>vignettes/Expression.Rmd</code></a></small>
      <div class="d-none name"><code>Expression.Rmd</code></div>
    </div>

    
    
<div class="section level2">
<h2 id="introduction">Introduction<a class="anchor" aria-label="anchor" href="#introduction"></a>
</h2>
<p>In the <a href="Introduction.html">Introduction</a> vignette we
introduced and discussed how phylotranscriptomics can be applied to
capture evolutionary signals in (developmental) transcriptomes.
Furthermore, in the <a href="Enrichment.html">Enrichment Analyses</a>
vignette we provide a use case to correlate specific groups or sets of
genes with their predicted evolutionary origin. Here, we aim to combine
previously introduced techniques with <em>classic</em> gene expression
analyses to detect possible functional causes for the observed
transcriptome conservation.</p>
<p>In other words, phylotranscriptomics allows us to detect stages or
periods of evolutionary conservation and is able to predict the
evolutionary origin of process or trait specific genes based on
enrichment analyses. By combining evolutionary enrichment analyses with
the functional annotation of process or trait specific genes (see <a href="https://github.com/ropensci/biomartr/blob/master/vignettes/Functional_Annotation.Rmd" class="external-link">Functional
Annotation</a> for details) the detection of evolutionary signals can be
correlated with functional processes. Then, performing gene expression
analyses on corresponding process or trait specific genes allows users
to detect potential causes of stage/period specific evolutionary
transcriptome conservation.</p>
<p>The following sections introduce main gene expression data analysis
techniques implemented in <code>myTAI</code>:</p>
<ul>
<li><p>Detection of Differentially Expressed Genes (DEGs)</p></li>
<li><p>Fold-Change</p></li>
<li><p>Welch t-test</p></li>
<li><p>Wilcoxon Rank Sum Test (Mann-Whitney U test)</p></li>
<li><p>Negative Binomial (Exact Tests)</p></li>
<li><p>Collapsing Replicate Samples</p></li>
<li><p>Filter for Expressed Genes</p></li>
<li><p>Compute the Statistical Significance of Each Replicate
Combination</p></li>
</ul>
</div>
<div class="section level2">
<h2 id="detection-of-differenentially-expressed-genes-degs">Detection of Differenentially Expressed Genes (DEGs)<a class="anchor" aria-label="anchor" href="#detection-of-differenentially-expressed-genes-degs"></a>
</h2>
<p>A variety of methods have been published to detect differentially
expressed genes. Some methods are based on non-statistical
quantification of expression differences (e.g. fold-change and
log-fold-change), but most methods are based on statistical tests to
quantify the significance of differences in gene expression between
samples. These statistical methods can furthermore be divided into two
methodological categories: parametric tests and non-parametric tests.
The <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> function available in <code>myTAI</code>
implements the most popular and useful methods to detect differentially
expressed genes. In the literature, different methods have been
introduced and discussed for microarray technologies versus RNA-Seq
technologies.</p>
<p>In this section we will introduce all methods implemented in
<code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> using small examples and will furthermore,
discuss published advantages and disadvantages of each method and each
mRNA quantification technology.</p>
<p><strong>Note that when using <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> it is assumed
that your input dataset has been normalized before passing it to
<code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>. For RNA-Seq data <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>
assumes that the libraries have been normalized to have the same size,
i.e., to have the same expected column sum under the null hypothesis (or
the lib.size argument in <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> is specified
accordingly).</strong></p>
</div>
<div class="section level2">
<h2 id="fold-changes">Fold-Changes<a class="anchor" aria-label="anchor" href="#fold-changes"></a>
</h2>
<p>A fold change in gene expression is simply the ratio of the gene
expression level of one sample against a second sample:
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mfrac><msub><mi>e</mi><mrow><mi>i</mi><mn>1</mn></mrow></msub><msub><mi>e</mi><mrow><mi>i</mi><mn>2</mn></mrow></msub></mfrac><annotation encoding="application/x-tex">\frac{e_{i1}}{e_{i2}}</annotation></semantics></math>,
where
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mi>e</mi><mrow><mi>i</mi><mn>1</mn></mrow></msub><annotation encoding="application/x-tex">e_{i1}</annotation></semantics></math>
is the expression level of gene
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math>
in sample one and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mi>e</mi><mrow><mi>i</mi><mn>2</mn></mrow></msub><annotation encoding="application/x-tex">e_{i2}</annotation></semantics></math>
is the expression level of gene
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math>
in sample two. In case replicate expression levels are present for each
sample the ratio of means of the corresponding replicates is computed:
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mfrac><msub><mover><mi>e</mi><mo accent="true">‾</mo></mover><mrow><mi>i</mi><mn>1</mn></mrow></msub><msub><mover><mi>e</mi><mo accent="true">‾</mo></mover><mrow><mi>i</mi><mn>2</mn></mrow></msub></mfrac><annotation encoding="application/x-tex">\frac{\bar{e}_{i1}}{\bar{e}_{i2}}</annotation></semantics></math>,
where
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mover><mi>e</mi><mo accent="true">‾</mo></mover><mrow><mi>i</mi><mn>1</mn></mrow></msub><annotation encoding="application/x-tex">\bar{e}_{i1}</annotation></semantics></math>
is the mean of replicate expression levels of gene
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math>
in sample one and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mover><mi>e</mi><mo accent="true">‾</mo></mover><mrow><mi>i</mi><mn>2</mn></mrow></msub><annotation encoding="application/x-tex">\bar{e}_{i2}</annotation></semantics></math>
is the mean of replicate expression levels of gene
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math>
in sample two.</p>
<ul>
<li><p><strong>Advantages:</strong> Given a small number of replicate
values the statistical evaluation of differentially expressed genes
might be biased (depending on the statistical test chosen) by underlying
sample distributions which are not fulfilled or because a small number
of replicate values is not sufficient enough to perform non-parametric
tests. Here, fold-changes provide a simple way to quantify gene
expression differences between samples by
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>n</mi><annotation encoding="application/x-tex">n</annotation></semantics></math>-fold
enrichment. In our opinion, although the process of choosing a threshold
for defining genes as being differentially expressed or not based on
fold-change values is purely subjective and relies on common sense, in
some cases this procedure will be more suitable than defining
differentially expressed genes based on p-values obtained from a test
statistic with violated test assumptions.</p></li>
<li><p><strong>Disadvantages:</strong> If used appropriately,
statistical tests not only systematically quantify the significance of
the observed gene-by-gene differences of expression, but furthermore,
accounts the variance of replicate expression levels when comparing the
mean difference of replicate expression levels between samples. Hence,
the gene specific variance between replicates is also quantified by the
p-value returned by the sufficient test statistic which is not
quantified by a simple fold-change measure.</p></li>
</ul>
<div class="section level3">
<h3 id="example-fold-change">Example: Fold-Change<a class="anchor" aria-label="anchor" href="#example-fold-change"></a>
</h3>
<p>For the following example we assume that
<code>PhyloExpressionSetExample[1:5,1:8]</code> stores 5 genes and 3
developmental stages with 2 replicate expression levels per stage.</p>
<div class="sourceCode" id="cb1"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the fold-change measure</span></span>
<span><span class="va">DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                  nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                  method        <span class="op">=</span> <span class="st">"foldchange"</span>,</span>
<span>                  stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/head.html" class="external-link">head</a></span><span class="op">(</span><span class="va">DEGs</span><span class="op">)</span></span></code></pre></div>
<pre><code>   Phylostratum      GeneID    S1-&gt;S2    S1-&gt;S3    S2-&gt;S1    S2-&gt;S3    S3-&gt;S1    S3-&gt;S2
1            1 at1g01040.2 1.6713881 2.0806706 0.5983051 1.2448758 0.4806143 0.8032930
2            1 at1g01050.1 1.0273222 1.2709185 0.9734045 1.2371177 0.7868325 0.8083305
3            1 at1g01070.1 1.3087379 1.4044799 0.7640949 1.0731560 0.7120073 0.9318310
4            1 at1g01080.2 0.7779572 0.7286769 1.2854177 0.9366542 1.3723503 1.0676299
5            1 at1g01090.1 0.3803866 0.2288961 2.6289042 0.6017460 4.3687939 1.6618307</code></pre>
<p>The resulting output shows all combinations of fold-changes between
samples (developmental stages). Here, <code>S1-&gt;S2</code> denotes
that the fold-change was computed for expression levels of stage
<code>S1</code> against stage <code>S2</code>.</p>
</div>
<div class="section level3">
<h3 id="example-log-fold-change">Example: Log-Fold-Change<a class="anchor" aria-label="anchor" href="#example-log-fold-change"></a>
</h3>
<p><strong>When selecting <code>method = "log-foldchange"</code> it is
assumed that the input <code>ExpressionSet</code> stores
<code>log2</code> expression levels. Here, we transform absolute
expression levels stored in <code>PhyloExpressionSetExample</code> to
<code>log2</code> expression levels using the <code><a href="../reference/tf.html">tf()</a></code> function
before log-fold-changes are computed.</strong></p>
<div class="sourceCode" id="cb3"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the logfold-change measure</span></span>
<span><span class="va">log.DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="fu"><a href="../reference/tf.html">tf</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,<span class="va">log2</span><span class="op">)</span>,</span>
<span>                      nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                      method        <span class="op">=</span> <span class="st">"log-foldchange"</span>,</span>
<span>                      stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/head.html" class="external-link">head</a></span><span class="op">(</span><span class="va">log.DEGs</span><span class="op">)</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID      S1-&gt;S2     S1-&gt;S3      S2-&gt;S1      S2-&gt;S3     S3-&gt;S1      S3-&gt;S2
1            1 at1g01040.2  0.74104679  1.0570486 -0.74104679  0.31600182 -1.0570486 -0.31600182
2            1 at1g01050.1  0.03888868  0.3458715 -0.03888868  0.30698280 -0.3458715 -0.30698280
3            1 at1g01070.1  0.38817621  0.4900360 -0.38817621  0.10185975 -0.4900360 -0.10185975
4            1 at1g01080.2 -0.36223724 -0.4566488  0.36223724 -0.09441158  0.4566488  0.09441158
5            1 at1g01090.1 -1.39446159 -2.1272350  1.39446159 -0.73277345  2.1272350  0.73277345</code></pre>
<p>The resulting output stores all combinations of log fold-changes
between samples (developmental stages).</p>
</div>
</div>
<div class="section level2">
<h2 id="welch-t-test">Welch t-test<a class="anchor" aria-label="anchor" href="#welch-t-test"></a>
</h2>
<p>The <code>Welch t-test</code> is a parametric test to statistically
quantify the difference of sample means in cases where the assumption of
homogeneity of variance (equal variances in the two populations) is
violated (Boslaugh, 2013). The <code>Welch t-test</code> is a sufficient
parameter test for small sample sizes and thus, has been used to detect
differentially expressed genes based on p-values returned by the test
statistic (Hahne et al., 2008).</p>
<p>In detail, the test statistic is computed as follows:</p>
<p><math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>t</mi><mo>=</mo><mfrac><mrow><msub><mover><mi>x</mi><mo accent="true">‾</mo></mover><mn>1</mn></msub><mo>−</mo><msub><mover><mi>x</mi><mo accent="true">‾</mo></mover><mn>2</mn></msub></mrow><msqrt><mrow><mfrac><msubsup><mi>s</mi><mn>1</mn><mn>2</mn></msubsup><msub><mi>n</mi><mn>1</mn></msub></mfrac><mo>+</mo><mfrac><msubsup><mi>s</mi><mn>2</mn><mn>2</mn></msubsup><msub><mi>n</mi><mn>2</mn></msub></mfrac></mrow></msqrt></mfrac></mrow><annotation encoding="application/x-tex">t = \frac{\bar{x}_1 - \bar{x}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}</annotation></semantics></math></p>
<p>where
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mover><mi>x</mi><mo accent="true">‾</mo></mover><mn>1</mn></msub><annotation encoding="application/x-tex">\bar{x}_1</annotation></semantics></math>
and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mover><mi>x</mi><mo accent="true">‾</mo></mover><mn>2</mn></msub><annotation encoding="application/x-tex">\bar{x}_2</annotation></semantics></math>
are sample means,
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>s</mi><mn>1</mn><mn>2</mn></msubsup><annotation encoding="application/x-tex">s_1^2</annotation></semantics></math>
and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msubsup><mi>s</mi><mn>2</mn><mn>2</mn></msubsup><annotation encoding="application/x-tex">s_2^2</annotation></semantics></math>
are the sample variances, and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mi>n</mi><mn>1</mn></msub><annotation encoding="application/x-tex">n_1</annotation></semantics></math>
and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mi>n</mi><mn>2</mn></msub><annotation encoding="application/x-tex">n_2</annotation></semantics></math>
are the sample sizes.</p>
<p>The degrees of freedom for Welch’s t-test are then computed as
follows:</p>
<p><math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>d</mi><mi>f</mi><mo>=</mo><mfrac><mrow><mo minsize="1.2" maxsize="1.2" stretchy="false" form="prefix">(</mo><mfrac><msubsup><mi>s</mi><mn>1</mn><mn>2</mn></msubsup><msub><mi>n</mi><mn>1</mn></msub></mfrac><mo>+</mo><mfrac><msubsup><mi>s</mi><mn>2</mn><mn>2</mn></msubsup><msub><mi>n</mi><mn>2</mn></msub></mfrac><msup><mo minsize="1.2" maxsize="1.2" stretchy="false" form="postfix">)</mo><mn>2</mn></msup></mrow><mrow><mfrac><msubsup><mi>s</mi><mn>1</mn><mn>4</mn></msubsup><mrow><msubsup><mi>n</mi><mn>1</mn><mn>2</mn></msubsup><mrow><mo stretchy="true" form="prefix">(</mo><msub><mi>n</mi><mn>1</mn></msub><mo>−</mo><mn>1</mn><mo stretchy="true" form="postfix">)</mo></mrow></mrow></mfrac><mo>+</mo><mfrac><msubsup><mi>s</mi><mn>2</mn><mn>4</mn></msubsup><mrow><msubsup><mi>n</mi><mn>2</mn><mn>2</mn></msubsup><mrow><mo stretchy="true" form="prefix">(</mo><msub><mi>n</mi><mn>2</mn></msub><mo>−</mo><mn>1</mn><mo stretchy="true" form="postfix">)</mo></mrow></mrow></mfrac></mrow></mfrac></mrow><annotation encoding="application/x-tex">df = \frac{\big(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}\big)^2}{\frac{s_1^4}{n_1^2 (n_1 - 1)} + \frac{s_2^4}{n_2^2 (n_2 - 1)}}</annotation></semantics></math></p>
<p>To perform a sufficient <code>Welch t-test</code> the following
assumptions about the input data need to be fulfilled to test whether
two samples come from populations with equal means:</p>
<p><strong>Assumptions about input data</strong></p>
<ul>
<li>independent samples</li>
<li>continuous data</li>
<li>(approximate) normality</li>
</ul>
<p>Nevertheless, although in most cases <code>log2</code> expression
levels are used to perform the <code>Welch t-test</code> assuming that
expression levels are log-normal distributed which approximates a normal
distribution in infinity, in most cases the small number of replicates
is not sufficient enough to fulfill the (approximate) normality
assumption of the <code>Welch t-test</code>.</p>
<p>Due to this fact, non-parametric, sampling based, or generalized
linear model based methods have been proposed to quantify p-values of
differential expression. Nevertheless, the <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>
function implements the <code>Welch t-test</code> for the detection of
differentially expressed genes, allowing users to compare the results
with more recent DEG detection methods/methodologies also implemented in
<code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>.</p>
<ul>
<li>
<strong>Advantages:</strong>
<ul>
<li>DEG detection based on statistical quantification</li>
<li>Parametric test resulting in a strong test statistic</li>
<li>Can handle small sample sizes</li>
</ul>
</li>
<li>
<strong>Disadvantages:</strong>
<ul>
<li>Test assumptions must be fulfilled to return sufficient
p-values<br>
</li>
<li>Can hardly assure normality with very sample sizes of
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>n</mi><mo>=</mo><mn>3</mn><mo>,</mo><mn>4</mn><mo>,</mo><mn>5</mn><mo>,</mo><mi>.</mi><mi>.</mi></mrow><annotation encoding="application/x-tex">n = 3,4,5,..</annotation></semantics></math>
(replicates)</li>
<li>Pairwise comparisons between different stages or experiments</li>
</ul>
</li>
</ul>
<div class="section level3">
<h3 id="example-welch-t-test">Example: Welch t-test<a class="anchor" aria-label="anchor" href="#example-welch-t-test"></a>
</h3>
<p>Performing <code>Welch t-test</code> with <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>
can be done by specifying <code>method = "t.test"</code>. Internally
<code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> performs a two sided <code>Welch t-test</code>.
This means that the <code>Welch t-test</code> quantifies only whether or
not a gene is significantly differentially expressed, but not the
direction of enrichment (over-expressed or under-expressed).</p>
<p>The <code>PhyloExpressionSetExample</code> we use in the following
example stores absolute expression levels. In case your
<code>ExpressionSet</code> also stores absolute expression levels (which
is likely due to the <code>ExpressionSet</code> standard for
Phylotranscriptomics analyses), you can use the <code><a href="../reference/tf.html">tf()</a></code>
function implemented in <code>myTAI</code> to transform absolute
expression levels to <code>log2</code> expression levels before
performing <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> with a <code>Welch t-test</code>,
e.g. <code>tf(PhyloExpressionSetExample[1:5,1:8],log2)</code>. In
general, using <code>log2</code> transformed expression levels as input
<code>ExpressionSet</code> of <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> allows us to (at
least) assume that samples (replicate expression levels) used to perform
the <code>Welch t-test</code> are log-normal distributed and therefore,
somewhat approximate normal distributed.</p>
<p>Please notice however, that RNA-Seq data can include count values of
0. So when transforming absolute counts to <code>log2</code> counts
infinity values of <code>log2(0) = -Inf</code> will be produced and
therefore, p-value computations will not be possible. To avoid this case
you could either remove RNA-Seq count values of 0 from the input dataset
using the <code><a href="../reference/Expressed.html">Expressed()</a></code> function (see section <em>Filter for
Expressed Genes</em>), e.g. pass
<code>tf(Expressed(PhyloExpressionSetExample[1:5,1:8], cut.off = 1),log2)</code>
as <code>ExpressionSet</code> argument to <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> or
shift all count values by a constant value, e.g. pass
<code>tf(PhyloExpressionSetExample[1:5,1:8], function(x) log2(x + 1))</code>
as <code>ExpressionSet</code> argument to <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>.</p>
<p>Internally, <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> will also check for 0 values in
input data and will automatically shift all expression levels by
<code>+1</code> in case 0 values are included.</p>
<div class="sourceCode" id="cb5"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by a Welch t-test</span></span>
<span><span class="va">ttest.DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="fu"><a href="../reference/tf.html">tf</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,<span class="va">log2</span><span class="op">)</span>,</span>
<span>                        nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                        method        <span class="op">=</span> <span class="st">"t.test"</span>,</span>
<span>                        stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">ttest.DEGs</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID     S1&lt;-&gt;S2    S1&lt;-&gt;S3    S2&lt;-&gt;S3
1            1 at1g01040.2 0.027832572 0.04020203 0.13481563
2            1 at1g01050.1 0.852379466 0.31471871 0.36326955
3            1 at1g01070.1 0.003200692 0.00113536 0.02236621
4            1 at1g01080.2 0.086426813 0.03092924 0.45999438
5            1 at1g01090.1 0.090387087 0.04638872 0.04978092</code></pre>
<p>The resulting <code>data.frame</code> stores the p-values of
stage-wise comparisons for each gene. To adjust p-values for multiple
testing of stage-wise comparisons you can specify the
<code>p.adjust.method</code> argument with one of the p-value adjustment
methods implemented in <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>.</p>
<p>In detail, correcting for multiple testing allows to appropriately
choose selection cut-offs for p-values fulfilling the differential
expression criteria. Hahne et al., 2008 (p. 87) give a nice example of
correcting for multiple testing to determine appropriate selection
cut-offs.</p>
<p>Please consult the documentation of <code><a href="https://rdrr.io/r/stats/p.adjust.html" class="external-link">?p.adjust</a></code> to see
which p-value adjustment methods are implemented in
<code><a href="../reference/DiffGenes.html">DiffGenes()</a></code>.</p>
<p>Please also consult these reviews <a href="http://www.stat.columbia.edu/~gelman/research/published/multiple2f.pdf" class="external-link">Gelman
et al., 2008</a>, and <a href="http://www.gs.washington.edu/academics/courses/akey/56008/lecture/lecture10.pdf" class="external-link">Slides</a>)
to decide whether or not to apply p-value adjustment to your own
dataset.</p>
<div class="sourceCode" id="cb7"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by a Welch t-test</span></span>
<span><span class="co"># and furthermore, adjust p-values for multiple comparison</span></span>
<span><span class="co"># using the Benjamini &amp; Hochberg (1995) method: method = "BH"</span></span>
<span><span class="va">ttest.DEGs.p_adjust</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet   <span class="op">=</span> <span class="fu"><a href="../reference/tf.html">tf</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,<span class="va">log2</span><span class="op">)</span>,</span>
<span>                                 nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                 method          <span class="op">=</span> <span class="st">"t.test"</span>,</span>
<span>                                 p.adjust.method <span class="op">=</span> <span class="st">"BH"</span>,</span>
<span>                                 stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span></span>
<span><span class="va">ttest.DEGs.p_adjust</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID    S1&lt;-&gt;S2   S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.06958143 0.0579859 0.2246927
2            1 at1g01050.1 0.85237947 0.3147187 0.4540869
3            1 at1g01070.1 0.01600346 0.0056768 0.1118311
4            1 at1g01080.2 0.11298386 0.0579859 0.4599944
5            1 at1g01090.1 0.11298386 0.0579859 0.1244523</code></pre>
<p>The resulting p-value adjusted <code>data.frame</code> can be used to
filter for differentially expressed genes. Here, specifying the
arguments: <code>comparison</code>, <code>alpha</code>, and
<code>filter.method</code> in <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> allows users to
obtain only significant differentially expressed genes.</p>
<div class="sourceCode" id="cb9"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># Detection of DEGs using the p-value returned by a Welch t-test</span></span>
<span><span class="co"># and furthermore, adjust p-values for multiple comparison</span></span>
<span><span class="co"># using the Benjamini &amp; Hochberg (1995) method: method = "BH"</span></span>
<span><span class="co"># and filter for significantly differentially expressed genes (alpha = 0.05) </span></span>
<span><span class="va">ttest.DEGs.p_adjust.filtered</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet   <span class="op">=</span> <span class="fu"><a href="../reference/tf.html">tf</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">10</span> ,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,<span class="va">log2</span><span class="op">)</span>,</span>
<span>                                          nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                          method          <span class="op">=</span> <span class="st">"t.test"</span>,</span>
<span>                                          p.adjust.method <span class="op">=</span> <span class="st">"BH"</span>,</span>
<span>                                          stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span>,</span>
<span>                                          comparison      <span class="op">=</span> <span class="st">"above"</span>,</span>
<span>                                          alpha           <span class="op">=</span> <span class="fl">0.05</span>,</span>
<span>                                          filter.method   <span class="op">=</span> <span class="st">"n-set"</span>,</span>
<span>                                          n               <span class="op">=</span> <span class="fl">1</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the genes fulfilling the filter criteria </span></span>
<span><span class="va">ttest.DEGs.p_adjust.filtered</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID    S1&lt;-&gt;S2   S1&lt;-&gt;S3   S2&lt;-&gt;S3
3            1 at1g01070.1 0.03200692 0.0113536 0.2192432</code></pre>
<p>In this example, only 1 out of 10 genes fulfills the p-value criteria
(<code>alpha = 0.05</code>) in at least one stage comparison.</p>
<div class="section level4">
<h4 id="rank-top-p-values">Rank top p-values<a class="anchor" aria-label="anchor" href="#rank-top-p-values"></a>
</h4>
<p>Finally, users can rank genes in increasing p-value order for each
stage comparison by typing:</p>
<div class="sourceCode" id="cb11"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="va">ttest.DEGs.p_adjust</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet   <span class="op">=</span> <span class="fu"><a href="../reference/tf.html">tf</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">500</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,<span class="va">log2</span><span class="op">)</span>,</span>
<span>                                 nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                 method          <span class="op">=</span> <span class="st">"t.test"</span>,</span>
<span>                                 p.adjust.method <span class="op">=</span> <span class="st">"BH"</span>,</span>
<span>                                 stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/head.html" class="external-link">head</a></span><span class="op">(</span><span class="va">ttest.DEGs.p_adjust</span><span class="op">[</span><span class="fu"><a href="https://rdrr.io/r/base/order.html" class="external-link">order</a></span><span class="op">(</span><span class="va">ttest.DEGs.p_adjust</span><span class="op">[</span> , <span class="st">"S1&lt;-&gt;S2"</span><span class="op">]</span>, decreasing <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span> , <span class="fl">1</span><span class="op">:</span><span class="fl">3</span><span class="op">]</span><span class="op">)</span></span></code></pre></div>
<pre><code>    Phylostratum      GeneID  S1&lt;-&gt;S2
54             1 at1g02400.1 0.151388
119            1 at1g03870.1 0.151388
137            1 at1g04380.1 0.151388
289            1 at1g08110.4 0.151388
383            1 at1g10360.1 0.151388
413            1 at1g11040.1 0.151388</code></pre>
<p>Here the line
<code>ttest.DEGs.p_adjust[order(ttest.DEGs.p_adjust[ , "S1&lt;-&gt;S2"], decreasing = FALSE) , 1:3]</code>
will sort p-values of stage comparison <code>"S1&lt;-&gt;S2"</code> in
increasing order.</p>
</div>
</div>
</div>
<div class="section level2">
<h2 id="wilcoxon-mann-whitney-test-mann-whitney-u-test">Wilcoxon-Mann-Whitney test (Mann-Whitney U test)<a class="anchor" aria-label="anchor" href="#wilcoxon-mann-whitney-test-mann-whitney-u-test"></a>
</h2>
<p>The Wilcoxon-Mann-Whitney test is a <em>nonparametric</em> test to
quantify the shift in empirical distribution parameters.
<em>Nonparametric</em> tests are useful when sample populations do not
meet the test assumptions of <em>parametric</em> tests.</p>
<div class="sourceCode" id="cb13"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by a Wilcoxon-Mann-Whitney test</span></span>
<span><span class="va">Wilcox.DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                        nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                        method        <span class="op">=</span> <span class="st">"wilcox.test"</span>,</span>
<span>                        stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">Wilcox.DEGs</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID   S1&lt;-&gt;S2   S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.3333333 0.3333333 0.3333333
2            1 at1g01050.1 1.0000000 0.3333333 0.3333333
3            1 at1g01070.1 0.3333333 0.3333333 0.3333333
4            1 at1g01080.2 0.3333333 0.3333333 0.6666667
5            1 at1g01090.1 0.3333333 0.3333333 0.3333333</code></pre>
<p>Again, users can adjust p-values by specifying the
<code>p.adjust.method</code> argument.</p>
<div class="sourceCode" id="cb15"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by a Wilcoxon-Mann-Whitney test</span></span>
<span><span class="co"># and furthermore, adjust p-values for multiple comparison</span></span>
<span><span class="co"># using the Benjamini &amp; Hochberg (1995) method: method = "BH"</span></span>
<span><span class="co"># and filter for significantly differentially expressed genes (alpha = 0.05)</span></span>
<span><span class="va">Wilcox.DEGs.adj</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet  <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                            nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                            method          <span class="op">=</span> <span class="st">"wilcox.test"</span>,</span>
<span>                            stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span>,</span>
<span>                            p.adjust.method <span class="op">=</span> <span class="st">"BH"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">Wilcox.DEGs.adj</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID   S1&lt;-&gt;S2   S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.4166667 0.3333333 0.4166667
2            1 at1g01050.1 1.0000000 0.3333333 0.4166667
3            1 at1g01070.1 0.4166667 0.3333333 0.4166667
4            1 at1g01080.2 0.4166667 0.3333333 0.6666667
5            1 at1g01090.1 0.4166667 0.3333333 0.4166667</code></pre>
</div>
<div class="section level2">
<h2 id="negative-binomial-exact-tests">Negative Binomial (Exact Tests)<a class="anchor" aria-label="anchor" href="#negative-binomial-exact-tests"></a>
</h2>
<p>Exact Tests for Differences between two groups of negative-binomial
counts implemented in <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> are based on the
<code>edgeR</code> function <code>exactTest()</code>. Please consult the
<a href="https://www.bioconductor.org/packages/release/bioc/html/edgeR.html" class="external-link">edgeR
Users Guide</a> for mathematical details.</p>
<div class="section level3">
<h3 id="install-edger-package">Install edgeR Package<a class="anchor" aria-label="anchor" href="#install-edger-package"></a>
</h3>
<p>The detection of DEGs using negative binomial models is based on the
powerful implementations provided by the <a href="https://www.bioconductor.org/packages/release/bioc/html/edgeR.html" class="external-link">edgeR</a>
package. Hence, before using the negative binomial models in
<code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> users need to install the edgeR package.</p>
<div class="sourceCode" id="cb17"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># install edgeR</span></span>
<span><span class="kw"><a href="https://rdrr.io/r/base/source.html" class="external-link">source</a></span><span class="op">(</span><span class="st">"http://bioconductor.org/biocLite.R"</span><span class="op">)</span></span>
<span><span class="fu">biocLite</span><span class="op">(</span><span class="st">"edgeR"</span><span class="op">)</span></span></code></pre></div>
</div>
<div class="section level3">
<h3 id="double-tail-method">Double Tail Method<a class="anchor" aria-label="anchor" href="#double-tail-method"></a>
</h3>
<p>This method computes two-sided p-values by doubling the smaller tail
probability (see <code>?exactTestByDeviance</code> for details). To
compute p-values for stagewise comparisons based on negative binomial
models, the <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> argument
<code>method = "doubletail"</code>, the number of replicates per stage
<code>nrep</code>, and <code>lib.size</code> quantifying the library
size to equalize sample library sizes by quantile-to-quantile
normalization need to be specified (see also
<code>?equalizeLibSizes</code>).</p>
<div class="sourceCode" id="cb18"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by the Double Tail Method</span></span>
<span><span class="va">DoubleTail.DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                        nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                        method        <span class="op">=</span> <span class="st">"doubletail"</span>,</span>
<span>                        lib.size      <span class="op">=</span> <span class="fl">1000</span>,</span>
<span>                        stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">DoubleTail.DEGs</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID    S1&lt;-&gt;S2     S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.26026604 0.110233012 0.6304508
2            1 at1g01050.1 0.95314428 0.598102712 0.6398757
3            1 at1g01070.1 0.55461941 0.456018563 0.8774231
4            1 at1g01080.2 0.58130025 0.487028051 0.8860005
5            1 at1g01090.1 0.03615134 0.001773543 0.2645537</code></pre>
<p>Again, users can adjust p-values by specifying the
<code>p.adjust.method</code> argument.</p>
<div class="sourceCode" id="cb20"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by the Double Tail Method</span></span>
<span><span class="co"># and furthermore, adjust p-values for multiple comparison</span></span>
<span><span class="co"># using the Benjamini &amp; Hochberg (1995) method: method = "BH"</span></span>
<span><span class="co"># and filter for significantly differentially expressed genes (alpha = 0.05)</span></span>
<span><span class="va">DoubleTail.DEGs.adj</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet  <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                                nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                method          <span class="op">=</span> <span class="st">"doubletail"</span>,</span>
<span>                                lib.size        <span class="op">=</span> <span class="fl">1000</span>,</span>
<span>                                stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span>,</span>
<span>                                p.adjust.method <span class="op">=</span> <span class="st">"BH"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">DoubleTail.DEGs.adj</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID   S1&lt;-&gt;S2     S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.6506651 0.275582530 0.8860005
2            1 at1g01050.1 0.9531443 0.598102712 0.8860005
3            1 at1g01070.1 0.7266253 0.598102712 0.8860005
4            1 at1g01080.2 0.7266253 0.598102712 0.8860005
5            1 at1g01090.1 0.1807567 0.008867715 0.8860005</code></pre>
</div>
<div class="section level3">
<h3 id="small-p-method">Small-P Method<a class="anchor" aria-label="anchor" href="#small-p-method"></a>
</h3>
<p>This method performs the method of small probabilities as proposed by
Robinson and Smyth (2008) (see <code>exactTestBySmallP</code> for
details). To compute p-values for stagewise comparisons based on
negative binomial models, the <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> argument
<code>method = "doubletail"</code>, the number of replicates per stage
<code>nrep</code>, and <code>lib.size</code> quantifying the library
size to equalize sample library sizes by quantile-to-quantile
normalization need to be specified (see also
<code>?equalizeLibSizes</code>).</p>
<div class="sourceCode" id="cb22"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by the Small-P Method</span></span>
<span><span class="va">SmallP.DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                        nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                        method        <span class="op">=</span> <span class="st">"smallp"</span>,</span>
<span>                        lib.size      <span class="op">=</span> <span class="fl">1000</span>,</span>
<span>                        stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">SmallP.DEGs</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID    S1&lt;-&gt;S2     S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.26026604 0.110233012 0.6304508
2            1 at1g01050.1 0.95314428 0.598102712 0.6398757
3            1 at1g01070.1 0.55461941 0.456018563 0.8774231
4            1 at1g01080.2 0.58130025 0.487028051 0.8860005
5            1 at1g01090.1 0.03615134 0.001773543 0.2645537</code></pre>
<p>Again, users can adjust p-values by specifying the
<code>p.adjust.method</code> argument.</p>
<div class="sourceCode" id="cb24"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by the Small-P Method</span></span>
<span><span class="co"># and furthermore, adjust p-values for multiple comparison</span></span>
<span><span class="co"># using the Benjamini &amp; Hochberg (1995) method: method = "BH"</span></span>
<span><span class="co"># and filter for significantly differentially expressed genes (alpha = 0.05)</span></span>
<span><span class="va">SmallP.DEGs.adj</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet  <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                                nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                method          <span class="op">=</span> <span class="st">"smallp"</span>,</span>
<span>                                lib.size        <span class="op">=</span> <span class="fl">1000</span>,</span>
<span>                                stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span>,</span>
<span>                                p.adjust.method <span class="op">=</span> <span class="st">"BH"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">SmallP.DEGs.adj</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID   S1&lt;-&gt;S2     S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.6506651 0.275582530 0.8860005
2            1 at1g01050.1 0.9531443 0.598102712 0.8860005
3            1 at1g01070.1 0.7266253 0.598102712 0.8860005
4            1 at1g01080.2 0.7266253 0.598102712 0.8860005
5            1 at1g01090.1 0.1807567 0.008867715 0.8860005</code></pre>
</div>
<div class="section level3">
<h3 id="deviance-method">Deviance Method<a class="anchor" aria-label="anchor" href="#deviance-method"></a>
</h3>
<p>This method uses the deviance goodness of fit statistics to define
the rejection region, and is therefore equivalent to a conditional
likelihood ratio test (see <code>exactTestByDeviance</code> for
details). To compute p-values for stagewise comparisons based on
negative binomial models, the <code><a href="../reference/DiffGenes.html">DiffGenes()</a></code> argument
<code>method = "doubletail"</code>, the number of replicates per stage
<code>nrep</code>, and <code>lib.size</code> quantifying the library
size to equalize sample library sizes by quantile-to-quantile
normalization need to be specified (see also
<code>?equalizeLibSizes</code>).</p>
<div class="sourceCode" id="cb26"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by the Deviance</span></span>
<span><span class="va">Deviance.DEGs</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                        nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                        method        <span class="op">=</span> <span class="st">"deviance"</span>,</span>
<span>                        lib.size      <span class="op">=</span> <span class="fl">1000</span>,</span>
<span>                        stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">Deviance.DEGs</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID    S1&lt;-&gt;S2     S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.26026604 0.110233012 0.6304508
2            1 at1g01050.1 0.95314428 0.598102712 0.6398757
3            1 at1g01070.1 0.55461941 0.456018563 0.8774231
4            1 at1g01080.2 0.58130025 0.487028051 0.8860005
5            1 at1g01090.1 0.03615134 0.001773543 0.2645537</code></pre>
<p>Again, users can adjust p-values by specifying the
<code>p.adjust.method</code> argument.</p>
<div class="sourceCode" id="cb28"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="st">"PhyloExpressionSetExample"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># Detection of DEGs using the p-value returned by the Deviance Method</span></span>
<span><span class="co"># and furthermore, adjust p-values for multiple comparison</span></span>
<span><span class="co"># using the Benjamini &amp; Hochberg (1995) method: method = "BH"</span></span>
<span><span class="co"># and filter for significantly differentially expressed genes (alpha = 0.05)</span></span>
<span><span class="va">Deviance.DEGs.adj</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/DiffGenes.html">DiffGenes</a></span><span class="op">(</span>ExpressionSet    <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">5</span>,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                                nrep            <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                method          <span class="op">=</span> <span class="st">"deviance"</span>,</span>
<span>                                lib.size        <span class="op">=</span> <span class="fl">1000</span>,</span>
<span>                                stage.names     <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span>,</span>
<span>                                p.adjust.method <span class="op">=</span> <span class="st">"BH"</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># look at the results</span></span>
<span><span class="va">Deviance.DEGs.adj</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID   S1&lt;-&gt;S2     S1&lt;-&gt;S3   S2&lt;-&gt;S3
1            1 at1g01040.2 0.6506651 0.275582530 0.8860005
2            1 at1g01050.1 0.9531443 0.598102712 0.8860005
3            1 at1g01070.1 0.7266253 0.598102712 0.8860005
4            1 at1g01080.2 0.7266253 0.598102712 0.8860005
5            1 at1g01090.1 0.1807567 0.008867715 0.8860005</code></pre>
</div>
</div>
<div class="section level2">
<h2 id="replicate-quality-check">Replicate Quality Check<a class="anchor" aria-label="anchor" href="#replicate-quality-check"></a>
</h2>
<p>Users can also perform replicate quality checks to quantify the
variability between replicate expression levels fo each stage
separately.</p>
<p>The <code><a href="../reference/PlotReplicateQuality.html">PlotReplicateQuality()</a></code> is designed to perform
customized replicate variability checks for any
<code>ExpressionSet</code> object storing replicates.</p>
<div class="sourceCode" id="cb30"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># visualize the sd() between replicates</span></span>
<span><span class="fu"><a href="../reference/PlotReplicateQuality.html">PlotReplicateQuality</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span> , <span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                     nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                     legend.pos   <span class="op">=</span> <span class="st">"topright"</span>,</span>
<span>                     ylim          <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">0</span>,<span class="fl">0.2</span><span class="op">)</span>,</span>
<span>                     lwd           <span class="op">=</span> <span class="fl">6</span><span class="op">)</span></span></code></pre></div>
<p>The resulting plot visualizes the kernel density estimates for the
variance (log variance) between replicates. Each curve represents the
density function for the replicate variation within one stage or
experiment. In this case the variance between replicates of
<code>Stage 1</code> to <code>Stage 3</code> (each including 2
replicates) seem to deviate from each other allowing the conclusion that
each stage has a different expression level variability between
replicates.</p>
<p>The <code>FUN</code> argument implemented in
<code><a href="../reference/PlotReplicateQuality.html">PlotReplicateQuality()</a></code> allows users to furthermore, specify
customized criteria quantifying replicate variability. Please notice
that the function specified in <code>FUN</code> will be performed
separately on each gene and stage.</p>
<p>In the following example the median absolute deviation function
<code><a href="https://rdrr.io/r/stats/mad.html" class="external-link">mad()</a></code> is used to quantify replicate variability.</p>
<div class="sourceCode" id="cb31"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># visualize the mad() between replicates</span></span>
<span><span class="fu"><a href="../reference/PlotReplicateQuality.html">PlotReplicateQuality</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span> , <span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span>,</span>
<span>                     nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                     FUN           <span class="op">=</span> <span class="va">mad</span>,</span>
<span>                     legend.pos    <span class="op">=</span> <span class="st">"topright"</span>,</span>
<span>                     ylim          <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">0</span>,<span class="fl">0.015</span><span class="op">)</span>,</span>
<span>                     lwd           <span class="op">=</span> <span class="fl">6</span><span class="op">)</span></span></code></pre></div>
<p>In general, users are not limited to specific functions implemented
in R. By writing customized functions such as
<code>FUN = function(x) return((x - mean(x))^2)</code> users can define
their own criteria to quantify replicate variability and can then apply
this criteria to <code><a href="../reference/PlotReplicateQuality.html">PlotReplicateQuality()</a></code> by specifying the
<code>FUN</code> argument.</p>
</div>
<div class="section level2">
<h2 id="collapsing-replicate-samples">Collapsing Replicate Samples<a class="anchor" aria-label="anchor" href="#collapsing-replicate-samples"></a>
</h2>
<p>After performing differential gene expression analyses, replicate
expression levels are collapsed to a single stage specific expression
level. For this purpose, <code>myTAI</code> implements the
<code><a href="../reference/CollapseReplicates.html">CollapseReplicates()</a></code> function, allowing users to combine
replicate expression levels stored in a standard
<code>PhyloExpressionSet</code> or <code>DivergenceExpressionSet</code>
object to a stage specific expression level using a specified window
function.</p>
<div class="sourceCode" id="cb32"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="kw"><a href="https://rdrr.io/r/base/library.html" class="external-link">library</a></span><span class="op">(</span><span class="va"><a href="https://drostlab.github.io/myTAI/">myTAI</a></span><span class="op">)</span></span>
<span></span>
<span><span class="co"># load example data</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># generate an example PhyloExpressionSet with replicates</span></span>
<span><span class="va">ExampleReplicateExpressionSet</span> <span class="op">&lt;-</span> <span class="va">PhyloExpressionSetExample</span><span class="op">[</span> ,<span class="fl">1</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span></span>
<span></span>
<span><span class="co"># rename stages</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/names.html" class="external-link">names</a></span><span class="op">(</span><span class="va">ExampleReplicateExpressionSet</span><span class="op">)</span><span class="op">[</span><span class="fl">3</span><span class="op">:</span><span class="fl">8</span><span class="op">]</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"Stage_1_Repl_1"</span>,<span class="st">"Stage_1_Repl_2"</span>,</span>
<span>                                               <span class="st">"Stage_2_Repl_1"</span>,<span class="st">"Stage_2_Repl_2"</span>,</span>
<span>                                               <span class="st">"Stage_3_Repl_1"</span>,<span class="st">"Stage_3_Repl_2"</span><span class="op">)</span></span>
<span><span class="co"># have a look at the example dataset</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/head.html" class="external-link">head</a></span><span class="op">(</span><span class="va">ExampleReplicateExpressionSet</span>, <span class="fl">5</span><span class="op">)</span></span></code></pre></div>
<pre><code>  Phylostratum      GeneID Stage_1_Repl_1 Stage_1_Repl_2 Stage_2_Repl_1 Stage_2_Repl_2 Stage_3_Repl_1 Stage_3_Repl_2
1            1 at1g01040.2       2173.635      1911.2001       1152.555      1291.4224       1000.253       962.9772
2            1 at1g01050.1       1501.014      1817.3086       1665.309      1564.7612       1496.321      1114.6435
3            1 at1g01070.1       1212.793      1233.0023        939.200       929.6195        864.218       877.2060
4            1 at1g01080.2       1016.920       936.3837       1181.338      1329.4734       1392.643      1287.9746
5            1 at1g01090.1      11424.567     16778.1685      34366.649     39775.6405      56231.569     66980.3673</code></pre>
<p>Now, assume that this example <code>PhyloExpressionSet</code> stores
three developmental stages and 2 biological replicates for each
developmental stage. Of course, we could now compute and visualize the
TAI profile by typing:</p>
<div class="sourceCode" id="cb34"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># visualize the TAI profile over 3 stages of development</span></span>
<span><span class="co"># and 2 replicates per stage</span></span>
<span><span class="fu"><a href="../reference/PlotPattern.html">PlotPattern</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">ExampleReplicateExpressionSet</span>,</span>
<span>            type          <span class="op">=</span> <span class="st">"l"</span>,</span>
<span>            lwd           <span class="op">=</span> <span class="fl">6</span><span class="op">)</span></span></code></pre></div>
<p>Usually, one would expect that variations in replicate values are
smaller than variations between developmental stages. In this example
however, we constructed replicate values that vary larger than
expression levels between developmental stages. For many applications it
might be useful to visualize TAI/TDI values of replicates as well, but
normally replicate values are collapsed to one gene and stage specific
value after differential gene expression analyses and replicate quality
control have been performed.</p>
<p>The following example illustrates how to collapse replicates with
<code><a href="../reference/CollapseReplicates.html">CollapseReplicates()</a></code>:</p>
<div class="sourceCode" id="cb35"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># combine the expression levels of the 2 replicates (const) per stage</span></span>
<span><span class="co"># using geom.mean as window function and rename new stages to: "S1","S2","S3"</span></span>
<span><span class="va">CollapssedPhyloExpressionSet</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/CollapseReplicates.html">CollapseReplicates</a></span><span class="op">(</span></span>
<span>                                       ExpressionSet <span class="op">=</span> <span class="va">ExampleReplicateExpressionSet</span>,</span>
<span>                                       nrep          <span class="op">=</span> <span class="fl">2</span>,</span>
<span>                                       FUN           <span class="op">=</span> <span class="va">geom.mean</span>,</span>
<span>                                       stage.names   <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="st">"S1"</span>,<span class="st">"S2"</span>,<span class="st">"S3"</span><span class="op">)</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># have a look at the collapsed PhyloExpressionSet</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/head.html" class="external-link">head</a></span><span class="op">(</span><span class="va">CollapssedPhyloExpressionSet</span><span class="op">)</span></span></code></pre></div>
<pre><code>   Phylostratum      GeneID         S1         S2         S3
1            1 at1g01040.2  2038.1982  1220.0147   981.4381
2            1 at1g01050.1  1651.6070  1614.2524  1291.4582
3            1 at1g01070.1  1222.8557   934.3975   870.6878
4            1 at1g01080.2   975.8215  1253.2189  1339.2866
5            1 at1g01090.1 13844.9740 36972.3612 61371.0937
6            1 at1g01120.1   815.3288   894.8987   905.8272</code></pre>
<p>The <code>nrep</code> argument specifies either a constant number of
replicates per stage or a numeric vector storing variable numbers of
replicates for each developmental stage. In our example, each
developmental stage had a constant (equal) number of replicates per
developmental stage (<code>nrep = 2</code>). In case a variable stage
specific number of replicates is present, one could specify
<code>nrep = c(2,3,2)</code> defining the case that developmental stage
1 stores 2 replicates, stage 2 stores 3 replicates, and stage 3 again,
stores 2 replicates.</p>
<p>The argument <code>FUN</code> specifies the window function to
collapse replicate expression levels to a single stage specific value.
In this example, we chose the <code><a href="../reference/geom.mean.html">geom.mean()</a></code> (geometric mean)
function implemented in <code>myTAI</code>, because our example
<code>PhyloExpressionSet</code> stores absolute expression levels.
Notice that the mathematical equivalent of performing arithmetic mean
(<code><a href="https://rdrr.io/r/base/mean.html" class="external-link">mean()</a></code>) computations on <code>log</code> expression levels
is to perform the geometric mean (<code><a href="../reference/geom.mean.html">geom.mean()</a></code>) on absolute
expression levels.</p>
<p>The <code>stage.names</code> argument then specifies the new names of
collapsed stages.</p>
</div>
<div class="section level2">
<h2 id="filter-for-expressed-genes">Filter for Expressed Genes<a class="anchor" aria-label="anchor" href="#filter-for-expressed-genes"></a>
</h2>
<p>After differential gene expression analyses and replicate aggregation
have been performed, some studies filter gene expression levels in
RNA-Seq count tables or microarray expression matrices for non-expressed
or outlier genes. For example, in most studies performing RNA-Seq
experiments FPKM/RPKM values &lt; 1 are remove from the processed
(final) count table.</p>
<p>For this purpose <code>myTAI</code> implements the
<code><a href="../reference/Expressed.html">Expressed()</a></code> function to filter (remove) expression levels
in RNA-Seq count tables or microarray expression matrices which do not
pass a defined expression threshold.</p>
<p>The <code><a href="../reference/Expressed.html">Expressed()</a></code> function takes a standard
<code>PhyloExpressionSet</code> or <code>DivergenceExpressionSet</code>
object storing a RNA-Seq count table (CT) or microarray gene expression
matrix and removes genes from this count table or gene expression matrix
that have an expression level below a defined <code>cut.off</code>
value.</p>
<p><code><a href="../reference/Expressed.html">Expressed()</a></code> allows users to choose from several gene
extraction methods (see <code><a href="../reference/Expressed.html">?Expressed</a></code> for details):</p>
<ul>
<li><p><code>const</code>: all genes that have at least one stage that
undercuts or exceeds the expression <code>cut.off</code> will be
excluded from the <code>ExpressionSet</code>. Hence, for a 7 stage
<code>ExpressionSet</code> genes passing the expression level
<code>cut.off</code> in 6 stages will be retained in the
<code>ExpressionSet</code>.</p></li>
<li><p><code>min-set</code>: genes passing the expression level
<code>cut.off</code> in <code>ceiling(n/2)</code> stages will be
retained in the <code>ExpressionSet</code>, where <code>n</code> is the
number of stages in the <code>ExpressionSet</code>.</p></li>
<li><p><code>n-set</code>: genes passing the expression level
<code>cut.off</code> in <code>n</code> stages will be retained in the
<code>ExpressionSet</code>. Here, the argument <code>n</code> is
defining the number of stages for which the threshold criteria should be
fulfilled.</p></li>
</ul>
<div class="sourceCode" id="cb37"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># check number of genes in PhyloExpressionSetExample</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span><span class="co">#&gt; [1] 25260</span></span>
<span></span>
<span><span class="co"># remove genes that have an expression level below 8000</span></span>
<span><span class="co"># in at least one developmental stage</span></span>
<span><span class="va">FilterConst</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/Expressed.html">Expressed</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span>,</span>
<span>                         cut.off       <span class="op">=</span> <span class="fl">8000</span>,</span>
<span>                         comparison    <span class="op">=</span> <span class="st">"below"</span>, </span>
<span>                         method        <span class="op">=</span> <span class="st">"const"</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">FilterConst</span><span class="op">)</span> <span class="co"># check number of retained genes</span></span>
<span><span class="co">#&gt; [1] 449</span></span></code></pre></div>
<p>Users will observe that only 449 out of 25260 genes in
<code>PhyloExpressionSetExample</code> have an absolute expression level
above <code>8000</code> when omitting genes using
<code>method = 'const'</code>. The argument <code>comparison</code>
specifies whether genes having expression levels below, above, or below
AND above (both) the <code>cut.off</code> value should be removed from
the dataset.</p>
<p>The following comparison methods can be selected:</p>
<ul>
<li>
<code>comparison = "below"</code>: define genes as not expressed
which undercut the <code>cut-off</code> threshold.</li>
<li>
<code>comparison = "above"</code>: define genes as outliers which
exceed the <code>cut-off</code> threshold.</li>
<li>
<code>comparison = "both"</code>: remove genes fulfilling the
<code>comparison = "below"</code> <strong>AND</strong>
<code>comparison = "above"</code> criteria.</li>
</ul>
<div class="sourceCode" id="cb38"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># again: check number of genes in PhyloExpressionSetExample</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span><span class="co">#&gt; [1] 25260</span></span>
<span></span>
<span><span class="co"># remove genes that have an expression level above 12000</span></span>
<span><span class="co"># in at least one developmental stage (outlier removal)</span></span>
<span><span class="va">FilterConst.above</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/Expressed.html">Expressed</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span>,</span>
<span>                               cut.off       <span class="op">=</span> <span class="fl">12000</span>,</span>
<span>                               comparison    <span class="op">=</span> <span class="st">"above"</span>, </span>
<span>                               method        <span class="op">=</span> <span class="st">"const"</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">FilterConst.above</span><span class="op">)</span> <span class="co"># check number of retained genes</span></span>
<span><span class="co">#&gt; [1] 23547</span></span></code></pre></div>
<p>For this example 25260 - 23547 = 1713 have been classified as
outliers (expression levels above 12000) and were removed from the
dataset.</p>
<div class="sourceCode" id="cb39"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># again: check number of genes in PhyloExpressionSetExample</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span><span class="co">#&gt; [1] 25260</span></span>
<span></span>
<span><span class="co"># remove genes that have an expression level below 8000 AND above 12000</span></span>
<span><span class="co"># in at least one developmental stage (non-expressed genes AND outlier removal)</span></span>
<span><span class="va">FilterConst.both</span> <span class="op">&lt;-</span>  <span class="fu"><a href="../reference/Expressed.html">Expressed</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span>,</span>
<span>                               cut.off       <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">8000</span>,<span class="fl">12000</span><span class="op">)</span>,</span>
<span>                               comparison    <span class="op">=</span> <span class="st">"both"</span>, </span>
<span>                               method        <span class="op">=</span> <span class="st">"const"</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">FilterConst.both</span><span class="op">)</span> <span class="co"># check number of retained genes</span></span>
<span><span class="co">#&gt; [1] 2</span></span></code></pre></div>
<p>When selecting <code>comparison = 'both'</code>, the
<code>cut.off</code> argument receives 2 threshold values: the
<em>below</em> <code>cut.off</code> as first element and the
<em>above</em> <code>cut.off</code> as second element. In this case
<code>cut.off = c(8000,12000)</code>. Here, only 2 genes fulfill these
criteria.</p>
<p>Analogously, users can specify the number of stages that should
fulfill the threshold criteria using the <code>n-set</code> method.</p>
<div class="sourceCode" id="cb40"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># remove genes that have an expression level below 8000</span></span>
<span><span class="co"># in at least 5 developmental stages (in this case: n = 2 stages fulfilling the criteria)</span></span>
<span><span class="va">FilterNSet</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/Expressed.html">Expressed</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span>,</span>
<span>                        cut.off       <span class="op">=</span> <span class="fl">8000</span>,</span>
<span>                        method        <span class="op">=</span> <span class="st">"n-set"</span>,</span>
<span>                        comparison    <span class="op">=</span> <span class="st">"below"</span>,</span>
<span>                        n             <span class="op">=</span> <span class="fl">2</span><span class="op">)</span></span>
<span></span>
<span><span class="fu"><a href="https://rdrr.io/r/base/nrow.html" class="external-link">nrow</a></span><span class="op">(</span><span class="va">FilterMinSet</span><span class="op">)</span> <span class="co"># check number of retained genes</span></span>
<span><span class="co">#&gt; [1] 20</span></span></code></pre></div>
<p>Here, 20 genes are fulfilling these criteria.</p>
</div>
<div class="section level2">
<h2 id="compute-the-statistical-significance-of-each-replicate-combination">Compute the Statistical Significance of Each Replicate
Combination<a class="anchor" aria-label="anchor" href="#compute-the-statistical-significance-of-each-replicate-combination"></a>
</h2>
<p>In some cases (high variability of replicates) it might be useful to
verify that there is no sequence of replicates (for all possible
combination of replicates) that results in a non-significant
<code>TAI</code> or <code>TDI</code> pattern, when the initial pattern
with combined replicates was shown to be significant.</p>
<p>The <code><a href="../reference/CombinatorialSignificance.html">CombinatorialSignificance()</a></code> function implemented in
<code>myTAI</code> allows users to compute the p-values quantifying the
statistical significance of the underlying pattern for all combinations
of replicates.</p>
<div class="section level3">
<h3 id="a-small-example">A small Example:<a class="anchor" aria-label="anchor" href="#a-small-example"></a>
</h3>
<p>Assume a <code>PhyloExpressionSet</code> stores 3 developmental
stages with 3 replicates measured for each stage. The 9 replicates in
total are denoted as:
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mn>1.1</mn><mo>,</mo><mn>1.2</mn><mo>,</mo><mn>1.3</mn><mo>,</mo><mn>2.1</mn><mo>,</mo><mn>2.2</mn><mo>,</mo><mn>2.3</mn><mo>,</mo><mn>3.1</mn><mo>,</mo><mn>3.2</mn><mo>,</mo><mn>3.3</mn></mrow><annotation encoding="application/x-tex">1.1, 1.2, 1.3, 2.1, 2.2, 2.3, 3.1, 3.2, 3.3</annotation></semantics></math>.
Now the function computes the statistical significance of each pattern
derived by the corresponding combination of replicates, e.g.</p>
<ul>
<li><p>1.1, 2.1, 3.1 : p-value for combination 1</p></li>
<li><p>1.1, 2.2, 3.1 : p-value for combination 2</p></li>
<li><p>1.1, 2.3, 3.1 : p-value for combination 3</p></li>
<li><p>1.2, 2.1, 3.1 : p-value for combination 4</p></li>
<li><p>1.2, 2.1, 3.1 : p-value for combination 5</p></li>
<li><p>1.2, 2.1, 3.1 : p-value for combination 6</p></li>
<li><p>1.3, 2.1, 3.1 : p-value for combination 7</p></li>
<li><p>1.3, 2.2, 3.1 : p-value for combination 8</p></li>
<li><p>1.3, 2.3, 3.1 : p-value for combination 9</p></li>
<li><p>…</p></li>
</ul>
<p>This procedure yields 27 p-values for the
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mn>3</mn><mn>3</mn></msup><annotation encoding="application/x-tex">3^3</annotation></semantics></math>
(<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>n</mi><mi>m</mi></msup><annotation encoding="application/x-tex">n^m</annotation></semantics></math>)
replicate combinations, where
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>n</mi><annotation encoding="application/x-tex">n</annotation></semantics></math>
denotes the number of developmental stages and
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>m</mi><annotation encoding="application/x-tex">m</annotation></semantics></math>
denotes the number of replicates per stage.</p>
<p>Note that in case users have a large amount of stages/experiments and
a large amount of replicates the computation time will increase by
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mi>n</mi><mi>m</mi></msup><annotation encoding="application/x-tex">n^m</annotation></semantics></math>.
For 11 stages and 4 replicates,
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msup><mn>4</mn><mn>11</mn></msup><annotation encoding="application/x-tex">4^{11}</annotation></semantics></math>
= 4194304 p-values have to be computed. Each p-value computation itself
is based on a permutation test running with
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mn>1</mn><mo>,</mo><mn>000</mn><mo>,</mo><mn>10</mn><mo>,</mo><mn>000</mn><mo>,</mo><mi>.</mi><mi>.</mi><mi>.</mi></mrow><annotation encoding="application/x-tex">1,000, 10,000, ...</annotation></semantics></math>
or more permutations. Be aware that this might take some time.</p>
<p>The p-value vector returned by this function can then be used to plot
the p-values to see whether an critical value
<math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>α</mi><annotation encoding="application/x-tex">\alpha</annotation></semantics></math>
is exceeded or not
(e.g. <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>α</mi><mo>=</mo><mn>0.05</mn></mrow><annotation encoding="application/x-tex">\alpha = 0.05</annotation></semantics></math>).</p>
<div class="sourceCode" id="cb41"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="co"># load a standard PhyloExpressionSet</span></span>
<span><span class="fu"><a href="https://rdrr.io/r/utils/data.html" class="external-link">data</a></span><span class="op">(</span><span class="va">PhyloExpressionSetExample</span><span class="op">)</span></span>
<span></span>
<span><span class="co"># we assume that the PhyloExpressionSetExample</span></span>
<span><span class="co"># consists of 3 developmental stages</span></span>
<span><span class="co"># and 2 replicates for stage 1, 3 replicates for stage 2,</span></span>
<span><span class="co"># and 2 replicates for stage 3</span></span>
<span><span class="co"># FOR REAL ANALYSES PLEASE USE: permutations = 1000 or 10000</span></span>
<span><span class="co"># BUT NOTE THAT THIS TAKES MUCH MORE COMPUTATION TIME</span></span>
<span><span class="va">p.vector</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/CombinatorialSignificance.html">CombinatorialSignificance</a></span><span class="op">(</span>ExpressionSet <span class="op">=</span> <span class="va">PhyloExpressionSetExample</span>,</span>
<span>                                      replicates    <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html" class="external-link">c</a></span><span class="op">(</span><span class="fl">2</span>,<span class="fl">3</span>,<span class="fl">2</span><span class="op">)</span>,</span>
<span>                                      TestStatistic <span class="op">=</span> <span class="st">"FlatLineTest"</span>,</span>
<span>                                      permutations  <span class="op">=</span> <span class="fl">100</span>,</span>
<span>                                      parallel      <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></span></code></pre></div>
<pre><code> [1] 2.436296e-03 2.288593e-02 1.608399e-03 1.185615e-02 1.835306e-06 1.077012e-05
 [7] 2.025515e-07 5.148342e-07 1.654885e-07 6.251145e-06 9.265520e-10 1.047479e-06</code></pre>
<p>Users will observe that none of the replicate combinations resulted
in p-values &gt; 0.05 and thus we can assume that the
phylotranscriptomic pattern computed based on collapsed replicates is
not biased by insignificant replicate combinations.</p>
<div class="sourceCode" id="cb43"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span><span class="fu"><a href="https://rdrr.io/r/base/any.html" class="external-link">any</a></span><span class="op">(</span><span class="va">p.vector</span> <span class="op">&gt;</span> <span class="fl">0.05</span><span class="op">)</span></span>
<span><span class="co">#&gt; FALSE</span></span></code></pre></div>
<p><code><a href="../reference/CombinatorialSignificance.html">CombinatorialSignificance()</a></code> can perform all significance
tests introduced in the <a href="Introduction.html">Introduction</a> and
<a href="Intermediate.html">Intermediate</a> vignettes.</p>
<p>Furthermore, the <code>parallel</code> argument allows users to
perform significance computations in parallel on a multicore machine.
This will speed up p-value computations for a large number of
combinations.</p>
</div>
</div>
  </main><aside class="col-md-3"><nav id="toc" aria-label="Table of contents"><h2>On this page</h2>
    </nav></aside>
</div>



    <footer><div class="pkgdown-footer-left">
  <p>Developed by Hajk-Georg Drost, Jaruwatana Sodai Lotharukpong, Nikola Kalábová, Kristian K Ullrich.</p>
</div>

<div class="pkgdown-footer-right">
  <p>Site built with <a href="https://pkgdown.r-lib.org/" class="external-link">pkgdown</a> 2.1.0.</p>
</div>

    </footer>
</div>





  </body>
</html>
