# This function is modified from import_qiime2.R in MicrobiotaProcess
# https://github.com/YuLab-SMU/MicrobiotaProcess/blob/master/R/import_qiime2.R
# and https://github.com/jbisanz/qiime2R
#' Import function to read the the output of dada2 as phyloseq object
#'
#' Import the qiime2 artifacts, including feature table, taxonomic table,
#' phylogenetic tree, representative sequence and sample metadata into
#' phyloseq object.
#'
#' @param otu_qza character, file path of the feature table from qiime2.
#' @param taxa_qza character, file path of the taxonomic table from qiime2,
#' default `NULL`.
#' @param sam_tab character, file path of the sample metadata in tsv format,
#' default `NULL`.
#' @param refseq_qza character, file path of the representative sequences from
#' qiime2, default `NULL`.
#' @param tree_qza character, file path of the phylogenetic tree from
#' qiime2, default `NULL`.
#' @export
#' @return [`phyloseq::phyloseq-class`] object.
#' @examples
#' otuqza_file <- system.file(
#' "extdata", "table.qza",
#' package = "microbiomeMarker"
#' )
#' taxaqza_file <- system.file(
#' "extdata", "taxonomy.qza",
#' package = "microbiomeMarker"
#' )
#' sample_file <- system.file(
#' "extdata", "sample-metadata.tsv",
#' package = "microbiomeMarker"
#' )
#' treeqza_file <- system.file(
#' "extdata", "tree.qza",
#' package = "microbiomeMarker"
#' )
#' ps <- import_qiime2(
#' otu_qza = otuqza_file, taxa_qza = taxaqza_file,
#' sam_tab = sample_file, tree_qza = treeqza_file
#' )
#' ps
import_qiime2 <- function(otu_qza,
taxa_qza = NULL,
sam_tab = NULL,
refseq_qza = NULL,
tree_qza = NULL) {
feature_tab <- read_qza(otu_qza)
if (!is.null(taxa_qza)) {
taxa_tab <- read_qza(taxa_qza)
taxa_tab <- subset_taxa_in_feature(taxa_tab, feature_tab)
taxa_tab <- parse_q2taxonomy(taxa_tab)
} else {
taxa_tab <- NULL
}
if (!is.null(sam_tab)) {
sam_tab <- read_q2sample_meta(sam_tab)
} else {
sam_tab <- NULL
}
if (!is.null(tree_qza)) {
tree <- read_qza(tree_qza)
} else {
tree <- NULL
}
# check the row.names of feature table: DNA sequence or not
# if is DNA, row.names(feature_tab) is set as refseq
if (is_dna_seq(row.names(feature_tab))) {
refseq <- row.names(feature_tab)
refseq_nm <- paste0("OTU", seq_along(refseq))
names(refseq) <- refseq_nm
# set the rownames of feature and taxa tab as OTU1, OTU2,...
if (!is.null(taxa_tab)) {
rownames(taxa_tab) <- refseq_nm
}
rownames(feature_tab) <- refseq_nm
} else {
refseq <- NULL
}
if (!is.null(refseq_qza)) {
refseq <- read_qza(refseq_qza)
}
if (!is.null(refseq)) {
refseq <- Biostrings::DNAStringSet(refseq)
} else {
refseq <- NULL
}
ps <- phyloseq(feature_tab, taxa_tab, sam_tab, tree, refseq)
if (inherits(ps, "otu_table")) {
warning("`otu_table` object is returned")
}
ps
}
#' Read the qza file output from qiime2
#'
#' Import the qiime2 artifacts to R.
#' @param file character, path of the input qza file. Only files in format of
#' `BIOMV210DirFmt` (feature table), `TSVTaxonomyDirectoryFormat` (taxonomic
#' table), `NewickDirectoryFormat` (phylogenetic tree ) and
#' `DNASequencesDirectoryFormat` (representative sequences) are supported
#' right now.
#' @param temp character, a temporary directory in which the qza file will be
#' decompressed to, default `tempdir()`.
#' @importFrom yaml read_yaml
#' @return [`phyloseq::otu_table-class`] object for feature table,
#' [`phyloseq::taxonomyTable-class`] object for taxonomic table,
#' [`ape::phylo`] object for phylogenetic tree,
#' [`Biostrings::DNAStringSet-class`] for representative sequences of taxa.
#' @noRd
read_qza <- function(file, temp = tempdir()) {
unzipped_file <- utils::unzip(file, exdir = temp)
meta_file <- grep("metadata.yaml", unzipped_file, value = TRUE)
metadata <- yaml::read_yaml(meta_file[1])
format <- metadata$format
uuid <- metadata$uuid
if (grepl("BIOMV", metadata$format)) {
biom_file <- file.path(temp, uuid, "data/feature-table.biom")
res <- read_q2biom(biom_file)
} else if (format == "TSVTaxonomyDirectoryFormat") {
taxa_file <- file.path(temp, uuid, "data/taxonomy.tsv")
res <- read_q2taxa(taxa_file)
} else if (format == "NewickDirectoryFormat") {
tree_file <- file.path(temp, uuid, "data/tree.nwk")
res <- read_tree(tree_file)
} else if (format == "DNASequencesDirectoryFormat") {
refseq_file <- file.path(temp, uuid, "data/dna-sequences.fasta")
res <- Biostrings::readDNAStringSet(refseq_file)
} else {
stop(
"Only files in format of 'BIOMV210DirFmt' ",
"'TSVTaxonomyDirectoryFormat', ",
"'NewickDirectoryFormat' and 'DNASequencesDirectoryFormat'",
" are supported."
)
}
res
}
#' Read qiime2 feature table
#'
#' Import qiime2 feature table into otu_table object
#'
#' @param file character, file name of the biom file.
#' @importFrom biomformat read_biom biom_data
#' @return [`phyloseq::otu_table-class`] object.
#' @noRd
read_q2biom <- function(file) {
biomobj <- suppressWarnings(read_biom(file))
feature_tab <- as(biom_data(biomobj), "matrix")
otu_table(feature_tab, taxa_are_rows = TRUE)
}
#' Read qiime2 taxa file
#' @keywords internal
#' @noRd
read_q2taxa <- function(file) {
taxa <- utils::read.table(file, sep = "\t", header = TRUE)
if ("Confidence" %in% names(taxa)) {
taxa$Confidence <- NULL
}
tax_table(as.matrix(taxa))
}
#' Read qiime2 sample meta data file
#' @keywords internal
#' @noRd
read_q2sample_meta <- function(file) {
phyloseq::import_qiime_sample_data(file)
}
#' Subset taxa according to the taxa in feature table
#' @keywords internal
#' @noRd
subset_taxa_in_feature <- function(taxa_tab, feature_tab) {
idx <- match(row.names(feature_tab), taxa_tab[, "Feature.ID"])
taxa_tab <- taxa_tab[idx, , drop = FALSE]
row.names(taxa_tab) <- taxa_tab[, "Feature.ID"]
taxa_tab
}
#' Parse qiime2 taxa in different taxonomic levels
#' @param taxa `tax_table` object.
#' @param sep character string containing a regular expression, separator
#' between different taxonomic levels, defaults to on compatible with both
#' GreenGenes and SILVA `; |;"`.
#' @param trim_rank_prefix logical whether remove leading characters from
#' taxonomic levels, e.g. k__ or D_0__, default `TRUE`.
#' @return [`phyloseq::taxonomyTable-class`] object.
#' @keywords internal
#' @noRd
parse_q2taxonomy <- function(taxa, sep = "; |;", trim_rank_prefix = TRUE) {
taxa <- data.frame(taxa)
if (trim_rank_prefix) {
# remove leading characters from GG
taxa$Taxon <- gsub("[kpcofgs]__", "", taxa$Taxon)
# remove leading characters from SILVA
taxa$Taxon <- gsub("D_\\d__", "", taxa$Taxon)
}
taxa <- tidyr::separate(taxa, .data$Taxon,
c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"),
sep = sep,
fill = "right"
)
taxa <- purrr::map_df(taxa, ~ ifelse(.x == "", NA_character_, .x)) %>%
as.data.frame()
rownames(taxa) <- taxa$Feature.ID
taxa$Feature.ID <- NULL
tax_table(as.matrix(taxa))
}
#' check the row.names of feature table is DNA sequence or not
#' This function is from
#' {github}YuLab-SMU/MicrobiotaProcess/blob/master/R/import_qiime2.R#L169-L177
#' @keywords internal
#' @noRd
is_dna_seq <- function(names) {
x <- unlist(strsplit(toupper(names[1]), split = ""))
freq <- mean(x %in% c("A", "G", "T", "C", "N", "X", "-"))
if (length(x) > 30 & freq > 0.9) {
return(TRUE)
} else {
return(FALSE)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.