MBQC-base mothur BL-10 methods

Methods

The mothur analysis was performed using the most recent stable version of mothur available at the time, mothur 1.39.3. We used the 64-bit 1.39.3 Linux build without readline support from the following GitHub URL: https://github.com/mothur/mothur/releases/download/v.1.39.3/Mothur.linux_64.noReadLine.zip

Demultiplexed reads in gzipped FASTQ format were downloaded from http://downloads.ihmpdcc.org/data/MBQCBD2/ for all 15 MBQC data subsets, HL-A through HL-O. The gzipped FASTQ files for subsets HL-B, HL-E, HL-H, and HL-K were then uncompressed and recompressed locally using gzip 1.3.5, and without making any other changes to the files. This unusual step was necessary in order to prevent mothur 1.39.3 from entering what appeared to be an infinite loop when trying to read certain of the original gzipped FASTQ files in those 4 subsets. The issue was reproducible on multiple operating systems and mothur versions and appeared only when the original gzipped files were used rather than the uncompressed or recompressed files.

A script was run to create a tab-delimited mothur manifest file for each of the 15 subsets. Each manifest file specifies the group/sample name for every demultiplexed sample (plus “unmatched” reads in some cases), along with the paths to the corresponding R1 and R2 gzipped FASTQ read files. All subsets except for HL-I contained paired-end reads. Initial attempts to process HL-I with mothur’s fastq.info command instead of make.contigs were unsuccessful due to the apparent difficulty of getting fastq.info to accept single-ended reads without losing the correct association of reads with samples. While there is almost certainly a way to do it, it did not appear to be possible with a single fastq.info call and, furthermore, it led to a bifurcation in the processing to accommodate only one of the 15 subsets. Therefore we decided to create “dummy” R2 files for HL-I in which each single-ended read was reverse-complemented and its quality string reversed to create an effectively identical R2 read for each input R1 read. The HL-I reads could then be processed with make.contigs in the same manner as all the rest. It is possible that this approach led to an inflation of the per-base quality scores in the assembled reads (mothur contigs), but we do not believe that any of the steps in the subsequent mothur pipeline were significantly impacted by such a change.

We then downloaded the following alignment and classification databases from the mothur web site for use in the pipeline:

1. The SILVA-based bacterial reference alignment database from the mothur MiSeq SOP,  https://www.mothur.org/w/images/9/98/Silva.bacteria.zip
2. The mothur-formatted v.9 RDP training set, also from the mothur MiSeq SOP, https://www.mothur.org/w/images/5/59/Trainset9_032012.pds.zip
3. The Greengenes 13.5.99 reference taxonomy database, http://www.mothur.org/w/images/9/9d/Gg_13_5_99.taxonomy.tgz

The first two databases were used as-is, but the last one, the Greengenes reference taxonomy, was modified by appending the Greengenes OTU id to the end of every taxonomic classification in the included gg_13_5_99.pds.tax file. So, for example, here is a line from the file before it was modified:

4484460 k__Bacteria;p__Firmicutes;c__Bacilli;o__Bacillales;f__Bacillaceae;g__Bacillus

And here is the same example line post-modification:

4484460 k__Bacteria;p__Firmicutes;c__Bacilli;o__Bacillales;f__Bacillaceae;g__Bacillus;4484460;

Making this change to the Greengenes database allowed us to determine when a contig had been classified down to the OTU level by mothur’s classify.seqs routine. However, it also added significantly to the runtime and memory usage of the naive Bayesian classifier.

After the databases were installed each MBQC subset (HL-A through HL-O) was run through mothur in batch mode using a mothur batch file that performs the following steps (this particular example is for the HL-H subset):

pcr.seqs(fasta=silva.bacteria.fasta, start=13862, end=23444, keepdots=F)
make.contigs(file=HLH_full.files, processors=1)
summary.seqs()
screen.seqs(fasta=current, group=current, maxambig=0)
unique.seqs()
count.seqs(name=current, group=current)
summary.seqs(count=current)
align.seqs(fasta=current, reference=silva.bacteria.pcr.fasta, flip=t)
summary.seqs(count=current)
screen.seqs(fasta=current, count=current, start=20, end=9562, maxhomop=8)
filter.seqs(fasta=current, vertical=T, trump=.)
unique.seqs(fasta=current, count=current)
summary.seqs(count=current)
pre.cluster(fasta=current, count=current, diffs=3)
chimera.vsearch(fasta=current, count=current, dereplicate=t)
remove.seqs(fasta=current, accnos=current)
classify.seqs(fasta=current, count=current, reference=trainset9_032012.pds.fasta, taxonomy=trainset9_032012.pds.tax, cutoff=80)
remove.lineage(fasta=current, count=current, taxonomy=current, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota)
remove.lineage(fasta=current, count=current, taxonomy=current, taxon=unknown)
summary.seqs(count=current)
split.abund(fasta=current, count=current, cutoff=10, accnos=true)
summary.seqs(fasta=HLH_full.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.abund.fasta, count=HLH_full.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.pick.pick.abund.count_table)
classify.seqs(fasta=HLH_full.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.abund.fasta, count=HLH_full.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.pick.pick.abund.count_table, reference=gg_13_5_99.fasta, taxonomy=gg_13_5_99.pds.tax, cutoff=80)
summary.seqs(count=current)

The preceding batch file is based on the recommended MiSeq SOP from https://www.mothur.org/wiki/MiSeq_SOP, with two major differences:

1. After SILVA alignment, filtering, and RDP classification the split.abund command was used with cutoff=10 to remove all unique representative sequences associated with fewer than 10 reads/contigs. After this step only the abundant sequences are processed further. This modification was made in order to reduce the number of unique sequences to be classified to a manageable number.
Instead of performing de novo OTU calling (i.e., some kind of clustering, the usual mothur default) on the filtered and aligned sequences, the classify.seqs command was used to assign abundant representative/unique sequences to Greengenes 13.5.99 OTUs.

2. To reduce the overall runtime the SILVA reference alignment was initially trimmed to just the region covered by the assembled paired end reads (in the very first pcr.seqs step.) The trim coordinates were determined empirically and except for HL-A and HL-D all other subsets were aligned to the 13862-23444 region of the SILVA reference alignment. The paired-end reads in all the other subsets assembled to ~250bp contigs by the mothur make.contigs command, but for HL-A and HL-D the paired-end reads assembled to contigs that were closer to 500bp and aligned with 6332-25436 and 5343-23440, respectively, in the SILVA reference. HL-A and HL-D were first processed the same way as the other subsets, using the 13862-23444 coordinates to effectively truncate them to ~250bp. They were also processed using the aforementioned larger coordinate ranges in order to avoid truncating the contigs at the SILVA alignment step. In addition, in the ~500bp runs the pre.cluster distance was also set to 5 instead of 3, as per mothur’s recommendation to allow a 1bp difference per 100bp of sequence.

After each of the mothur batch runs completed the final classification results were parsed with a script to determine the (total and unique) number of sequences classified down to the Greengenes OTU level (at the default 80% confidence threshold for the classifier’s bootstrapping procedure.) Those per-OTU read counts were mapped back to the original MBQC sample ids and merged with the available metadata to create an OTU table to compare with those produced by the other bioinformatic analyses.