Antimicrobial resistance (AMR) poses a threat to public health. Clinical microbiology laboratories typically rely on culturing bacteria for antimicrobial-susceptibility testing (AST). As the implementation costs and technical barriers fall, whole-genome sequencing (WGS) has emerged as a ‘one-stop’ test for epidemiological and predictive AST results. Few published comparisons exist for the myriad analytical pipelines used for predicting AMR. To address this, we performed an inter-laboratory study providing sets of participating researchers with identical short-read WGS data from clinical isolates, allowing us to assess the reproducibility of the bioinformatic prediction of AMR between participants, and identify problem cases and factors that lead to discordant results. We produced ten WGS datasets of varying quality from cultured carbapenem-resistant organisms obtained from clinical samples sequenced on either an Illumina NextSeq or HiSeq instrument. Nine participating teams (‘participants’) were provided these sequence data without any other contextual information. Each participant used their choice of pipeline to determine the species, the presence of resistance-associated genes, and to predict susceptibility or resistance to amikacin, gentamicin, ciprofloxacin and cefotaxime. We found participants predicted different numbers of AMR-associated genes and different gene variants from the same clinical samples. The quality of the sequence data, choice of bioinformatic pipeline and interpretation of the results all contributed to discordance between participants. Although much of the inaccurate gene variant annotation did not affect genotypic resistance predictions, we observed low specificity when compared to phenotypic AST results, but this improved in samples with higher read depths. Had the results been used to predict AST and guide treatment, a different antibiotic would have been recommended for each isolate by at least one participant. These challenges, at the final analytical stage of using WGS to predict AMR, suggest the need for refinements when using this technology in clinical settings. Comprehensive public resistance sequence databases, full recommendations on sequence data quality and standardization in the comparisons between genotype and resistance phenotypes will all play a fundamental role in the successful implementation of AST prediction using WGS in clinical microbiology laboratories.

Background Human microbiome research is interdisciplinary, making concise organizing and reporting of results across the different styles of epidemiology, biology, bioinformatics, and statistics a challenge. Commonly used reporting guidelines for observational or genetic studies lack key aspects specific to microbiome studies. Methods A multidisciplinary group of microbiome researchers reviewed elements of available reporting guidelines for observational and genetic studies, and adapted these for application to human microbiome studies. New reporting elements were developed for laboratory, bioinformatic, and statistical analysis specific to microbiome studies, and other parts of these checklists were streamlined to keep reporting manageable. Results STORMS is a 18-item checklist for reporting on human microbiome studies, organized into six sections covering all sections of a scientific publication, presented as a table with space for author-provided details and intended for inclusion in supplementary materials. Conclusions STORMS provides guidance for authors and standardization for interdisciplinary microbiome studies, facilitating complete and concise reporting. Availability STORMS is downloadable as a versioned spreadsheet from [<a href="http://storms.waldronlab.io">storms.waldronlab.io</a>][1]. ### Competing Interest Statement The authors have declared no competing interest. [1]: <http://storms.waldronlab.io>

A relevant parameter when looking at sequencing and alignment quality is the edit distance to the reference genome, roughly equivalent to the&nbsp;accuracy of the reads. (Roughly, because we ignore true variants between the sample and reference).&nbsp;The edit distance or Levenshtein distance can be defined as the number of single letter (nucleotide) changes that have to be made to one string (read) for it to be equal to another string (reference genome). Since the error profile of Nanopore sequencing is dominated by insertions and deletions, the edit distance isn&rsquo;t just the number of single nucleotide mismatches. This technology results in a wide range of read lengths, therefore I scale the edit distance to the length of the aligned fragment. Longer reads shouldn&rsquo;t get penalised more than shorter reads. It&rsquo;s important to take the alignment length and not the read length since the ends of reads can be clipped substantially by the aligner, sometimes tens of kb. &nbsp;The experiments below were done using MinION data from the human genome in a bottle sample NA12878. My scripts are in Python, so I&rsquo;ll add some code snippets to this post. For parsing bam files and extracting the relevant bits of information, I use pysam&nbsp;which is pretty convenient and well-documented. Those snippets are not the full script, but just a minimal example to get the job done. A full example of the code is at the bottom of the post. Getting this information from an alignment done using bwa mem&nbsp; is trivial since bwa sets the NM bam tag which is an integer with precisely what I&nbsp;need: the edit distance. A function to extract the NM tag using a list comprehension is added below. .gist table { margin-bottom: 0; } import pysam def extractNMFromBam(bam): ''' loop over a bam file and get the edit distance to the reference genome stored in the NM tag scale by aligned read length ''' samfile = pysam.AlignmentFile(bam, "rb") return [read.get_tag("NM")/read.query_alignment_length for read in samfile.fetch()] view rawgetNMtag.py hosted with &#10084; by GitHub For aligners such as GraphMap&nbsp;this is less trivial since the NM tag is not set. However, another tag comes to the rescue: the MD tag, which stores a string containing matching numbers of nucleotides and non-matching nucleotides. Interesting information about the MD tag can be found here. I found it a quite tough representation of the read to wrap my head around, which resulted in some wrong interpretations. My first naive implementation counted the number of matching nucleotides (the integers in the MD string) and subtracted those matches from the total alignment length to get the number of mismatched nucleotides. The list comprehension is quite long, but essentially I split the MD string on all occurrences of A, C, T, G or ^ (indicating a deletion) and sum the obtained integers. This sum is subtracted from the aligned read length and divided by the same. .gist table { margin-bottom: 0; } import pysam import re def extractMDFromBam(bam): ''' loop over a bam file and get the edit distance to the reference genome mismatches are stored in the MD tag scale by aligned read length ''' samfile = pysam.AlignmentFile(bam, "rb") return [(read.query_alignment_length &ndash; sum([int(item) for item in re.split('[ACTG^]', read.get_tag("MD")) if not item == '']))/read.query_alignment_length for read in samfile.fetch()] view rawMD-editDistance-fromBam_1.py hosted with &#10084; by GitHub As a sanity check I looped over a bwa mem aligned bam file to extract both the scaled NM tag and scaled MD-derived edit distance and plot those against each other, for which you can see the code below followed by the images. Between gathering the information from the bam file in lists and plotting the data I convert my lists to a numpy array and create a pandas DataFrame (see bottom of post). Since the first plot is rather overcrowded I also added a kernel density estimation to get a better idea of the density of the dots. .gist table { margin-bottom: 0; } import seaborn as sns import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy import stats def makePlot(datadf): plot = sns.jointplot( x='editDistancesNM', y='editDistancesMD', data=datadf, kind="scatter", color="#4CB391", stat_func=stats.pearsonr, space=0, joint_kws={"s": 1}, size=10) plot.savefig('EditDistancesCompared_scatter.png', format='png', dpi=1000) plot = sns.jointplot( x='editDistancesNM', y='editDistancesMD', data=datadf, kind="kde", color="#4CB391", stat_func=stats.pearsonr, space=0, size=10) plot.savefig('EditDistancesCompared_kde.png', format='png', dpi=1000) view rawplotNMvsMD-implementation.py hosted with &#10084; by GitHub The Pearson correlation coefficient is not too bad, but obviously, we want to get exactly the same as the NM tag. It&rsquo;s clear that my implementation of getting the edit distance from the MD tag returns an underestimation of the edit distance. After thinking a while I decided to create a question on biostars&nbsp;where&nbsp;Santosh Anand&nbsp;joined my quest. He suggested counting the mismatches in the MD string, rather than the matches, which I implemented below. So this time I split the MD string on all numbers and sum the length of the mismatched nucleotides. The plots obtained by this approach are shown below the code. .gist table { margin-bottom: 0; } import pysam import re def extractMDFromBam(bam): ''' loop over a bam file and get the edit distance to the reference genome mismatches are stored in the MD tag scale by aligned read length ''' samfile = pysam.AlignmentFile(bam, "rb") return [sum([len(item) for item in re.split('[0-9^]', read.get_tag("MD"))]) / read.query_alignment_length for read in samfile.fetch()] view rawMD-editDistance-fromBam_2.py hosted with &#10084; by GitHub That&rsquo;s an improvement, but we&rsquo;re not yet there. At this point&nbsp;Santosh asked for a few examples of reads which showed a clear deviation in NM and MD-derived edit distance. I wrote a function to do just that, which prints out: The NM tag The MD tag derived edit distance The MD tag The CIGAR string .gist table { margin-bottom: 0; } import sys import re import pysam def extractDisagreement(bam): samfile = pysam.AlignmentFile(bam, "rb") for read in samfile.fetch(): NMdef = read.get_tag("NM")/read.query_alignment_length MDdef = sum([len(item) for item in re.split('[0-9^]', read.get_tag("MD"))])/read.query_alignment_length if NMdef &ndash; MDdef &gt; 0.2: print('\t'.join( [ str(read.get_tag("NM")), str(sum([len(item) for item in re.split('[0-9^]', read.get_tag("MD"))])), read.get_tag("MD"), read.cigarstring, ]) ) view rawgetDisagreementMD-NM-tagEditDistance.py hosted with &#10084; by GitHub This is an example of a read showing a clear mismatch between the NM tag and the MD-derived edit distance. 73 34 21G1C3A1T1A2G4A8T2G4A2G0G0G19^G2G10T6C1A2A0G0G0A5T6^GA5A0G0C0C1G3A5A6C1C2 1684H20M3I7M3I10M1I5M2I1M1I1M3I6M1I4M1I13M2I5M4I3M7I6M1D6M2I1M2I26M1I3M1I2M1I3M2D16M1I6M3I10M10942H The hero on my quest then found that we were still missing the insertions, which are not present in the MD tag but can only be counted by parsing the CIGAR string. I copied the image below from his post. So my successful implementation also parses the CIGAR string to get the insertions and add those to the mismatches. In my final code piece below I show the full code of my eva…