In this document, we will go through a few steps used in my study of the potato microbiota as explained during the BioData Basel meetup on April, 9th 2026. A specificity of this project is that we isolated bacterial strains from the same plants we took microbiota samples, increasing the chance that the observed Amplicon Sequence Variant (ASV) correspond to strains. This gave us possibility of a testing hypotheses on the role of some strains based on the changes observed in the plant’s microbiota.

Mixing metagenomics and culturomics is challenging and sometimes messy. Unfortunately the “messy” nature of the task is rarely reported and deserves more attention. The goal of this presentation is to show you some concrete examples of the challenges we meet with the current technology.

For clarity reasons, I show only here a small fraction of the work we did. If you want the whole story, the preprint and the link to all the data and codes are available on BioRxiv

In the first part I show you how we detected differentially abundant strains between infected and uninfected plants using the package edgeR, a package normally used to detect differentially expressed genes in RNAseq studies.

In the second part I show you an example of how to represent the differential abundance in the whole rhizosphere. I will also show you the limits of the aggregation per Phylum/Order/Genus that we see frequently in microbiome papers.

In the third part we see how to create a custom BLAST database and query it with the goal of detecting the ASV corresponding to each bacterial strain we isolated.

This whole document is supposed to be read along with the attached pptx file.

1 - Differential abundance strains detection

I show you how we detected differentially abundant strains between infected and uninfected plants using the package edgeR, a package normally used to detect differentially expressed genes in RNAseq studies.

Import data

## Import ASV table
asvtab_filtered <- read.csv( "data/asvmat_br.csv",
                             check.names = FALSE,
                             row.names = 1)

## Turn it into a matrix
asvmat_BR <- as.matrix(asvtab_filtered)

asvtab_filtered |> 
  head() |> 
  kable() |> 
  kable_styling() |> 
  scroll_box(width = "750px")
2_BR15_1C3 2_BR15_1C4 2_BR15_1C5 2_BR15_1V1 2_BR15_1V3 2_BR15_1V5 2_BR15_2C2 2_BR15_2C4 2_BR15_2C5 2_BR15_2V1 2_BR15_2V4 2_BR15_2V5
asv000003 189 584 181 397 194 472 709 312 348 475 817 428
asv000004 114 85 51 0 0 226 80 115 237 6 279 12
asv000005 12 23 23 0 0 40 5 21 40 3 23 5
asv000006 199 662 175 371 170 454 627 277 312 458 794 416
asv000007 114 100 60 0 0 183 89 117 265 8 214 15
asv000008 24 32 12 0 0 27 4 15 31 0 27 10 
## Import phylogeny tab
phylo_filtered <- read.csv( "data/phylogeny_br.csv",
                            check.names = FALSE,
                            row.names = 1)

## Turn it into a matrix
phylo_BR <- as.matrix(phylo_filtered)

phylo_filtered |> 
  head() |> 
  kable() |> 
  kable_styling()
Kingdom Phylum Class Order Family Genus
asv000003 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingobium
asv000004 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae Streptomyces
asv000005 Bacteria Actinobacteriota Actinobacteria Streptosporangiales Streptosporangiaceae Nonomuraea
asv000006 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingobium
asv000007 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae Streptomyces
asv000008 Bacteria Actinobacteriota Actinobacteria Streptosporangiales Streptosporangiaceae Nonomuraea 
## Import sample metadata
meta_BR <- read.csv("data/metadata_br.csv",
                     check.names = FALSE,
                     row.names = 1)

meta_BR |> 
  head() |> 
  kable() |> 
  kable_styling() |> 
  scroll_box(width = "750px")
gen cultiv comp time treat_batch treat rep plant sample_harvest dna_extract_batch tube_nb dna_extract_date qubit1 qubit2 qubit_mean pcr_batch
2_BR15_1C3 2 B R 15 1 C 3 BC06 13/04/2020 1 51 16/12/2020 8.00 8.48 8.24 1
2_BR15_1C4 2 B R 15 1 C 4 BC07 13/04/2020 2 52 23/12/2020 6.04 6.12 6.08 2
2_BR15_1C5 2 B R 15 1 C 5 BC08 13/04/2020 3 53 28/12/2020 9.88 10.10 9.99 3
2_BR15_1V1 2 B R 15 1 V 1 BV02 13/04/2020 1 54 16/12/2020 9.24 9.40 9.32 1
2_BR15_1V3 2 B R 15 1 V 3 BV22 13/04/2020 2 55 23/12/2020 8.32 8.44 8.38 2
2_BR15_1V5 2 B R 15 1 V 5 BV30 13/04/2020 4 56 29/12/2020 9.00 9.20 9.10 4 
## Import reference sequences
sequences <- read.csv("data/ref_sequences_br.csv",
                     check.names = FALSE)

Create DGE object

Here we combine the matrices into a DGE object to be analysed by edgeR.

## Create DGEList objects
dge_br <- DGEList(counts = asvmat_BR,
                 genes = phylo_BR,
                 samples = meta_BR)

## Convert type of the "treatment" column to factor
dge_br$samples$treat <- factor(dge_br$samples$treat,
                             levels = c("C", "V"))

## TMM normalization (Trimmed Mean of the M-values)
dge_br <- edgeR::calcNormFactors(dge_br)

Find differentially abundant taxa

## Create Multidimensional Scaling plot
plotMDS(dge_br
        , labels = rownames(dge_br$samples)
        , col = as.numeric(dge_br$samples$treat)
        , cex = 0.5)

The color represents the infection status (red is infected).

Most of the variability is linked to the infection status.

## Design the model's formula for the Generalized Linear Model
design <- model.matrix(~ treat , data = dge_br$samples )

## Estimate the variability of counts of each ASV (necessary for negative binomial GLM)
glmest <- estimateGLMRobustDisp(dge_br, design = design)

## Fits negative binomial GLM. prior.count adds a pseudocount to stabilize estimates for low-occurring ASVs.
fit_glm <- glmFit(glmest, design = design, prior.count = 0.125)

## GLM Likelihood Ratio test, checks if the infection status has an effect on the counts of each ASV.
## Calculates AND ADJUSTS p-value.
lrt_treat <- glmLRT(fit_glm, coef = 2)

## Sorts ASV per significance
tt_treat <- as.data.frame(topTags(lrt_treat, n = Inf ))

## Show tt_treat
tt_treat |> 
  head() |> 
  kable() |> 
  kable_styling() |> 
  scroll_box(width = "750px")
Kingdom Phylum Class Order Family Genus logFC logCPM LR PValue FDR
asv000561 Bacteria Bacteroidota Bacteroidia Chitinophagales Saprospiraceae Phaeodactylibacter -9.575301 10.85905 135.92952 0 0
asv000547 Bacteria Bacteroidota Bacteroidia Chitinophagales Saprospiraceae Phaeodactylibacter -8.939168 10.87846 113.19125 0 0
asv001179 Bacteria Bacteroidota Bacteroidia Flavobacteriales Flavobacteriaceae Flavobacterium -9.623207 10.90245 112.89278 0 0
asv001158 Bacteria Bacteroidota Bacteroidia Flavobacteriales Flavobacteriaceae Flavobacterium -9.498248 10.78260 111.82150 0 0
asv001588 Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales NA NA -8.260053 9.61971 85.29664 0 0
asv000839 Bacteria Bacteroidota Bacteroidia Chitinophagales Saprospiraceae NA 9.190066 10.50159 69.46067 0 0

This table shows how differentially abundant each ASV is between the infected and uninfected plants. They are ordered accoding to the p-value.

## Export differentially expressed ASVs counts 
write.csv(tt_treat, file = "data/all_ASV.csv")

We can then plot the results as a volcano plot.

## Summarize differential abundance (used for the volcano plot below)
treat_summary <- t(
    summary( de <- decideTests(lrt_treat)) )

## Create volcano plot
plotMD( lrt_treat, main = "Differential abundance in the rhizosphere", cex=0.5 )
abline( h=c(-1,1), col = "blue" )
text( x = 12, y = -4, labels = paste0("Up: ",treat_summary[3], 
                                      "\n NotSign: ", treat_summary[2], 
                                      "\n Down: ", treat_summary[1]) )

In the plot above, each point represents one ASV. The x axis represents the average abundance (log scale), and the y axis represents the log-fold change between the infected and uninfected plants. The red dots represent the ASV that became more abundance upon infection, the blue dots are the ASV that became less abundant, the black dots did not change significantly.

2 - Community shifts visualization

Here I show you an example of how to represent the differential abundance in the whole rhizosphere.

Some more data wrangling is necessary between DA ASV detection and this step:

  • Creation of a giant phylogenetic tree with all the ASVs (in Unix on HPC, aligned with MAFFT, tree created with FastTree)
  • Aggregate the taxa at the Order phylogenetic level, and calculate 3 indicators:
    • the total abundance of each Order.
    • Proportion of differentially abundant ASVs, weighted by ASV abundance, scaled with the order’s abundance.
    • Average log-fold change of differentially abundant ASVs, weighted by ASV abundance, scaled with the total abundance.

we will skip to focus on tree building and comments.

Import and filter objects

## phylogenetic trees
tree <- readRDS("data/tree_objects/ASV_phylo.rds")

This is a phylogenetic tree, aggregated at the order level.

## Phylogenetic information
tree_info <- readRDS("data/tree_objects/tree_info.rds") |> 
  select(-Family, -Genus)

tree_info |> 
  head() |> 
  kable() |> 
  kable_styling(full_width = FALSE)
Kingdom Phylum Class Order
Micropepsales Bacteria Proteobacteria Alphaproteobacteria Micropepsales
Rhizobiales Bacteria Proteobacteria Alphaproteobacteria Rhizobiales
Parvibaculales Bacteria Proteobacteria Alphaproteobacteria Parvibaculales
Rhodospirillales Bacteria Proteobacteria Alphaproteobacteria Rhodospirillales
Ferrovibrionales Bacteria Proteobacteria Alphaproteobacteria Ferrovibrionales
Dongiales Bacteria Proteobacteria Alphaproteobacteria Dongiales 
## List of orders belonging to each phylum (used for colors)
taxalist <- readRDS("data/tree_objects/taxa_list.rds")

map(taxalist, head, 3)[1:5]
## $Proteobacteria
## [1] "Micropepsales"  "Rhizobiales"    "Parvibaculales"
## 
## $Bdellovibrionota
## [1] "Bdellovibrionales" "na_Bdellovibrion." "Bacteriovoracales"
## 
## $Dependentiae
## [1] "Babeliales"
## 
## $Myxococcota
## [1] "Blfdi19"      "Myxococcales" "mle1-27"     
## 
## $Acidobacteriota
## [1] "Subgroup 7"     "PAUC26f"        "Bryobacterales"

This is a list of the Orders included in each Phylum.

## Abundance matrix (used for annotations)
abund_matrix <- readRDS("data/tree_objects/abund_matrix.rds")

abund_matrix  |> 
  head() |> 
  kable() |> 
  kable_styling(full_width = FALSE)
B I
Micropepsales 369782.747 438629.875
Rhizobiales 2671773.800 2441710.203
Parvibaculales 1978.973 1320.667
Rhodospirillales 112384.401 130058.304
Ferrovibrionales 8435.939 3225.957
Dongiales 52407.014 49983.107

This is the matrix of abundance.

## Local indicator (described above,used for more annotations)
DA_proportion_matrix <- readRDS("data/tree_objects/abund_ratio_matrix.rds")

DA_proportion_matrix  |> 
  head() |> 
  kable() |> 
  kable_styling(full_width = FALSE)
B I
Micropepsales 0.1506374 0.0126441
Rhizobiales 0.1733440 0.0242824
Parvibaculales NA NA
Rhodospirillales 0.0542960 0.1678773
Ferrovibrionales NA NA
Dongiales NA NA

This is the matrix showing the ratio of differential abundance present in each Order.

## Global indicator (described above, used for even more annotations)
logfc_direction_matrix <- readRDS("data/tree_objects/logfc_matrix.rds")

logfc_direction_matrix  |> 
  head() |> 
  kable() |> 
  kable_styling(full_width = FALSE)
B I
Rhizobiales 96703.519 -11361.976
Micropepsales -7315.006 -1324.272
na_Proteobacteria 2625.931 -2210.329
Parvibaculales NA NA
Caulobacterales 45797.812 -76891.229
Rhodospirillales -2419.906 4781.282

This is the matrix showing the average log-fold change of the significantly differentially abundant ASVs.

Create base tree

## Make basic tree
p <- ggtree(tree, 
            ladderize = TRUE, 
            layout="fan",
            open.angle = 2.5) %<+% tree_info +
  geom_tiplab(size = 2.5, 
              align = TRUE, 
              linesize = 0.5,
              linetype = "dotted",
              show.legend = FALSE)+
  labs(title = "Generation 2, rhizosphere")+
  theme(legend.title = element_text(size = 16),
        legend.text = element_text(size = 14))

## rotate tree and add colors
q <- p  |>  rotate(230) |> rotate(225)|> rotate(228)|> rotate(224)|> rotate(215)|> rotate(217)|> rotate(219)|> rotate(205)|> rotate(201) |> rotate(213)|> rotate(226)|> rotate(214) |> rotate(194) |> rotate(189) |> rotate(190) |> rotate(163) |> rotate(164) |>
  groupOTU(taxalist, "Phylum") +
  aes(color = Phylum) +
  scale_color_manual(values = c("deeppink2", "orangered", "orange", "olivedrab3", "limegreen","springgreen4", "lightseagreen", "steelblue2", "royalblue1", "grey60","purple1", "orchid2", "hotpink1"))

Add informative matrices

r <- gheatmap(q,
         abund_matrix,
         offset=0.88,
         width =0.15,
         color = "white",
         colnames_angle = 45,
         font.size = 2.5,
         colnames_offset_y = 0,
         colnames_offset_x = 0
         ) +
  scale_fill_gradientn(
    colors=c("#FFFFFF", "#FDE725FF", 
             "#94D840FF", "#3CBC75FF", "#1F968BFF", "#2D718EFF", 
             "#3F4788FF", "#440154FF"),
    limits = c(0, 6e6),
                      name = "TMM-normalized \nAbundance",
                      na.value = "grey"
                      ) +
  new_scale_fill()

## Add density heatmap
s <- gheatmap(r,
              DA_proportion_matrix,
              offset=0.66,
              width = 0.15,
              color = "white",
              colnames_angle = 45,
              font.size = 2.5,
              colnames_offset_y = 0,
              colnames_offset_x = 0
              )+

  scale_fill_gradientn(colors = c("white", "#3F007D"),
                      limits = c(0, 1),
                      name = "Ratio",
                      na.value = "white"
                      ) +
  new_scale_fill()

## Add averaged logfc heatmap
t <- gheatmap(s,
         logfc_direction_matrix,
         offset=0.45,
         width =0.15,
         color = "white",
         colnames_angle = 45,
         font.size = 2.5,
         colnames_offset_y = 0,
         colnames_offset_x = 0
         )+
  scale_fill_gradientn(colors=c("#053061", "#2166AC", "#FFFFFF", "#B2182B", "#530019"),
                      values= scales::rescale(c(-3.5e5, -1e5, 0, 1e5, 2.95e5)),
                      limits = c(-3.5e5, 2.95e5),
                      name = "Weighted Averaged \nLog-fold Change",
                      na.value = "grey78"
                      )

t

To read the tree properly, please open the generated pdf file below.

This tree and the information present in it is complicated, but it is the best I can do to represent such a complex dataset. For us, this tree has been a great source of information and it takes some time to get into it and understand each information. Once it’s done though, a glance is enough to compare our results with results presented in many articles.

One of the reasons the ggtree package is great is because of the fantastic tutorials available on internet. I learned everything from this page.

3 - Blast database creation and query

Database Creation and query

–> See the script blast_db_creation_query.sh

Get list of matching strains

## Import output
alignment <- read.csv("save_me/ASV_isolates_alignment.csv")

## Get the list of strains of interest
strain_list <- alignment |>
  ## Remove short alignments
  filter(length > 400) |>
  ## Pull strain names only
  pull(stitle) |>
  ## remove duplicate strains
  unique()

## Print complete list of strain of interest
cat(sprintf("In the end, we have %d strains of interest to phenotype:\n",
            length(strain_list)))
## In the end, we have 458 strains of interest to phenotype:
## Get list of perfect matches
perfect_match <- alignment |> 
  ## Filter
  filter(length > 400,      # Remove short alignments
         pident == 100      # Keep only perfect matches
         ) |>
  ## Keep only one match for each strain
  distinct(saccver, .keep_all = TRUE) |> 
  ## Select columns
  select(strain = saccver, ASV = qaccver, length)

## Print
perfect_match |> 
  kable(label = "List of perfect matches") |> 
  kable_styling(full_width = FALSE) |> 
  scroll_box(height = "400px")
strain ASV length
BR014A asv000047 414
BS236 asv000047 414
BS222 asv000047 414
BS177 asv000047 414
BS139 asv000047 414
BS200 asv000047 414
BS232 asv000047 414
BR066 asv000047 414
BS199 asv000047 414
BS150 asv000047 414
BL096 asv000047 414
BL074 asv000047 414
AS035 asv000047 414
AS001 asv000047 414
AR039 asv000047 414
AS013 asv000047 414
BR101 asv000049 413
BR010 asv000049 413
BS192 asv000049 413
BS181 asv000049 413
AR150 asv000049 413
BS130 asv000051 415
BS119 asv000051 415
AS029 asv000068 416
BR107 asv000068 416
BL049 asv000068 416
BS196 asv000068 416
AS040 asv000068 416
AS022 asv000068 416
AR057 asv000068 416
AR042 asv000068 416
AR005 asv000068 416
AR002 asv000068 416
AR001 asv000068 416
BS132 asv000070 416
BS180 asv000077 416
AR037 asv000077 416
AL004 asv000077 416
BR134 asv000077 416
AR055 asv000077 416
AR206 asv000077 416
AS085 asv000086 413
BR157 asv000086 413
AS050 asv000086 413
AR097 asv000086 413
BS107 asv000134 416
BS079 asv000216 416
BS179 asv000216 416
BS187 asv000216 416
AR095 asv000216 416
AS146 asv000244 416
BS239 asv000330 416
AL025 asv000330 416
AL010 asv000330 416
AL040 asv000330 416
BS140 asv000330 416
BR077 asv000330 416
AR117 asv000330 416
BR121 asv000372 416
AR021 asv000372 416
BS231 asv000375 409
BS103 asv000375 409
BR094 asv000375 409
BR144 asv000375 409
AS108 asv000395 412
BS088 asv000395 412
AR010 asv000395 412
AS010 asv000395 412
AR051 asv000395 412
AR046 asv000395 412
AS121b asv000395 412
AS106 asv000395 412
BR081 asv000395 412
AL005 asv000438 413
BR133 asv000438 413
AR009 asv000438 413
BR112 asv000438 413
AS056 asv000475 417
AS006 asv000475 417
AR008 asv000475 417
AR068 asv000512 414
AR093 asv000512 414
AS024 asv000668 413
AR076 asv000668 413
BR034 asv000668 413
AR144 asv000668 413
BS038 asv000668 413
BS037 asv000668 413
BS158 asv000668 413
AS182 asv000668 413
AS181 asv000668 413
BS061 asv000668 413
AS026 asv000668 413
AR190 asv000668 413
AR178 asv000668 413
AR176 asv000668 413
AR143 asv000668 413
BR109 asv000670 409
AR071 asv000701 415
AR168 asv000701 415
AR125 asv000701 415
AR088 asv000701 415
AS123 asv000843 416
AS093 asv001155 413
AR201 asv001155 413
AR132 asv001155 413
BS131 asv001219 414
BS118 asv001219 414
AL091 asv001219 414
AL090 asv001219 414
BS146 asv001219 414
BL080 asv001219 414
BS134 asv001219 414
BR088 asv001219 414
BL052 asv001219 414
BL046 asv001219 414
AS096 asv001219 414
AS094 asv001219 414
AS132 asv001219 414
AS144 asv001219 414
AS021 asv001219 414
BL047a asv001872 412
AS133 asv001872 412
BR013 asv001884 409
BS071 asv002046 417
BR099 asv002046 417
BR015 asv002046 417
AS178 asv002046 417
AR179 asv002046 417
AR177 asv002046 417
BS237 asv002208 417
BS027 asv002231 414
BS229 asv002231 414
BR012 asv002231 414
BR046 asv002231 414
AR203b asv002231 414
BS113 asv002388 416
BR136 asv002388 416
AS027 asv002573 413
AS180 asv002573 413
AR032 asv002573 413
AS080 asv002573 413
BS157 asv002615 412
AL028 asv002748 414
AS092 asv002849 416
AS057 asv002849 416
BS172 asv002982 415
AR073 asv003008 415
AR089 asv003008 415
AR024 asv003008 415
AR202 asv003008 415
BS141 asv003045 417
AL039 asv003204 417
AL034 asv003204 417
AR061 asv003204 417
BL047b asv003269 414
BL047c asv003269 414
BR050 asv003333 409
BR049 asv003333 409
BS123 asv003333 409
BR156 asv003333 409
AR118 asv003333 409
AR149 asv003603 416
BS080 asv003647 416
BS138 asv004665 417
AS058 asv005149 417
BL071 asv005149 417
BL067 asv005149 417
BR153 asv005149 417
BR150 asv005149 417
BR126 asv005149 417
AS064 asv005149 417
AS054 asv005149 417
AR083 asv005149 417
AL064 asv005149 417
AS070 asv005149 417
AR116 asv005149 417
AR081 asv005149 417
AR077 asv005149 417
AR079 asv005979 414
AR078 asv005979 414
BR127 asv007669 416 
## Report on perfect matches
cat(sprintf("%d ASVs are perfectly matched by at least a strain, with a total of %d strains.\n",
            length(unique(perfect_match$ASV)),
            nrow(perfect_match)))
## 45 ASVs are perfectly matched by at least a strain, with a total of 182 strains.
## Get list with 1 mismatch (simplified)
mismatch <- alignment |> 
  ## Filter
  filter(length > 400,      # Remove short alignments
         pident != 100,    # Keep only imperfect matches
         !(saccver %in% perfect_match$strain),   # Remove strains that have a perfect match
         ) |>
  ## Keep only one match for each strain
  distinct(saccver, .keep_all = TRUE) |> 
  ## Select columns
  select(strain = saccver, ASV = qaccver, mismatch)

## Print
mismatch |> 
  kable(label = "List of strains with a mismatch") |> 
  kable_styling(full_width = FALSE) |> 
  scroll_box(height = "400px")
strain ASV mismatch
BS227 asv000004 1
BS195a asv000004 1
BS101 asv000004 1
BS217 asv000004 1
BS195B asv000004 1
BS166 asv000004 1
BS089 asv000004 1
BR163 asv000004 1
BR111 asv000004 1
BR104 asv000004 1
BR097 asv000004 1
BR048 asv000004 1
BL087 asv000004 1
AS136 asv000004 1
AS113 asv000004 1
AS111 asv000004 1
AS104 asv000004 1
AS099 asv000004 1
AS076 asv000004 1
AS073 asv000004 1
AS041 asv000004 1
AS037 asv000004 1
AS032 asv000004 1
AS004 asv000004 1
AR083D asv000004 1
AR083C asv000004 1
AR072 asv000004 1
AR069 asv000004 1
AR064 asv000004 1
AR062 asv000004 1
AR036 asv000004 1
AR029 asv000004 1
AR027 asv000004 1
AR018 asv000004 1
AR013 asv000004 1
AR012 asv000004 1
BR070 asv000004 1
AR065 asv000004 1
BS220 asv000004 1
BS185 asv000004 1
BS093c asv000004 1
BS093b asv000004 1
BS086 asv000004 1
BR117 asv000004 1
AS131 asv000004 1
AS114 asv000004 1
AS078 asv000004 1
AS020 asv000004 1
AR109 asv000004 1
AR108 asv000004 1
AR003 asv000004 1
BS197b asv000004 1
BS188 asv000004 1
BS167 asv000004 1
BS112 asv000004 1
BS096 asv000004 1
BS090 asv000004 1
BS084 asv000004 1
BS082 asv000004 1
BS238 asv000004 1
BS226 asv000004 1
BS207 asv000004 1
BS205 asv000004 1
BS160 asv000004 1
BS159 asv000004 1
BS152 asv000004 1
BS149 asv000004 1
BS128 asv000004 1
BR079 asv000004 1
BR068 asv000004 1
BR061 asv000004 1
BR045 asv000004 1
BR044 asv000004 1
BR165 asv000004 1
BR164 asv000004 1
BR159 asv000004 1
BR143 asv000004 1
BR142 asv000004 1
BR132 asv000004 1
BL038 asv000004 1
AS079 asv000004 1
AS053 asv000004 1
AS002 asv000004 1
AS116 asv000004 1
AS115 asv000004 1
AR094 asv000004 1
AR091 asv000004 1
AR082 asv000004 1
AR006 asv000004 1
AR044 asv000004 1
AR033 asv000004 1
AR030 asv000004 1
AR019 asv000004 1
AR113 asv000004 1
AR101 asv000004 1
AL071 asv000004 1
AL045 asv000004 1
BS078 asv000004 1
BR113 asv000004 1
AS105 asv000004 1
AS098a asv000004 1
AS008 asv000004 1
AR124 asv000004 1
AR115 asv000004 1
AR085 asv000004 1
AR043 asv000004 1
AR038b asv000004 1
AL083 asv000004 1
BS211 asv000004 1
AS100 asv000004 1
AR092 asv000004 1
BR085 asv000004 1
BR074 asv000004 1
AS098b asv000004 1
AR050 asv000004 1
BS215a asv000004 1
BS202 asv000004 1
BS155 asv000004 1
BS093a asv000004 1
BR131 asv000004 1
AS121c asv000004 1
AS121a asv000004 1
AS014 asv000004 1
AS007 asv000004 1
AR092b asv000004 1
AR083b asv000004 1
AR015 asv000004 1
AR004 asv000004 1
AL074 asv000004 1
BR122 asv000049 1
AS044 asv000051 1
AR060 asv000051 1
BR138 asv000051 1
AS176 asv000051 1
AS063 asv000051 1
AS062 asv000051 1
BR154 asv000051 1
BR098 asv000068 1
AR025 asv000077 1
AR017 asv000086 1
BS213 asv000149 1
AS112 asv000149 1
BL040B asv000149 1
BL040A asv000149 1
AS118 asv000149 1
AS097 asv000149 1
AR020 asv000149 1
BR092a asv000149 1
BR092b asv000149 1
BS234 asv000149 1
BS145 asv000149 1
BR060 asv000149 1
BR114 asv000149 1
AL067 asv000149 1
BR057 asv000149 1
BR051 asv000149 1
BS214 asv000149 1
AR170 asv000149 1
BS106 asv000216 1
BS105 asv000216 1
BS144 asv000216 1
BR130 asv000216 1
AL076 asv000216 1
BS184 asv000248 1
BS102 asv000285 1
AR070 asv000575 1
AR066 asv000575 1
BS240 asv000668 1
AS071 asv000701 1
BR087 asv000701 1
BR078 asv000701 1
BR075 asv000701 1
BR137 asv000701 1
BR135 asv000701 1
BR093 asv000701 1
AR114 asv000701 1
BR058 asv000701 1
AR040 asv001199 1
AS048 asv001199 1
BS216c asv001199 1
BS216a asv001199 1
BS197a asv001199 1
AS005 asv001199 1
BS133 asv001219 1
AS018 asv001219 1
AR098 asv001277 1
BR065 asv001327 1
BS115 asv001631 1
BS235 asv001977 1
AS086 asv002271 1
AR067 asv002271 1
AL048 asv002271 1
AL043 asv002271 1
AL016 asv002271 1
AL006 asv002271 1
BS163 asv002271 1
BS110 asv002271 1
BR145 asv002271 1
AS011 asv002271 1
AR111 asv002271 1
BS203 asv002271 1
BS189 asv002271 1
BS143 asv002271 1
BR120 asv002271 1
BR103 asv002271 1
BR067 asv002271 1
BR128 asv002271 1
AS047 asv002271 1
AS038 asv002271 1
AS030 asv002271 1
AS019 asv002271 1
AR105 asv002271 1
AR049 asv002271 1
AR048 asv002271 1
BR155 asv002271 1
BR089 asv002271 1
AS095 asv002271 1
AS081 asv002271 1
AS077 asv002271 1
AS075 asv002271 1
AS015 asv002271 1
AS012 asv002271 1
AR103 asv002271 1
AR084 asv002271 1
BR064 asv002346 1
AR016 asv002346 1
BR119 asv002346 1
BR095 asv002346 1
BR129 asv002346 1
BS161 asv002346 1
BS083 asv002346 1
BR106 asv002346 1
BR100 asv002346 1
AR063 asv002448 1
BS148 asv002488 1
AR122 asv003008 1
AS109 asv003126 1
BS233 asv003607 1
BS228 asv003607 1
AS127 asv003607 1
AS107 asv003607 1
BR162 asv003607 1
AS049 asv003625 1
BR071 asv003647 1
BS124 asv003647 1
BS219c asv003681 1
BS219b asv003681 1
BS129b asv003681 1
BS170 asv003790 1
BS092 asv003945 1
AS028 asv004156 1
AR075 asv004400 1
AS061 asv004651 1
AR056 asv004651 1
BR052 asv004651 1
AR135 asv004651 1
BS230 asv004743 1
AS052 asv004743 1
BS135 asv004823 1
BS137b asv005080 1
AS145 asv005080 1
AS140A asv005080 1
BS136 asv005080 1
AS143 asv005080 1
AS142 asv005080 1
AS141 asv005080 1
AS139a asv005080 1
AR034 asv005080 1
BS218 asv005080 1
BS137a asv005080 1
BS140b asv005080 1
AS140b asv005080 1
AS060 asv005242 1
AS087 asv005242 1
AR022 asv006178 1
AL013 asv006178 1 
## Report on perfect matches
cat(sprintf("%d ASVs are matched by at least a strain with a mismatch, with a total of %d strains.\n",
            length(unique(mismatch$ASV)),
            nrow(mismatch)))
## 39 ASVs are matched by at least a strain with a mismatch, with a total of 276 strains.

We did not have the means to screen all the strains of interest, so we chosed based on the following criteria:

  • We pick all the strains that have a perfect 16S/ASV match,
  • For each ASV that does not have a perfect match, we pick 1 strain with 1 mismatch.
  • We did not test the filamentous bacteria because the 16SrRNA gene is a mediocre identity marker for them, and detected very low changes upon infection. They deserve attention though because they are the soil’s champions.

If you want the whole story, the preprint and the link to all the data and codes are available on BioRxiv.

The End

sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Paris
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggnewscale_0.5.2   ggtreeExtra_1.20.1 ggtree_4.0.4       pheatmap_1.0.13   
##  [5] RColorBrewer_1.1-3 viridis_0.6.5      viridisLite_0.4.3  edgeR_4.8.2       
##  [9] limma_3.66.0       kableExtra_1.4.0   knitr_1.51         reshape2_1.4.5    
## [13] lubridate_1.9.5    forcats_1.0.1      stringr_1.6.0      dplyr_1.2.0       
## [17] purrr_1.2.1        readr_2.2.0        tidyr_1.3.2        tibble_3.3.1      
## [21] ggplot2_4.0.2      tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] ggiraph_0.9.6           tidyselect_1.2.1        farver_2.1.2           
##  [4] S7_0.2.1                fastmap_1.2.0           lazyeval_0.2.2         
##  [7] fontquiver_0.2.1        digest_0.6.39           timechange_0.4.0       
## [10] lifecycle_1.0.5         statmod_1.5.1           tidytree_0.4.7         
## [13] magrittr_2.0.4          compiler_4.5.2          rlang_1.1.7            
## [16] sass_0.4.10             tools_4.5.2             yaml_2.3.12            
## [19] labeling_0.4.3          htmlwidgets_1.6.4       plyr_1.8.9             
## [22] xml2_1.5.2              aplot_0.2.9             withr_3.0.2            
## [25] grid_4.5.2              gdtools_0.5.0           scales_1.4.0           
## [28] MASS_7.3-65             cli_3.6.5               rmarkdown_2.31         
## [31] treeio_1.34.0           generics_0.1.4          rstudioapi_0.18.0      
## [34] tzdb_0.5.0              ape_5.8-1               cachem_1.1.0           
## [37] parallel_4.5.2          ggplotify_0.1.3         yulab.utils_0.2.4      
## [40] vctrs_0.7.1             jsonlite_2.0.0          fontBitstreamVera_0.1.1
## [43] gridGraphics_0.5-1      hms_1.1.4               patchwork_1.3.2        
## [46] systemfonts_1.3.1       locfit_1.5-9.12         jquerylib_0.1.4        
## [49] glue_1.8.0              stringi_1.8.7           gtable_0.3.6           
## [52] pillar_1.11.1           rappdirs_0.3.4          htmltools_0.5.9        
## [55] R6_2.6.1                textshaping_1.0.5       evaluate_1.0.5         
## [58] lattice_0.22-7          ggfun_0.2.0             fontLiberation_0.1.0   
## [61] bslib_0.10.0            Rcpp_1.1.1              svglite_2.2.2          
## [64] gridExtra_2.3           nlme_3.1-168            xfun_0.56              
## [67] fs_1.6.6                pkgconfig_2.0.3

Annex: Exporting ASVs sequences as a fasta file

We know what we want about the ASVs: Are they differentially abundant or not? We will now export their sequences as a FASTA file for the next step of the analysis: Finding the corresponding strains.

## Get the list of all the ASVs
all_asv <- rownames(tt_treat)

## Get their DNA sequence
all_seq <- sequences |> 
  filter(asv %in% all_asv) |> 
  column_to_rownames("asv")

all_seq  |> 
  head() |> 
  kable() |> 
  kable_styling()|> 
  scroll_box(width = "750px")
## Create a function to export the sequences as fasta
export_to_fasta <- function(df, output_file) {
  # Extract row names and sequences
  # We replace spaces with underscores because BLAST+ reads sequence IDs up to the first space
  strain_names <- gsub(" ", "_", rownames(df)) 
  
  # Assuming the sequence is in the first column
  sequences <- as.character(df[, 1]) 
  
  # Create the >Header lines
  headers <- paste0(">", strain_names)
  
  # Interleave headers and sequences (Row 1: Header, Row 2: Seq, Row 3: Header...)
  fasta_lines <- c(rbind(headers, sequences))
  
  # Write the lines to the specified file
  writeLines(fasta_lines, output_file)
  
  cat("Successfully exported", length(strain_names), "sequences to", output_file, "\n")
}

## Save the sequences of all ASV as fasta
export_to_fasta(all_seq, "data/all_ASV.fasta")