{'fraction_non_alphanumeric': 0.03318333007239288, 'fraction_numerical': 0.009704558794756409, 'mean_word_length': 5.001878816345703, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CAMAMED: a pipeline for composition-aware mapping-based analysis of metagenomic data

Abstract Metagenomics is the study of genomic DNA recovered from a microbial community. Both assembly-based and mapping-based methods have been used to analyze metagenomic data. When appropriate gene catalogs are available, mapping-based methods are preferred over assembly based approaches, especially for analyzing the data at the functional level. In this study, we introduce CAMAMED as a composition-aware mapping-based metagenomic data analysis pipeline. This pipeline can analyze metagenomic samples at both taxonomic and functional profiling levels. Using this pipeline, metagenome sequences can be mapped to non-redundant gene catalogs and the gene frequency in the samples are obtained. Due to the highly compositional nature of metagenomic data, the cumulative sum-scaling method is used at both taxa and gene levels for compositional data analysis in our pipeline. Additionally, by mapping the genes to the KEGG database, annotations related to each gene can be extracted at different functional levels such as KEGG ortholog groups, enzyme commission numbers and reactions. Furthermore, the pipeline enables the user to identify potential biomarkers in case-control metagenomic samples by investigating functional differences. The source code for this software is available from https://github.com/mhnb/camamed. Also, the ready to use Docker images are available at https://hub.docker.com.


INTRODUCTION
Metagenomics is an interdisciplinary research field which lies in the intersection of molecular genomics, microbial ecology and data analysis. The main subject of this field is metagenome, which refers to as the total genomic content of microorganisms present in a certain environment. Metagenomics is based on microbial culture-independent methods, including high-throughput genome-sequencing techniques. In this way, DNA from the inhabiting microorganisms (i.e. the microbiome) are extracted from a certain environment, e.g. intestine, and are studied using various computational techniques (1).
There are two main approaches for analyzing metagenomic data: (i) taxonomic profiling, which describes the phylogenetic diversity of microorganisms in samples; and (ii) functional profiling, which includes computational strategies for mapping genomic sequences to 'functional' groups, to gain an insight on the functional capacities of microbiota in samples (2).
Comprehensive functional analysis of microbiome can significantly improve our understanding of biochemical capabilities of the microbial community (3). A possible strategy for functional profiling of microbiome is to assemble reads to larger contigs, and then, predict the gene content for functional profiling (4). For poorly unexplored microbiomes, assembly based strategies are inevitable, despite the fact that most conventional assemblies are of relatively low accuracy (5). In contrast, for those environments that have already been studied extensively, like the human gut, good microbiome gene catalogs have previously been compiled (6). In such cases, mapping-based analysis can be used for functional analysis of metagenomes.
In shotgun metagenomics studies, the collection of microorganisms is studied through direct DNA sequencing without any culture and isolation. By comparing the frequency of genes mapped to the gene catalog, functional dif-ferences between samples can be studied. Thus, gene abundance data are affected by high levels of systematic variability, which can greatly reduce statistical power and increase false positives (7). There are many reasons for the systematic variation in metagenome data that can affect the observed abundance of genes and microorganisms. One of the important reasons is the difference in the depth of sequencing so that each sample has a different number of sequencing reads (8). Other reasons for systematic variability include inconsistencies in sampling methods, DNA extraction, variation in the quality of sequencing runs, errors in read mapping and incomplete reference gene catalogs (9). In addition, the systematic variability can be due to differences in the average genome size of microorganisms, species richness and GC-content related to reads, which can affect the observed gene abundance (7,10).
Next-Generation Sequencing (NGS) data are also inherently compositional. Compositional means that the relative abundance of each nucleotide fragments is dependent to the abundance of other fragments. This property is related to the sequencing equipment and underlying methodology, and the resulted sequences are affected by the bias involved in amplification and subsequent nucleotide sampling (11). Hence, the composition is a result of this ambiguity in measurements that are an unclear part of the whole (e.g. metagenomic count data generated by NGS sequencing). The compositional data analysis (CoDA) refers to handling and resolving this bias (12).
Metagenomic count data also faces more severe challenges compared to the other NGS data. One of these challenges is the highly variable number of sequenced reads or sequencing depth in different samples. The second challenge is the very high percentage of zeros in metagenomic count data referred to as zero-inflation (roughly between 50 and 90%) (13). Also, metagenomic data are very high dimensional in comparison with the other NGS data (e.g. in a sample of the gut microbiome gene catalog, there are ∼10 M gene sequences or features (6)). On the other hand, due to the low frequency of DNA sampling, the very rare taxa are not recorded, which is called technical zeros. Also, some taxa may not be captured through their missing population, known as structural zeros. Another challenges are the size of the study and large variance in taxa distributions (overdispersion) (14).
Normalization processes can identify and eliminate systematic variability and compositional bias, so it is an essential step in data preprocessing and analysis. Many normalization methods have been proposed for high-dimensional count data, but for most of them, their performance has not been evaluated on metagenomic data (7). Various approaches have been proposed yet to address the challenges involving the compositional data. For example, to solve the problem of uneven sequencing depth, two approaches are introduced. First, rescaling the read counts in different samples to achieve a fixed value for their library size, and second, re-sampling reads to achieve the number of fixed reads for all samples (2). Also, many CoDA methods use transformation instead of normalization. These methods map the data to real space using log-ratio transformation, which makes it possible to use conventional statistical methods for further analysis. These methods try to use a reference such as the geometric mean of the subset feature for data transformation (15).
As explained, one of the major challenges in metagenomic compositional data is sparsity or zero-inflation, which becomes more acute for gene-centric count data. Several packages and tools have been developed to handle this problem and improve the accuracy of comparative gene abundance studies. Some packages have been developed specifically to deal with the zeros in the metagenomic data. The metagenomeSeq uses a zero-inflated log-normal model for gene abundance data (16,17). This method assumes that the zero-inflation is sample-specific and depends on the depth of sequencing (18). Ratio approach for identifying differential abundance (RAIDA) used a statistical model that first converts counts into relative frequencies, which are described by a log-normal distribution. RAIDA assumes that most features are not differentially abundant which makes it suitable to analyze metagenomic data at the taxa and gene levels. Also, this tool was developed to comparative analysis of metagenomic data samples in two different conditions, which can be generalized to more than two conditions (19). Also, there are several zero-inflated statistical models for metagenomic data, including zeroinflated negative-binomial and zero-inflated beta regression models (20,21).
However, other packages and tools for CoDA are provided, including ANCOM (22), ZIBSeq (21), CPL (23), DESeq2 (24) and edgeR (25). These methods exploit different statistical methods and try to handle the compositional bias and zero-inflation in the high-throughput sequencing data. It is important to note that most of the packages described above, including DESeq2 and edgeR, were originally developed for RNA-seq data analysis. Some of the most widely used CoDA packages and their properties are organized in Table 1.
In addition, a large number of normalization methods have been proposed for compositional data. Some of these methods are provided for RNA-seq data (26), some for operational taxonomic units (OTUs) data generated by amplicon sequencing (27) and some for metagenome data. But comparative studies of these methods show that there is a large dependency between performance and data characteristics. For example, the methods that have better performance for RNA-seq data will not necessarily be suitable for metagenome data (7). Some of these normalization methods are described below. Total sum scaling (TSS) is a standard normalization method for count data that is obtained by dividing the individual counts by the total counts in the sample such that the sum of the normalized values is 1 (28). TSS with a fixed scaling factor may harm OTU counts due to technological sequencing biases. Cumulative sum scaling (CSS) (16) re-scales samples based on the low-frequency (relatively constant and independent) quartiles, and does not eliminate the effect of high-frequency samples. Additionally, CSS as a highly cited log-normal model, has been used at taxa/gene level for metagenomic CoDA in many studies (2,18,29,30). Trimmed mean of M-values (TMM) estimates scale factors between samples for use in statistical analysis to identify differential expression. TMM normalization assumes that most genes are not differentially expressed between samples (31). Another family of meth- Metagenomic data (taxa and gene level) suitable for handling sparsity in very high dimensional data. It is also highly efficient for detecting rare samples in metagenomic data. (16) ods for normalizing compositional data is based on the logratio transformation. The mutual dependence between features in a composition implies that the analysis of individual features is performed to a reference baseline that transforms each sample into a new space and the statistical analysis will be done in this new space. Based on the choice of reference, different log-ratio based methods were developed. The centered log-ratio (CLR), additive log-ratio (ALR) and relative log expression (RLE) transformations use different strategies for selecting the reference (11,32). Table 1 details some of the packages, related normalization methods, properties and their advantages for CoDA. For functional analysis of metagenomic data, it is necessary to use a pipeline that considers the compositional nature of metagenomic data. In the present paper, we introduce CAMAMED, a mapping-based software pipeline to perform taxonomic and functional analysis of metagenomic data. Therefore, the proper normalization method is very important for metagenomic data analysis. Considering the properties previously described for this data (e.g. sparsity, high-dimensionality, rarefaction, undersampling and vast differences in sequencing depths, etc.), the metagenomeSeq and CSS as its normalization method is one of the most widely used methods in metagenomic studies (2,17,18,33). However, one of the most important drawbacks reported for metagenomeSeq in different studies is the high false-positive rate for differential abundance analysis when the sample size is small (34,35). CAMAMED uses the metagenomeSeq due to its high capability of handling compositional bias in sparse metagenomic data in taxa, gene, KO, EC number and reaction levels. CAMAMED is implemented using Python3 and Shell programming based on the Linux operating system. It is designed for nonprofessional users and is relatively easy to execute. Also, for easier use of CAMAMED, two Docker images have been constructed that enables users to run the CAMAMED pipeline without involving in installation details and dependencies. These images are available at www.hub.docker.com.

Metagenomic dataset
In this study, we used 80 metagenomic shotgun-sequenced fecal samples. This dataset (named originally 'cohort1') includes samples from 24 healthy (control), 27 colorectal adenoma and 29 colorectal carcinoma individuals. The samples had been sequenced using the Illumina platform and paired-end sequencing methods (36). Also, to obtain the abundance of genes in the samples, we used a previouslyreported gene catalog with 9.88 × 10 6 non-redundant gene sequences. This catalog, called the integrated gene catalog (IGC), contains the nucleotide and protein gene sequences of the human gut microbiome (6). For more information on the 'Materials and Methods', see Supplementary File 1.

Pipeline overview
CAMAMED is an automatic and easy-to-use computational pipeline for taxonomic and functional profiling of metagenomic data. Figure 1 shows the overall workflow of this pipeline. Pre-processes. CAMAMED can get the input data in FASTQ, FASTA or SRA file formats, in both paired-end and single-end modes. If the sequence format is SRA, it will be converted to FASTQ or FASTA format using the SRA Toolkit v2.8.2 (http://ncbi.github.io/sra-tools). Then, the quality control of the sequence dataset is performed using FastQC v0.11.5 (https://www.bioinformatics.babraham. ac.uk/projects/fastqc). The statistical characteristics of samples are extracted using SeqKit v0.10.1 (37).
Mapping processes. To find the bacterial frequency of the samples at different taxonomic levels, we use MetaPhlAn v2.0. MetaPhlAn (Metagenomic Phylogenetic Analysis) is a computational tool for profiling the composition of microbial communities from metagenomic shotgun sequencing data (38). CAMAMED uses a mapping-based strategy for metagenome data processing. To obtain the best mapping results, an appropriate gene catalog should be provided (see Ref. (6) as an example). If an in-house gene catalog is created, one can use CD-HIT v4.7 to remove potential duplicate sequences from the gene catalog. CD-HIT is a widely used program for clustering and comparing protein or nucleotide sequences (39). At this step, metagenomic read sequences are mapped to the gene catalog using MO-SAIK v2.2.3 software. This software uses the hash-based algorithm to map nucleotide sequences very quickly to the nucleotide gene catalog (40).

KEGG annotations.
After mapping reads to the gene catalog, one can compute the genes frequencies in each sample. For functional profiling, we use the KEGG database, which includes a comprehensive list of genes and genomes, together with their biochemical annotations (41). To extract the KEGG orthology (KO) groups associated with each of the genes, KEGG Automatic Annotation Server (KAAS, Ver. 2.1) (42) can be used. By using this tool, nucleotide and amino acid sequences can be mapped to KEGG and the associated KOs are retrieved. Alternatively, one may use GhostKOALA Ver. 2.2 (43), which maps amino acid sequences to KEGG. After extracting the KOs associated with each gene, for each KO, we obtain the enzyme commission (EC) numbers and reaction IDs. CAMAMED currently includes the KO-EC-reaction relations associated with KEGG. It is always possible for the user to retrieve the latest update of these data from KEGG, but this might be a time-consuming task. Then, for each sample, CAMAMED Data normalization. After extracting the frequencies of genes, these data should be corrected for the compositional bias. Before this step, we remove the genes to which less than five reads in all of the samples have been mapped (44). We then, divide the frequency of each gene to the length of the gene, in order to compute a normalized abundance of the genes with different length to calculate KO, EC number and reaction frequencies. Then, we use the metagenomeSeq package with CSS method for compositional correction and normalization (16). We also use this method to correct compositional bias in species-level in taxonomic data. Now, to calculate the normalized frequencies of KO, EC number, reaction, we use the sum of normalized frequency for their subset genes. In the following a brief explanation of the CSS method is provided. Suppose raw data is represented as a count matrix M (m, n) where m and n represent the number of features and samples, respectively. The raw data in the count matrix with c ij represents the number of times that feature i was observed in sample j. Also, the sum of counts for sample j is calculated as s j = i c i j . A usual method for normalizing feature value isc i j = c i j /s j , which is called total-sum normalization and has its own important drawbacks (improper handling of compositional bias). To avoid these drawbacks, CSS considers the lth quantile of sample j to be q l j , meaning that sample j has l features with counts less than q l j . For l = 0.95 * m . q l j represents the 95th percentile of the count distribution for sample j. Also, s l j = i |c i j ≤q l j c i j is the sum of feature counts up to lth quantile in sample j (16). CSS selectsl ≤ m to calculate the normalization scaling factor in each sample and defines normalized counts as c i j = (c i j / sˆl j ) (N), where N is the normalization constant that is selected equally for all samples. It is recommended to select N as the median of scaling factors sˆl j of all samples. The counts for samples with a scaling factor close to N can be considered as reference samples, and the counts for other samples can be calculated relative to the reference samples (16). Selecting an appropriate quantile based onl is critical to the normalization of data, and its value is project-specific and is chosen based on data driven methods that use experimental details such as sample preparation and sequencing (16,26). Also, CSS-normalized sample abundance can be well approximated with zero-inflated log-normal model in studies with a large number of samples. Therefore, logarith-mic transformation is used on normalized count data. This transformation controls the diversity of features measured in the samples (16).

Statistical analyses.
A possible application of metagenome analysis is when a comparative study is performed, to detect those 'biomarkers' which are underrepresented in case versus control samples. Using the non-parametric Kruskal-Wallis H test, one can identify the components that are changed significantly in different sample groups. Kruskal-Wallis test is used for comparing the value of a variable in two or more groups. The one-way analysis of variance (ANOVA) test can also be used in this pipeline assuming the data distribution is normal. We then use the Benjamini-Hochberg (BH) method to correct for false discovery rates.
By this stage, we have normalized data for the frequency of taxa, genes, KOs, EC numbers and reactions. We also specify the group labels to which each sample belongs. All the code and tools used in CAMAMED are collected in two Docker images that the user can easily run CAMAMED without having to worry about the installation details.

Test case: colorectal adenoma and carcinoma
To test the applicability of CAMAMED, we used 80 metagenomic shotgun-sequenced fecal samples obtained from 24 healthy (control), 27 colorectal adenoma and 29 colorectal carcinoma individuals (36). The results of the mapping analysis can depend on the reference gene catalog used. We used IGC, a previously-reported gene catalog of the human gut (6), to evaluate CAMAMED.
By applying CAMAMED for mapping the reads to the gene catalog, we found 3 354 281 genes to which at least one read was mapped. Table 2 shows the results of the Kruskal-Wallis test for levels of the species, gene, KO, EC number and reaction (for a significance level of P-value ≤ 0.01).
The results in Table 2 show that the ratio of significantly changed genes (that is, 1.92%) is not different from what is expected by chance (as p was assumed to be significant at the level of 0.01). This observation suggests that functional analysis cannot be ideally performed at the level of genes, in contrast to what has been previously proposed (45). In contrast, the ratio of significantly changed species is 2.67%, which explains why taxonomic profiling of microbiome data is widely used in the literature. However, one should note that 2.67% of 374 species means only 10 species, which might not provide us with enough number of features to be used in a successful classification of samples to control, adenoma and carcinoma groups. In case of other functional features, with the same level of significance (P ≤ 0.01) a greater ratio of KO groups, EC numbers and reactions might be detected as significant (2.98, 2.55 and 2.45%, respectively). Also, the last column of Table 2 shows the percentage of significantly changed entities after P-value adjustment, using the BH method associated with the fourth column. Therefore, we recommend the metagenome functional analysis to be performed at these levels, rather than the gene level (30). CAMAMED is currently the only available pipeline for this kind of analysis. Note that CAMAMED is an integrated pipeline of more than ten well-known bioinformatics tools that were previously used in other studies. To ensure the correctness of our implementation, we also applied CAMAMED on another datasets (30), and observed that all the results are reproducible.

DISCUSSION
Handling metagenomic data is a time-consuming and elaborate task. A number of pipelines are currently available for facilitating metagenomic analysis. Different aspects of the available metagenome analysis pipelines were compared in Table 3. These tools use assembly or mapping-based methods to process sequences. The level of annotation (taxon, gene, etc.) reported to the user is a vital aspect when tools are compared. The level of annotation in these tools is commonly taxon or gene, while CAMAMED annotates the samples at five levels, including taxon, gene, KO, EC number and reaction. This exclusive feature enables the user to analyze the samples at the functional level, which is reported to be more robust compared to taxonomic or genomic changes (30,46).
In recent years, in most of the metagenomic analyses, taxonomic profiles have been used as markers in case-control groups (49). In the functional analysis of metagenomes, on the other hand, case-control differences are studied at the gene level (50). Using the CAMAMED pipeline, not only one can easily analyze metagenome data at taxonomical (taxon) and functional (gene) level, but also it is possible to go further by analyzing the potential functional differences at other functional levels, that is, KO, EC number and reaction.
Methods for analyzing microbiome data sometimes assume, although implicitly, that sequencing data can be used equivalently in place of environmental data. However, microbiome sequencing data is often compositional and may not represent the original distribution of the samples in the environment (51). Due to the sparse nature of metagenomic data, different methods have been proposed for compositional bias correction. To this end, the CSS method is one of the most popular and powerful methods to handle this challenge (17,52). To demonstrate the composition-awareness of CSS and hence CAMAMED, Norouzi-Beirami et al. used two independent gut metagenome datasets on colorectal cancer (30). In this study, taxa and gene abundance data were extracted and compositional bias removal was performed independently on the data using the CSS method. The feature set extracted as a colorectal cancer marker from the first dataset has been accurately evaluated on the second dataset. However, without applying the CSS, the results from these datasets were not consistent (30). Furthermore, CAMAMED considers the correction for zeroinflation and compositional bias in the metagenomic data, which is, to the best of our knowledge, largely neglected in other pipelines.
The CAMAMED pipeline performs all the steps involved in the analysis of metagenomic data in a (semi-)automatic and step-by-step manner. Most of the tools used in this pipeline do not need to be installed separately by the user, which liberates the nonprofessional user from being engaged in potential software installation obstacles. It is necessary to emphasize that CAMAMED is a mappingbased pipeline for analyzing metagenomic data at the taxonomic and functional level. After running CAMAMED on metagenomic samples, normalized datasets are extracted at five levels of taxon, gene, KO, EC number and reaction. Such output can be then exploited by the user for additional machine learning and statistical studies.
Also, preparing the Docker images for CAMAMED make it possible that the user can employ it without involving in installation details and dependencies. These images also make using CAMAMED easier and increase the life of the software.

DATA AVAILABILITY
The software manual and the software package are available from https://github.com/mhnb/camamed and the Docker images with titled camamed pipeline and camamed pipeline db at www.hub.docker.com.