dbcAmplicons pipeline:Bioinformatics

Input files: barcode, primer, samples

Barcode Table

Requires 3 columns: BarcodeID [a name for the pair], Index1 (Read2 in RC), Index2 (Read3) in a plain tab-delimited text file. Orientation is important, but you can change in the preprocess arguments. First line is a comment and just help me remembers.

Barcode table for the workshop.

Primer Table

Requires 4 columns: the read in which the primer should be checked for (allowable are P5/P7, R1/R2, READ1/READ2, F/R, FORWARD/REVERSE, Primer Pair ID describes which should be found ‘together’, Primer ID individual id, and sequence (IUPAC ambiguity characters are allowed).

Primer table for the workshop.

Samples Table

Requires 4 columns and a header: SampleID samples name, PrimerPairID same as in primer file, barcodeID same as in barcode file, and ProjectID which represents the file prefix for the output and can include a path. SampleID, PrimerPairID, BarcodeID pairs must be unique. In addition for PrimerPairID, can be comma separated, * (match any primer), or ‘-’ should match no primer.

Additional columns are allowed and will be added to the biom file in dbcAmplicons abundances.

Samples table for the workshop.

Sequence Reads

Typically you receive fastq file(s) from the sequencing provider.

• Fastq files are actually not raw data from the provider, “raw” data is actually bcl files.
• Sequencing provider will run an application bcl2fastq with a sample sheet to produce demultiplexed (by barcode) fastq files.
• For dbcAmplicons you want to request from your sequencing provider non-demultiplexed fastq (so one set of fastqs for the entire run) with the index reads.

fastq

fastq files combine the sequence and quality scores into 1 file. Each sequence here has 4 lines (should be enforced strictly), header, sequence, historical ‘+’, and quality.

CASAVA 1.8 Read IDs

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG

• EAS139 the unique instrument name
• 136 the run id
• FC706VJ the flowcell id
• 2 flowcell lane
• 2104 tile number within the flowcell lane
• 15343 ’x’-coordinate of the cluster within the tile
• 197393 ’y’-coordinate of the cluster within the tile
• 1 the member of a pair, 1 or 2 (paired-end or mate-pair reads only)
• Y Y if the read fails filter (read is bad), N otherwise
• 18 0 when none of the control bits are on, otherwise it is an even number
• ATCACG index sequence

Quality scores

Quality scores are paired 1 to 1 with sequence characters.

Each quality character has a numerical value associated with it (ASCII value). In Illumina 1.8+ you subtract 33 from the ascii value associated with the quality character to get the quality score.

dbcAmplicons

Preprocessing reads

Application specific downstream

Population Community Profiling ( i.e. microbial, bacterial, fungal, etc. )

dbcAmplicons Python Application

• Screen - Using Bowtie2, screen targets against a reference fasta file, separating reads by those that produce matches and those that do not match sequences in the reference database.
• Classify - Wrapper around the MSU Ribosomal Database Project (RDP) Classifier for Bacterial and Archaeal 16S rRNA sequences, Fungal 28S rRNA, fungal ITS regions
• Abundance - Reduce RDP classifier results to abundance tables (or biom file format), rows are taxa and columns are samples ready for additional community analysis.
• Extract - Extract the reads associated with as taxonomic group for separate processing.

Targeted Re-sequencing

Set up R functions

• Consensus - Reduce reads to consensus sequence for each sample and amplicon.
• Most Common – Reduce reads to the most commonly occurring read in the sample and amplicon ( that is present in at least 5% and 5 reads, by default )
• Haplotypes – Impute the different haplotypes in the sample and amplicon

Supplemental Scripts

• convert2Readto4Read.py - For when samples are processed by someone else
• splitReadsBySample.py - To facilitate upload to the SRA
• preprocPair_with_inlineBC.py - Cut out inline BC and create 4 reads for standard input processing Will work with ”Mills lab” protocol
• dbcVersionReport.sh - Print out version numbers of all tools

Validate

Validate the barcode, primer and sample sheet (Must have all 3). Can be performed before data is available to check and be ready for preprocessing once available.

1. Read in the metadata input tables: Barcodes, Primers (optional), Samples (optional).
2. Validate that there are no duplicate sample/barcode/primer combos.
3. Validate that all barcodes listed in the sample sheet are in the barcode table.
4. Validate that all primers listed in the sample sheet are in the primer table.

Report any errors.

Preprocessing

Preprocess reads for barcode and primer, match to sample sheet, only “legitimate” (contain matching barcode/primer) will be output.

1. Read in the metadata input tables: Barcodes, Primers (optional), Samples (optional)
2. Read in a batch of reads (default 100,000), for each read.
3. Compare index barcodes to the barcode table, note best matching barcode.
4. Compare 5’ end of reads to the primer table, note best matching primer.
5. Compare to barcode:primer pair to the sample table, note sampleID and projectID.
6. If its a legitimate reads (contains matching barcode,primer,sample) output the read pair to the output file.
7. Output Identified_Barcodes.txt file.

Output: Preprocessed reads
Identified_Barcodes.txt

Barcode Comparison - compares each barcode to all possible barcodes and returns the best match < desired edit distance.

Primer Comparison - compares the beginning (primer region) of each read to all possible primers and returns the best match < specified maximimum Levenshtein distance + final 4 exact match

The new read header

Join

Uses Flash2 to merge reads that overlap to produce a longer (or sometimes shorter read).

Modification to original Flash algorithm include:

• Performs complete overlaps with adapter trimming Allows for different sized reads (after cutting primer off)
• Discards reads with > 50% Q of 10 or less, which are indicative of adapter/primer dimers

Output: prefix.notCombined_1.fastq.gz prefix.notCombined_2.fastq.gz
prefix.extendedFrags.fastq.gz
prefix.hist
prefix.histogram

Overlapping and Adapter trimming by overlapping reads.

Consider the three scenarios below

Insert size > length of the number of cycles

Insert size < length of the number of cycles (10bp min)

Insert size < length of the read length

Overlapping produces histograms of the overlapped Reads

Classify

For each read determine which organism is likely came from.

• Uses the RDP (Ribosomal Database Project) classifier for bacterial and archaeal 16S, fungal LSU, ITS warcup/unite databases. You can provide your own training database

• Classifies sequences to the closest taxonomic reference provides a bootstrap score for reliability

• Concatenates Paired-end reads

  • Can trim off low quality ends, to some value Q

Output: fixrank file

Ribosomal Database Project (RDP) - naïve Bayesian Classifier

• Compares each read to a database
• Database is updates periodically
• Compares by k-mers (15 mers)
• 100 bootstraps to establish confidence in result

Order does not matter, no 3% !

• Drawbacks
  • Accepts only fasta (though website implies fastq) files
  • Can be slow
  • Down to genus only (for 16s, species for ITS)
  • Kmer database are based on whole 16s
  • Cannot group together unknown OTUs that represent unique taxa

Clustering

Clustering – “Because of the increasing sizes of today’s amplicon datasets, fast and greedy de novo clustering heuristics are the preferred and only practical approach to produce OTUs”. I DISAGREE

Shared steps in these algorithms are:

1. An amplicon is drawn out of the amplicon pool and becomes the center of a new OTU (centroid selection)
2. This centroid is then compared to all other amplicons remaining in the pool.
3. Amplicons for which the distance is within a global clustering threshold, t (e.g. 3%), to the centroid are moved from the pool to the OUT
4. The OTU is then closed. These steps are repeated as long as amplicons remain in the pool.

Why I’m not a fan

• Little to no biological rational to any of the clustering parameters, modify the parameters to get a result you like.
• Dependent on ordering, reorder our reads you can get different set of OTUs. Often not repeatable from run to run.
• 3% (or any other cutoff) is BS. Most clustering algorithms do not consider sequencing errors.
• If you generate more data you have to start the clustering process all over again as population of sequences matters.
• I’m sure there is more

A Comparison study

Abundance

Takes fixrank file(s) produced from classify and outputs abundance tables and taxa_info table

Abundance tables

• Rows are taxa
• Columns are samples
• Counts of the number of amplicons for each taxa/samples

Proportions tables

• Same as abundance but each cell is the proportion of amplicons (so counts in cell divided by the columns sum)

Biom file (Biological Observation Matrix)

• JSON/hd5 file format for microbiome files
• http://biom-format.org
• Abundance tables are 0 heavy, a biom file removes the 0’s as well as stores extra metadata

The BIOM file format (canonically pronounced biome) is designed to be a general-use format for representing biological sample by observation contingency tables. BIOM is a recognized standard for the Earth Microbiome Project and is a Genomics Standards Consortium supported project. Contains the abundance counts, the sample names, full taxonomic string [domain through genus/species], and any sample metadata in the sample sheet.

Abundance file

Proportions file

Taxa Info file

Supplies extra information about the tax identified in the experiment as well as the full taxonomic path.

Future Directions

• dbcAmplicons is a data reduction pipeline, produces abundance/biome files, post processing most typically done in R.
• Update to python 3
• Change to allow for bcl2fastq demultiplexed files.
• Replace primer algorithm and FLASH2 with htstream.
• Include screening of diversity sample in preprocessing to get an idea of actual proportion in the pool. Available in htstream.
• Replace RDP classification with another scheme
  • Dada2 implementation of RDP classifier
• Use amplicon length in classification
• Add in ability to generate a phylo tree.
• Output data for rarefaction curves

Post Processing

Pretty much do all of my post analysis (abundance table, Biome) in R

• Common Packages
  • Vegan
  • Vegetarian
  • Phyloseq (uses vegan, ade4, ape, picante)

Ecological Diversity Analysis

• how does the structure of taxa across samples/groups compare.

Ordination Analysis (multivariate analysis)

• Visualize the relative similarity/dissimilarity across samples, test for taxa/environment relationships

Differential Abundance Analysis (univariate analysis)

• Uses tools from RNAseq (limma, edgeR)

Visualization

• temporal
• heatmaps
• ’trees’
• more

Multi-community Analysis

Because dbcAmplicons is suitable for multiple amplicons across multiple communities, cross community analysis is possible.