OBITools V4

Eric Coissac

1/17/23

Preface

The first version of OBITools started to be developed in 2005. This was at the beginning of the DNA metabarcoding story at the Laboratoire d’Ecologie Alpine (LECA) in Grenoble. At that time, with Pierre Taberlet and François Pompanon, we were thinking about the potential of this new methodology under development. PIerre and François developed more the laboratory methods, while I was thinking more about the tools for analysing the sequences produced. Two ideas were behind this development. I wanted something modular, and something easy to extend. To achieve the first goal, I decided to implement obitools as a suite of unix commands mimicking the classic unix commands but dedicated to sequence files. The basic unix commands are very useful for automatically manipulating, parsing and editing text files. They work in flow, line by line on the input text. The result is a new text file that can be used as input for the next command. Such a design makes it possible to quickly develop a text processing pipeline by chaining simple elementary operations. The OBITools are the exact counterpart of these basic Unix commands, but the basic information they process is a sequence (potentially spanning several lines of text), not a single line of text. Most OBITools consume sequence files and produce sequence files. Thus, the principles of chaining and modularity are respected. In order to be able to easily extend the OBITools to keep up with our evolving ideas about processing DNA metabarcoding data, it was decided to develop them using an interpreted language: Python. Python 2, the version available at the time, allowed us to develop the OBITools efficiently. When parts of the algorithms were computationally demanding, they were implemented in C and linked to the Python code. Even though Python is not the most efficient language available, even though computers were not as powerful as they are today, the size of the data we could produce using 454 sequencers or early solexa machines was small enough to be processed in a reasonable time.

The first public version of obitools was OBITools2 (Boyer et al. 2016), this was actually a cleaned up and documented version of OBITools that had been running at LECA for years and was not really distributed except to a few collaborators. This is where OBITools started its public life from then on. The DNA metabarcoding spring schools provided and still provide user training every year. But OBITools2 soon suffered from two limitations: it was developed in Python2, which was increasingly abandoned in favour of Python3, and the data size kept increasing with the new illumina machines. Python’s intrinsic slowness coupled with the increasing size of the datasets made OBITools computation times increasingly long. The abandonment of all maintenance of Python2 by its developers also imposed the need for a new version of OBITools.

OBITools3 was the first response to this crisis. Developed and maintained by Céline Mercier, OBITools3 attempted to address several limitations of OBITools2. It is a complete new code, mainly developed in Python3, with most of the lower layer code written in C for efficiency. OBITools3 has also abandoned text files for binary files for the same reason of efficiency. They have been replaced by a database structure that keeps track of every operation performed on the data.

Here we present OBITools4 which can be seen as a return to the origins of OBITools. While OBITools3 offered traceability of analyses, which is in line with the concept of open science, and faster execution, OBITools2 was more versatile and not only usable for the analysis of DNA metabarcoding data. OBITools4 is the third full implementation of OBITools. The idea behind this new version is to go back to the original design of OBITools which ran on text files containing sequences, like the classic Unix commands, but running at least as fast as OBITools3 and taking advantage of the multicore architecture of all modern laptops. For this, the idea of relying on an interpreted language was abandoned. The OBITools4 are now fully implemented in the GO language with the exception of a few small pieces of specific code already implemented very efficiently in C. OBITools4 also implement a new format for the annotations inserted in the header of every sequences. Rather tha relying on a format specific to OBITools, by default OBITools4 use the JSON format. This simplifies the writing of parsers in any languages, and thus allows obitools to easiestly interact with other software.

1 Installation of the OBITools

1.1 Availability of the OBITools

The OBITools are open source and protected by the CeCILL 2.1 license.

All the sources of the OBITools4 can be downloaded from the metabarcoding git server (https://git.metabarcoding.org).

1.2 Prerequisites

The OBITools4 are developped using the GO programming language, we stick to the latest version of the language, today the 1.21.41.21.4. If you want to download and compile the sources yourself, you first need to install the corresponding compiler on your system. Some parts of the soft are also written in C, therefore a recent C compiler is also requested, GCC on Linux or Windows, the Developer Tools on Mac.

Whatever the installation you decide for, you will have to ensure that a C compiler is available on your system.

1.3 Installation with the install script

An installation script that compiles the new OBITools on your Unix-like system is available online. The easiest way to run it is to copy and paste the following command into your terminal

curl -L https://metabarcoding.org/obitools4/install.sh | bash

By default, the script installs the OBITools commands and other associated files into the /usr/local directory. The names of the commands in the new OBITools4 are mostly identical to those in OBITools2. Therefore, installing the new OBITools may hide or delete the old ones. If you want both versions to be available on your system, the installation script offers two options:

-i, –install-dir Directory where OBITools are installed (as example use /usr/local not /usr/local/bin).

-p, –obitools-prefix Prefix added to the OBITools command names if you want to have several versions of obitools at the same time on your system (as example -p g will produce gobigrep command instead of obigrep).

You can use these options by following the installation command:

curl -L https://metabarcoding.org/obitools4/install.sh | \
      bash -s -- --install-dir test_install --obitools-prefix k

In this case, the binaries will be installed in the test_install directory and all command names will be prefixed with the letter k. Thus obigrep will be named kobigrep.

1.4 Compilation from sources

2 File formats usable with OBITools

OBITools manipulate have to manipulate DNA sequence data and taxonomical data. They can use some supplentary metadata describing the experiment and produce some stats about the processed DNA data. All the manipulated data are stored in text files, following standard data format.

2.1 The DNA sequence data

Sequences can be stored following various format. OBITools knows some of them. The central formats for sequence files manipulated by OBITools scripts are the FASTA and FASTQ format. OBITools extends the both these formats by specifying a syntax to include in the definition line data qualifying the sequence. All file formats use the IUPAC code for encoding nucleotides.

Moreover these two formats that can be used as input and output formats, OBITools4 can read the following format :

EBML flat file format (use by ENA)

Genbank flat file format

ecoPCR output files

2.1.1 The IUPAC Code

The International Union of Pure and Applied Chemistry (IUPAC) defined the standard code for representing protein or DNA sequences.

Code | Nucleotide
A | Adenine
C | Cytosine
G | Guanine
T | Thymine
U | Uracil
R | Purine (A or G)
Y | Pyrimidine (C, T, or U)
M | C or A
K | T, U, or G
W | T, U, or A
S | C or G
B | C, T, U, or G (not A)
D | A, T, U, or G (not C)
H | A, T, U, or C (not G)
V | A, C, or G (not T, not U)
N | Any base (A, C, G, T, or U)

2.1.2 The FASTA sequence format

The FASTA format is certainly the most widely used sequence file format. This is certainly due to its great simplicity. It was originally created for the Lipman and Pearson FASTA program. OBITools use in more of the classical FASTA format an extended version of this format where structured data are included in the title line.

In FASTA format a sequence is represented by a title line beginning with a > character and the sequences by itself following the IUPAC code. The sequence is usually split other severals lines of the same length (expect for the last one)

>my_sequence this is my pretty sequence
ACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGT
GTGCTGACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTGTTT
AACGACGTTGCAGTACGTTGCAGT

This is no special format for the title line excepting that this line should be unique. Usually the first word following the > character is considered as the sequence identifier. The end of the title line corresponding to a description of the sequence. Several sequences can be concatenated in a same file. The description of the next sequence is just pasted at the end of the record of the previous one

>sequence_A this is my first pretty sequence
ACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGT
GTGCTGACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTGTTT
AACGACGTTGCAGTACGTTGCAGT
>sequence_B this is my second pretty sequence
ACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGT
GTGCTGACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTGTTT
AACGACGTTGCAGTACGTTGCAGT
>sequence_C this is my third pretty sequence
ACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGT
GTGCTGACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTACGTTGCAGTGTTT
AACGACGTTGCAGTACGTTGCAGT

2.1.2.1 File extensions

There is no standard file extension for a FASTA file, but .fa and .fasta, are commonly used.

2.1.3 The FASTQ sequence format1

The FASTQ format is a text file format for storing both biological sequences (only nucleic acid sequences) and the associated sequence quality scores. Every nucleotide of the sequence and its associated quality score are each encoded by a single ASCII character. This format was originally developed by the Wellcome Trust Sanger Institute to link a FASTA sequence file to the corresponding quality data, but is now became the de facto standard for storing results from high-throughput sequencers (Cock et al. 2010).

OBITools considers that a FASTQ file uses four lines to encode a sequence record.

Line 1 begins with a ‘@’ character and is followed by a sequence identifier and an optional description (like a FASTA title line).

Line 2 is the sequence letters, in upper or lower case, but OBITools only write sequences as lower cases.

Line 3 begins with a ‘+’ character and is optionally followed by the same sequence identifier (and any description) again.

Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.

A FASTQ file looks like this:

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

The character ‘!’ represents the lowest quality while ‘~’ is the highest. Here are the quality value characters in left-to-right increasing order of quality (ASCII):

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]
^_`abcdefghijklmnopqrstuvwxyz{|}~

If the original Sanger FASTQ files also allowed the sequence and quality strings to be wrapped (split over multiple lines), it is not supported by OBITools as it make parsing complicated due to the unfortunate choice of “@” and “+” as markers (these characters can also occur in the quality string).

Sequence quality scores

The Phred quality value Q is an integer mapping of p (i.e., the probability that the corresponding base call is incorrect). Two different equations have been in use. The first is the standard Sanger variant to assess reliability of a base call, otherwise known as Phred quality score:

Qsanger=−10log10p Q_\text{sanger} = -10 \, \log_{10} p

The Solexa pipeline (i.e., the software delivered with the Illumina Genome Analyzer) earlier used a different mapping, encoding the odds 𝐩/(1−𝐩)\mathbf{p}/(1-\mathbf{p}) instead of the probability 𝐩\mathbf{p}:

Qsolexa-prior to v.1.3=−10log10p1−p Q_\text{solexa-prior to v.1.3} = -10 \; \log_{10} \frac{p}{1-p}

Although both mappings are asymptotically identical at higher quality values, they differ at lower quality levels (i.e., approximately 𝐩>0.05\mathbf{p} > 0.05, or equivalently, 𝐐<13\mathbf{Q} < 13).

Figure 2.1: Relationship between Q and p using the Sanger (red) and Solexa (black) equations (described above). The vertical dotted line indicates 𝐩=0.05\mathbf{p}= 0.05, or equivalently, Q=13Q = 13.

Encoding

The FASTQ format had differente way of encoding the Phred quality score along the time. Here a breif history of these changes is presented. OBITools

Sanger format can encode a Phred quality score from 0 to 93 using ASCII 33 to 126 (although in raw read data the Phred quality score rarely exceeds 60, higher scores are possible in assemblies or read maps).

Solexa/Illumina 1.0 format can encode a Solexa/Illumina quality score from -5 to 62 using ASCII 59 to 126 (although in raw read data Solexa scores from -5 to 40 only are expected)

Starting with Illumina 1.3 and before Illumina 1.8, the format encoded a Phred quality score from 0 to 62 using ASCII 64 to 126 (although in raw read data Phred scores from 0 to 40 only are expected).

Starting in Illumina 1.5 and before Illumina 1.8, the Phred scores 0 to 2 have a slightly different meaning. The values 0 and 1 are no longer used and the value 2, encoded by ASCII 66 “B”.

Sequencing Control Software, Version 2.6, (Catalog # SY-960-2601, Part # 15009921 Rev. A, November 2009, page 30) states the following: If a read ends with a segment of mostly low quality (Q15 or below), then all of the quality values in the segment are replaced with a value of 2 (encoded as the letter B in Illumina’s text-based encoding of quality scores)… This Q2 indicator does not predict a specific error rate, but rather indicates that a specific final portion of the read should not be used in further analyses. Also, the quality score encoded as “B” letter may occur internally within reads at least as late as pipeline version 1.6, as shown in the following example:

@HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGA
TTTGTTGGGGGAGACATTTTTGTGATTGCCTTGAT
+HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
efcfffffcfeefffcffffffddf`feed]`]_Ba_^__[YBBBBBBBBBBRTT\]][ dddd`
ddd^dddadd^BBBBBBBBBBBBBBBBBBBBBBBB

An alternative interpretation of this ASCII encoding has been proposed. Also, in Illumina runs using PhiX controls, the character ‘B’ was observed to represent an “unknown quality score”. The error rate of ‘B’ reads was roughly 3 phred scores lower the mean observed score of a given run.

Starting in Illumina 1.8, the quality scores have basically returned to the use of the Sanger format (Phred+33).

OBITools follows the Sanger format. Nevertheless, It is possible to read files encoded following the Solexa/Illumina format by applying a shift of 62 (see the option –solexa of the OBITools commands).

2.1.3.1 File extensions

There is no standard file extension for a FASTQ file, but .fq and .fastq, are commonly used.

2.2 The taxonomy files

Many OBITools are able to take into account taxonomic data. This is done by specifying a directory containing all :doc:NCBI taxonomy dump files <./taxdump>.

2.3 The sample description file

A key file for OBITools4 is the file describing all samples (PCR) analyzed in the processed sequencing library file. This file, often called the ngsfilter file, is a tab separated values (TSV) file. The format of this file is exactly identical to that used in OBITools2 and OBITools4.

{tsv, .smaller} wolf_diet 13a_F730603 aattaac TTAGATACCCCACTATGC TAGAACAGGCTCCTCTAG F wolf_diet 15a_F730814 gaagtag TTAGATACCCCACTATGC TAGAACAGGCTCCTCTAG F wolf_diet 26a_F040644 gaatatc TTAGATACCCCACTATGC TAGAACAGGCTCCTCTAG F wolf_diet 29a_F260619 gcctcct TTAGATACCCCACTATGC TAGAACAGGCTCCTCTAG F

At least six columns must be present in every line of the file.

The first column contains the name of the experience:

An experiment name groups a set of sample together. Sequences belonging to the experiment are tagged with an attribute experiment containing the name of the experiment in their annotation.

The second column contains the sample identifier in the experiment

The sample identifier must be unique in the experiment. The obimultiplex and obitagpcr commands add to all the sequences bellonging to the same sample an attribute sample containing the sample identifier

The third column contains description of the tag used to identify sequences corresponding to this sample

The fourth column contains the forward primer sequence

The fifth column contains the reverse primer sequence

The sixth column must always contain the character F (full length)

This article uses material from the Wikipedia article FASTQ format which is released under the Creative Commons Attribution-Share-Alike License 3.0↩︎

3 OBITools V4 Tutorial

Here is a short tutorial on how to analyze DNA metabarcoding data produced on Illumina sequencers using:

the OBITools

some basic Unix commands

3.1 Wolves’ diet based on DNA metabarcoding

The data used in this tutorial correspond to the analysis of four wolf scats, using the protocol published in Shehzad et al. (2012) for assessing carnivore diet. After extracting DNA from the faeces, the DNA amplifications were carried out using the primers TTAGATACCCCACTATGC and TAGAACAGGCTCCTCTAG amplifiying the 12S-V5 region (Riaz et al. 2011), together with a wolf blocking oligonucleotide.

The complete data set can be downloaded here: the tutorial dataset

Once the data file is downloaded, using a UNIX terminal unarchive the data from the tgz file.

tar zxvf wolf_diet.tgz

That command create a new directory named wolf_data containing every required data files:

fastq <fastq> files resulting of aGA IIx (Illumina) paired-end (2 x 108 bp) sequencing assay of DNA extracted and amplified from four wolf faeces:

wolf_F.fastq

wolf_R.fastq

the file describing the primers and tags used for all samples sequenced:

wolf_diet_ngsfilter.txt The tags correspond to short and specific sequences added on the 5' end of each primer to distinguish the different samples

the file containing the reference database in a fasta format:

db_v05_r117.fasta This reference database has been extracted from the release 117 of EMBL using obipcr

To not mix raw data and processed data a new directory called results is created.

mkdir results

3.2 Step by step analysis

3.2.1 Recover full sequence reads from forward and reverse partial reads

When using the result of a paired-end sequencing assay with supposedly overlapping forward and reverse reads, the first step is to recover the assembled sequence.

The forward and reverse reads of the same fragment are at the same line position in the two fastq files obtained after sequencing. Based on these two files, the assembly of the forward and reverse reads is done with the obipairing utility that aligns the two reads and returns the reconstructed sequence.

In our case, the command is:

obipairing --min-identity=0.8 \
           --min-overlap=10 \
           -F wolf_data/wolf_F.fastq \
           -R wolf_data/wolf_R.fastq \
           > results/wolf.fastq

The --min-identity and --min-overlap options allow discarding sequences with low alignment quality. If after the aligment, the overlaping parts of the reads is shorter than 10 base pairs or the similarity over this aligned region is below 80% of identity, in the output file, the forward and reverse reads are not aligned but concatenated, and the value of the mode attribute in the sequence header is set to joined instead of alignment.

3.2.2 Remove unaligned sequence records

Unaligned sequences (:pymode=joined) cannot be used. The following command allows removing them from the dataset:

obigrep -p 'annotations.mode != "join"' \
        results/wolf.fastq > results/wolf.ali.fastq

The -p requires a go like expression. annotations.mode != "join" means that if the value of the mode annotation of a sequence is different from join, the corresponding sequence record will be kept.

The first sequence record of wolf.ali.fastq can be obtained using the following command line:

head -n 4 results/wolf.ali.fastq

The folling piece of code appears on thew window of tour terminal.

@HELIUM_000100422_612GNAAXX:7:108:5640:3823#0/1 {"ali_dir":"left","ali_length":62,"mode":"alignment","pairing_mismatches":{"(T:26)->(G:13)":62,"(T:34)->(G:18)":48},"score":484,"score_norm":0.968,"seq_a_single":46,"seq_ab_match":60,"seq_b_single":46}
ccgcctcctttagataccccactatgcttagccctaaacacaagtaattaatataacaaaattgttcgccagagtactaccggcaatagcttaaaactcaaaggacttggcggtgctttatacccttctagaggagcctgttctaaggaggcgg
+
CCCCCCCBCCCCCCCCCCCCCCCCCCCCCCBCCCCCBCCCCCCC<CcCccbe[`F`accXV<TA\RYU\\ee_e[XZ[XEEEEEEEEEE?EEEEEEEEEEDEEEEEEECCCCCCCCCCCCCCCCCCCCCCCACCCCCACCCCCCCCCCCCCCCC

3.2.3 Assign each sequence record to the corresponding sample/marker combination

Each sequence record is assigned to its corresponding sample and marker using the data provided in a text file (here wolf_diet_ngsfilter.txt). This text file contains one line per sample, with the name of the experiment (several experiments can be included in the same file), the name of the tags (for example: aattaac if the same tag has been used on each extremity of the PCR products, or aattaac:gaagtag if the tags were different), the sequence of the forward primer, the sequence of the reverse primer, the letter T or F for sample identification using the forward primer and tag only or using both primers and both tags, respectively (see obimultiplex for details).

obimultiplex -t wolf_data/wolf_diet_ngsfilter.txt \
             -u results/unidentified.fastq \
             results/wolf.ali.fastq \
             > results/wolf.ali.assigned.fastq

This command creates two files:

unidentified.fastq containing all the sequence records that were not assigned to a sample/marker combination

wolf.ali.assigned.fastq containing all the sequence records that were properly assigned to a sample/marker combination

Note that each sequence record of the wolf.ali.assigned.fastq file contains only the barcode sequence as the sequences of primers and tags are removed by the obimultiplex program. Information concerning the experiment, sample, primers and tags is added as attributes in the sequence header.

For instance, the first sequence record of wolf.ali.assigned.fastq is:

@HELIUM_000100422_612GNAAXX:7:108:5640:3823#0/1_sub[28..127] {"ali_dir":"left","ali_length":62,"direction":"direct","experiment":"wolf_diet","forward_match":"ttagataccccactatgc","forward_mismatches":0,"forward_primer":"ttagataccccactatgc","forward_tag":"gcctcct","mode":"alignment","pairing_mismatches":{"(T:26)->(G:13)":35,"(T:34)->(G:18)":21},"reverse_match":"tagaacaggctcctctag","reverse_mismatches":0,"reverse_primer":"tagaacaggctcctctag","reverse_tag":"gcctcct","sample":"29a_F260619","score":484,"score_norm":0.968,"seq_a_single":46,"seq_ab_match":60,"seq_b_single":46}
ttagccctaaacacaagtaattaatataacaaaattgttcgccagagtactaccggcaatagcttaaaactcaaaggacttggcggtgctttataccctt
+
CCCBCCCCCBCCCCCCC<CcCccbe[`F`accXV<TA\RYU\\ee_e[XZ[XEEEEEEEEEE?EEEEEEEEEEDEEEEEEECCCCCCCCCCCCCCCCCCC

3.2.4 Dereplicate reads into uniq sequences

The same DNA molecule can be sequenced several times. In order to reduce both file size and computations time, and to get easier interpretable results, it is convenient to work with unique sequences instead of reads. To dereplicate such reads into unique sequences, we use the obiuniq command.

Definition: Dereplicate reads into unique sequences

compare all the reads in a data set to each other

group strictly identical reads together

output the sequence for each group and its count in the original dataset (in this way, all duplicated reads are removed)

Definition adapted from Seguritan and Rohwer (2001)

For dereplication, we use the obiuniq command with the -m sample. The -m sample option is used to keep the information of the samples of origin for each uniquesequence.

obiuniq -m sample \
        results/wolf.ali.assigned.fastq \
        > results/wolf.ali.assigned.uniq.fasta

Note that obiuniq returns a fasta file.

The first sequence record of wolf.ali.assigned.uniq.fasta is:

>HELIUM_000100422_612GNAAXX:7:93:6991:1942#0/1_sub[28..126] {"ali_dir":"left","ali_length":63,"count":1,"direction":"reverse","experiment":"wolf_diet","forward_match":"ttagataccccactatgc","forward_mismatches":0,"forward_primer":"ttagataccccactatgc","forward_tag":"gaatatc","merged_sample":{"26a_F040644":1},"mode":"alignment","pairing_mismatches":{"(A:10)->(G:34)":76,"(C:06)->(A:34)":58},"reverse_match":"tagaacaggctcctctag","reverse_mismatches":0,"reverse_primer":"tagaacaggctcctctag","reverse_tag":"gaatatc","score":730,"score_norm":0.968,"seq_a_single":45,"seq_ab_match":61,"seq_b_single":45}
ttagccctaaacataaacattcaataaacaagaatgttcgccagagaactactagcaaca
gcctgaaactcaaaggacttggcggtgctttatatccct

The run of obiuniq has added two key=values entries in the header of the fasta sequence:

"merged_sample":{"29a_F260619":1}: this sequence have been found once in a single sample called 29a_F260619

"count":1 : the total count for this sequence is 11

To keep only these two attributes, we can use the obiannotate command:

obiannotate -k count -k merged_sample \
  results/wolf.ali.assigned.uniq.fasta \
  > results/wolf.ali.assigned.simple.fasta

The first five sequence records of wolf.ali.assigned.simple.fasta become:

>HELIUM_000100422_612GNAAXX:7:26:18930:11105#0/1_sub[28..127] {"count":1,"merged_sample":{"29a_F260619":1}}
ttagccctaaacacaagtaattaatataacaaaatwattcgcyagagtactacmggcaat
agctyaaarctcamagrwcttggcggtgctttataccctt
>HELIUM_000100422_612GNAAXX:7:58:5711:11399#0/1_sub[28..127] {"count":1,"merged_sample":{"29a_F260619":1}}
ttagccctaaacacaagtaattaatataacaaaattattcgccagagtwctaccgssaat
agcttaaaactcaaaggactgggcggtgctttataccctt
>HELIUM_000100422_612GNAAXX:7:100:15836:9304#0/1_sub[28..127] {"count":1,"merged_sample":{"29a_F260619":1}}
ttagccctaaacatagataattacacaaacaaaattgttcaccagagtactagcggcaac
agcttaaaactcaaaggacttggcggtgctttataccctt
>HELIUM_000100422_612GNAAXX:7:55:13242:9085#0/1_sub[28..126] {"count":4,"merged_sample":{"26a_F040644":4}}
ttagccctaaacataaacattcaataaacaagagtgttcgccagagtactactagcaaca
gcctgaaactcaaaggacttggcggtgctttacatccct
>HELIUM_000100422_612GNAAXX:7:86:8429:13723#0/1_sub[28..127] {"count":7,"merged_sample":{"15a_F730814":5,"29a_F260619":2}}
ttagccctaaacacaagtaattaatataacaaaattattcgccagagtactaccggcaat
agcttaaaactcaaaggactcggcggtgctttataccctt

3.2.5 Denoise the sequence dataset

To have a set of sequences assigned to their corresponding samples does not mean that all sequences are biologically meaningful i.e. some of these sequences can contains PCR and/or sequencing errors, or chimeras.

Tag the sequences for PCR errors (sequence variants)

The obiclean program tags sequence variants as potential error generated during PCR amplification. We ask it to keep the head sequences (-H option) that are sequences which are not variants of another sequence with a count greater than 5% of their own count (-r 0.05 option).

obiclean -s sample -r 0.05 -H \
  results/wolf.ali.assigned.simple.fasta \
      > results/wolf.ali.assigned.simple.clean.fasta

One of the sequence records of wolf.ali.assigned.simple.clean.fasta is:

>HELIUM_000100422_612GNAAXX:7:66:4039:8016#0/1_sub[28..127] {"count":17,"merged_sample":{"13a_F730603":17},"obiclean_head":true,"obiclean_headcount":1,"obiclean_internalcount":0,"obi
clean_samplecount":1,"obiclean_singletoncount":0,"obiclean_status":{"13a_F730603":"h"},"obiclean_weight":{"13a_F730603":25}}
ctagccttaaacacaaatagttatgcaaacaaaactattcgccagagtactaccggcaac
agcccaaaactcaaaggacttggcggtgcttcacaccctt

To remove such sequences as much as possible, we first discard rare sequences and then rsequence variants that likely correspond to artifacts.

Get some statistics about sequence counts

obicount results/wolf.ali.assigned.simple.clean.fasta

time="2023-02-23T18:43:46+01:00" level=info msg="Appending results/wolf.ali.assigned.simple.clean.fasta file\n"
 2749 36409 273387

The dataset contains 43134313 sequences variant corresponding to 42452 sequence reads. Most of the variants occur only a single time in the complete dataset and are usualy named singletons

obigrep -p 'sequence.Count() == 1' results/wolf.ali.assigned.simple.clean.fasta \
    | obicount

time="2023-02-23T18:43:46+01:00" level=info msg="Reading sequences from stdin in guessed\n"
time="2023-02-23T18:43:46+01:00" level=info msg="Appending results/wolf.ali.assigned.simple.clean.fasta file\n"
time="2023-02-23T18:43:46+01:00" level=info msg="On output use JSON headers"
 2623 2623 261217

In that dataset sigletons corresponds to 35113511 variants.

Using R and the ROBIFastread package able to read headers of the fasta files produced by OBITools, we can get more complete statistics on the distribution of occurrencies.

library(ROBIFastread)
library(ggplot2)

seqs <- read_obifasta("results/wolf.ali.assigned.simple.clean.fasta",keys="count")

ggplot(data = seqs,  mapping=aes(x = count)) +
  geom_histogram(bins=100) +
  scale_y_sqrt() +
  scale_x_sqrt() +
  geom_vline(xintercept = 10, col="red", lty=2) +
  xlab("number of occurrencies of a variant")

In a similar way it is also possible to plot the distribution of the sequence length.

ggplot(data = seqs,  mapping=aes(x = nchar(sequence))) +
  geom_histogram() +
  scale_y_log10() +
  geom_vline(xintercept = 80, col="red", lty=2) +
  xlab("sequence lengths in base pair")

Keep only the sequences having a count greater or equal to 10 and a length shorter than 80 bp

Based on the previous observation, we set the cut-off for keeping sequences for further analysis to a count of 10. To do this, we use the obigrep <scripts/obigrep> command. The -p 'count>=10' option means that the python expression :pycount>=10 must be evaluated to :pyTrue for each sequence to be kept. Based on previous knowledge we also remove sequences with a length shorter than 80 bp (option -l) as we know that the amplified 12S-V5 barcode for vertebrates must have a length around 100bp.

obigrep -l 80 -p 'sequence.Count() >= 10' results/wolf.ali.assigned.simple.clean.fasta \
    > results/wolf.ali.assigned.simple.clean.c10.l80.fasta

The first sequence record of results/wolf.ali.assigned.simple.clean.c10.l80.fasta is:

>HELIUM_000100422_612GNAAXX:7:22:2603:18023#0/1_sub[28..127] {"count":12182,"merged_sample":{"15a_F730814":7559,"29a_F260619":4623},"obiclean_head":true,"obiclean_headcount":2,"obiclean_internalcount":0,"obiclean_samplecount":2,"obiclean_singletoncount":0,"obiclean_status":{"15a_F730814":"h","29a_F260619":"h"},"obiclean_weight":{"15a_F730814":9165,"29a_F260619":6275}}
ttagccctaaacacaagtaattaatataacaaaattattcgccagagtactaccggcaat
agcttaaaactcaaaggacttggcggtgctttataccctt

At that time in the data cleanning we have conserved :

obicount results/wolf.ali.assigned.simple.clean.c10.l80.fasta

time="2023-02-23T18:43:47+01:00" level=info msg="Appending results/wolf.ali.assigned.simple.clean.c10.l80.fasta file\n"
 26 31337 2585

3.2.6 Taxonomic assignment of sequences

Once denoising has been done, the next step in diet analysis is to assign the barcodes to the corresponding species in order to get the complete list of species associated to each sample.

Taxonomic assignment of sequences requires a reference database compiling all possible species to be identified in the sample. Assignment is then done based on sequence comparison between sample sequences and reference sequences.

Download the taxonomy

It is always possible to download the complete taxonomy from NCBI using the following commands.

mkdir TAXO
cd TAXO
curl http://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz \
   | tar -zxvf -
cd ..

For people have a low speed internet connection, a copy of the taxdump.tar.gz file is provided in the wolf_data directory. The NCBI taxonomy is dayly updated, but the one provided here is ok for running this tutorial.

To build the TAXO directory from the provided taxdump.tar.gz, you need to execute the following commands

mkdir TAXO
cd TAXO
tar zxvf wolf_data/taxdump.tar.gz
cd ..

Build a reference database

One way to build the reference database is to use the obipcr program to simulate a PCR and extract all sequences from a general purpose DNA database such as genbank or EMBL that can be amplified in silico by the two primers (here TTAGATACCCCACTATGC and TAGAACAGGCTCCTCTAG) used for PCR amplification.

The two steps to build this reference database would then be

Today, the easiest database to download is Genbank. But this will take you more than a day and occupy more than half a terabyte on your hard drive. In the wolf_data directory, a shell script called download_gb.sh is provided to perform this task. It requires that the programs wget2 and curl are available on your computer.

Use obipcr to simulate amplification and build a reference database based on the putatively amplified barcodes and their recorded taxonomic information.

As these steps can take a long time (about a day for the download and an hour for the PCR), we already provide the reference database produced by the following commands so you can skip its construction. Note that as the Genbank and taxonomic database evolve frequently, if you run the following commands you may get different results.

Download the sequences

mkdir genbank
cd genbank
../wolf_data/install_gb.sh
cd ..

DO NOT RUN THIS COMMAND EXCEPT IF YOU ARE REALLY CONSIENT OF THE TIME AND DISK SPACE REQUIRED.

Use obipcr to simulate an in silico` PCR

obipcr -t TAXO -e 3 -l 50 -L 150 \
       --forward TTAGATACCCCACTATGC \
       --reverse TAGAACAGGCTCCTCTAG \
       --no-order \
       genbank/Release-251/gb*.seq.gz
       > results/v05.pcr.fasta

Note that the primers must be in the same order both in wolf_diet_ngsfilter.txt and in the obipcr command. The part of the path indicating the Genbank release can change. Please check in your genbank directory the exact name of your release.

Clean the database

filter sequences so that they have a good taxonomic description at the species, genus, and family levels (obigrep command command below).

remove redundant sequences (obiuniq command below).

ensure that the dereplicated sequences have a taxid at the family level (obigrep command below).

ensure that sequences each have a unique identification (obiannotate command below)

obigrep -t TAXO \
          --require-rank species \
          --require-rank genus \
          --require-rank family \
          results/v05.ecopcr > results/v05_clean.fasta

obiuniq -c taxid \
        results/v05_clean.fasta \
        > results/v05_clean_uniq.fasta

obirefidx -t TAXO results/v05_clean_uniq.fasta \
        > results/v05_clean_uniq.indexed.fasta

Warning

From now on, for the sake of clarity, the following commands will use the filenames of the files provided with the tutorial. If you decided to run the last steps and use the files you have produced, you'll have to use results/v05_clean_uniq.indexed.fasta instead of wolf_data/db_v05_r117.indexed.fasta.

3.2.7 Assign each sequence to a taxon

Once the reference database is built, taxonomic assignment can be carried out using the obitag command.

obitag -t TAXO -R wolf_data/db_v05_r117.indexed.fasta \
       results/wolf.ali.assigned.simple.clean.c10.l80.fasta \
       > results/wolf.ali.assigned.simple.clean.c10.l80.taxo.fasta

The obitag adds several attributes in the sequence record header, among them:

obitag_bestmatch=ACCESSION where ACCESSION is the id of hte sequence in the reference database that best aligns to the query sequence;

obitag_bestid=FLOAT where FLOAT*100 is the percentage of identity between the best match sequence and the query sequence;

taxid=TAXID where TAXID is the final assignation of the sequence by obitag

scientific_name=NAME where NAME is the scientific name of the assigned taxid.

The first sequence record of wolf.ali.assigned.simple.clean.c10.l80.taxo.fasta is:

>HELIUM_000100422_612GNAAXX:7:81:18704:12346#0/1_sub[28..126] {"count":88,"merged_sample":{"26a_F040644":88},"obiclean_head":true,"obiclean_headcount":1,"obiclean_internalcount":0,"obiclean_samplecount":1,"obiclean_singletoncount":0,"obiclean_status":{"26a_F040644":"h"},"obiclean_weight":{"26a_F040644":208},"obitag_bestid":0.9207920792079208,"obitag_bestmatch":"AY769263","obitag_difference":8,"obitag_match_count":1,"obitag_rank":"clade","scientific_name":"Boreoeutheria","taxid":1437010}
ttagccctaaacataaacattcaataaacaagaatgttcgccagaggactactagcaata
gcttaaaactcaaaggacttggcggtgctttatatccct

3.2.8 Generate the final result table

Some unuseful attributes can be removed at this stage.

obiclean_head

obiclean_headcount

obiclean_internalcount

obiclean_samplecount

obiclean_singletoncount

obiannotate  --delete-tag=obiclean_head \
             --delete-tag=obiclean_headcount \
             --delete-tag=obiclean_internalcount \
             --delete-tag=obiclean_samplecount \
             --delete-tag=obiclean_singletoncount \
  results/wolf.ali.assigned.simple.clean.c10.l80.taxo.fasta \
  > results/wolf.ali.assigned.simple.clean.c10.l80.taxo.ann.fasta

The first sequence record of wolf.ali.assigned.simple.c10.l80.clean.taxo.ann.fasta is then:

>HELIUM_000100422_612GNAAXX:7:84:16335:5083#0/1_sub[28..126] {"count":96,"merged_sample":{"26a_F040644":11,"29a_F260619":85},"obiclean_status":{"26a_F040644":"s","29a_F260619":"h"},"obiclean_weight":{"26a_F040644":14,"29a_F260619":110},"obitag_bestid":0.9595959595959596,"obitag_bestmatch":"AC187326","obitag_difference":4,"obitag_match_count":1,"obitag_rank":"subspecies","scientific_name":"Canis lupus familiaris","taxid":9615}
ttagccctaaacataagctattccataacaaaataattcgccagagaactactagcaaca
gattaaacctcaaaggacttggcagtgctttatacccct

3.2.9 Looking at the data in R

library(ROBIFastread)
library(vegan)

Le chargement a nécessité le package : permute

Le chargement a nécessité le package : lattice

This is vegan 2.6-4

library(magrittr)

diet_data <- read_obifasta("results/wolf.ali.assigned.simple.clean.c10.l80.taxo.fasta")
diet_data %<>% extract_features("obitag_bestmatch","obitag_rank","scientific_name",'taxid')

diet_tab <- extract_readcount(diet_data,key="obiclean_weight")
diet_tab

4 x 26 sparse Matrix of class "dgCMatrix"

  [[ suppressing 26 column names 'HELIUM_000100422_612GNAAXX:7:9:9623:15395#0/1_sub[28..127]', 'HELIUM_000100422_612GNAAXX:7:20:15729:20493#0/1_sub[28..126]', 'HELIUM_000100422_612GNAAXX:7:25:11714:14251#0/1_sub[28..126]' ... ]]

29a_F260619 353 16  . 25 13     .   1 .  .  .  . 391  . 110  .    .  .  .  .  .
26a_F040644   .  . 14  .  . 12830 481 . 31  . 15   . 52  14 17    . 88  .  . 72
13a_F730603   .  .  .  .  1     .  19 9  . 25  .   .  .   .  . 8409  . 15 22  .
15a_F730814   .  .  .  .  .     .   5 4  .  .  .   .  .   .  .    .  .  .  .  .

29a_F260619  . 44  . 6139   .  .
26a_F040644  .  . 18    . 208 43
13a_F730603 20  .  .    .   .  .
15a_F730814  .  .  . 9165   .  .

This file contains 26 sequences. You can deduce the diet of each sample:

13a_F730603: Cervus elaphus

15a_F730814: Capreolus capreolus

26a_F040644: Marmota sp. (according to the location, it is Marmota marmota)

29a_F260619: Capreolus capreolus

Note that we also obtained a few wolf sequences although a wolf-blocking oligonucleotide was used.

4 Specifying the data input to OBITools commands

4.1 Specifying input format

Five sequence formats are accepted for input files. Fasta (Section 2.1.2) and Fastq (Section 2.1.3) are the main ones, EMBL and Genbank allow the use of flat files produced by these two international databases. The last one, ecoPCR, is maintained for compatibility with previous OBITools and allows to read ecoPCR outputs as sequence files.

--ecopcr : Read data following the ecoPCR output format.

--embl Read data following the EMBL flatfile format.

--genbank Read data following the Genbank flatfile format.

Several encoding schemes have been proposed for quality scores in Fastq format. Currently, OBITools considers Sanger encoding as the standard. For reasons of compatibility with older datasets produced with Solexa sequencers, it is possible, by using the following option, to force the use of the corresponding quality encoding scheme when reading these older files.

--solexa Decodes quality string according to the Solexa specification. (default: false)

5 Controling OBITools outputs

5.1 Specifying output format

Only two output sequence formats are supported by OBITools, Fasta and Fastq. Fastq is used when output sequences are associated with quality information. Otherwise, Fasta is the default format. However, it is possible to force the output format by using one of the following two options. Forcing the use of Fasta results in the loss of quality information. Conversely, when the Fastq format is forced with sequences that have no quality data, dummy qualities set to 40 for each nucleotide are added.

--fasta-output Read data following the ecoPCR output format.

--fastq-output Read data following the EMBL flatfile format.

OBITools allows multiple input files to be specified for a single command.

--no-order When several input files are provided, indicates that there is no order among them. (default: false). Using such option can increase a lot the processing of the data.

5.2 The Fasta and Fastq annotations format

OBITools extend the Fasta and Fastq formats by introducing a format for the title lines of these formats allowing to annotate every sequence. While the previous version of OBITools used an ad-hoc format for these annotation, this new version introduce the usage of the standard JSON format to store them.

On input, OBITools automatically recognize the format of the annotations, but two options allows to force the parsing following one of them. You should normally not need to use these options.

--input-OBI-header FASTA/FASTQ title line annotations follow OBI format. (default: false)

--input-json-header FASTA/FASTQ title line annotations follow json format. (default: false)

On output, by default annotation are formatted using the new JSON format. For compatibility with previous version of OBITools and with external scripts and software, it is possible to force the usage of the previous OBITools format.

--output-OBI-header|-O output FASTA/FASTQ title line annotations follow OBI format. (default: false)

--output-json-header output FASTA/FASTQ title line annotations follow json format. (default: false)

6 Options common to most of the OBITools commands

6.1 Helpful options

--help, -h

Display a friendly help message.

--no-progressbar

6.2 System related options

Managing parallel execution of tasks

A new feature of OBITools V4 is the ability to run multiple tasks in parallel, reading files, calculating on the data, formatting and writing the results. Each of these tasks can itself be parallelized by dividing the data into batches and running the calculation on several batches in parallel. This allows the overall calculation time of an OBITools command to be reduced considerably. The parameters organizing the parallel calculation are determined automatically to use the maximum capacity of your computer. But in some circumstances, it is necessary to override these default settings either to try to optimize the computation on a given machine, or to limit the OBITools to using only a part of the computational capacity. There are two options for doing this.

--max-cpu

OBITools V4 are able to run in parallel on all the CPU cores available on the computer. It is sometime required to limit the computation to a smaller number of cores. That option specify the maximum number of cores that the OBITools command can use. This behaviour can also be set up using the OBIMAXCPU environment variable.

--workers, -w

If your computer has 8 cores, but you want to limit OBITools to use only two of them you have several solution:

If you want to set the limit for a single execution you can use the –max-cpu option

obiconvert --max-cpu 2 --fasta-output data.fastq > data.fasta

or you can precede the command by setting the environment variable OBIMAXCPU

OBIMAXCPU=2 obiconvert --fasta-output data.fastq > data.fasta

If you want to set the limit to your complete session, you have to export OBIMAXCPU

export OBIMAXCPU=2

all the following OBITools commands will be limited to use at max 2 CPU cores.

If all the time you want to impose this limit, you must include the above export command in your .bashrc file.

OBITools debuging related options

--debug

7 OBITools expression language

Several OBITools (e.g. obigrep, obiannotate) allow the user to specify some simple expressions to compute values or define predicates. This expressions are parsed and evaluated using the gval go package, which allows for evaluating go-Like expression.

7.1 Variables usable in the expression

sequence is the sequence object on which the expression is evaluated.

annotationsis a map object containing every annotations associated to the currently processed sequence.

7.2 Function defined in the language

Instrospection functions

len(x)

It is a generic function allowing to retreive the size of a object. It returns the length of a sequences, the number of element in a map like annotations, the number of elements in an array. The reurned value is an int.

contains(map,key)

Tests if the map contains a value assciated to key

Cast functions

int(x)

Converts if possible the x value to an integer value. The function returns an int.

numeric(x)

Converts if possible the x value to a float value. The function returns a float.

bool(x)

Converts if possible the x value to a boolean value. The function returns a bool.

String related functions

printf(format,...)

Allows to combine several values to build a string. format follows the classical C printf syntax. The function returns a string.

subspc(x)

substitutes every space in the x string by the underscore (_) character. The function returns a string.

Condition function

ifelse(condition,val1,val2)

The condition value has to be a bool value. If it is true the function returns val1, otherwise, it is returning val2.

7.2.1 Sequence analysis related function

composition(sequence)

The nucleotide composition of the sequence is returned as as map indexed by a, c, g, or t and each value is the number of occurrences of that nucleotide. A fifth key others accounts for all others symboles.

gcskew(sequence)

Computes the excess of g compare to c of the sequence, known as the GC skew.

SkewGC=G−CG+C Skew_{GC}=\frac{G-C}{G+C}

7.3 Accessing to the sequence annotations

The annotations variable is a map object containing all the annotations associated to the currently processed sequence. Index of the map are the attribute names. It exists to possibillities to retreive an annotation. It is possible to use the classical [] indexing operator, putting the attribute name quoted by double quotes between them.

annotations["direction"]

The above code retreives the direction annotation. A second notation using the dot (.) is often more convenient.

annotations.direction

Special attributes of the sequence are accessible only by dedicated methods of the sequence object.

The sequence identifier : Id()

THe sequence definition : Definition()

sequence.Id()

8 Metabarcode design and quality assessment

8.1 obipcr

Replace the ecoPCR original OBITools

9 File format conversions

9.1 obiconvert

10 Sequence annotations

10.1 obiannotate

10.2 obitag

10.3 obitagpcr

11 Computations on sequences

11.1 obipairing

Replace the illuminapairedends original OBITools

Alignment procedure

obipairing is introducing a new alignment algorithm compared to the illuminapairedend command of the OBITools V2. Nethertheless this new algorithm has been design to produce the same results than the previous, except in very few cases.

The new algorithm is a two-step procedure. First, a FASTN-type algorithm (Lipman and Pearson 1985) identifies the best offset between the two matched readings. This identifies the region of overlap.

In the second step, the matching regions of the two reads are extracted along with a flanking sequence of Δ\Delta base pairs. The two subsequences are then aligned using a “one side free end-gap” dynamic programming algorithm. This latter step is only called if at least one mismatch is detected by the FASTP step.

Unless the similarity between the two reads at their overlap region is very low, the addition of the flanking regions in the second step of the alignment ensures the same alignment as if the dynamic programming alignment was performed on the full reads.

The scoring system

In the dynamic programming step, the match and mismatch scores take into account the quality scores of the two aligned nucleotides. By taking these into account, the probability of a true match can be calculated for each aligned base pair.

If we consider a nucleotide read with a quality score QQ, the probability of misreading this base (PEP_E) is : PE=10−Q10 P_E = 10^{-\frac{Q}{10}}

Thus, when a given nucleotide XX is observed with the quality score QQ. The probability that XX is really an XX is :

P(X=X)=1−PE P(X=X) = 1 - P_E

Otherwise, XX is actually one of the three other possible nucleotides (XE1X_{E1}, XE2X_{E2} or XE3X_{E3}). If we suppose that the three reading error have the same probability :

P(X=XE1)=P(X=XE3)=P(X=XE3)=PE3 P(X=X_{E1}) = P(X=X_{E3}) = P(X=X_{E3}) = \frac{P_E}{3}

At each position in an alignment where the two nucleotides X1X_1 and X2X_2 face each other (not a gapped position), the probability of a true match varies depending on whether X1=X2X_1=X_2, an observed match, or X1≠X2X_1 \neq X_2, an observed mismatch.

Probability of a true match when X1=X2X_1=X_2

That probability can be divided in two parts. First X1X_1 and X2X_2 have been correctly read. The corresponding probability is :

PTM=(1−PE1)(1−PE2)=(1−10−Q110)(1−10−Q210) \begin{aligned} P_{TM} &= (1- PE_1)(1-PE_2)\\ &=(1 - 10^{-\frac{Q_1}{10} } )(1 - 10^{-\frac{Q_2}{10}} ) \end{aligned}

Secondly, a match can occure if the true nucleotides read as X1X_1 and X2X_2 are not X1X_1 and X2X_2 but identical.

P(X1==XE1)∩P(X2==XE1)=PE1PE29P(X1==XEx)∩P(X2==XEx)=PE1PE23 \begin{aligned} P(X_1==X_{E1}) \cap P(X_2==X_{E1}) &= \frac{P_{E1} P_{E2}}{9} \\ P(X_1==X_{Ex}) \cap P(X_2==X_{Ex}) & = \frac{P_{E1} P_{E2}}{3} \end{aligned}

The probability of a true match between X1X_1 and X2X_2 when X1=X2X_1 = X_2 an observed match :

P(MATCH|X1=X2)=(1−PE1)(1−PE2)+PE1PE23 \begin{aligned} P(MATCH | X_1 = X_2) = (1- PE_1)(1-PE_2) + \frac{P_{E1} P_{E2}}{3} \end{aligned}

Probability of a true match when X1≠X2X_1 \neq X_2

That probability can be divided in three parts.

X1X_1 has been correctly read and X2X_2 is a sequencing error and is actually equal to X1X_1. Pa=(1−PE1)PE23 P_a = (1-P_{E1})\frac{P_{E2}}{3}

X2X_2 has been correctly read and X1X_1 is a sequencing error and is actually equal to X2X_2. Pb=(1−PE2)PE13 P_b = (1-P_{E2})\frac{P_{E1}}{3}

X1X_1 and X2X_2 corresponds to sequencing error but are actually the same base XExX_{Ex} Pc=2PE1PE29 P_c = 2\frac{P_{E1} P_{E2}}{9}

Consequently : P(MATCH|X1≠X2)=(1−PE1)PE23+(1−PE2)PE13+2PE1PE29 \begin{aligned} P(MATCH | X_1 \neq X_2) = (1-P_{E1})\frac{P_{E2}}{3} + (1-P_{E2})\frac{P_{E1}}{3} + 2\frac{P_{E1} P_{E2}}{9} \end{aligned}

Probability of a match under the random model

The second considered model is a pure random model where every base is equiprobable, hence having a probability of occurrence of a nucleotide equals 0.250.25. Under that hypothesis

P(MATCH|Random model)=0.25 P(MATCH | \text{Random model}) = 0.25

The score is a log ration of likelyhood

Score is define as the logarithm of the ratio between the likelyhood of the observations considering the sequencer error model over tha likelyhood u

Evolution of the match and mismatch scores when the quality of base is 20 while the second range from 10 to 40.

11.2 obimultiplex

Replace the ngsfilter original OBITools

11.3 obicomplement

11.4 obiclean

11.5 obiuniq

12 Sequence sampling and filtering

12.1 obigrep – filters sequence files according to numerous conditions

The obigrep command is somewhat analogous to the standard Unix grep command. It selects a subset of sequence records from a sequence file. A sequence record is a complex object consisting of an identifier, a set of attributes (a key, defined by its name, associated with a value), a definition, and the sequence itself. Instead of working text line by text line like the standard Unix tool, obigrep selection is done sequence record by sequence record. A large number of options allow you to refine the selection on any element of the sequence. obigrep allows you to specify multiple conditions simultaneously (which take on the value TRUE or FALSE) and only those sequence records which meet all conditions (all conditions are TRUE) are selected. obigrep is able to work on two paired read files. The selection criteria apply to one or the other of the readings in each pair depending on the mode chosen (--paired-mode option). In all cases the selection is applied in the same way to both files, thus maintaining their consistency.

12.1.1 The options usable with obigrep

12.1.1.1 Selecting sequences based on their caracteristics

Sequences can be selected on several of their caracteristics, their length, their id, their sequence. Options allow for specifying the condition if selection.

Selection based on the sequence

Sequence records can be selected according if they match or not with a pattern. The simplest pattern is as short sequence (e.g AACCTT). But the usage of regular patterns allows for looking for more complex pattern. As example, A[TG]C+G matches a A, followed by a T or a G, then one or several C and endly a G.

--sequence|-s PATTERN

Regular expression pattern to be tested against the sequence itself. The pattern is case insensitive. A complete description of the regular pattern grammar is available here.

Examples:

Selects only the sequence records that contain an EcoRI restriction site.

obigrep -s 'GAATTC' seq1.fasta > seq2.fasta

: Selects only the sequence records that contain a stretch of at least 10 A.

obigrep -s 'A{10,}' seq1.fasta > seq2.fasta

: Selects only the sequence records that do not contain ambiguous nucleotides.

obigrep -s '^[ACGT]+$' seq1.fasta > seq2.fasta

--min-count | -c COUNT

only sequences reprensenting at least COUNT reads will be selected. That option rely on the count attribute. If the count attribute is not defined for a sequence record, it is assumed equal to 11.

--max-count | -C COUNT

only sequences reprensenting no more than COUNT reads will be selected. That option rely on the count attribute. If the count attribute is not defined for a sequence record, it is assumed equal to 11.

Examples

Selecting sequence records representing at least five reads in the dataset.

obigrep -c 5 data_SPER01.fasta > data_norare_SPER01.fasta

13 Utilities

13.1 obicount

obicount counts the number of sequence records, the sum of the count attributes, and the sum of the length of all the sequences.

Example:

obicount seq.fasta

Prints the number of sequence records contained in the seq.fasta file and the sum of their count attributes.

Options specific to the command

--reads|-r Prints read counts.

--symbols|-s Prints symbol counts.

--variants|-v Prints variant counts.

13.2 obidistribute

13.3 obifind

Replace the ecofind original OBITools.

Appendix A — Annexes

A.1 Sequence attributes

ali_dir (string)

Set by the obipairing tool

The attribute can contain 2 string values left or right.

The alignment generated by obipairing is a 3’-end gap free algorithm. Two cases can occur when aligning the forward and reverse reads. If the barcode is long enough, both the reads overlap only on their 3’ ends. In such case, the alignment direction ali_dir is set to left. If the barcode is shorter than the read length, the paired reads overlap by their 5’ ends, and the complete barcode is sequenced by both the reads. In that later case, ali_dir is set to right.

ali_length (int)

Set by the obipairing tool

Length of the aligned parts when merging forward and reverse reads

count (int)

Set by the obiuniq tool

Getter : method Count()

Setter : method SetCount(int)

The count attribute indicates how-many strictly identical reads have been merged in a single record. It contains an integer value. If it is absent this means that the sequence record represents a single occurrence of the sequence.

The Count() method allows to access to the count attribute as an integer value. If the count attribute is not defined for the given sequence, the value 1 is returned

merged_* (map[string]int)

Set by the obiuniq tool

The -m option of the obiuniq tools allows for keeping track of the distribution of the values stored in given attribute of interest. Often this option is used to summarise distribution of a sequence variant accross samples when obiuniq is run after running obimultiplex. The actual name of the attribute depends on the name of the monitored attribute. If -m option is used with the attribute sample, then this attribute names merged_sample.

mode (string)

Set by the obipairing tool

The attribute can contain 2 string values join or alignment.

obitag_ref_index (map[string]string)

Set by the obirefidx tool.

It resumes to which taxonomic annotation a match to that sequence must lead according to the number of differences existing between the query sequence and the reference sequence having that tag.

   {"0":"9606@Homo sapiens@species",
    "2":"207598@Homininae@subfamily",
    "3":"9604@Hominidae@family",
    "8":"314295@Hominoidea@superfamily",
    "10":"9526@Catarrhini@parvorder",
    "12":"1437010@Boreoeutheria@clade",
    "16":"9347@Eutheria@clade",
    "17":"40674@Mammalia@class",
    "22":"117571@Euteleostomi@clade",
    "25":"7776@Gnathostomata@clade",
    "29":"33213@Bilateria@clade",
    "30":"6072@Eumetazoa@clade"}

pairing_mismatches (map[string]string)

Set by the obipairing tool

seq_a_single (int)

Set by the obipairing tool

seq_ab_match (int)

Set by the obipairing tool

seq_b_single (int)

Set by the obipairing tool

score (int)

Set by the obipairing tool

score_norm (float)

Set by the obipairing tool

The value ranges between 0 and 1.

Score of the alignment between forward and reverse reads expressed as a fraction of identity.

References

Boyer, Frédéric, Céline Mercier, Aurélie Bonin, Yvan Le Bras, Pierre Taberlet, and Eric Coissac. 2016. “obitools: a unix-inspired software package for DNA metabarcoding.” Molecular Ecology Resources 16 (1): 176–82. https://doi.org/10.1111/1755-0998.12428.

Cock, Peter JA, Christopher J Fields, Naohisa Goto, Michael L Heuer, and Peter M Rice. 2010. “The Sanger FASTQ File Format for Sequences with Quality Scores, and the Solexa/Illumina FASTQ Variants.” Nucleic Acids Research 38 (6): 1767–71.

Lipman, D J, and W R Pearson. 1985. “Rapid and sensitive protein similarity searches.” Science 227 (4693): 1435–41. http://www.ncbi.nlm.nih.gov/pubmed/2983426.

Riaz, Tiayyba, Wasim Shehzad, Alain Viari, François Pompanon, Pierre Taberlet, and Eric Coissac. 2011. “ecoPrimers: inference of new DNA barcode markers from whole genome sequence analysis.” Nucleic Acids Research 39 (21): e145. https://doi.org/10.1093/nar/gkr732.

Seguritan, V, and F Rohwer. 2001. “FastGroup: a program to dereplicate libraries of 16S rDNA sequences.” BMC Bioinformatics 2 (October): 9. https://doi.org/10.1186/1471-2105-2-9.

Shehzad, Wasim, Tiayyba Riaz, Muhammad A Nawaz, Christian Miquel, Carole Poillot, Safdar A Shah, Francois Pompanon, Eric Coissac, and Pierre Taberlet. 2012. “Carnivore diet analysis based on next-generation sequencing: Application to the leopard cat (Prionailurus bengalensis) in Pakistan.” Molecular Ecology 21 (8): 1951–65. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-294X.2011.05424.x.