From 8f8239a16590c11428ebce6695830807a17bc910 Mon Sep 17 00:00:00 2001 From: mariaAVC Date: Thu, 16 Oct 2025 14:06:21 +0200 Subject: [PATCH 1/4] Some edits to ebame10_regression_pt1 and half of ebame10_regression_pt2 --- vignettes/ebame10_regression_pt1.Rmd | 6 ++- vignettes/ebame10_regression_pt2.Rmd | 65 +++++++++++++++++++++------- 2 files changed, 54 insertions(+), 17 deletions(-) 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..f822709 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,32 @@ 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 set the working directory to ...). 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")) ``` -# Asking MV to clean up the mess that is below this :) +Note: If you get the message `Error: '/.../...layer-additional-data.txt' does not exist.`, double-check you replaced the directory address correctly. -load raoBust. You may need to install it if you haven't. +# Asking MV to clean up the mess that is below this :) -``` -remotes::install.packages("statdivlab/raoBust") -library(raoBust) -``` +Our first step is to extract and format our data to something R likes. We will do this in a few steps. -Generate presence-absence data +In the following code we: +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 +74,34 @@ 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,16 +113,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")) - ``` +What do we observe in the output? Explain the output. Why does it make sense that the p-value is large? ``` From 67f90ba1d19c806cd859cc9bbdb010ca1d1e2615 Mon Sep 17 00:00:00 2001 From: mariaAVC Date: Thu, 16 Oct 2025 14:17:41 +0200 Subject: [PATCH 2/4] Rendering both labs I have worked on --- docs/articles/ebame10_regression_pt1.html | 6 +- docs/articles/ebame10_regression_pt2.html | 86 +++++++++++++++++++---- 2 files changed, 75 insertions(+), 17 deletions(-) diff --git a/docs/articles/ebame10_regression_pt1.html b/docs/articles/ebame10_regression_pt1.html index 65d177a..aa16baf 100644 --- a/docs/articles/ebame10_regression_pt1.html +++ b/docs/articles/ebame10_regression_pt1.html @@ -123,7 +123,7 @@

Setup
  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 +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. @@ -172,7 +172,9 @@

    Att I want to copy the column light and call it group. give me a cmd line one-liner thx

    (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!

    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

    “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..bbb14c6 100644 --- a/docs/articles/ebame10_regression_pt2.html +++ b/docs/articles/ebame10_regression_pt2.html @@ -105,25 +105,43 @@

    Setup 
    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:

    +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
     base=$(basename "$file" .db)
     anvi-export-functions -c "$file" -o "genes/${base}.txt"
 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 set the working directory to …). 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.

    -
    remotes::install.packages("statdivlab/raoBust")
-library(raoBust)
    -

    Generate presence-absence data

    +

    Our first step is to extract and format our data to something R +likes. We will do this in a few steps.

    +

    In the following code we: 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 usingdistinct` +(repeated rows would have been deleted)

    what_is_there <- list.files(paste0(dir, "genes"), pattern = "*.txt", full.names = T) %>%
   lapply(function(f) {
     read_tsv(f) %>%
@@ -131,15 +149,39 @@ 
    Asking MV to clean up
   }) %>%
   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 %>%
   mutate(presence = 1) %>%
   dplyr::select(sample,accession,presence) %>%
@@ -147,13 +189,27 @@ 
    Asking MV to clean up
   pivot_wider(values_from = presence, values_fill = 0, names_from = accession) %>%
   pivot_longer(!sample,names_to = "cog", values_to = "presence") %>% 
   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? Explain the output. Why does it +make sense that the p-value is large?

    genes_long %>%
   filter(cog == "COG0592") %>% 
   count(presence,light)
    From ebe82e7d9429d60de8c20f66eb2e0ef44c20aac0 Mon Sep 17 00:00:00 2001 From: mariaAVC Date: Thu, 16 Oct 2025 14:26:03 +0200 Subject: [PATCH 3/4] Fixing a small formatting mistake --- docs/articles/ebame10_regression_pt2.html | 56 ++++++++++++++--------- vignettes/ebame10_regression_pt2.Rmd | 14 +++--- 2 files changed, 42 insertions(+), 28 deletions(-) diff --git a/docs/articles/ebame10_regression_pt2.html b/docs/articles/ebame10_regression_pt2.html index bbb14c6..662d95c 100644 --- a/docs/articles/ebame10_regression_pt2.html +++ b/docs/articles/ebame10_regression_pt2.html @@ -130,18 +130,23 @@

    Asking MV to clean up

    Our first step is to extract and format our data to something R likes. We will do this in a few steps.

    -

    In the following code we: 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 usingdistinct` -(repeated rows would have been deleted)

    +

    In the following code we:

    +
      +
    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. +
    2. Joined all these tables by their rows with +bind_rows().
      3. +
    3. Kept only those whose annotation source equals +COG14_FUNCTION with filter.
      4. +
    4. Remove the not-longer-needed columns “source” and “e-value” with +select.
      5. +
    5. Change the name of the column “function” with “cog_function” with +`rename``.
      6. +
    6. Made sure there are not identical rows using distinct +(repeated rows would have been deleted).
      7. +
    what_is_there <- list.files(paste0(dir, "genes"), pattern = "*.txt", full.names = T) %>%
   lapply(function(f) {
     read_tsv(f) %>%
@@ -166,22 +171,29 @@ 
    Asking MV to clean up
 
    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
+“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. 
+
    2. Removing not-longer-useful columns and keeping only “sample”,
+“accession” and “presence” with select
+
      3. 
+
    3. Remove repeated rows with distinct
+
      4. 
+
    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.
      
+1.
      5. 
+
    5. Return to the long format with pivot_longer, keeping
+all the zeros filled up in the previous step.
      6. 
+
    6. We add covariate information contained in the table “meta”, adding
+the columns “clade” and “light” with inner_join.
      7. 
+
    
 
    genes_long <- what_is_there %>%
   mutate(presence = 1) %>%
   dplyr::select(sample,accession,presence) %>%
diff --git a/vignettes/ebame10_regression_pt2.Rmd b/vignettes/ebame10_regression_pt2.Rmd
index f822709..3ecbade 100644
--- a/vignettes/ebame10_regression_pt2.Rmd
+++ b/vignettes/ebame10_regression_pt2.Rmd
@@ -59,12 +59,13 @@ Note: If you get the message `Error: '/.../...layer-additional-data.txt' does no
 Our first step is to extract and format our data to something R likes. We will do this in a few steps.
 
 In the following code we:
-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)
+
+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) %>%
@@ -95,6 +96,7 @@ 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`

From 1e3827c56f6495ac9f2e8119a5f9d73c38c96152 Mon Sep 17 00:00:00 2001
From: mariaAVC 
Date: Thu, 16 Oct 2025 15:10:58 +0200
Subject: [PATCH 4/4] Last changes to part 2 lab

---
 docs/articles/ebame10_regression_pt1.html |  3 +-
 docs/articles/ebame10_regression_pt2.html | 35 +++++++++++++++--------
 vignettes/ebame10_regression_pt2.Rmd      | 25 +++++++++-------
 3 files changed, 40 insertions(+), 23 deletions(-)

diff --git a/docs/articles/ebame10_regression_pt1.html b/docs/articles/ebame10_regression_pt1.html
index aa16baf..7fba0a2 100644
--- a/docs/articles/ebame10_regression_pt1.html
+++ b/docs/articles/ebame10_regression_pt1.html
@@ -33,7 +33,8 @@
 
  3. 
   
   
    4. 
 
    Ok! Time to open it in R. Open your favourite interactive R editor
-(and set the working directory to …). 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.
    
+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.
@@ -220,8 +221,8 @@ 
    Asking MV to clean up
   filter(cog == "COG0592") %>%
   glm_test(presence ~ light, data = .,
            family = binomial(link = "logit"))
    
-
    What do we observe in the output? 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 %>%
   filter(cog == "COG0592") %>% 
   count(presence,light) 
    
@@ -229,27 +230,37 @@ 
    Asking MV to clean up
 
    cog_fns %>%
   filter(accession == "COG0592") 
    
 
    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") |>
              mutate(cog = .y),
            .id = NULL) %>%
   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. 
+
    2. Plotting the estimated difference in odds, along with the robust
+standard error.
      3. 
+
    3. We can order the estimated coefficients from those indicating higher
+differences overall (some repetitions in our estimates happened. Why
+could this be?)
      4. 
+
    
 
    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)))
    
diff --git a/vignettes/ebame10_regression_pt2.Rmd b/vignettes/ebame10_regression_pt2.Rmd
index 3ecbade..067ae81 100644
--- a/vignettes/ebame10_regression_pt2.Rmd
+++ b/vignettes/ebame10_regression_pt2.Rmd
@@ -40,7 +40,7 @@ done
 ```
 
 
-Ok! Time to open it in R. Open your favourite interactive R editor (and set the working directory to ...). 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.  
+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.  
 
 
 ```
@@ -138,8 +138,7 @@ genes_long %>%
            family = binomial(link = "logit"))
 ```
 
-What do we observe in the output? 
-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 %>%
@@ -154,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") |>
@@ -168,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)))