Microbiology, Metagenomics and Bioinformatics

Johan Bengtsson-Palme, University of Gothenburg | Wisconsin Institute for Discovery

Browsing Posts tagged Species classification

I am very happy to announce that a first public beta version of Metaxa2 version 2.2 has been released today! This new version brings two big and a number of small improvements to the Metaxa2 software (1). The first major addition is the introduction of the Metaxa2 Database Builder, which allows the user to create custom databases for virtually any genetic barcoding region. The second addition, which is related to the first, is that the classifier has been rewritten to have a more solid mathematical foundation. I have been promising that these updates were coming “soon” for one and a half years, but finally the end-product is good enough to see some real world testing. Bear in mind though that this is still a beta version that could contain obscure bugs. Here follows a list of new features (with further elaboration on a few below):

  • The Metaxa2 Database Builder
  • Support for additional barcoding genes, virtually any genetic region can now be used for taxonomic classification in Metaxa2
  • The Metaxa2 database repository, which can be accessed through the new metaxa2_install_database tool
  • Improved classification scoring model for better clarity and sensitivity
  • A bundled COI database for athropods, showing off the capabilities of the database builder
  • Support for compressed input files (gzip, zip, bzip, dsrc)
  • Support for auto-detection of database locations
  • Added output of probable taxonomic origin for sequences with reliability scores at each rank, made possible by the updated classifier
  • Added the -x option for running only the extraction without the classification step
  • Improved memory handling for very large rRNA datasets in the classifier (millions of sequences)
  • This update also fixes a bug in the metaxa2_rf tool that could cause bias in very skewed datasets with small numbers of taxa

The new version of Metaxa2 can be downloaded here, and for those interested I will spend the rest of this post outlining the Metaxa2 Database Builder. The information below is also available in a slightly extended version in the software manual.

The major enhancement in Metaxa2 version 2.2 is the ability to use custom databases for classification. This means that the user can now make their own database for their own barcoding region of choice, or download additional databases from the Metaxa2 Database Repository. The selection of other databases is made through the “-g” option already existing in Metaxa2. As part of these changes, we have also updated the classification scoring model for better stringency and sensitivity across multiple databases and different genes. The old scoring system can still be used by specifying the –scoring_model option to “old”.

There are two different main operating modes of the Metaxa2 Database Builder, as well as a hybrid mode combining the features of the two other modes. The divergent and conserved modes work in almost completely different ways and deal with two different types of barcoding regions. The divergent mode is designed to deal with barcoding regions that exhibit fairly large variation between taxa within the same taxonomic domain. Such regions include, e.g., the eukaryotic ITS region, or the trnL gene used for plant barcoding. In the other mode – the conserved mode – a highly conserved barcoding region is expected (at least within the different taxonomic domains). Genes that fall into this category would be, e.g., the 16S SSU rRNA, and the bacterial rpoB gene. This option would most likely also be suitable for barcoding within certain groups of e.g. plants, where similarity of the barcoding regions can be expected to be high. There is also a third mode – the hybrid mode – that incorporates features of both the other. The hybrid mode is more experimental in nature, but could be useful in situations where both the other modes perform poorer than desired.

In the divergent (default) mode, the database builder will start by clustering the input sequences at 20% identity using USEARCH (2). All clusters generated from this process are then individually aligned using MAFFT (3). Those alignments are split into two regions, which are used to build two hidden Markov models for each cluster of sequences. These models will be less precise, but more sensitive than those generated in the conserved mode. In the divergent mode, the database builder will attempt to extract full-length sequences from the input data, but this may be less successful than in the conserved mode.

In the conserved mode, on the other hand, the database builder will first extract the barcoding region from the input sequences using models built from a reference sequence provided (see above) and the Metaxa2 extractor (1). It will then align all the extracted sequences using MAFFT and determine the conservation of each position in the alignment. When the criteria for degree of conservation are met, all conserved regions are extracted individually and are then re-aligned separately using MAFFT. The re-aligned sequences are used to build hidden Markov models representing the conserved regions with HMMER (4). In this mode, the classification database will only consist of the extracted full-length sequences.

In the hybrid mode, finally, the database builder will cluster the input sequences at 20% identity using USEARCH, and then proceed with the conserved mode approach on each cluster separately .

The actual taxonomic classification in Metaxa2 is done using a sequence database. It was shown in the original Metaxa2 paper that replacing the built-in database with a generic non-processed one was detrimental to performance in terms of accuracy (1). In the database builder, we have tried to incorporate some of the aspects of the manual database curation we did for the built-in database that can be automated. By default, all these filtration steps are turned off, but enabling them might drastically increase the accuracy of classifications based on the database.

To assess the accuracy of the constructed database, the Metaxa2 Database Builder allows for testing the detection ability and classification accuracy of the constructed database. This is done by sub-dividing the database sequences into subsets and rebuilding the database using a smaller (by default 90%), randomly selected, set of the sequence data (5). The remaining sequences (10% by default) are then classified using Metaxa2 with the subset database. The number of detections, and the numbers of correctly or incorrectly classified entries are recorded and averaged over a number of iterations (10 by default). This allows for obtaining a picture of the lower end of the accuracy of the database. However, since the evaluation only uses a subset of all sequences included in the full database, the performance of the full database actually constructed is likely to be slightly better. The evaluation can be turned on using the “–evaluate T” option.

Metaxa2 2.2 also introduces the database repository, from which the user can download additional databases for Metaxa2. To download new databases from the repository, the metaxa2_install_database command is used. This is a simple piece of software but requires internet access to function.

References

  1. Bengtsson-Palme J, Hartmann M, Eriksson KM, Pal C, Thorell K, Larsson DGJ, Nilsson RH: Metaxa2: Improved Identification and Taxonomic Classification of Small and Large Subunit rRNA in Metagenomic Data. Molecular Ecology Resources (2015). doi: 10.1111/1755-0998.12399 [Paper link]
  2. Edgar RC: Search and clustering orders of magnitude faster than BLAST. Bioinformatics, 26, 2460–2461 (2010).
  3. Katoh K, Standley DM: MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular Biology and Evolution, 30, 772–780 (2013).
  4. Eddy SR: Accelerated profile HMM searches. PLoS Computational Biology, 7, e1002195 (2011).
  5. Richardson RT, Bengtsson-Palme J, Johnson RM: Evaluating and Optimizing the Performance of Software Commonly Used for the Taxonomic Classification of DNA Sequence Data. Molecular Ecology Resources, 17, 4, 760–769 (2017). doi: 10.1111/1755-0998.12628

Yesterday, Molecular Ecology Resources put online an unedited version of a recent paper which I co-authored. This time, Rodney Richardson at Ohio State University has made a tremendous work of evaluating three taxonomic classification software – the RDP Naïve Bayesian Classifier, RTAX and UTAX – on a set of DNA barcoding regions commonly used for plants, namely the ITS2, and the matK, rbcL, trnL and trnH genes.

In the paper (1), we discuss the results, merits and limitations of the classifiers. In brief, we found that:

  • There is a considerable trade-off between accuracy and sensitivity for the classifiers tested, which indicates a need for improved sequence classification tools (2)
  • UTAX was superior with respect to error rate, but was exceedingly stringent and thus suffered from a low assignment rate
  • The RDP Naïve Bayesian Classifier displayed high sensitivity and low error at the family and order levels, but had a genus-level error rate of 9.6 percent
  • RTAX showed high sensitivity at all taxonomic ranks, but at the same time consistently produced the high error rates
  • The choice of locus has significant effects on the classification sensitivity of all tested tools
  • All classifiers showed strong relationships between database completeness, classification sensitivity and classification accuracy

We believe that the methods of comparison we have used are simple and robust, and thereby provides a methodological and conceptual foundation for future software evaluations. On a personal note, I will thoroughly enjoy working with Rodney and Reed again; I had a great time discussing the ins and outs of taxonomic classification with them! The paper can be found here.

References and notes

  1. Richardson RT, Bengtsson-Palme J, Johnson RM: Evaluating and Optimizing the Performance of Software Commonly Used for the Taxonomic Classification of DNA Sequence Data. Molecular Ecology Resources, Early view (2016). doi: 10.1111/1755-0998.12628 [Paper link]
  2. This is something that several classifiers also showed in the evaluation we did for the Metaxa2 paper (3). Interestingly enough, Metaxa2 is better at maintaining high accuracy also as sensitivity is increased.
  3. Bengtsson-Palme J, Hartmann M, Eriksson KM, Pal C, Thorell K, Larsson DGJ, Nilsson RH: Metaxa2: Improved identification and taxonomic classification of small and large subunit rRNA in metagenomic data. Molecular Ecology Resources, 15, 6, 1403–1414 (2015). doi: 10.1111/1755-0998.12399 [Paper link]

I am very happy to announce that Metaxa2 version 2.1 has been released today. This new version brings a lot of important improvements to the Metaxa2 software (1), in particular by the introduction of the Metaxa2 Diversity Tools. This is the list of new features (further elaboration follows below):

  • The Metaxa2 Diversity Tools:
    • metaxa2_dc – a tool for collecting several .taxonomy.txt output files into one large abundance matrix, suitable for analysis in, e.g., R
    • metaxa2_rf – generates rarefaction curves based on the .taxonomy.txt output
    • metaxa2_si – species inference based on guessing species data from the other species present in the .taxonomy.txt output file
    • metaxa2_uc – a tool for determining if the community composition of a sample is significantly different from others through resampling analysis
  • Added a new detection mode for detection of multiple rRNA in the same sequence, e.g. a genome
  • Added the --reference option to improve the use of Metaxa2 as a tool to sort out host rRNA sequences from a dataset
  • Added the --split_pairs option causing Metaxa2 to output paired-end sequences into two separate files, which is nice for further analysis of rRNA reads
  • The default setting for the --align option has been changed to ‘none
  • Automatic detection of which BLAST package that is installed
  • Fixed a bug causing the last read of paired-end FASTA input to be ignored
  • Fixed an occasionally occurring BLAST+ related warning message
  • Fixed a bug that could cause the classifier to crash on highly divergent BLAST matches

The new version of Metaxa2 can be downloaded here, and for those interested I will spend the rest of this post outlining the new features.

Metaxa2 Diversity Tools
One often requested feature of Metaxa2 is the ability to further make simple analysis from the data after classification. The Metaxa2 Diversity Tools included in Metaxa2 2.1 is a seed for such an effort (although not close to a full-fledge community analysis package compared to QIIME (2) or Mothur (3)). The set currently consist of four tools

The Metaxa2 Data Collector (metaxa2_dc) is the simplest of them (but probably the most requested), designed to merge the output of several *.level_X.txt files from the Metaxa2 Taxonomic Traversal Tool into one large abundance matrix, suitable for further analysis in, for example, R. The Metaxa2 Species Inference tool (metaxa2_si) can be used to further infer taxon information on, for example, the species level at a lower reliability than what would be permitted by the Metaxa2 classifier, using a complementary algorithm. The idea is that is if only a single species is present in, e.g., a family and a read is assigned to this family, but not classified to the species level, that sequence will be inferred to the same species as the other reads, given that it has more than 97% sequence identity to its best reference match. This can be useful if the user really needs species or genus classifications but many organisms in the studied species group have similar rRNA sequences, making it hard for the Metaxa2 classifier to classify sequences to the species level.

The Metaxa2 Rarefaction analysis tool (metaxa2_rf) performs a rarefaction analysis based on the output from the Metaxa2 classifier, taking into account also the unclassified portion of rRNAs. The Metaxa2 Uniqueness of Community analyzer (metaxa2_uc), finally, allows analysis of whether the community composition of two or more samples or groups is significantly different. Using resampling of the community data, the null hypothesis that the taxonomic content of two communities is drawn from the same set of taxa (given certain abundances) is tested. All these tools are further described in the manual.

The genome mode
Metaxa2 has long been said not to be useful for predicting rRNA in longer sequences, such as full genomes or chromosomes, since it has traditionally only looked for a single rRNA hit. With Metaxa2 2.1, it is now possible to use Metaxa2 on longer sequences to detect multiple rRNA occurrences. To do this, you need to change the operating mode using the new --mode option to either ‘auto‘ or ‘genome‘. The auto mode will treat sequences longer than 2500 bp as “genome” sequences and look for multiple matches in these.

The reference mode
Another feature request that has been addressed in the new Metaxa2 version is the ability to filter out certain sequences from the data set. For example, you may want to exclude all rRNA sequences that are derived from to host organism, but keep the analysis of all other rRNA reads. This is now possible by supplying a file of reference rRNA sequences to exclude in FASTA format to the --reference option.

Experimental Usearch support
Finally, we have toyed around with support for Usearch (4) instead of BLAST (5) as the search algorithm for the classification step. However, this is far from fine-tuned and it is included as an experimental feature that you may use on your own risk! We recommend that you not use it for classification of data for publication yet. However, we are interested in how this works for you, so if you like you may test to run the Usearch algorithm in parallel with your BLAST-based analysis and compare the results and send me your input on how it works. You can read more about using Usearch at the end of the Metaxa2 manual.

References

  1. Bengtsson-Palme J, Hartmann M, Eriksson KM, Pal C, Thorell K, Larsson DGJ, Nilsson RH: Metaxa2: Improved Identification and Taxonomic Classification of Small and Large Subunit rRNA in Metagenomic Data. Molecular Ecology Resources (2015). doi: 10.1111/1755-0998.12399 [Paper link]
  2. Caporaso JG, Kuczynski J, Stombaugh J et al.: QIIME allows analysis of high-throughput community sequencing data. Nature Methods, 7, 335–336 (2010).
  3. Schloss PD, Westcott SL, Ryabin T et al.: Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Applied and Environmental Microbiology, 75, 7537–7541 (2009).
  4. Edgar RC: Search and clustering orders of magnitude faster than BLAST. Bioinformatics, 26, 2460–2461 (2010).
  5. Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ: Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res, 25, 3389–3402 (1997).

Metaxa2 has been updated to version 2.0.2 and can be downloaded from the Metaxa2 web site. The 2.0.2 update fixes two minor bugs; one causing the “.graph” file to display incorrect or no names for the regions of the LSU regions, and one causing misreporting of the number of sequences in single-end FASTQ files (paired-end files were reported correctly). The update also brings a slightly improved classifier. Thanks to Marco Severgnini for reporting the FASTQ file issue! The update is available here.

After almost a year in different stages of review and revision, in which the paper (but not the software) saw a total transformation, I am happy to announce that the paper describing Metaxa2 has been accepted in Molecular Ecology Resources and is available in a rudimentary online early form. The figures in this version are not that pretty, but those who wants to read the paper asap, you have the possibility to do so.

This means that if you have been using Metaxa2 for a publication, there is now a new preferred way of citing this, namely:

Bengtsson-Palme J, Hartmann M, Eriksson KM, Pal C, Thorell K, Larsson DGJ, Nilsson RH: Metaxa2: Improved Identification and Taxonomic Classification of Small and Large Subunit rRNA in Metagenomic Data. Molecular Ecology Resources (2015). doi: 10.1111/1755-0998.12399

The paper (1), apart from describing the new Metaxa version, also brings a very thorough evaluation of the software, compared to other tools for taxonomic classification implemented in QIIME (2). In short, we show that:

  • Metaxa2 can make trustworthy taxonomic classifications even with reads as short as 100 bp
  • Generally, the performance is reliable across the entire SSU rRNA gene, regardless of which V-region a read is derived from
  • Metaxa2 can reliably recapture species composition from short-read metagenomic data, comparable with results of amplicon sequencing
  • Metaxa2 outperforms other popular tools such as Mothur (3), the RDP Classifier (4), Rtax (5) and the QIIME implementation of Uclust (6) in terms of proportion of correctly classified reads from metagenomic data
  • The false positive rate of Metaxa2 is very close to zero; far superior to many of the above mentioned tools, many of which assume that reads must derive from the rRNA gene

Metaxa2 can be downloaded here. We have already used it for around two years internally, and it forms the base of the taxonomic classifications in e.g. our recently published paper on antibiotic resistance in a polluted Indian lake (7).

References

  1. Bengtsson-Palme J, Hartmann M, Eriksson KM, Pal C, Thorell K, Larsson DGJ, Nilsson RH: Metaxa2: Improved Identification and Taxonomic Classification of Small and Large Subunit rRNA in Metagenomic Data. Molecular Ecology Resources (2015). doi: 10.1111/1755-0998.12399 [Paper link]
  2. Caporaso JG, Kuczynski J, Stombaugh J et al.: QIIME allows analysis of high-throughput community sequencing data. Nature Methods, 7, 335–336 (2010).
  3. Schloss PD, Westcott SL, Ryabin T et al.: Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Applied and Environmental Microbiology, 75, 7537–7541 (2009).
  4. Wang Q, Garrity GM, Tiedje JM, Cole JR: Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Applied and Environmental Microbiology, 73, 5261–5267 (2007).
  5. Soergel DAW, Dey N, Knight R, Brenner SE: Selection of primers for optimal taxonomic classification of environmental 16S rRNA gene sequences. The ISME Journal, 6, 1440–1444 (2012).
  6. Edgar RC: Search and clustering orders of magnitude faster than BLAST. Bioinformatics, 26, 2460–2461 (2010).
  7. Bengtsson-Palme J, Boulund F, Fick J, Kristiansson E, Larsson DGJ: Shotgun metagenomics reveals a wide array of antibiotic resistance genes and mobile elements in a polluted lake in India. Frontiers in Microbiology, 5, 648 (2014).

My colleague Henrik Nilsson has been interviewed by the ResearchGate news team about the recent effort to better annotate ITS data for plant pathogenic fungi. It’s an interesting read, and I think Henrik nicely underscores why large-scale efforts for improving and correcting sequence annotations are important. You can read the interview here, and the paper they talk about is referenced below.

Nilsson RH, Hyde KD, Pawlowska J, Ryberg M, Tedersoo L, Aas AB, Alias SA, Alves A, Anderson CL, Antonelli A, Arnold AE, Bahnmann B, Bahram M, Bengtsson-Palme J, Berlin A, Branco S, Chomnunti P, Dissanayake A, Drenkhan R, Friberg H, Frøslev TG, Halwachs B, Hartmann M, Henricot B, Jayawardena R, Jumpponen A, Kauserud H, Koskela S, Kulik T, Liimatainen K, Lindahl B, Lindner D, Liu J-K, Maharachchikumbura S, Manamgoda D, Martinsson S, Neves MA, Niskanen T, Nylinder S, Pereira OL, Pinho DB, Porter TM, Queloz V, Riit T, Sanchez-García M, de Sousa F, Stefaczyk E, Tadych M, Takamatsu S, Tian Q, Udayanga D, Unterseher M, Wang Z, Wikee S, Yan J, Larsson E, Larsson K-H, Kõljalg U, Abarenkov K: Improving ITS sequence data for identification of plant pathogenic fungi. Fungal Diversity, Volume 67, Issue 1 (2014), 11–19. doi: 10.1007/s13225-014-0291-8 [Paper link]

Another paper I have co-authored related to the UNITE database for fungal rDNA ITS sequences is now published as an Online Early article in Fungal Diversity. The paper describes an effort to improve the annotation of ITS sequences from fungal plant pathogens. Why is this important? Well, apart from fungal plant pathogens being responsible for great economic losses in agriculture, the paper is also conceptually important as it shows that together we can accomplish a substantial improvement to the metadata in sequence databases. In this work we have hunted down high-quality reference sequences for various plant pathogenic fungi, and re-annotated incorrectly or insufficiently annotated ITS sequences from the same fungal lineages. In total, the 59 authors have made 31,954 changes to UNITE database data, on average 540 changes per author. While one, or a few, persons could not feasibly have made this effort alone, this work shows that in larger consortia vast improvements can be made to the quality of databases, by distributing the work among many scientists. In many ways, this relates to proposals to “wikify” GenBank, and after Rfam and Pfam it might now be time to take the user-contribution model to, at least, the RefSeq portion of GenBank, which despite its description as being “comprehensive, integrated, non-redundant, [and] well-annotated” still contains errors and examples of non-usable annotation. More on that at a later point…

Paper reference:

Nilsson RH, Hyde KD, Pawlowska J, Ryberg M, Tedersoo L, Aas AB, Alias SA, Alves A, Anderson CL, Antonelli A, Arnold AE, Bahnmann B, Bahram M, Bengtsson-Palme J, Berlin A, Branco S, Chomnunti P, Dissanayake A, Drenkhan R, Friberg H, Frøslev TG, Halwachs B, Hartmann M, Henricot B, Jayawardena R, Jumpponen A, Kauserud H, Koskela S, Kulik T, Liimatainen K, Lindahl B, Lindner D, Liu J-K, Maharachchikumbura S, Manamgoda D, Martinsson S, Neves MA, Niskanen T, Nylinder S, Pereira OL, Pinho DB, Porter TM, Queloz V, Riit T, Sanchez-García M, de Sousa F, Stefaczyk E, Tadych M, Takamatsu S, Tian Q, Udayanga D, Unterseher M, Wang Z, Wikee S, Yan J, Larsson E, Larsson K-H, Kõljalg U, Abarenkov K: Improving ITS sequence data for identification of plant pathogenic fungi. Fungal Diversity Online early (2014). doi: 10.1007/s13225-014-0291-8 [Paper link]

Our paper on the most recent developments of the UNITE database for fungal rDNA ITS sequences has just been published as an Early view article in Molecular Ecology. In this paper, we aim to ease two of the major problems facing the identification of newly generated fungal ITS sequences: the lack of a sufficiently goof reference dataset, and the lack of a way to refer to fungal species without a Latin name. As part of a solution, we have introduced the term species hypothesis for all fungal species represented by at least two ITS sequences. The UNITE database has an easy-to-use web-based sequence management system, and we encourage everybody that can improve on the annotations or metadata of a fungal lineage to do so.

My main contribution on this paper has been to tailor ITSx functionality for the UNITE database, so that ITS data could be more easily processed for the Species Hypotheses.

Paper reference:
Kõljalg U, Nilsson RH, Abarenkov K, Tedersoo L, Taylor AFS, Bahram M, Bates ST, Bruns TT, Bengtsson-Palme J, Callaghan TM, Douglas B, Drenkhan T, Eberhardt U, Dueñas M, Grebenc T, Griffith GW, Hartmann M, Kirk PM, Kohout P, Larsson E, Lindahl BD, Lücking R, Martín MP, Matheny PB, Nguyen NH, Niskanen T, Oja J, Peay KG, Peintner U, Peterson M, Põldmaa K, Saag L, Saar I, Schüßler A, Senés C, Smith ME, Suija A, Taylor DE, Telleria MT, Weiß M, Larsson KH: Towards a unified paradigm for sequence-based identification of Fungi. Accepted in Molecular Ecology. doi: 10.1111/mec.12481 [Paper link]

One potential use for Metaxa (paper) is to include it in a pipeline for classification of SSU rRNA in metagenomic data (or other environmental sequencing sets). However, as Metaxa is provided from this site, it only classifies SSUs to the domain level (archaea, bacteria and eukaryotes, with the addition of chloroplasts and mitochondria). It is also able to do some (pretty rough) species guesses using the “--guess_species T” option. An easy solution to implement would be to pass the Metaxa output, e.g. “metaxa_output.bacteria.fasta” to BLAST, and compare all these sequences to the sequences in e.g. the SILVA or GreenGenes database. There is, however, a way to improve this, which uses Metaxa’s ability to compares sequences to custom databases. In this tutorial, I will show you how to achieve this.

Before we start, you will of course need to download and install Metaxa, and its required software packages (BLAST, HMMER, MAFFT). When you have done this, we can get going with the database customization. I will in this tutorial use the SILVA database for SSU classification. However, the basic idea for the tutorial should be easily applicable to GreenGenes and other rRNA databases as well.

  1. Visit SILVA through this link, and download the file named “SSURef_106_tax_silva.fasta.tgz”. The file is pretty big so it may take a while to download it. If you’re running Metaxa on a server, you’ll have to get the SILVA-file to the server somehow.
  2. Unzip and untar the file (Mac OS X makes this neatly by doubleclicking the file, on linux you can do it on the command line by typing “tar -xvzf SSURef_106_tax_silva.fasta.tgz“). This will give you a FASTA-file.
  3. The FASTA-file needs to be prepared a bit for Metaxa usage. First, we need to give Metaxa identifiers it can understand. Metaxa identifies sequences’ origins by the last character in their identifier, e.g. “>A16379.1.1496.B”. Here, “.B” indicates that this is a bacterial sequence. We are now going to use the unix command sed to process the file and insert the appropriate identifiers.
    1. We begin with the archaeal sequences. To get those straight, we type:
      sed "s/ Archaea;/.A - Archaea;/" SSURef_106_tax_silva.fasta > temp1
      Notice that we direct the output to a temporary file. It is bad practice to replace the input file with the output file, so we work with two temp-files instead.
    2. The next step is also easy, now we find all eukaryote sequences and add E:s to the identifiers:
      sed "s/ Eukaryota;/.E - Eukaryota;/" temp1 > temp2
    3. Now it becomes a little more complicated, as SILVA classes mitochondrial and chloroplast SSU sequences as subclasses of bacteria. However, there is a neat little trick we can use. First we do the same with the bacterial sequences as with the archaeal and eukaryote:
      sed "s/ Bacteria;/.B - Bacteria;/" temp2 > temp1
    4. Now, we can use two a little more complicated commands to annotate the mitochondrial and chloroplast sequences:
      sed "s/\.B - \(Bacteria;.*;[Mm]itochondria;\)/.M - \1/" temp1 > temp2
      sed "s/\.B - \(Bacteria;.*;[Cc]hloroplast;\)/.C - \1/" temp2 > temp1
    5. We also need to get “rid” of the unclassified sequences, by assigning them to the “other” origin (O):
      sed "s/ Unclassified;/.O - Unclassified;/" temp1 > temp2
  4. That wasn’t too complicated, was it? We can now check the number of different sequences in the file by typing the pretty complicated command:
    grep ">" temp2 | cut -f 1 -d " " | rev | cut -f 1 -d "." | sort | uniq -c
    If you have been working with the same files as me, you should now see the following numbers:
    23172 A
    471949 B
    3712 C
    55937 E
    534 M
    226 O
  5. At this stage, we need to remove the full taxonomy from the FASTA headers, as Metaxa cannot handle species names of this length. We do this by typing:
    sed "s/ - .*;/ - /" temp2 > temp1
  6. We can now change the temp-file into a FASTA file, and delete the other temp-file:
    mv temp1 SSURef.fasta
    rm temp2
  7. We now need to configure Metaxa to use the database. First, we format a BLAST-database from the FASTA-file we just created:
    formatdb -i SSURef.fasta -t "SSURef Metaxa DB" -o T -p F
  8. With that done, we can now run Metaxa using this database instead of the classification database that comes with the program. By specifying that we want to guess the species origin of sequences, we can get (as accurate as SILVA lets us be) which species each sequence in our set come from. We do this by using the -d and the --guess_species options:
    metaxa -i test.fasta -d SSURef.fasta -o TEST --guess_species T --cpu 2
    The input in this case was the test file that comes with Metaxa. Note also that we’re using two CPUs to get multithreaded speeds. Remember that you must provide the full (or relative) path to the database files we just created, if you are not running Metaxa from the same directory as the database resides in.
  9. The output should now look like this (taken from the bacterial file):
    >coryGlut_Bielefeld_dna Bacterial 16S SSU rRNA, best species guess: Corynebacterium glutamicum
    CGAACGCTG...
    >gi|116668568:792344-793860 Bacterial 16S SSU rRNA, best species guess: Arthrobacter sp. J3.40
    TGAACGCTG...
    >gi|117927211:c1399163-1397655 Bacterial 16S SSU rRNA, best species guess: Acidothermus cellulolyticus
    >CGAACGCTG...

    And so on. As you can see the species names are now located at the end of each definition line, and can easily be extracted using sed, e.g. “grep ">" TEST.bacteria.fasta | sed "s/.*: //"“.

And that’s it. It’s pretty simple, and can easily be scripted. In fact, I have already made the bash script for you. That means that the short version is, download the script, download the sequence file from SILVA, move into the directory you have downloaded the file to and run the script by typing: ./prepare_silva_for_metaxa.sh

A few notes at the end. The benefit of using this approach is that we maintain the sorting capabilities, marking of uncertain sequences and error checking of Metaxa, but we don’t have to add another BLAST step after Metaxa has finished. However, as this database we create is a lot bigger than the database that comes with Metaxa, the running time of the classification step will be substantially longer. This is in most cases acceptable, as that time is the same as the time it would have taken to run BLAST on the Metaxa output. It should also be noted that this approach limits Metaxa’s ability of classifying 12S sequences, as there are no such sequences in SILVA. Good luck with classifying your metagenome SSUs (and if you use Metaxa in your research, remember to cite the paper)!