The complete import workflow is shown below. To keep the installed vignette quick to build, the executed examples use a deterministic subset containing four samples from each treatment group and the 100 most abundant OTUs. These results illustrate filtering behavior and should not be interpreted as estimates from the full study.
# Get the data files from package
input_files <- pkg.data$input_files
biom_file <- input_files$biom_files$silva # Path to silva biom file
tree_file <- input_files$tree_files$silva # Path to silva tree file
metadata_file <- input_files$metadata$two_groups # Path to Nephele metadata
parse_func <- parse_taxonomy_silva_128 # A custom phyloseq parsing function for silva annotations
# Get the phyloseq object
phy_obj <- create_phyloseq(biom_file = biom_file,
tree_file = tree_file,
metadata_file = metadata_file,
parse_func = parse_func)
# Get the taxmap object in the raw format
raw_metacoder <- as_MicrobiomeR_format(obj = phy_obj, format = "raw_format")Metacoder is a highly useful package that comes with tons of features
right out of the box. Functions such as
metacoder::filter_ambiguous_taxa(),
metacoder::filter_taxa(), and
taxa::filter_obs() for instance can almost always be used
in your workflow.
# Remove Archaea from the taxmap object
metacoder_obj <- metacoder::filter_taxa(
obj = raw_metacoder,
taxon_names == "Archaea",
subtaxa = TRUE,
invert = TRUE
)
# Ambiguous Annotation Filter - Remove taxonomies with ambiguous names
metacoder_obj <- metacoder::filter_ambiguous_taxa(metacoder_obj,
subtaxa = TRUE)Three functions are provided by MicrobiomeR that do some basic filtering:
sample_id_filter() for filtering samples.taxon_id_filter() for filtering by taxon_id, which
includes intermediate taxa.otu_id_filter() for filtering by otu_id, which only
includes leaf taxa.These functions all take the same parameters, most notably a transformation function (.f_transform), a filtering function (.f_filter), and a conditional function (.f_condition).
# Low Sample Filter - Remove the low samples
# The sample filter should generally be implemented first
metacoder_obj <- sample_id_filter(obj = metacoder_obj,
.f_filter = ~sum(.),
.f_condition = ~.>= 20,
validated = TRUE)The advanced filtering functions available with MicrobiomeR do several things. They wrap the basic filtering functions mentioned above, they wrap common metacoder and taxa functions, and they mimic the tools found in the phyloseq package. Below I’ve mentioned some of them and how they relate to the phyloseq package:
agglomerate_taxmap() is equivalent to
phyloseq::tax_glom().otu_proportion_filter() is seen in the first step in
phyloseq’s preprocessing
vignette.Prevalence Filtering functions filter observations by
their prevalence across samples.
otu_prevelance_filter() is seen in
phyloseq’s prevalence
filtering vignette.taxa_prevalence_filter() is seen in
phyloseq’s taxonomic
filtering vignette.cov_filter() is seen in the ninth step in phyloseq’s preprocessing
vignette.
# Master Threshold Filter - Add the otu_proportions table and then filter OTUs based on min %
metacoder_obj <- otu_proportion_filter(obj = metacoder_obj,
otu_percentage = 0.00001)
# Taxon Prevalence Filter - Add taxa_abundance and taxa_proportions and then filter OTUs that do not
# appear more than a certian amount of times in a certain percentage of samples at the specified
# agglomerated rank. This is considered a supervised method, because it relies on intermediate
# taxonomies to filter the data.
# The default minimum abundance is 5 and the sample percentage is 0.5 (5%).
# Phylum
metacoder_obj <- taxa_prevalence_filter(obj = metacoder_obj,
rank = "Phylum")
# Class
metacoder_obj <- taxa_prevalence_filter(obj = metacoder_obj,
rank = "Class",
validated = TRUE)
# Order
metacoder_obj <- taxa_prevalence_filter(obj = metacoder_obj,
rank = "Order",
validated = TRUE)
# OTU Prevalence Filter - Filter OTUs that do not appear more than a certian amount of times in a
# certain percentage of samples. This is considered an unsupervised method, because it relies only
# on the leaf OTU ids to filter the data.
metacoder_obj <- otu_prevalence_filter(obj = metacoder_obj,
validated = TRUE)
# Coefficient of Variation Filter - Filter OTUs based on the coefficient of variation
metacoder_obj <- cov_filter(obj = metacoder_obj,
coefficient_of_variation = 3,
validated = TRUE)As mentioned previously taxmap filtering can be done in any way that
fits your needs with the taxa and metacoder
packages. However, MicrobiomeR also provides some utility based function
for filtering/manipulating your observation data by hand.
Observation data can be accessed within the taxmap object. Make sure you
don’t manipulate existing data inside of your taxmap object unless
you’re absolutely sure you know what you’re doing.
# Get the abundance tables from the taxmap object
taxa_abund <- metacoder_obj$data$taxa_abundance
otu_abund <- metacoder_obj$data$otu_abundance
# Transposing with one ID (taxon_id)
taxa_abund %>% transposer(ids = "taxon_id", header_name = "samples")
#> # A tibble: 8 × 56
#> samples ab ad ae af ag aj ak am an ao aq
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Sample_1 50238 17480 2687 25968 2339 11520 5757 25968 2339 203 956
#> 2 Sample_10 33168 19658 122 13312 76 18146 1227 13312 76 285 95
#> 3 Sample_11 26325 9905 21 13811 2588 5873 3903 13811 2588 129 14
#> 4 Sample_12 28593 10827 25 16727 1014 9209 1428 16727 1014 190 22
#> 5 Sample_35 59356 8640 364 50168 184 8171 292 50168 184 177 287
#> 6 Sample_36 63042 3989 18052 40999 2 3656 4 40999 2 329 0
#> 7 Sample_40 27870 8191 65 16213 3401 7815 259 16213 3401 117 58
#> 8 Sample_41 28389 2582 17 25778 12 2325 241 25778 12 16 0
#> # ℹ 44 more variables: ar <dbl>, at <dbl>, au <dbl>, aw <dbl>, ax <dbl>,
#> # ay <dbl>, ba <dbl>, bb <dbl>, bd <dbl>, be <dbl>, bg <dbl>, bh <dbl>,
#> # bi <dbl>, bj <dbl>, bk <dbl>, bm <dbl>, bn <dbl>, bo <dbl>, bp <dbl>,
#> # bq <dbl>, bs <dbl>, bt <dbl>, bu <dbl>, bw <dbl>, cb <dbl>, cc <dbl>,
#> # cd <dbl>, ce <dbl>, cf <dbl>, cg <dbl>, ci <dbl>, cj <dbl>, cl <dbl>,
#> # cm <dbl>, co <dbl>, cr <dbl>, cs <dbl>, ct <dbl>, cu <dbl>, cv <dbl>,
#> # cw <dbl>, cx <dbl>, dn <dbl>, dy <dbl>
# Re-Transposing with one ID (taxon_id)
taxa_abund %>%
transposer("taxon_id", "samples") %>%
transposer("samples", "taxon_id")
#> # A tibble: 55 × 9
#> taxon_id Sample_1 Sample_10 Sample_11 Sample_12 Sample_35 Sample_36 Sample_40
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 ab 50238 33168 26325 28593 59356 63042 27870
#> 2 ad 17480 19658 9905 10827 8640 3989 8191
#> 3 ae 2687 122 21 25 364 18052 65
#> 4 af 25968 13312 13811 16727 50168 40999 16213
#> 5 ag 2339 76 2588 1014 184 2 3401
#> 6 aj 11520 18146 5873 9209 8171 3656 7815
#> 7 ak 5757 1227 3903 1428 292 4 259
#> 8 am 25968 13312 13811 16727 50168 40999 16213
#> 9 an 2339 76 2588 1014 184 2 3401
#> 10 ao 203 285 129 190 177 329 117
#> # ℹ 45 more rows
#> # ℹ 1 more variable: Sample_41 <dbl>
# Transposing with two IDs (taxon_id, otu_id)
# The column headers will be a combination of the categorical data that will
# be parsed and split back into individual columns if retransposed.
otu_abund %>%
transposer("otu_id", "samples")
#> # A tibble: 8 × 88
#> samples `bd<_>GQ135595.1.1382` `bg<_>EU462500.1.1401` `bn<_>EU475206.1.1402`
#> <chr> <dbl> <dbl> <dbl>
#> 1 Sample_1 681 3260 702
#> 2 Sample_10 10 330 24
#> 3 Sample_11 85 2248 66
#> 4 Sample_12 298 18 197
#> 5 Sample_35 1678 500 567
#> 6 Sample_36 4 2 13
#> 7 Sample_40 3 0 5
#> 8 Sample_41 25 22 16
#> # ℹ 84 more variables: `bn<_>KF842208.1.1404` <dbl>,
#> # `bn<_>KF843564.1.1378` <dbl>, `bo<_>FJ679534.1.1449` <dbl>,
#> # `bo<_>JQ190140.1.1370` <dbl>, `bq<_>JQ191146.1.1321` <dbl>,
#> # `bs<_>EU462462.1.1381` <dbl>, `bs<_>JQ190288.1.1362` <dbl>,
#> # `bt<_>EU472593.1.1401` <dbl>, `bu<_>HQ785970.1.1442` <dbl>,
#> # `bw<_>EU461194.1.1423` <dbl>, `cb<_>DQ804726.1.1392` <dbl>,
#> # `cc<_>EU472179.1.1375` <dbl>, `cd<_>EU467376.1.1388` <dbl>, …
# Re-Transposing with two IDs (taxon_id, otu_id)
# When you transpose categorical data, the column headers are lost.
# To re-transpose you have to supply these headers
otu_abund %>%
transposer("otu_id", "samples") %>%
transposer("samples", "otu_id",
separated_categories = c("taxon_id", "otu_id"))
#> # A tibble: 87 × 10
#> taxon_id otu_id Sample_1 Sample_10 Sample_11 Sample_12 Sample_35 Sample_36
#> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 bd GQ135595… 681 10 85 298 1678 4
#> 2 bg EU462500… 3260 330 2248 18 500 2
#> 3 bn EU475206… 702 24 66 197 567 13
#> 4 bn KF842208… 700 60 132 81 466 12
#> 5 bn KF843564… 2229 75 83 487 805 7
#> 6 bo FJ679534… 191 504 23 264 159 0
#> 7 bo JQ190140… 649 827 1390 393 484 2
#> 8 bq JQ191146… 1731 27 7 3 77 0
#> 9 bs EU462462… 451 196 144 165 228 51
#> 10 bs JQ190288… 488 17 59 105 199 20
#> # ℹ 77 more rows
#> # ℹ 2 more variables: Sample_40 <dbl>, Sample_41 <dbl>
# Transforming to proportions/percentages by COLUMN
taxa_abund %>%
transformer(func = ~./sum(.))
#> # A tibble: 55 × 9
#> taxon_id Sample_1 Sample_10 Sample_11 Sample_12 Sample_35 Sample_40 Sample_41
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 ab 0.177 0.168 0.171 0.168 0.169 0.169 0.167
#> 2 ad 0.0616 0.0995 0.0643 0.0637 0.0246 0.0496 0.0152
#> 3 ae 0.00947 0.000618 0.000136 0.000147 0.00103 0.000394 0.0000999
#> 4 af 0.0916 0.0674 0.0896 0.0984 0.143 0.0982 0.151
#> 5 ag 0.00825 0.000385 0.0168 0.00596 0.000523 0.0206 0.0000705
#> 6 aj 0.0406 0.0919 0.0381 0.0542 0.0232 0.0473 0.0137
#> 7 ak 0.0203 0.00621 0.0253 0.00840 0.000830 0.00157 0.00142
#> 8 am 0.0916 0.0674 0.0896 0.0984 0.143 0.0982 0.151
#> 9 an 0.00825 0.000385 0.0168 0.00596 0.000523 0.0206 0.0000705
#> 10 ao 0.000716 0.00144 0.000837 0.00112 0.000503 0.000709 0.0000940
#> # ℹ 45 more rows
#> # ℹ 1 more variable: Sample_36 <dbl>
# Transforming by ROW is also possible
# This function will transpose/retranspose so other information is needed
taxa_abund %>%
transformer(by = "row",
func = ~./sum(.),
ids = "taxon_id",
header_name = "samples",
separated_categories = c("taxon_id"))
#> # A tibble: 55 × 9
#> taxon_id Sample_1 Sample_10 Sample_11 Sample_12 Sample_35 Sample_36 Sample_40
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 ab 0.158 0.105 0.0830 0.0902 0.187 0.199 0.0879
#> 2 ad 0.215 0.242 0.122 0.133 0.106 0.0491 0.101
#> 3 ae 0.126 0.00571 0.000983 0.00117 0.0170 0.845 0.00304
#> 4 af 0.128 0.0656 0.0680 0.0824 0.247 0.202 0.0799
#> 5 ag 0.243 0.00790 0.269 0.105 0.0191 0.000208 0.354
#> 6 aj 0.173 0.272 0.0880 0.138 0.122 0.0548 0.117
#> 7 ak 0.439 0.0936 0.298 0.109 0.0223 0.000305 0.0198
#> 8 am 0.128 0.0656 0.0680 0.0824 0.247 0.202 0.0799
#> 9 an 0.243 0.00790 0.269 0.105 0.0191 0.000208 0.354
#> 10 ao 0.140 0.197 0.0892 0.131 0.122 0.228 0.0809
#> # ℹ 45 more rows
#> # ℹ 1 more variable: Sample_41 <dbl>This function was # borrowed from .
# Get analyzed Data
metacoder_obj <- as_MicrobiomeR_format(obj = metacoder_obj, format = "analyzed_format")
# Create agglomerated taxmaps for phylum/class
phylum_mo <- agglomerate_taxmap(obj = metacoder_obj, rank = "Phylum")
class_mo <- agglomerate_taxmap(obj = metacoder_obj, rank = "Class")
# Get some observation data
phylum_data <- phylum_mo$data$stats_tax_data
class_data <- class_mo$data$stats_tax_data
# Take the Phylum in the "phylum_data" and cross reference these in the "class_data".
# Reutrn the "wilcox_p_value" of the class_data.
class_p_value <- vlookup(lookup_vector = phylum_data$Phylum,
df = class_data,
match_var = "Phylum",
return_var = "wilcox_p_value")
# Create a new column in the phylum_data.
new_data <- phylum_data %>% dplyr::mutate(class_p_value = class_p_value)
new_data[c("taxon_id", "Phylum", "Class", "class_p_value", "wilcox_p_value")]
#> # A tibble: 5 × 5
#> taxon_id Phylum Class class_p_value wilcox_p_value
#> <chr> <chr> <chr> <dbl> <dbl>
#> 1 ab NA NA 1 1
#> 2 ad Firmicutes NA 0.0286 0.0286
#> 3 ae Proteobacteria NA 0.886 0.886
#> 4 af Bacteroidetes NA 0.0571 0.0571
#> 5 ag Spirochaetae NA 0.486 0.486