diff --git a/docs/articles/ebame10_regression_pt1.html b/docs/articles/ebame10_regression_pt1.html
index 65d177a..7fba0a2 100644
--- a/docs/articles/ebame10_regression_pt1.html
+++ b/docs/articles/ebame10_regression_pt1.html
@@ -33,7 +33,8 @@
 <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/ebame10_regression_pt1.html">EBAME10: Functional enrichment in anvi'o</a></li>
+<li><a class="dropdown-item" href="../articles/ebame10_pt3.html">EBAME10: Differential Abundance with anvi'o ouput</a></li>
+    <li><a class="dropdown-item" href="../articles/ebame10_regression_pt1.html">EBAME10: Functional enrichment in anvi'o</a></li>
     <li><a class="dropdown-item" href="../articles/ebame10_regression_pt2.html">EBAME10: Functional enrichment in R</a></li>
     <li><a class="dropdown-item" href="../articles/pnf_extn_spaces.html">PNF: Extension spaces for phylogenetic trees</a></li>
     <li><a class="dropdown-item" href="../articles/pnf_radEmu.html">PNF: radEmu for differential abundance</a></li>
@@ -123,7 +124,7 @@ <h2 id="setup">Setup<a class="anchor" aria-label="anchor" href="#setup"></a>
 <ol style="list-style-type: decimal">
 <li>A list of the genomes that we want to include in the analysis. Here,
 we’re going to use the external genomes file
-<code>external-genomes.txt</code>. (Check it out its format with
+<code>external-genomes.txt</code>. (Check out its format with
 <code>head external-genomes.txt</code>.)</li>
 <li>The name of the file we want to save the results. I’m going to call
 it <code>enrichment-output.txt</code>.</li>
@@ -172,7 +173,9 @@ <h2 id="attempting-to-run-anvi-compute-functional-enrichment-across-genomes">Att
 I want to copy the column light and call it group. give me a cmd line one-liner thx</code></pre>
 <p>(You don’t have to say thanks, but maybe we’ll be less likely to turn
 into paperclips)</p>
-<pre><code>awk 'BEGIN{OFS="\t"} NR==1{$(NF+1)="group"} NR&gt;1{$(NF+1)=$3}1' layer-additional-data.txt &gt; layer-additional-data-with-group.txt\n</code></pre>
+<pre><code>awk 'BEGIN{OFS="\t"} NR==1{$(NF+1)="group"} NR&gt;1{$(NF+1)=$3}1' layer-additional-data.txt &gt; layer-additional-data-with-group.txt</code></pre>
+<p>It created a new document with the additional column (Check with
+<code>head layer-additional-data-with-group.txt</code>)</p>
 <p>Ok! Let’s try again!</p>
 <pre><code>anvi-compute-functional-enrichment-across-genomes -e external-genomes.txt -o enrichment-output.txt --annotation-source COG14_FUNCTION -G layer-additional-data-with-group.txt</code></pre>
 <p>“Config Error: The database at ‘Prochlorococcus_31_genomes/AS9601.db’
diff --git a/docs/articles/ebame10_regression_pt2.html b/docs/articles/ebame10_regression_pt2.html
index d4393fc..de9414a 100644
--- a/docs/articles/ebame10_regression_pt2.html
+++ b/docs/articles/ebame10_regression_pt2.html
@@ -33,7 +33,8 @@
 <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/ebame10_regression_pt1.html">EBAME10: Functional enrichment in anvi'o</a></li>
+<li><a class="dropdown-item" href="../articles/ebame10_pt3.html">EBAME10: Differential Abundance with anvi'o ouput</a></li>
+    <li><a class="dropdown-item" href="../articles/ebame10_regression_pt1.html">EBAME10: Functional enrichment in anvi'o</a></li>
     <li><a class="dropdown-item" href="../articles/ebame10_regression_pt2.html">EBAME10: Functional enrichment in R</a></li>
     <li><a class="dropdown-item" href="../articles/pnf_extn_spaces.html">PNF: Extension spaces for phylogenetic trees</a></li>
     <li><a class="dropdown-item" href="../articles/pnf_radEmu.html">PNF: radEmu for differential abundance</a></li>
@@ -105,25 +106,48 @@ <h2 id="setup">Setup<a class="anchor" aria-label="anchor" href="#setup"></a>
 <pre><code>mkdir genes </code></pre>
 <p>If you’ve setup your annotation databases correctly, you can pull out
 the gene annotations and save them in the folder <code>genes</code> by
-copying the following loop into your terminal:</p>
+copying the following loop into your terminal (may take a few minutes,
+so continue with setting R in the mean time):</p>
 <pre><code>for file in *.db; do
     base=$(basename "$file" .db)
     anvi-export-functions -c "$file" -o "genes/${base}.txt"
 done</code></pre>
-<p>Ok! Time to open it in R. Open your favourite interactive R editor.
-Create a variable with the file path where the genomes DBs are. Load in
-the high/low light data.</p>
+<p>Ok! Time to open it in R. Open your favourite interactive R editor
+and create a variable with the file path where the genomes DBs are (if
+you already forgot, you can find it with <code>echo "$PWD"</code>). Load
+in the high/low light data.</p>
 <pre><code><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://tidyverse.tidyverse.org" class="external-link">tidyverse</a></span><span class="op">)</span></span>
+<span></span>
+<span><span class="co">#Change the following dir with what you obtain with echo above.</span></span>
+<span></span>
 <span><span class="va">dir</span> <span class="op">&lt;-</span> <span class="st">"presentations/2510-ebame/Prochlorococcus_31_genomes/"</span></span>
 <span><span class="va">meta</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://readr.tidyverse.org/reference/read_delim.html" class="external-link">read_tsv</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/paste.html" class="external-link">paste0</a></span><span class="op">(</span><span class="va">dir</span>, <span class="st">"layer-additional-data.txt"</span><span class="op">)</span><span class="op">)</span></span></code></pre>
+<p>Note: If you get the message
+<code>Error: '/.../...layer-additional-data.txt' does not exist.</code>,
+double-check you replaced the directory address correctly.</p>
 </div>
 <div class="section level2">
 <h2 id="asking-mv-to-clean-up-the-mess-that-is-below-this">Asking MV to clean up the mess that is below this :)<a class="anchor" aria-label="anchor" href="#asking-mv-to-clean-up-the-mess-that-is-below-this"></a>
 </h2>
-<p>load raoBust. You may need to install it if you haven’t.</p>
-<pre><code><span><span class="fu">remotes</span><span class="fu">::</span><span class="fu">install.packages</span><span class="op">(</span><span class="st">"statdivlab/raoBust"</span><span class="op">)</span></span>
-<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">raoBust</span><span class="op">)</span></span></code></pre>
-<p>Generate presence-absence data</p>
+<p>Our first step is to extract and format our data to something R
+likes. We will do this in a few steps.</p>
+<p>In the following code we:</p>
+<ol style="list-style-type: decimal">
+<li>List all files in the genes folder (extracted from anvi’o above)
+with <code>list.files</code>, and for each file, we read as a table
+(with <code>read_tsv</code>) and added a column named “sample” with the
+name of the source file (with <code>mutate</code>).</li>
+<li>Joined all these tables by their rows with
+<code><a href="https://dplyr.tidyverse.org/reference/bind_rows.html" class="external-link">bind_rows()</a></code>.</li>
+<li>Kept only those whose annotation source equals
+<code>COG14_FUNCTION</code> with <code>filter</code>.</li>
+<li>Remove the not-longer-needed columns “source” and “e-value” with
+<code>select</code>.</li>
+<li>Change the name of the column “function” with “cog_function” with
+`rename``.</li>
+<li>Made sure there are not identical rows using <code>distinct</code>
+(repeated rows would have been deleted).</li>
+</ol>
 <pre><code><span><span class="va">what_is_there</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/list.files.html" class="external-link">list.files</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/paste.html" class="external-link">paste0</a></span><span class="op">(</span><span class="va">dir</span>, <span class="st">"genes"</span><span class="op">)</span>, pattern <span class="op">=</span> <span class="st">"*.txt"</span>, full.names <span class="op">=</span> <span class="cn">T</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://rdrr.io/r/base/lapply.html" class="external-link">lapply</a></span><span class="op">(</span><span class="kw">function</span><span class="op">(</span><span class="va">f</span><span class="op">)</span> <span class="op">{</span></span>
 <span>    <span class="fu"><a href="https://readr.tidyverse.org/reference/read_delim.html" class="external-link">read_tsv</a></span><span class="op">(</span><span class="va">f</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
@@ -131,15 +155,46 @@ <h2 id="asking-mv-to-clean-up-the-mess-that-is-below-this">Asking MV to clean up
 <span>  <span class="op">}</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/bind_rows.html" class="external-link">bind_rows</a></span><span class="op">(</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">source</span> <span class="op">==</span> <span class="st">"COG14_FUNCTION"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
-<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="op">-</span><span class="va">source</span>, <span class="op">-</span><span class="va">`e_value`</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
+<span>  <span class="fu">dplyr</span><span class="fu">::</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="op">-</span><span class="va">source</span>, <span class="op">-</span><span class="va">`e_value`</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/rename.html" class="external-link">rename</a></span><span class="op">(</span>cog_function <span class="op">=</span> <span class="va">`function`</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
 <span>  <span class="va">distinct</span></span></code></pre>
-<p>Create a look-up table to connect accessions to functions.</p>
+<p>“what_is_there” is now a table listing all genes found in each of the
+sample genomes, along with their accession code and the functions
+annotated to them by COG14_function.</p>
+<p>Now, we will create a look-up table to connect accessions to
+functions.</p>
 <pre><code><span><span class="va">cog_fns</span> <span class="op">&lt;-</span> <span class="va">what_is_there</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
-<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="va">accession</span>, <span class="va">cog_function</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
+<span>  <span class="fu">dplyr</span><span class="fu">::</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="va">accession</span>, <span class="va">cog_function</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
 <span>  <span class="va">distinct</span></span>
 <span><span class="va">cog_fns</span></span></code></pre>
-<p>Now, make the zero rows:</p>
+<p>How many distinct accession-function combinations we obtained?
+<code>dim(cog_fns)</code></p>
+<p>Now, we finish by constructing the table of gene presence/absence
+with a useful format for the regression model we will fit later.</p>
+<p>In the following code we construct our table “genes_long” from
+“what_is_there” by:</p>
+<ol style="list-style-type: decimal">
+<li>Add a column called “presence” full of 1’s with <code>mutate</code>
+(this eventually will indicate the gene in that row is present in the
+sample on the same row)</li>
+<li>Removing not-longer-useful columns and keeping only “sample”,
+“accession” and “presence” with <code>select</code>
+</li>
+<li>Remove repeated rows with <code>distinct</code>
+</li>
+<li>Change the format of the table from long to wide, using “sample” as
+the name of the rows, “accession” as the name of the columns and filling
+up the table with “presence” with <code>pivot_longer</code>. Note that
+<code>values_fill = 0</code> indicates that whenever a combination of
+“sample”-“accession” was not a part of the table (ie we are missing the
+value 1 in presence), to fill it up with 0. As a result, missing genes
+in samples will have presence 0, while present genes will have presence
+1.</li>
+<li>Return to the long format with <code>pivot_longer</code>, keeping
+all the zeros filled up in the previous step.</li>
+<li>We add covariate information contained in the table “meta”, adding
+the columns “clade” and “light” with <code>inner_join</code>.</li>
+</ol>
 <pre><code><span><span class="va">genes_long</span> <span class="op">&lt;-</span> <span class="va">what_is_there</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span>presence <span class="op">=</span> <span class="fl">1</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu">dplyr</span><span class="fu">::</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="va">sample</span>,<span class="va">accession</span>,<span class="va">presence</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
@@ -147,13 +202,27 @@ <h2 id="asking-mv-to-clean-up-the-mess-that-is-below-this">Asking MV to clean up
 <span>  <span class="fu"><a href="https://tidyr.tidyverse.org/reference/pivot_wider.html" class="external-link">pivot_wider</a></span><span class="op">(</span>values_from <span class="op">=</span> <span class="va">presence</span>, values_fill <span class="op">=</span> <span class="fl">0</span>, names_from <span class="op">=</span> <span class="va">accession</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://tidyr.tidyverse.org/reference/pivot_longer.html" class="external-link">pivot_longer</a></span><span class="op">(</span><span class="op">!</span><span class="va">sample</span>,names_to <span class="op">=</span> <span class="st">"cog"</span>, values_to <span class="op">=</span> <span class="st">"presence"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate-joins.html" class="external-link">inner_join</a></span><span class="op">(</span><span class="va">meta</span>, by <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">"sample"</span> <span class="op">=</span> <span class="st">"isolate"</span><span class="op">)</span><span class="op">)</span> </span></code></pre>
-<p>Run the model once:</p>
+<p>“genes_long” is now a table that contains a list of all the samples
+and all the genes, with a 1 in “presence” if the gene is in the sample,
+and 0 if it is missing, and with information on the sample by having
+“clade” and “light”.</p>
+<p>We are ready for some modeling!!</p>
+<p>We will be employing the library raoBust, which implements
+generalized linear models, along with robust score test. You may need to
+install it if you haven’t.</p>
+<pre><code><span><span class="fu">remotes</span><span class="fu">::</span><span class="fu"><a href="https://remotes.r-lib.org/reference/install_github.html" class="external-link">install_github</a></span><span class="op">(</span><span class="st">"statdivlab/raoBust"</span><span class="op">)</span></span>
+<span></span>
+<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">raoBust</span><span class="op">)</span></span></code></pre>
+<p>Let’s start with a small case. How can we learn if there is any
+difference in the odds of the gene with accession “COG0592” being
+present between samples with high or low light?</p>
+<p>The following code fits a model to estimate that very difference.</p>
 <pre><code><span><span class="va">genes_long</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">cog</span> <span class="op">==</span> <span class="st">"COG0592"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://rdrr.io/pkg/raoBust/man/glm_test.html" class="external-link">glm_test</a></span><span class="op">(</span><span class="va">presence</span> <span class="op">~</span> <span class="va">light</span>, data <span class="op">=</span> <span class="va">.</span>,</span>
 <span>           family <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/family.html" class="external-link">binomial</a></span><span class="op">(</span>link <span class="op">=</span> <span class="st">"logit"</span><span class="op">)</span><span class="op">)</span></span></code></pre>
-<p>Explain the output. Why does it make sense that the p-value is
-large?</p>
+<p>What do we observe in the output? Why does it make sense that the
+p-value is large?</p>
 <pre><code><span><span class="va">genes_long</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">cog</span> <span class="op">==</span> <span class="st">"COG0592"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/count.html" class="external-link">count</a></span><span class="op">(</span><span class="va">presence</span>,<span class="va">light</span><span class="op">)</span> </span></code></pre>
@@ -161,27 +230,37 @@ <h2 id="asking-mv-to-clean-up-the-mess-that-is-below-this">Asking MV to clean up
 <pre><code><span><span class="va">cog_fns</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">accession</span> <span class="op">==</span> <span class="st">"COG0592"</span><span class="op">)</span> </span></code></pre>
 <p>Ok, how do we run this at scale?</p>
+<p>For this we can create a table for each gene (using
+<code>group_split</code> based in the column cog), and running the model
+above for each. This may take a few minutes!</p>
 <pre><code><span><span class="va">coefs</span> <span class="op">&lt;-</span> <span class="va">genes_long</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/group_split.html" class="external-link">group_split</a></span><span class="op">(</span><span class="va">cog</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
-<span>  <span class="va">head</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span> </span>
 <span>  <span class="fu"><a href="https://rlang.r-lib.org/reference/set_names.html" class="external-link">set_names</a></span><span class="op">(</span><span class="fu"><a href="https://purrr.tidyverse.org/reference/map.html" class="external-link">map_chr</a></span><span class="op">(</span><span class="va">.</span>, <span class="op">~</span> <span class="fu"><a href="https://rdrr.io/r/base/unique.html" class="external-link">unique</a></span><span class="op">(</span><span class="va">.x</span><span class="op">$</span><span class="va">cog</span><span class="op">)</span><span class="op">)</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://purrr.tidyverse.org/reference/map.html" class="external-link">map</a></span><span class="op">(</span><span class="op">~</span> <span class="fu"><a href="https://rdrr.io/pkg/raoBust/man/glm_test.html" class="external-link">glm_test</a></span><span class="op">(</span><span class="va">presence</span> <span class="op">~</span> <span class="va">light</span>, data <span class="op">=</span> <span class="va">.x</span>, family <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/stats/family.html" class="external-link">binomial</a></span><span class="op">(</span>link <span class="op">=</span> <span class="st">"logit"</span><span class="op">)</span><span class="op">)</span><span class="op">$</span><span class="va">coef</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://purrr.tidyverse.org/reference/map_dfr.html" class="external-link">imap_dfr</a></span><span class="op">(</span><span class="op">~</span> <span class="fu"><a href="https://tibble.tidyverse.org/reference/as_tibble.html" class="external-link">as_tibble</a></span><span class="op">(</span><span class="va">.x</span>, rownames <span class="op">=</span> <span class="st">"term"</span><span class="op">)</span> <span class="op">|&gt;</span></span>
 <span>             <span class="fu"><a href="https://dplyr.tidyverse.org/reference/mutate.html" class="external-link">mutate</a></span><span class="op">(</span>cog <span class="op">=</span> <span class="va">.y</span><span class="op">)</span>,</span>
 <span>           .id <span class="op">=</span> <span class="cn">NULL</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/relocate.html" class="external-link">relocate</a></span><span class="op">(</span><span class="va">cog</span><span class="op">)</span></span></code></pre>
-<p>Now, look at the results!</p>
+<p>Now, look at the results! We can do this in several ways like:</p>
+<ol style="list-style-type: decimal">
+<li>Directly looking at the estimated numbers!</li>
+<li>Plotting the estimated difference in odds, along with the robust
+standard error.</li>
+<li>We can order the estimated coefficients from those indicating higher
+differences overall (some repetitions in our estimates happened. Why
+could this be?)</li>
+</ol>
 <pre><code><span><span class="va">coefs</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">term</span> <span class="op">==</span> <span class="st">"lightLL"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
-<span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="op">-</span><span class="va">`Non-robust Std Error`</span>, <span class="op">-</span><span class="va">`Non-robust Wald p`</span>, <span class="op">-</span><span class="va">`Robust Wald p`</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
+<span>  <span class="fu">dplyr</span><span class="fu">::</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/select.html" class="external-link">select</a></span><span class="op">(</span><span class="op">-</span><span class="va">`Non-robust Std Error`</span>, <span class="op">-</span><span class="va">`Non-robust Wald p`</span>, <span class="op">-</span><span class="va">`Robust Wald p`</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/arrange.html" class="external-link">arrange</a></span><span class="op">(</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/desc.html" class="external-link">desc</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/MathFun.html" class="external-link">abs</a></span><span class="op">(</span><span class="va">Estimate</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span>
-<span></span>
-<span><span class="va">out</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
+<span>  </span>
+<span><span class="va">coefs</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">term</span> <span class="op">==</span> <span class="st">"lightLL"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/ggplot.html" class="external-link">ggplot</a></span><span class="op">(</span><span class="fu"><a href="https://ggplot2.tidyverse.org/reference/aes.html" class="external-link">aes</a></span><span class="op">(</span>x <span class="op">=</span> <span class="va">Estimate</span>, y <span class="op">=</span> <span class="va">`Robust Std Error`</span><span class="op">)</span><span class="op">)</span> <span class="op">+</span></span>
 <span>  <span class="fu"><a href="https://ggplot2.tidyverse.org/reference/geom_jitter.html" class="external-link">geom_jitter</a></span><span class="op">(</span><span class="op">)</span></span>
 <span></span>
-<span><span class="va">out</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
+<span><span class="va">coefs</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/filter.html" class="external-link">filter</a></span><span class="op">(</span><span class="va">term</span> <span class="op">==</span> <span class="st">"lightLL"</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/count.html" class="external-link">count</a></span><span class="op">(</span><span class="va">Estimate</span><span class="op">)</span> <span class="op"><a href="https://magrittr.tidyverse.org/reference/pipe.html" class="external-link">%&gt;%</a></span></span>
 <span>  <span class="fu"><a href="https://dplyr.tidyverse.org/reference/arrange.html" class="external-link">arrange</a></span><span class="op">(</span><span class="fu"><a href="https://dplyr.tidyverse.org/reference/desc.html" class="external-link">desc</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/MathFun.html" class="external-link">abs</a></span><span class="op">(</span><span class="va">Estimate</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></span></code></pre>
diff --git a/vignettes/ebame10_regression_pt1.Rmd b/vignettes/ebame10_regression_pt1.Rmd
index 53d1c0b..bfc262b 100644
--- a/vignettes/ebame10_regression_pt1.Rmd
+++ b/vignettes/ebame10_regression_pt1.Rmd
@@ -46,7 +46,7 @@ We want to identify which COGs were present more often in one type of genomes. A
 
 We actually don't need to build a pangenome to answer this question*. All we need is
 
-1. A list of the genomes that we want to include in the analysis. Here, we're going to use the external genomes file `external-genomes.txt`. (Check it out its format with `head external-genomes.txt`.)
+1. A list of the genomes that we want to include in the analysis. Here, we're going to use the external genomes file `external-genomes.txt`. (Check out its format with `head external-genomes.txt`.)
 2. The name of the file we want to save the results. I'm going to call it `enrichment-output.txt`. 
 3. The annotation source. I'm going to choose `COG14_FUNCTION`. 
 4. The groups we want to compare across. This info is contained in `layer-additional-data.txt`. 
@@ -101,9 +101,11 @@ I want to copy the column light and call it group. give me a cmd line one-liner
 (You don't have to say thanks, but maybe we'll be less likely to turn into paperclips)
 
 ```
-awk 'BEGIN{OFS="\t"} NR==1{$(NF+1)="group"} NR>1{$(NF+1)=$3}1' layer-additional-data.txt > layer-additional-data-with-group.txt\n
+awk 'BEGIN{OFS="\t"} NR==1{$(NF+1)="group"} NR>1{$(NF+1)=$3}1' layer-additional-data.txt > layer-additional-data-with-group.txt
 ```
 
+It created a new document with the additional column (Check with `head layer-additional-data-with-group.txt`)
+
 Ok! Let's try again!
 
 ```
diff --git a/vignettes/ebame10_regression_pt2.Rmd b/vignettes/ebame10_regression_pt2.Rmd
index 3756b8c..067ae81 100644
--- a/vignettes/ebame10_regression_pt2.Rmd
+++ b/vignettes/ebame10_regression_pt2.Rmd
@@ -30,7 +30,7 @@ Let's make a new directory to keep this info in.
 mkdir genes 
 ```
 
-If you've setup your annotation databases correctly, you can pull out the gene annotations and save them in the folder `genes` by copying the following loop into your terminal:
+If you've setup your annotation databases correctly, you can pull out the gene annotations and save them in the folder `genes` by copying the following loop into your terminal (may take a few minutes, so continue with setting R in the mean time):
 
 ```
 for file in *.db; do
@@ -39,24 +39,33 @@ for file in *.db; do
 done
 ```
 
-Ok! Time to open it in R. Open your favourite interactive R editor. Create a variable with the file path where the genomes DBs are. Load in the high/low light data.  
+
+Ok! Time to open it in R. Open your favourite interactive R editor and create a variable with the file path where the genomes DBs are (if you already forgot, you can find it with `echo "$PWD"`). Load in the high/low light data.  
+
 
 ```
 library(tidyverse)
+
+#Change the following dir with what you obtain with echo above.
+
 dir <- "presentations/2510-ebame/Prochlorococcus_31_genomes/"
 meta <- read_tsv(paste0(dir, "layer-additional-data.txt"))
 ```
 
+Note: If you get the message `Error: '/.../...layer-additional-data.txt' does not exist.`, double-check you replaced the directory address correctly.
+
 # Asking MV to clean up the mess that is below this :)
 
-load raoBust. You may need to install it if you haven't. 
+Our first step is to extract and format our data to something R likes. We will do this in a few steps.
 
-```
-remotes::install.packages("statdivlab/raoBust")
-library(raoBust)
-```
+In the following code we:
 
-Generate presence-absence data
+1. List all files in the genes folder (extracted from anvi'o above) with `list.files`, and for each file, we read as a table (with `read_tsv`) and added a column named "sample" with the name of the source file (with `mutate`).
+2. Joined all these tables by their rows with `bind_rows()`.
+3. Kept only those whose annotation source equals `COG14_FUNCTION` with `filter`.
+4. Remove the not-longer-needed columns "source" and "e-value" with `select`.
+5. Change the name of the column "function" with "cog_function" with `rename``.
+6. Made sure there are not identical rows using `distinct` (repeated rows would have been deleted).
 
 ```
 what_is_there <- list.files(paste0(dir, "genes"), pattern = "*.txt", full.names = T) %>%
@@ -66,21 +75,35 @@ what_is_there <- list.files(paste0(dir, "genes"), pattern = "*.txt", full.names
   }) %>%
   bind_rows() %>%
   filter(source == "COG14_FUNCTION") %>% 
-  select(-source, -`e_value`) %>% 
+  dplyr::select(-source, -`e_value`) %>% 
   rename(cog_function = `function`) %>% 
   distinct
 ```
 
-Create a look-up table to connect accessions to functions.
+"what_is_there" is now a table listing all genes found in each of the sample genomes, along with their accession code and the functions annotated to them by COG14_function. 
+
+Now, we will create a look-up table to connect accessions to functions.
 
 ```
 cog_fns <- what_is_there %>% 
-  select(accession, cog_function) %>% 
+  dplyr::select(accession, cog_function) %>% 
   distinct
 cog_fns
 ```
 
-Now, make the zero rows:
+How many distinct accession-function combinations we obtained? `dim(cog_fns)`
+
+Now, we finish by constructing the table of gene presence/absence with a useful format for the regression model we will fit later. 
+
+In the following code we construct our table "genes_long" from "what_is_there" by:
+
+1. Add a column called "presence" full of 1's with `mutate` (this eventually will indicate the gene in that row is present in the sample on the same row)
+2. Removing not-longer-useful columns and keeping only "sample", "accession" and "presence" with `select`
+3. Remove repeated rows with `distinct`
+4. Change the format of the table from long to wide, using "sample" as the name of the rows, "accession" as the name of the columns and filling up the table with "presence" with `pivot_longer`. Note that `values_fill = 0` indicates that whenever a combination of "sample"-"accession" was not a part of the table (ie we are missing the value 1 in presence), to fill it up with 0. As a result, missing genes in samples will have presence 0, while present genes will have presence 1. 
+5. Return to the long format with `pivot_longer`, keeping all the zeros filled up in the previous step. 
+6. We add covariate information contained in the table "meta", adding the columns "clade" and "light" with `inner_join`.
+
 
 ```
 genes_long <- what_is_there %>%
@@ -92,17 +115,30 @@ genes_long <- what_is_there %>%
   inner_join(meta, by = c("sample" = "isolate")) 
 ```
 
-Run the model once:
+"genes_long" is now a table that contains a list of all the samples and all the genes, with a 1 in "presence" if the gene is in the sample, and 0 if it is missing, and with information on the sample by having "clade" and "light". 
+
+We are ready for some modeling!! 
+
+We will be employing the library raoBust, which implements generalized linear models, along with robust score test. You may need to install it if you haven't. 
+
+```
+remotes::install_github("statdivlab/raoBust")
+
+library(raoBust)
+```
+
+Let's start with a small case. How can we learn if there is any difference in the odds of the gene with accession "COG0592" being present between samples with high or low light? 
+
+The following code fits a model to estimate that very difference.
 
 ```
 genes_long %>%
   filter(cog == "COG0592") %>%
   glm_test(presence ~ light, data = .,
            family = binomial(link = "logit"))
-
 ```
 
-Explain the output. Why does it make sense that the p-value is large?
+What do we observe in the output? Why does it make sense that the p-value is large?
 
 ```
 genes_long %>%
@@ -117,12 +153,14 @@ cog_fns %>%
   filter(accession == "COG0592") 
 ```
 
-Ok, how do we run this at scale?
+Ok, how do we run this at scale? 
+
+For this we can create a table for each gene (using `group_split` based in the column cog),
+and running the model above for each. This may take a few minutes! 
 
 ```
 coefs <- genes_long %>%
   group_split(cog) %>%
-  head %>% 
   set_names(map_chr(., ~ unique(.x$cog))) %>%
   map(~ glm_test(presence ~ light, data = .x, family = binomial(link = "logit"))$coef) %>%
   imap_dfr(~ as_tibble(.x, rownames = "term") |>
@@ -131,20 +169,24 @@ coefs <- genes_long %>%
   relocate(cog)
 ```
 
-Now, look at the results!
+Now, look at the results! We can do this in several ways like: 
+
+1. Directly looking at the estimated numbers! 
+2. Plotting the estimated difference in odds, along with the robust standard error.
+3. We can order the estimated coefficients from those indicating higher differences overall (some repetitions in our estimates happened. Why could this be?)
 
 ```
 coefs %>%
   filter(term == "lightLL") %>%
-  select(-`Non-robust Std Error`, -`Non-robust Wald p`, -`Robust Wald p`) %>%
+  dplyr::select(-`Non-robust Std Error`, -`Non-robust Wald p`, -`Robust Wald p`) %>%
   arrange(desc(abs(Estimate)))
-
-out %>%
+  
+coefs %>%
   filter(term == "lightLL") %>%
   ggplot(aes(x = Estimate, y = `Robust Std Error`)) +
   geom_jitter()
 
-out %>%
+coefs %>%
   filter(term == "lightLL") %>%
   count(Estimate) %>%
   arrange(desc(abs(Estimate)))