How to download taxa profile of multiple samples

Data derived from curatedMetagenomicData

Table Of Contents

Introduction

R package curatedMetagenomicData provides more than 20,000 metaphlan’s results of samples. However, to download all the results seems not easy when using simple R codes. Here, we give users codes with step by step to grap all data.

Since phyloseq object contains abundance table, taxa table and metadata etc, we convert the final metadata and profile including absolute and relative abundance into phyloseq object.

Loading R packages

suppressMessages(library(dplyr))
suppressMessages(library(tibble))
suppressMessages(library(curatedMetagenomicData))
library(phyloseq)

Studies in curatedMetagenomicData

All the studies’ information were stored in sampleMetadata. Therefore, we obtain the studies’ names by extracting the information from it.

data("sampleMetadata")
availablediseases <- pull(sampleMetadata, study_condition) %>%
  table() %>%
  sort(decreasing = TRUE)

studies <- lapply(names(availablediseases), function(x) {
  dplyr::filter(sampleMetadata, study_condition %in% x) %>%
    dplyr::pull(study_name) %>%
    unique()
})
names(studies) <- names(availablediseases)

Acquiring absolute and relative abundance per study

  • Merging metadata without duplicated samples

  • Whether there is absolute or relative abundance profile per study

  • Whether the samples are identical between metadata and profile

metadata <- data.frame()
profile_counts <- data.frame()
profile_rb <- data.frame()

for (i in 1:length(studies)) {
  disease_dir <- paste0("./dataset/", names(studies)[i])
  
  if (!dir.exists(disease_dir)) {
    dir.create(disease_dir, recursive = TRUE)
  }
  
  for (disease in studies[[i]]) {
    rb_pattern <- paste0(disease, ".relative_abundance")
    
    # absolute abundance
    dat_counts_list <- curatedMetagenomicData(
      pattern = rb_pattern,
      dryrun = FALSE,
      counts = TRUE,
      rownames = "long")
      
    # relative_abundance
    dat_RB_list <- curatedMetagenomicData(
      pattern = rb_pattern,
      dryrun = FALSE,
      counts = FALSE,
      rownames = "long")
    
    # whether there is a taxa profile
    if (length(dat_counts_list) > 0) {
      dat_counts <- dat_counts_list[[1]]
      disease_path_counts <- paste0(disease_dir, "/", paste0(disease, "_counts.RData"))
      save(dat_counts, file = disease_path_counts)
      
      # metadata & profile
      temp_metadata <- colData(dat_counts) %>%
        data.frame()
      temp_profile_counts <- assay(dat_counts) %>%
        data.frame() %>%
        tibble::rownames_to_column("TaxaID")
      
      if (length(profile_counts) == 0) {
        profile_counts <- temp_profile_counts
      } else {
        
        colnames(profile_counts) <- gsub("-", ".", colnames(profile_counts))
        colnames(temp_profile_counts) <- gsub("-", ".", colnames(temp_profile_counts))
        
        # the overlap samples
        sid_count <- intersect(colnames(temp_profile_counts), colnames(profile_counts))
        sid_count <- sid_count[-which(sid_count == "TaxaID")]
        if (length(sid_count) > 0) {
          profile_counts <- profile_counts[, !colnames(profile_counts)%in%sid_count]
        } 
        profile_counts <- dplyr::full_join(temp_profile_counts, 
                                           profile_counts, 
                                           by = "TaxaID")
        
        colnames(profile_counts) <- gsub("-", ".", colnames(profile_counts))
      }
    }
    # whether there is a taxa profile
    if (length(dat_RB_list) > 0) {
      dat_RB <- dat_RB_list[[1]]
      disease_path_RB <- paste0(disease_dir, "/", paste0(disease, "_RB.RData"))
      save(dat_RB, file = disease_path_RB) 
      
      # metadata & profile
      temp_metadata <- colData(dat_RB) %>%
        data.frame()
      temp_profile_rb <- assay(dat_RB) %>%
        data.frame() %>%
        tibble::rownames_to_column("TaxaID")
      
      if (length(profile_rb) == 0) {
        profile_rb <- temp_profile_rb
      } else {
        
        colnames(profile_rb) <- gsub("-", ".", colnames(profile_rb))
        colnames(temp_profile_rb) <- gsub("-", ".", colnames(temp_profile_rb))
        
        # the overlap samples
        sid_rb <- intersect(colnames(temp_profile_rb), colnames(profile_rb))
        sid_rb <- sid_rb[-which(sid_rb == "TaxaID")]
        if (length(sid_rb) > 0) {
          profile_rb <- profile_rb[, !colnames(profile_rb)%in%sid_rb]
        }        
        
        profile_rb <- dplyr::full_join(temp_profile_rb, 
                                       profile_rb, 
                                       by = "TaxaID")
        
        colnames(profile_rb) <- gsub("-", ".", colnames(profile_rb))
      }
      
    }    
    
    if (any(length(dat_counts_list) > 0, length(dat_RB_list) > 0)) {
      
      if (length(metadata) == 0) {
        metadata <- temp_metadata
      } else {
        # diff: metadata vs temp_metadata
        diff1 <- setdiff(colnames(metadata), colnames(temp_metadata))
        diff1_df <- data.frame(matrix(NA, nrow = nrow(temp_metadata), ncol = length(diff1)),
                               row.names = rownames(temp_metadata))
        colnames(diff1_df) <- diff1
        temp_metadata_new <- cbind(temp_metadata, diff1_df)
        
        # diff: temp_metadata vs metadata
        diff2 <- setdiff(colnames(temp_metadata), colnames(metadata))
        diff2_df <- data.frame(matrix(NA, nrow = nrow(metadata), ncol = length(diff2)),
                               row.names = rownames(metadata))
        colnames(diff2_df) <- diff2
        metadata_new <- cbind(metadata, diff2_df) 
        
        rownames(metadata) <- gsub("-", ".", rownames(metadata))
        rownames(temp_metadata) <- gsub("-", ".", rownames(temp_metadata))        
        
        # the overlap samples
        sid <- intersect(rownames(temp_metadata), rownames(metadata))
        if (length(sid) > 0) {
          metadata <- metadata[!rownames(metadata)%in%sid, ]
        }
        
        metadata <- rbind(temp_metadata_new, metadata_new) %>%
          dplyr::distinct()
        
        rownames(metadata) <- gsub("-", ".", rownames(metadata))
      }      
    }
  }
}

Converting datasets into phyloseq

  • Extracting taxonomic levels from 1st column of profile

  • Removing prevalence of features less than 0

# extract taxnomic levels
import_metaphlan_taxa <- function(data_metaphlan2,
                                  taxa_level = NULL) {

  taxa_rank <- c("Kingdom", "Phylum", "Class",
                 "Order", "Family", "Genus",
                 "Species", "Strain")

  # rename 1st column into "ID"
  colnames(data_metaphlan2)[which(colnames(data_metaphlan2) == "ID" |
                                    colnames(data_metaphlan2) == "clade_name")] <- "ID"

  # remove the "NCBI_tax_id" column
  if (any(colnames(data_metaphlan2) %in% "NCBI_tax_id")) {
    data_metaphlan2 <- data_metaphlan2 %>% dplyr::select(-NCBI_tax_id)
  }

  # remove sample names with "metaphlan_bugs_list" suffix
  colnames(data_metaphlan2) <- gsub("_metaphlan_bugs_list$", "", colnames(data_metaphlan2))

  if (is.null(taxa_level)) {
    taxa_level <- "Species"
  }
  ind_number <- match(taxa_level, taxa_rank)

  # remove k__Archaea & k__Viruses & k__Eukaryota & Others which are not bacteria
  if (length(grep("k__Bacteria", data_metaphlan2$ID)) > 0) {
    taxa_bacteria <- data_metaphlan2 %>%
      dplyr::slice(grep("k__Bacteria", ID)) %>%
      tibble::column_to_rownames("ID")
  } else {
    stop("There are no bacteria identified in this profile")
  }

  # dividing by 100 to normalize the sum into 1
  taxa_bac_per <- (taxa_bacteria / 100)  %>%
    tibble::rownames_to_column("ID")

  taxa <- taxa_bac_per %>%
    dplyr::group_by(ID) %>%
    dplyr::mutate(Number=length(unlist(strsplit(ID, "|", fixed = TRUE)))) %>%
    dplyr::select(ID, Number) %>%
    dplyr::filter(Number == ind_number)

  abundance_table <- taxa_bac_per %>%
    dplyr::filter(ID%in%taxa$ID)

  taxa_table <- lapply(abundance_table$ID, strsplit, split = "|", fixed = TRUE)
  taxa_table <- lapply(taxa_table, unlist)
  taxa_table <- do.call(rbind, taxa_table)
  taxa_table <- data.frame(taxa_table)
  names(taxa_table) <- taxa_rank[1:ind_number]
  rownames(taxa_table) <- taxa_table[, ind_number]

  abundance_table_res <- abundance_table %>%
    tibble::column_to_rownames("ID")
  rownames(abundance_table_res) <- rownames(taxa_table)
  abundance_table_res <- abundance_table_res[match(rownames(taxa_table),
                                                   row.names(abundance_table_res)), ,
                                             drop = FALSE]

  res <- list(tax_tab=taxa_table,
              abu_tab=abundance_table_res)

  return(res)
}

# Import function to read the the output of metaphlan as phyloseq object
get_metaphlan_phyloseq <- function(otu_tab,
                                   sam_tab,
                                   tax_tab = NULL) {

  if (!is.null(tax_tab)) {
    tax_tab <- phyloseq::tax_table(as(tax_tab, "matrix"))
  }

  if (!is.null(sam_tab)) {
    sam_tab <- phyloseq::sample_data(sam_tab %>% data.frame())
  }

  asv_tab <- phyloseq::otu_table(as(otu_tab, "matrix"),
                                 taxa_are_rows = TRUE)

  ps <- phyloseq::phyloseq(asv_tab, tax_tab, sam_tab)

  return(ps)
}

# Trim samples or taxa in `otu_table` by Occurrences
run_trim <- function(object,
                     taxa_level = NULL,
                     cutoff = 0.1,
                     trim = c("none", "both", "feature", "sample")){

  trim <- match.arg(trim, c("none", "both", "feature", "sample"))
  if (inherits(object, "phyloseq")) {
    # ps <- preprocess_ps(object)
    if (!is.null(taxa_level)) {
      ps_taxa <- summarize_taxa(object, taxa_level = taxa_level)
    } else {
      ps_taxa <- object
    }
    prf <- as(phyloseq::otu_table(ps_taxa), "matrix")
  } else if (inherits(object, "environment")) {
    prf <- as(object$.Data, "matrix")
  } else {
    prf <- object
  }

  if (trim == "feature") {
    tmp1 <- trim_FeatureOrSample(prf, 1, cutoff)
    remain_features <- rownames(tmp1)
    remain_samples <- colnames(prf)
  } else if (trim == "sample") {
    tmp2 <- trim_FeatureOrSample(prf, 2, cutoff)
    remain_features <- rownames(prf)
    remain_samples <- rownames(tmp2)
  } else if(trim == "both") {
    tmp1 <- trim_FeatureOrSample(prf, 1, cutoff)
    tmp2 <- trim_FeatureOrSample(prf, 2, cutoff)
    remain_features <- rownames(tmp1)
    remain_samples <- rownames(tmp2)
  } else if(trim == "none") {
    return(object)
  }

  if (length(remain_features) > 1 & length(remain_samples) > 1) {
    prf_remain <- prf[remain_features, remain_samples]
  } else if (length(remain_features) > 1 & length(remain_samples) == 1) {
    prf_remain <- prf[remain_features, remain_samples] %>%
      data.frame()
    colnames(prf_remain) <- remain_samples
  } else if (length(remain_features) == 1 & length(remain_samples) > 1) {
    prf_remain <- prf[remain_features, remain_samples] %>%
      data.frame()
    rownames(prf_remain) <- remain_features
  } else if (length(remain_features) == 1 & length(remain_samples) == 1) {
    prf_remain <- prf[remain_features, remain_samples] %>%
      data.frame()

    colnames(prf_remain) <- remain_samples
    rownames(prf_remain) <- remain_features
  } else {
    stop("No sample or taxa were remained, please reinput your cutoff")
  }

  if (inherits(object, "phyloseq")) {
    otu_table(ps_taxa) <- phyloseq::otu_table(prf_remain, taxa_are_rows = taxa_are_rows(ps_taxa))
    object <- ps_taxa
  } else if (inherits(object, "environment")) {
    object <- phyloseq::otu_table(prf_remain, taxa_are_rows = taxa_are_rows(object))
  } else {
    object <- prf_remain
  }

  return(object)
}

# the data is trimmed by threshold
trim_FeatureOrSample <- function(x, nRow, threshold) {

  df_occ <- apply(x, nRow, function(x) {
    length(x[c(which(!is.na(x) & x!=0))]) / length(x)
  }) %>%
    data.frame() %>% stats::setNames("Occ") %>%
    tibble::rownames_to_column("type")
  if(nRow == 1){
    rownames(df_occ) <- rownames(x)
  }else{
    rownames(df_occ) <- colnames(x)
  }
  df_KEEP <- apply(df_occ > threshold, 1, all) %>%
    data.frame() %>% stats::setNames("Status") %>%
    dplyr::filter(Status)

  return(df_KEEP)
}

# generate phyloseq
get_phyloseq <- function(y, x=metadata) {
  
  # x = metadata
  # y = profile_counts

  phenotype <- x %>%
    # dplyr::mutate(V1=gsub("\\.*", "", V1)) %>%
    dplyr::distinct() %>%
    tibble::column_to_rownames("V1")
  profile <- y %>%
    dplyr::select(-V1) %>%
    dplyr::rename(ID=TaxaID)
  # colnames(profile) <- gsub("\\.*", "", colnames(profile))
  
  sid <- intersect(rownames(phenotype),
                   colnames(profile))
  phen <- phenotype[rownames(phenotype) %in% sid, ]
  prof <- profile %>%
    dplyr::select(all_of(c("ID", rownames(phen))))
  
  dat_res <- import_metaphlan_taxa(data_metaphlan2 = prof)
  
  phy <- get_metaphlan_phyloseq(otu_tab = dat_res$abu_tab,
                         sam_tab = phen,
                         tax_tab = dat_res$tax_tab)
  
  res <- run_trim(object = phy, 
                  trim = "feature",
                  cutoff = 0)
  
  return(res)
}
  • Acquiring phyloseq object
## Counts 
phy_counts <- get_phyloseq(y = profile_counts)

## Relative abundance
phy_rb <- get_phyloseq(y = profile_rb)


if (!dir.exists("./phyloseq/")) {
  dir.create("./phyloseq/", recursive = T)
}

file_name_counts <- paste0("./phyloseq/curatedMetagenomic_", Sys.Date(), "_count.RDS")
saveRDS(phy_counts, file_name_counts, compress = TRUE)

file_name_rb <- paste0("./phyloseq/curatedMetagenomic_", Sys.Date(), "_relative.RDS")
saveRDS(phy_rb, file_name_rb, compress = TRUE)

Finally, taxonomic profile of 20,089 samples were downloaded from curatedMetagenomicData.

Summary

Download taxa profile of multiple samples curatedMetagenomicData using R.

  1. Extracting Studies’ information.
  2. Downloading results in one batch.
  3. Converting metadata and profile into phyloseq object.