SUMMARY: gofasta comprises a set of command-line utilities for handling alignments of short assembled genomes in a genomic epidemiology context. It was developed for processing large numbers of closely related SARS-CoV-2 viral genomes and should be useful with other densely sampled pathogen genomic datasets. It provides functions to convert sam-format pairwise alignments between assembled genomes to fasta format; to annotate mutations in multiple sequence alignments, and to extract sets of sequences by genetic distance measures for use in outbreak investigations. AVAILABILITY AND IMPLEMENTATION: gofasta is an open-source project distributed under the MIT license. Binaries are available at https://github.com/virus-evolution/gofasta, from Bioconda, and through the Go programming language’s package management system. Source code and further documentation, including walkthroughs for common use cases, are available on the GitHub repository

Jackson, Ben

English

PubMed

PubMed Central

gofasta: command-line utilities for genomic epidemiology research

Summarygofasta comprises a set of command-line utilities for handling alignments of short assembled genomes in a genomic epidemiology context. It was developed for processing large numbers of closely related SARS-CoV-2 viral genomes, and should be useful with other densely sampled pathogen genomic datasets. It provides functions to convert sam-format pairwise alignments between assembled genomes to fasta format; to annotate mutations in multiple sequence alignments, and to extract sets of sequences by genetic distance measures for use in outbreak investigations.Availability and Implementationgofasta is an open-source project distributed under the MIT license. Binaries are available at https://github.com/virus-evolution/gofasta, from Bioconda, and through the Go programming language’s package management system. Source code and further documentation, including walkthroughs for common use cases, are available on the GitHub repository

Edinburgh Research Explorer

     Edinburgh Research Explorer                                      gofastaCitation for published version:Jackson, B & Borgwardt, K (ed.) 2022, 'gofasta: Command-line utilities for genomic epidemiology research',Bioinformatics. https://doi.org/10.1093/bioinformatics/btac424Digital Object Identifier (DOI):10.1093/bioinformatics/btac424Link:Link to publication record in Edinburgh Research ExplorerDocument Version:Publisher's PDF, also known as Version of recordPublished In:BioinformaticsGeneral rightsCopyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)and / or other copyright owners and it is a condition of accessing these publications that users recognise andabide by the legal requirements associated with these rights.Take down policyThe University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorercontent complies with UK legislation. If you believe that the public display of this file breaches copyright pleasecontact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately andinvestigate your claim.Download date: 24. Sep. 2022Genome analysisgofasta: command-line utilities for genomicepidemiology researchBen Jackson *Institute of Ecology and Evolution, University of Edinburgh, Edinburgh EH9 3FL, UK*To whom correspondence should be addressed.Associate Editor: Karsten BorgwardtReceived on January 27, 2022; revised on May 9, 2022; editorial decision on June 20, 2022; accepted on June 20, 2022AbstractSummary: gofasta comprises a set of command-line utilities for handling alignments of short assembled genomes in agenomic epidemiology context. It was developed for processing large numbers of closely related SARS-CoV-2 viralgenomes and should be useful with other densely sampled pathogen genomic datasets. It provides functions to convertsam-format pairwise alignments between assembled genomes to fasta format; to annotate mutations in multiple se-quence alignments, and to extract sets of sequences by genetic distance measures for use in outbreak investigations.Availability and implementation: gofasta is an open-source project distributed under the MIT license. Binaries areavailable at https://github.com/virus-evolution/gofasta, from Bioconda, and through the Go programming lan-guage’s package management system. Source code and further documentation, including walkthroughs for com-mon use cases, are available on the GitHub repository.Contact: ben.jackson@ed.ac.uk1. Introductiongofasta provides utilities that handle fasta-format alignments of viralconsensus genomes in a genomic epidemiology context. Its develop-ment was motivated by the scale and nature of SARS-CoV-2 genomedatasets, which are orders of magnitudes larger than any previouspathogen sequence dataset, which necessitates efficient and scalableanalysis methods (Hodcroft et al., 2021). gofasta is a command-linetool written in Go (https://go.dev) and is freely available under theMIT license as source code and in binary form from https://github.com/virus-evolution/gofasta, as well as in binary form from Bioconda(Grüning et al., 2018). It is designed for use on Unix-like (Linux andMac OSX) operating systems. Where appropriate, it defaults to read-ing input from standard in (stdin) and writing output to standard out(stdout), so that commands can be incorporated into shell pipelines; inorder to deal with large datasets, data are preferentially streamed, ra-ther than loaded into memory. gofasta incorporates bıogo (Kortschaket al., 2017) and makes use of the bitwise coding scheme for nucleoti-des by Emmanuel Paradis (Paradis et al., 2004). Its main functionalityis outlined below, including some brief use cases. Further documenta-tion is available at https://github.com/virus-evolution/gofasta.2. Functionality2.1 Sam to fasta format conversionPhylogenetic analysis, as is commonly performed in genomic epidemio-logical contexts (reviewed in Gardy and Loman, 2018), requires thatsequences be aligned to each other. The accepted gold standard ismultiple sequence alignment, but these methods are computationallyexpensive (reviewed in Edgar and Batzoglou, 2006), so are intractablefor very large numbers of sequences and/or when fast turnaround timesare required. When sequences are closely related, as is the case current-ly with SARS-CoV-2 for example, multiple pairwise alignments to acommon reference sequence offer an efficient approximation. Thistechnique is widely used in current analysis pipelines and software(Aksamentov et al., 2021; Moshiri, 2021; Moustafa and Planet, 2021;O’Toole et al., 2021a). gofasta provides functions to convert sequencealignment/map (sam) format files (Li et al., 2009) representing pairwisealignments between assembled virus genomes including those outputby minimap2 (Li, 2018) to fasta format (https://www.ncbi.nlm.nih.gov/genbank/fastaformat/). This is the pipeline that has been used byPangolin (O’Toole et al., 2021a) and the UK’s daily SARS-CoV-2processing pipelines (https://github.com/COG-UK) during the SARS-CoV-2 pandemic to date. As an example, the current Pangolin pipelinefor aligning SARS-CoV-2 sequences is equivalent to:minimap2 –a –x asm20 --score–N¼0 --sam–hit–only--secondary¼no reference.fasta unaligned.consen-sus.fasta –o aligned.samgofasta sam toMultiAlign –s aligned.sam --start266 --end 29674 --pad –o aligned.fastaIn this routine, the alignment is trimmed to set reference coordi-nates using the --start and --end flags and the trimmed regionsare replaced with Ns to preserve reference length with the --padVC The Author(s) 2022. Published by Oxford University Press. 1This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unre-stricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.Bioinformatics, 2022, 1–3https://doi.org/10.1093/bioinformatics/btac424Advance Access Publication Date: 5 July 2022Applications NoteDownloaded from https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btac424/6631223 by guest on 04 August 2022flag. These options may be useful, for example, to mask out untrans-lated regions in the case that they are more prone to errors due tolower sequencing coverage when using tiled amplicon schemes(Quick et al., 2017), or to restrict downstream analyses to specificportions of the genome.Equivalently, one can issue:minimap2 –t37 –a –x asm20 --score–N¼0 --sam–hit–only--secondary¼no reference.fasta unaligned.consen-sus.fasta j gofasta sam toma –t3 --start 266 --end 29674--pad > aligned.fastawhich has the advantage that the fasta file is written concurrent-ly and no intermediate sam file is written to disk, saving time andspace. Timings for file format conversion of 1 million SARS-CoV-2genomes are presented in Figure 1.In the routine above insertions relative to the reference are omit-ted from the converted pairwise alignments so that all sequences inthe fasta-format representation are the same (reference) length.gofasta also provides a routine (toPairAlign) that will convert asam file to pairwise alignments in fasta format (one file per se-quence) which include insertions relative to the reference and thereference sequence itself by default.2.2. Genetic distances and searching for nearestneighboursgofasta provides two routines which are designed to aid outbreakinvestigations by searching for neighbouring sequences by geneticdistance. Both are optimized for comparing a smaller set of querysequences against a larger set of target sequences: the focal querysequences are loaded into memory and the background dataset oftarget sequences is streamed from disk, so the target file can be arbi-trarily large.The first routine, closest, takes alignments in fasta format asinput and implements three standard distance measures for pairwisecomparisons, two of which are uncorrected for multiple hits (thenumber of nucleotide differences, and the number of nucleotidedifferences per site) and one of which is a corrected evolutionary dis-tance (TN93; Tamura and Nei, 1993). The user has the option to re-turn the closest n neighbours or all neighbours under a set distance dfor each query sequence. Catchments are either returned as lists oftarget neighbours, sorted by distance, with ties broken by genomecompleteness, or can be printed in long form with the distance be-tween each query–target neighbour pair written.The second routine, updown, implements an uncorrected dis-tance measure (the number of nucleotide differences). Sequencesmust be pre-aligned, then are reduced to the set of ATGC single nu-cleotide polymorphisms (SNPs) by which they differ by from a com-mon reference sequence, which is interpreted as being ancestral.Patterns of derived SNPs are used to infer the phylogenetic relation-ships between sequences (Fig. 2 and Fig. S1 in O’Toole et al.,2021b). For low-diversity datasets, this routine has the potential tobe faster than traditional distance measures, because it avoids thecomputational cost of comparing every alignment position. Findingthe 1000 closest genetic neighbours for a single query among a tar-get of 10.5 million SARS-CoV-2 sequences takes 2 min and 40 s on adesktop computer with Intel processors and a solid-state hard drive.The equivalent command using closest and TN93 distance takes8 min and 10 s. Peak memory usage for both routines was less than10 MB for these examples.2.3 Mutation callinggofasta provides a routine (snps) for describing all individual nu-cleotide differences between a set of aligned sequences and a singlereference sequence. It also provides routines to describe insertion,deletion, amino acid and nucleotide mutations from pairwise align-ments in sam format (sam variants) or a multiple sequence align-ment in fasta format (variants), given a reference sequence (or itsidentity in the multiple sequence alignment) and an annotation ofthe protein-coding content of the reference genome. Mutations areeither annotated per-genome or can be aggregated, to provide fre-quencies per input file. In the latter case, an optional threshold canbe applied to return only mutations that are present above a givenfrequency.As a brief example, it would be possible to find the frequenciesof all the amino acid changes at residue 681 in the Spike gene, andthe nucleotide changes underlying them, from the sample of SARS-CoV-2 sequences in aligned.fasta like:gofasta variants --msa aligned.fasta --annotationMN908947.gb --aggregate --append-snps j grep"^aa:S:P681"aa:S:P681H(nuc:C23604A),0.004000000aa:S:P681R(nuc:C23604G),0.983000000AcknowledgementsThank you to Áine O’Toole and Andrew Rambaut for conceptual input, andJT McCrone and three anonymous reviewers for helpful comments on themanuscript.FundingThis work was supported by funding from the Medical Research Council(MRC) part of UK Research & Innovation (UKRI), the National Institute ofHealth Research (NIHR) [MC_PC_19027] and Genome Research Limited,operating as the Wellcome Sanger Institute.Conflict of Interest: none declared.Data availability statementNo new data were generated or analysed in support of this research.Fig. 1. The additional time taken to convert a sam-format alignment of 1 millionSARS-CoV-2 genomes to fasta format. Times are given for reading and writingfrom/to a rotating hard disk drive (dot-dashed line), from/to a solid state drive (solidline) and the additional time above minimap2’s runtime needed to write the fastafile to a hard disk drive while reading the sam file from standard in (dashed line). Ineach case, the sam file was generated on a server with 40 logical CPUs by runningminimap2 with 32 threads and was 28 GB in size. gofasta was run with differentnumbers of threads (represented on the x-axis)2 B.JacksonDownloaded from https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btac424/6631223 by guest on 04 August 2022ReferencesAksamentov,I. et al. (2021) Nextclade: clade assignment, mutation calling andquality control for viral genomes. J. Open Source Softw., 6, 3773.Edgar,R.C. and Batzoglou,S. (2006) Multiple sequence alignment. Curr. Opin.Struct. Biol., 16, 368–373.Gardy,J.L. and Loman,N.J. (2018) Towards a genomics-informed, real-time,global pathogen surveillance system. Nat. Rev. Genet., 19, 9–20.Grüning,B. et al.; Bioconda Team. (2018) Bioconda: sustainable and comprehen-sive software distribution for the life sciences. Nat. Methods, 15, 475–476.Hodcroft,E.B. et al. (2021) Want to track pandemic variants faster? Fix thebioinformatics bottleneck. Nature, 591, 30–33.Kortschak,R.D. et al. (2017) Bıogo: a simple high-performance bioinformaticstoolkit for the go language. J. Open Source Softw., 2, 167.Li,H. et al.; 1000 Genome Project Data Processing Subgroup. (2009) The sequencealignment/map format and SAMtools. Bioinformatics, 25, 2078–2079.Li,H. (2018) Minimap2: pairwise alignment for nucleotide sequences.Bioinformatics, 34, 3094–3100.Moshiri,N. (2021) ViralMSA: massively scalable reference-guided multiple se-quence alignment of viral genomes. Bioinformatics, 37, 714–716.Moustafa,A.M. and Planet,P.J. (2021) Emerging SARS-CoV-2 diversityrevealed by rapid whole-genome sequence typing. Genome Biol. Evol., 13,evab197. https://doi.org/10.1093/gbe/evab197.O’Toole,Á. et al. (2021a) Assignment of epidemiological lineages in an emerg-ing pandemic using the pangolin tool. Virus Evol., 7, veab064.O’Toole,Á. et al. (2021b) Genomics-informed outbreak investigations ofSARS-CoV-2 using civet. medRxiv, https://doi.org/10.1101/2021.12.13.21267267.Paradis,E. et al. (2004) APE: analyses of phylogenetics and evolution in R lan-guage. Bioinformatics, 20, 289–290.Quick,J. et al. (2017) Multiplex PCR method for MinION and illuminasequencing of zika and other virus genomes directly from clinical samples.Nat. Protoc., 12, 1261–1276.Tamura,K. and Nei,M. (1993) Estimation of the number of nucleotide substi-tutions in the control region of mitochondrial DNA in humans and chim-panzees. Mol. Biol. Evol., 10, 512–526.gofasta 3Downloaded from https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btac424/6631223 by guest on 04 August 2022

gofasta:Command-line utilities for genomic epidemiology research

https://www.research.ed.ac.uk/files/289298103/btac424.pdf

gofasta:Command-line utilities for genomic epidemiology research

Abstract

Similar works

Full text

Available Versions

PubMed Central

Edinburgh Research Explorer