metascope_id: Identify which genomes are represented in a processed sample

In compbiomed/MetaScope: Tools and functions for preprocessing 16S and metagenomic sequencing microbiome data

Identify which genomes are represented in a processed sample

Description

This function will read in a .bam or .csv.gz file, annotate the taxonomy and genome names, reduce the mapping ambiguity using a mixture model, and output a .csv file with the results. Currently, it assumes that the genome library/.bam files use NCBI accession names for reference names (rnames in .bam file).

Usage

metascope_id(
  input_file,
  input_type = "csv.gz",
  aligner = "bowtie2",
  db = "ncbi",
  db_feature_table = NULL,
  accession_path = NULL,
  tmp_dir = dirname(input_file),
  out_dir = dirname(input_file),
  convEM = 1/10000,
  maxitsEM = 25,
  update_bam = FALSE,
  num_species_plot = NULL,
  out_fastas = FALSE,
  num_genomes = 100,
  num_reads = 50,
  quiet = TRUE
)

Arguments

input_file	The .bam or .csv.gz file of sample reads to be identified.
input_type	Extension of file input. Should be either "bam" or "csv.gz". Default is "csv.gz".
aligner	The aligner which was used to create the bam file. Default is "bowtie2" but can also be set to "subread" or "other".
db	Currently accepts one of c("ncbi", "silva", "other") Default is "ncbi", appropriate for samples aligned against indices compiled from NCBI whole genome databases. Alternatively, usage of an alternate database (like Greengenes2) should be specified with "other".
db_feature_table	If db = "other", a data.frame must be supplied with two columns, "Feature ID" matching the names of the alignment indices, and a second character column supplying the taxon identifying information.
accession_path	(character) Filepath to NCBI accessions SQL database. See taxonomzr::prepareDatabase().
tmp_dir	Path to a directory to which bam and updated bam files can be saved. Required.
out_dir	The directory to which the .csv output file will be output. Defaults to dirname(input_file).
convEM	The convergence parameter of the EM algorithm. Default set at 1/10000.
maxitsEM	The maximum number of EM iterations, regardless of whether the convEM is below the threshhold. Default set at 50. If set at 0, the algorithm skips the EM step and summarizes the .bam file 'as is'.
update_bam	Whether to update BAM file with new read assignments. Default is FALSE. If TRUE, requires input_type = "bam" such that a BAM file is the input to the function.
num_species_plot	The number of genome coverage plots to be saved. Default is NULL, which saves coverage plots for the ten most highly abundant species.
out_fastas	Logical, whether or not to output fasta files of reads. Default is FALSE.
num_genomes	Number of genomes to output fasta files for out_fastas. Default is 100.
num_reads	Number of reads per genome per fasta file for out_fastas. Default is 50.
quiet	Turns off most messages. Default is TRUE.

Value

This function returns a .csv file with annotated read counts to genomes with mapped reads. The function itself returns the output .csv file name. Depending on the parameters specified, can also output an updated BAM file, and fasta files for additional analysis downstream.

Examples

#### Align reads to reference library and then apply metascope_id()
## Assuming filtered bam files already exist

## Create temporary directory
file_temp <- tempfile()
dir.create(file_temp)

## Get temporary accessions database
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

#### Subread aligned bam file

## Create object with path to filtered subread csv.gz file
filt_file <- "subread_target.filtered.csv.gz"
bamPath <- system.file("extdata", filt_file, package = "MetaScope")
file.copy(bamPath, file_temp)

## Run metascope id with the aligner option set to subread
metascope_id(input_file = file.path(file_temp, filt_file),
             aligner = "subread", num_species_plot = 0,
             input_type = "csv.gz", accession_path = tmp_accession)

#### Bowtie 2 aligned .csv.gz file

## Create object with path to filtered bowtie2 bam file
bowtie_file <- "bowtie_target.filtered.csv.gz"
bamPath <- system.file("extdata", bowtie_file, package = "MetaScope")
file.copy(bamPath, file_temp)

## Run metascope id with the aligner option set to bowtie2
metascope_id(file.path(file_temp, bowtie_file), aligner = "bowtie2",
             num_species_plot = 0, input_type = "csv.gz",
             accession_path = tmp_accession)

## Remove temporary directory
unlink(file_temp, recursive = TRUE)

compbiomed/MetaScope documentation built on June 10, 2025, 11:52 a.m.