Skip to contents
library(MiscMetabar)
#> Warning: replacing previous import 'S4Arrays::makeNindexFromArrayViewport' by
#> 'DelayedArray::makeNindexFromArrayViewport' when loading 'SummarizedExperiment'
library(formattable)

Filter samples

Phyloseq package already propose a function to select samples (subset_samples()), but in some case, the subset internal function is painful. MiscMetabar propose a complementary function ([subset_samples_pq()]) which is more versatile but need to be used with caution because the order of the condition must match the orders of the samples.

data(data_fungi)
cond_samp <- grepl("A1", data_fungi@sam_data[["Sample_names"]])
subset_samples_pq(data_fungi, cond_samp)
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 1420 taxa and 9 samples ]
#> sample_data() Sample Data:       [ 9 samples by 7 sample variables ]
#> tax_table()   Taxonomy Table:    [ 1420 taxa by 12 taxonomic ranks ]
#> refseq()      DNAStringSet:      [ 1420 reference sequences ]

Filter Taxa

Filter taxa using condition(s)

Phyloseq package already propose a function to select samples (subset_taxa()), but in some case, the subset internal function is painful. MiscMetabar propose a complementary function ([subset_taxa_pq()]) which is more versatile and is based on the names of the taxa to match the condition and the taxa in the phyloseq object. In the example code above, we filter fungi from the Phylum “Basidiomycota” using phyloseq::subset_taxa() and then select only taxa with more than 1000 nb_sequences using subset_taxa_pq().

df_basidio <- subset_taxa(data_fungi, Phylum == "Basidiomycota")
df_basidio_abundant <- subset_taxa_pq(
  df_basidio,
  colSums(df_basidio@otu_table) > 1000
)

Filter taxa using blast score against a database

In some cases, we want to select only a given clade of taxon. One solution is to select only taxa with information given by taxonomic assignment (e.g. using the function dada2::assignTaxonomy()). However, in some cases, this information lead to false positive and false negative selection. MiscMetabar implement an other method relying on blast software (Altschul et al. 1990 & 1997) and database. The idea is to set a cutoff in four parameters to select only taxa which are close enough to sequences in the database :

  • id_filter (default: 90) cut of in identity percent to keep result.
  • bit_score_filter (default: 50) cut of in bit score to keep result. The higher the bit-score, the better the sequence similarity. The bit-score is the requires size of a sequence database in which the current match could be found just by chance. The bit-score is a log2 scaled and normalized raw-score. Each increase by one doubles the required database size (2bit-score).
  • min_cover_filter (default: 50) cut of in query cover (%) to keep result.
  • e_value_filter (default: 1e-30) cut of in e-value (%) to keep result. The BLAST E-value is the number of expected hits of similar quality (score) that could be found just by chance.
path_db <- system.file("extdata",
  "100_sp_UNITE_sh_general_release_dynamic.fasta",
  package = "MiscMetabar", mustWork = TRUE
)

suppressWarnings(blast_error_or_not <-
  try(system("blastn 2>&1", intern = TRUE), silent = TRUE))

if (!is(blast_error_or_not, "try-error")) {
  df_blast_80 <- filter_asv_blast(df_basidio, fasta_for_db = path_db)
  df_blast_50 <- filter_asv_blast(df_basidio,
    fasta_for_db = path_db,
    id_filter = 50, e_value_filter = 10,
    bit_score_filter = 20, min_cover_filter = 20
  )

  track_formattable <-
    track_wkflow(
      list(
        "raw data" = df_basidio,
        "id_filter = 80" = df_blast_80,
        "id_filter = 50" = df_blast_50
      )
    )
}
if (!is(blast_error_or_not, "try-error")) {
  formattable(
    track_formattable,
    list(
      area(col = nb_sequences) ~ color_bar("cyan", na.rm = TRUE),
      area(col = nb_clusters) ~ normalize_bar("yellowgreen",
        na.rm = TRUE, min = 0.3
      ),
      area(col = nb_samples) ~ normalize_bar("lightpink",
        na.rm = TRUE, min = 0.3
      )
    )
  )
}
nb_sequences nb_clusters nb_samples
raw data 784054 345 185
id_filter = 80 249494 45 127
id_filter = 50 780242 326 172

Filter taxa using a known taxa for control

To filter out contamination, one solution is to add a proportion of a known taxa which is not present in the environment of the study. In that case we can define some threshold for each sample to discard taxon based on pseudo-abundance. In the example code above, we select taxon using the ASV_50 as control through 6 different algorithms.

res_seq <-
  suppressWarnings(
    subset_taxa_tax_control(data_fungi, as.numeric(data_fungi@otu_table[, 50]),
      method = "cutoff_seq"
    )
  )
res_mixt <-
  suppressWarnings(
    subset_taxa_tax_control(data_fungi, as.numeric(data_fungi@otu_table[, 50]),
      method = "cutoff_mixt"
    )
  )
res_diff <- suppressWarnings(
  subset_taxa_tax_control(
    data_fungi,
    as.numeric(data_fungi@otu_table[, 50]),
    method = "cutoff_diff",
    min_diff_for_cutoff = 2
  )
)
res_min <-
  suppressWarnings(
    subset_taxa_tax_control(
      data_fungi,
      as.numeric(data_fungi@otu_table[, 50]),
      method = "min",
      min_diff_for_cutoff = 2
    )
  )
res_max <-
  suppressWarnings(
    subset_taxa_tax_control(
      data_fungi,
      as.numeric(data_fungi@otu_table[, 50]),
      method = "max",
      min_diff_for_cutoff = 2
    )
  )
res_mean <-
  suppressWarnings(
    subset_taxa_tax_control(
      data_fungi,
      as.numeric(data_fungi@otu_table[, 50]),
      method = "mean",
      min_diff_for_cutoff = 2
    )
  )
track_formattable <- track_wkflow(list(
  "raw data" = data_fungi,
  "cutoff_seq" = res_seq,
  "cutoff_mixt" = res_mixt,
  "cutoff_diff" = res_diff,
  "min" = res_min,
  "max" = res_max,
  "mean" = res_mean
))

formattable(
  track_formattable,
  list(
    area(col = nb_sequences) ~ color_bar("cyan"),
    area(col = nb_clusters) ~ normalize_bar("yellowgreen",
      na.rm = TRUE,
      min = 0.3
    ),
    area(col = nb_samples) ~ normalize_bar("lightpink", na.rm = TRUE)
  )
)
nb_sequences nb_clusters nb_samples
raw data 1839124 1420 185
cutoff_seq 1817543 1402 185
cutoff_mixt 1768917 1273 185
cutoff_diff 1805274 1378 185
min 1836038 1417 185
max 1754990 1222 185
mean 1789186 1317 185

Session information

sessionInfo()
#> R version 4.3.2 (2023-10-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] formattable_0.2.1 MiscMetabar_0.7.8 dplyr_1.1.4       dada2_1.30.0     
#> [5] Rcpp_1.0.12       ggplot2_3.4.4     phyloseq_1.46.0  
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-7                pbapply_1.7-2              
#>   [3] deldir_2.0-2                permute_0.9-7              
#>   [5] rlang_1.1.3                 magrittr_2.0.3             
#>   [7] ade4_1.7-22                 matrixStats_1.2.0          
#>   [9] compiler_4.3.2              mgcv_1.9-0                 
#>  [11] png_0.1-8                   systemfonts_1.0.5          
#>  [13] vctrs_0.6.5                 reshape2_1.4.4             
#>  [15] stringr_1.5.1               pkgconfig_2.0.3            
#>  [17] crayon_1.5.2                fastmap_1.1.1              
#>  [19] XVector_0.42.0              utf8_1.2.4                 
#>  [21] Rsamtools_2.18.0            rmarkdown_2.25             
#>  [23] ragg_1.2.7                  purrr_1.0.2                
#>  [25] xfun_0.42                   zlibbioc_1.48.0            
#>  [27] cachem_1.0.8                GenomeInfoDb_1.38.6        
#>  [29] jsonlite_1.8.8              biomformat_1.30.0          
#>  [31] rhdf5filters_1.14.1         DelayedArray_0.28.0        
#>  [33] Rhdf5lib_1.24.2             BiocParallel_1.36.0        
#>  [35] jpeg_0.1-10                 parallel_4.3.2             
#>  [37] cluster_2.1.4               R6_2.5.1                   
#>  [39] bslib_0.6.1                 stringi_1.8.3              
#>  [41] RColorBrewer_1.1-3          GenomicRanges_1.54.1       
#>  [43] jquerylib_0.1.4             SummarizedExperiment_1.32.0
#>  [45] iterators_1.0.14            knitr_1.45                 
#>  [47] mixtools_2.0.0              IRanges_2.36.0             
#>  [49] Matrix_1.6-1.1              splines_4.3.2              
#>  [51] igraph_2.0.1.1              tidyselect_1.2.0           
#>  [53] abind_1.4-5                 yaml_2.3.8                 
#>  [55] vegan_2.6-4                 codetools_0.2-19           
#>  [57] hwriter_1.3.2.1             lattice_0.21-9             
#>  [59] tibble_3.2.1                plyr_1.8.9                 
#>  [61] Biobase_2.62.0              withr_3.0.0                
#>  [63] ShortRead_1.60.0            evaluate_0.23              
#>  [65] desc_1.4.3                  survival_3.5-7             
#>  [67] RcppParallel_5.1.7          kernlab_0.9-32             
#>  [69] Biostrings_2.70.2           pillar_1.9.0               
#>  [71] MatrixGenerics_1.14.0       foreach_1.5.2              
#>  [73] stats4_4.3.2                plotly_4.10.4              
#>  [75] generics_0.1.3              RCurl_1.98-1.14            
#>  [77] S4Vectors_0.40.2            munsell_0.5.0              
#>  [79] scales_1.3.0                glue_1.7.0                 
#>  [81] lazyeval_0.2.2              tools_4.3.2                
#>  [83] interp_1.1-6                data.table_1.15.0          
#>  [85] GenomicAlignments_1.38.2    fs_1.6.3                   
#>  [87] rhdf5_2.46.1                grid_4.3.2                 
#>  [89] tidyr_1.3.1                 ape_5.7-1                  
#>  [91] latticeExtra_0.6-30         colorspace_2.1-0           
#>  [93] nlme_3.1-163                GenomeInfoDbData_1.2.11    
#>  [95] cli_3.6.2                   textshaping_0.3.7          
#>  [97] fansi_1.0.6                 segmented_2.0-2            
#>  [99] viridisLite_0.4.2           S4Arrays_1.2.0             
#> [101] gtable_0.3.4                sass_0.4.8                 
#> [103] digest_0.6.34               BiocGenerics_0.48.1        
#> [105] SparseArray_1.2.4           htmlwidgets_1.6.4          
#> [107] memoise_2.0.1               htmltools_0.5.7            
#> [109] pkgdown_2.0.7               multtest_2.58.0            
#> [111] lifecycle_1.0.4             httr_1.4.7                 
#> [113] MASS_7.3-60