Commit 2f7ac4b6 by Frédéric Boyer

Relecture d'Aurélie

parent a87b0f73
Wolves' diet based on DNA metabarcoding
==================================================================================
=======================================
How to analyze DNA metabarcoding data produced on Illumina sequencers using:
Here is a tutorial on how to analyze DNA metabarcoding data produced on Illumina
sequencers using:
- the *OBITools*
- some basic *Unix* commands
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.
+-------------------------------------------------------------+
| Good to remember: I am working with tons of sequences |
+-------------------------------------------------------------+
| It is always a good idea to have a look at the intermediate |
| results or to evaluate the best parameter of each step. |
| results or to evaluate the best parameter for each step. |
| Some commands are designed for that purpose, for example |
| use : |
| you can use : |
| |
| - :doc:`obicount <scripts/obicount>` to count for the number|
| - :doc:`obicount <scripts/obicount>` to count the number |
| of sequence records in a file |
| - :doc:`obihead <scripts/obihead>` and |
| - :doc:`obihead <scripts/obihead>` and |
| :doc:`obitail <scripts/obitail>` to view the first |
| or last sequence records of a file |
| - :doc:`obistat <scripts/obistat>` to get some basic |
| statistics (count, mean, standard deviation) on the |
| attributes (key=value couples) in the fasta header of each|
| attributes (key=value combinations) in the header of each |
| sequence record (see The `extended OBITools fasta format` |
| in the :doc:`fasta format <fasta>` description) |
| - any *Unix* command such as ``less``, ``awk``, ``sort``, |
......@@ -37,7 +44,7 @@ Data
The data needed to run the tutorial are the following:
- the :doc:`fastq <fastq>` files resulting of a GA IIx (Illumina) paired-end (2 x 108 bp)
- :doc:`fastq <fastq>` files resulting of a GA IIx (Illumina) paired-end (2 x 108 bp)
sequencing assay of DNA extracted and amplified from
four wolf faeces:
......@@ -47,25 +54,22 @@ The data needed to run the tutorial are the following:
- 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 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 fasta format:
- 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 :doc:`ecoPCR <scripts/ecoPCR>`
This reference database has been extracted from the release 117 of EMBL using
:doc:`ecoPCR <scripts/ecoPCR>`
- the NCBI taxonomy formatted in the :doc:`ecoPCR <scripts/ecoPCR>` format (see the :doc:`obiconvert <scripts/obiconvert>` utility for details) :
- the NCBI taxonomy formatted in the :doc:`ecoPCR <scripts/ecoPCR>` format (see the
:doc:`obiconvert <scripts/obiconvert>` utility for details) :
* ``embl_r117.ndx``
* ``embl_r117.rdx``
* ``embl_r117.tdx``
These data 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; Riaz et al. 2011),
and with a wolf blocking oligonucleotide.
Step by step analysis
---------------------
......@@ -74,11 +78,11 @@ Step by step analysis
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 from these two forward
and reverse 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 fragments are *at the same line position* in the two files.
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
:doc:`illuminapairedend <scripts/illuminapairedend>` utility that aligns the two reads
and returns the reconstructed sequence.
......@@ -89,90 +93,100 @@ In our case, the command is:
> illuminapairedend --score-min=40 -r wolf_R.fastq wolf_F.fastq > wolf.fastq
The :py:mod:`--score-min` option allows to avoid returning badly aligned sequence. If the alignment score is below 40, the
forward and reverse reads are not aligned but concatenated, and the value of the :py:mod:`mode` attribute is set to :py:mod:`joined`
instead of :py:mod:`alignment`
The :py:mod:`--score-min` option allows discarding sequences with low alignment quality.
If the alignment score is below 40, the forward and reverse reads are not aligned but
concatenated, and the value of the :py:mod:`mode` attribute in the sequence header is set
to :py:mod:`joined` instead of :py:mod:`alignment`
Remove not aligned sequence records
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Remove unaligned sequence records
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The not aligned sequence records (:py:mod:`mode=joined`) cannot be used. The following command allows to
remove them from the dataset:
Unaligned sequences (:py:mod:`mode=joined`) cannot be used. The following command allows
removing them from the dataset:
.. code-block:: bash
> obigrep -p 'mode!="joined"' wolf.fastq > wolf.ali.fastq
The :py:mod:`-p` requires a *python* expression. :py:mod:`mode!="joined"` means that if the value of
the :py:mod:`mode` attribute is different from :py:mod:`joined`, the corresponding sequence record will be kept.
The :py:mod:`-p` requires a *python* expression. :py:mod:`mode!="joined"` means that if
the value of the :py:mod:`mode` attribute is different from :py:mod:`joined`, the
corresponding sequence record will be kept.
The first sequence record of ``wolf.ali.fastq`` can be obtained using the following command line:
The first sequence record of ``wolf.ali.fastq`` can be obtained using the following
command line:
.. code-block:: bash
> obihead --without-progress-bar -n 1 wolf.ali.fastq
> obihead --without-progress-bar -n 1 wolf.ali.fastq
And the result is:
.. code-block:: bash
@HELIUM_000100422_612GNAAXX:7:119:14871:19157#0/1_CONS ali_length=61;
direction=left; seq_ab_match=47; sminR=40.0; seq_a_mismatch=7; seq_b_deletion=1;
seq_b_mismatch=7; seq_a_deletion=1; score_norm=1.89772607661;
score=115.761290673; seq_a_insertion=0; mode=alignment; sminL=40.0;
seq_a_single=46; seq_b_single=46; seq_b_insertion=0;
ccgcctcctttagataccccactatgcttagccctaaacacaagtaattattataacaaaatcattcgccagagtgtagc
gggagtaggttaaaactcaaaggacttggcggtgctttatacccttctagaggagcctgttctaaggaggcgg
+
ddddddddddddddddddddddcddddcacdddddddddddddc\d~b~~~b~~~~~~b`ryK~|uxyXk`}~ccBccBc
ccBcBcccBcBccccccc~~~~b|~~xdbaddaaWcccdaaddddadacddddddcddadbbddddddddddd
@HELIUM_000100422_612GNAAXX:7:119:14871:19157#0/1_CONS ali_length=61;
direction=left; seq_ab_match=47; sminR=40.0; seq_a_mismatch=7; seq_b_deletion=1;
seq_b_mismatch=7; seq_a_deletion=1; score_norm=1.89772607661;
score=115.761290673; seq_a_insertion=0; mode=alignment; sminL=40.0;
seq_a_single=46; seq_b_single=46; seq_b_insertion=0;
ccgcctcctttagataccccactatgcttagccctaaacacaagtaattattataacaaaatcattcgccagagtgtagc
gggagtaggttaaaactcaaaggacttggcggtgctttatacccttctagaggagcctgttctaaggaggcgg
+
ddddddddddddddddddddddcddddcacdddddddddddddc\d~b~~~b~~~~~~b`ryK~|uxyXk`}~ccBccBc
ccBcBcccBcBccccccc~~~~b|~~xdbaddaaWcccdaaddddadacddddddcddadbbddddddddddd
Assign each sequence record to its sample
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Assign each sequence record to the corresponding sample/marker combination
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Each sequence record is assigned to its original sample and to its experiment by using the information
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 experiment can be indicated 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 are different), the sequence of the
forward primer, the sequence of the reverse primer, the letter ``F`` or ``T`` for sample identification using the forward primer and tag
only or using both primers and both tags, respectively (see :doc:`ngsfilter <scripts/ngsfilter>` for details).
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 :doc:`ngsfilter <scripts/ngsfilter>`
for details).
.. code-block:: bash
> ngsfilter -t wolf_diet_ngsfilter.txt -u unidentified.fastq wolf.ali.fastq > wolf.ali.assigned.fastq
> ngsfilter -t wolf_diet_ngsfilter.txt -u unidentified.fastq wolf.ali.fastq > \
wolf.ali.assigned.fastq
This command creates two files:
- ``unidentified.fastq`` containing all the sequence records that were not assigned to a sample
- ``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
- ``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 were removed. The information concerning the experiment, the sample,
primers and the tags are added as several attributes in the sequence heading.
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
:doc:`ngsfilter <scripts/ngsfilter>` program. Information concerning the experiment,
sample, primers and tags is added as attributes in the sequence header.
The first sequence record of ``wolf.ali.assigned.fastq`` is:
For instance, the first sequence record of ``wolf.ali.assigned.fastq`` is:
.. code-block:: bash
@HELIUM_000100422_612GNAAXX:7:119:14871:19157#0/1_CONS_SUB_SUB status=full;
seq_ab_match=47; sminR=40.0; ali_length=61; tail_quality=67.0;
reverse_match=tagaacaggctcctctag; seq_a_deletion=1; sample=29a_F260619;
forward_match=ttagataccccactatgc; forward_primer=ttagataccccactatgc;
reverse_primer=tagaacaggctcctctag; sminL=40.0; forward_score=72.0;
score=115.761290673; seq_a_mismatch=7; forward_tag=gcctcct; seq_b_mismatch=7;
experiment=wolf_diet; mid_quality=69.4210526316; avg_quality=69.1045751634;
seq_a_single=46; score_norm=1.89772607661; reverse_score=72.0;
direction=forward; seq_b_insertion=0; seq_b_deletion=1; seq_a_insertion=0;
seq_length_ori=153; reverse_tag=gcctcct; seq_length=99; mode=alignment;
head_quality=67.0; seq_b_single=46;
ttagccctaaacacaagtaattattataacaaaatcattcgccagagtgtagcgggagtaggttaaaactcaaaggact
tggcggtgctttataccctt
+
cacdddddddddddddc\d~b~~~b~~~~~~b`ryK~|uxyXk`}~ccBccBcccBcBcccBcBccccccc~~~~b|~~
xdbaddaaWcccdaadddda
@HELIUM_000100422_612GNAAXX:7:119:14871:19157#0/1_CONS_SUB_SUB status=full;
seq_ab_match=47; sminR=40.0; ali_length=61; tail_quality=67.0;
reverse_match=tagaacaggctcctctag; seq_a_deletion=1; sample=29a_F260619;
forward_match=ttagataccccactatgc; forward_primer=ttagataccccactatgc;
reverse_primer=tagaacaggctcctctag; sminL=40.0; forward_score=72.0;
score=115.761290673; seq_a_mismatch=7; forward_tag=gcctcct; seq_b_mismatch=7;
experiment=wolf_diet; mid_quality=69.4210526316; avg_quality=69.1045751634;
seq_a_single=46; score_norm=1.89772607661; reverse_score=72.0;
direction=forward; seq_b_insertion=0; seq_b_deletion=1; seq_a_insertion=0;
seq_length_ori=153; reverse_tag=gcctcct; seq_length=99; mode=alignment;
head_quality=67.0; seq_b_single=46;
ttagccctaaacacaagtaattattataacaaaatcattcgccagagtgtagcgggagtaggttaaaactcaaaggact
tggcggtgctttataccctt
+
cacdddddddddddddc\d~b~~~b~~~~~~b`ryK~|uxyXk`}~ccBccBcccBcBcccBcBccccccc~~~~b|~~
xdbaddaaWcccdaadddda
......@@ -181,11 +195,13 @@ The first sequence record of ``wolf.ali.assigned.fastq`` is:
Dereplicate reads into uniq sequences
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The same DNA molecule can have been 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 uniq *sequences* instead of *reads*. To *dereplicate* such *reads* into uniq *sequences*, we use the :doc:`obiuniq <scripts/obiuniq>` command.
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 :doc:`obiuniq <scripts/obiuniq>` command.
+-------------------------------------------------------------+
| Definition: Dereplicate reads into uniq sequences |
| Definition: Dereplicate reads into unique sequences |
+-------------------------------------------------------------+
| 1. compare all the reads in a data set to each other |
| 2. group strictly identical reads together |
......@@ -197,33 +213,42 @@ it is convenient to work with uniq *sequences* instead of *reads*. To *dereplica
+-------------------------------------------------------------+
We use the :doc:`obiuniq <scripts/obiuniq>` command with the `-m sample`. The `-m sample` option is used to keep the information of the samples of origin for each uniq sequence.
For dereplication, we use the :doc:`obiuniq <scripts/obiuniq>` command with the `-m
sample`. The `-m sample` option is used to keep the information of the samples of origin
for each unique sequence.
.. code-block:: bash
> obiuniq -m sample wolf.ali.assigned.fastq > wolf.ali.assigned.uniq.fasta
Note that :doc:`obiuniq <scripts/obiuniq>` returns a fasta file.
The first sequence record of ``wolf.ali.assigned.uniq.fasta`` is:
.. code-block:: bash
>HELIUM_000100422_612GNAAXX:7:119:14871:19157#0/1_CONS_SUB_SUB_CMP ali_length=61; seq_ab_match=47;
sminR=40.0; tail_quality=67.0; reverse_match=ttagataccccactatgc; seq_a_deletion=1;
forward_match=tagaacaggctcctctag; forward_primer=tagaacaggctcctctag; reverse_primer=ttagataccccactatgc;
sminL=40.0; merged_sample={'29a_F260619': 1}; forward_score=72.0; seq_a_mismatch=7; forward_tag=gcctcct;
seq_b_mismatch=7; score=115.761290673; mid_quality=69.4210526316; avg_quality=69.1045751634;
seq_a_single=46; score_norm=1.89772607661; reverse_score=72.0; direction=reverse; seq_b_insertion=0;
experiment=wolf_diet; seq_b_deletion=1; seq_a_insertion=0; seq_length_ori=153; reverse_tag=gcctcct;
count=1; seq_length=99; status=full; mode=alignment; head_quality=67.0; seq_b_single=46;
>HELIUM_000100422_612GNAAXX:7:119:14871:19157#0/1_CONS_SUB_SUB_CMP ali_length=61;
seq_ab_match=47; sminR=40.0; tail_quality=67.0; reverse_match=ttagataccccactatgc;
seq_a_deletion=1; forward_match=tagaacaggctcctctag; forward_primer=tagaacaggctcctctag;
reverse_primer=ttagataccccactatgc; sminL=40.0; merged_sample={'29a_F260619': 1};
forward_score=72.0; seq_a_mismatch=7; forward_tag=gcctcct; seq_b_mismatch=7;
score=115.761290673; mid_quality=69.4210526316; avg_quality=69.1045751634;
seq_a_single=46; score_norm=1.89772607661; reverse_score=72.0; direction=reverse;
seq_b_insertion=0; experiment=wolf_diet; seq_b_deletion=1; seq_a_insertion=0;
seq_length_ori=153; reverse_tag=gcctcct; count=1; seq_length=99; status=full;
mode=alignment; head_quality=67.0; seq_b_single=46;
aagggtataaagcaccgccaagtcctttgagttttaacctactcccgctacactctggcg
aatgattttgttataataattacttgtgtttagggctaa
The run of :doc:`obiuniq <scripts/obiuniq>` has added two key=values entries in the header of the fasta sequence :
- :py:mod:`merged_sample={'29a_F260619': 1}` : this sequence have been found once in a single sample
- :py:mod:`count=1` : the total number of counts for this sequence is 1
To keep only these two ``key=value`` informations, we can use the :doc:`obiannotate <scripts/obiannotate>` command:
The run of :doc:`obiuniq <scripts/obiuniq>` has added two key=values entries in the header
of the fasta sequence:
- :py:mod:`merged_sample={'29a_F260619': 1}`: this sequence have been found once in a
single sample called 29a_F260619
- :py:mod:`count=1` : the total count for this sequence is 1
To keep only these two ``key=value`` attributes, we can use the
:doc:`obiannotate <scripts/obiannotate>` command:
.. code-block:: bash
......@@ -232,7 +257,7 @@ To keep only these two ``key=value`` informations, we can use the :doc:`obiannot
wolf.ali.assigned.uniq.fasta > $$ ; mv $$ wolf.ali.assigned.uniq.fasta
The first five sequence records of ``wolf.ali.assigned.uniq.fasta`` becomes:
The first five sequence records of ``wolf.ali.assigned.uniq.fasta`` become:
.. code-block:: bash
......@@ -253,21 +278,25 @@ The first five sequence records of ``wolf.ali.assigned.uniq.fasta`` becomes:
gaacattcttgtttattgaatgtttatgtttagggctaa
Denoising the sequence dataset
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Denoise the sequence dataset
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
To have a set of sequences assigned to their original 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. To remove such sequences as much as possible, we first remove rare sequences
and then remove sequences variants that likely correspond to artifacts.
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. To remove such sequences as much as possible, we
first discard rare sequences and then rsequence variants that likely correspond to
artifacts.
Get the counts statistics
~~~~~~~~~~~~~~~~~~~~~~~~~
Get the count statistics
~~~~~~~~~~~~~~~~~~~~~~~~
In that case, we use :doc:`obistat <scripts/obistat>` to get the counting statistics on the 'count' attribute (the count attribute has been set by the :doc:`obiuniq <scripts/obiuniq>` command). By piping
the result in the *Unix* commands ``sort`` and ``head`` we keep only the counting statistics for the 20 lowest values of the 'count' attributes.
In that case, we use :doc:`obistat <scripts/obistat>` to get the counting statistics on
the 'count' attribute (the count attribute has been added by the :doc:`obiuniq
<scripts/obiuniq>` command). By piping the result in the *Unix* commands ``sort`` and
``head``, we keep only the count statistics for the 20 lowest values of the 'count'
attribute.
.. code-block:: bash
......@@ -306,34 +335,35 @@ The dataset contains 3504 sequences occurring only once.
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 :doc:`obigrep <scripts/obigrep>` command.
The ``-p 'count>=10'`` option means that the ``python`` expression :py:mod:`count>=10` must be
evaluated to :py:mod:`True` 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 vertebrate must have a length arround 100bp.
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 :doc:`obigrep <scripts/obigrep>`
command.
The ``-p 'count>=10'`` option means that the ``python`` expression :py:mod:`count>=10`
must be evaluated to :py:mod:`True` 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.
.. code-block:: bash
> obigrep -l 80 -p 'count>=10' wolf.ali.assigned.uniq.fasta \
> wolf.ali.assigned.uniq.c10.l80.fasta
> wolf.ali.assigned.uniq.c10.l80.fasta
The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.fasta`` is:
.. code-block:: bash
.. code-block:: bash
>HELIUM_000100422_612GNAAXX:7:22:8540:14708#0/1_CONS_SUB_SUB count=12335; merged_sample={'29a_F260619': 4697, '15a_F730814': 7638};
aagggtataaagcaccgccaagtcctttgagttttaagctattgccggtagtactctggc
gaataattttgttatattaattacttgtgtttagggctaa
Clean the sequences for PCR/sequencing errors (sequence variants)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
As a final step of denoising, using the :doc:`obiclean <scripts/obiclean>` we keep the `Head` sequences
(``-H`` option) that are sequences with no variants with greater count or sequences with no variants
with 20-fold greater (``-r 0.05`` option).
As a final denoising step, using the :doc:`obiclean <scripts/obiclean>` program, we keep
the `head` sequences (``-H`` option) that are sequences with no variants with a count
greater than 5% of their own count (``-r 0.05`` option).
.. code-block:: bash
......@@ -347,7 +377,9 @@ The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.fasta`` is:
>HELIUM_000100422_612GNAAXX:7:22:8540:14708#0/1_CONS_SUB_SUB
merged_sample={'29a_F260619': 4697, '15a_F730814': 7638};
obiclean_count={'29a_F260619': 5438, '15a_F730814': 8642}; obiclean_head=True;
obiclean_cluster={'29a_F260619': 'HELIUM_000100422_612GNAAXX:7:22:8540:14708#0/1_CONS_SUB_SUB', '15a_F730814': 'HELIUM_000100422_612GNAAXX:7:22:8540:14708#0/1_CONS_SUB_SUB'};
obiclean_cluster={'29a_F260619':
'HELIUM_000100422_612GNAAXX:7:22:8540:14708#0/1_CONS_SUB_SUB', '15a_F730814':
'HELIUM_000100422_612GNAAXX:7:22:8540:14708#0/1_CONS_SUB_SUB'};
count=12335; obiclean_internalcount=0; obiclean_status={'29a_F260619': 'h', '15a_F730814': 'h'};
obiclean_samplecount=2; obiclean_headcount=2; obiclean_singletoncount=0;
aagggtataaagcaccgccaagtcctttgagttttaagctattgccggtagtactctggc
......@@ -356,39 +388,42 @@ The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.fasta`` is:
Taxonomic assignment of sequences
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Once denoising has been done, the next step in diet analysis is to relate the barcodes to their respective
species in order to get the list of species associated to each sample.
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.
The taxonomic assignement of sequences requires a reference database compiling all possible species to be
identified in the sample. The assignment is then done based on sequence comparisons between the sample
sequences and the reference sequences.
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.
Build a reference database
~~~~~~~~~~~~~~~~~~~~~~~~~~
As a rough way to build the reference database, we use the :doc:`ecoPCR <scripts/ecoPCR>` program to simulate
a PCR and to extract all sequences from the EMBL that may be amplified in silico by the two primers
(`TTAGATACCCCACTATGC` and `TAGAACAGGCTCCTCTAG`) extracted from the samples description used to assign each
read to its sample (file ``wolf_diet_ngsfilter.txt``).
One way to build the reference database is to use the :doc:`ecoPCR <scripts/ecoPCR>`
program to simulate a PCR and to extract all sequences from the EMBL that may be amplified
`in silico` by the two primers (`TTAGATACCCCACTATGC` and `TAGAACAGGCTCCTCTAG`) used for
PCR amplification.
The full list of steps to do in order to build this reference database would then be:
The full list of steps for building this reference database would then be:
1. Download the whole set of EMBL sequences (available from: ftp://ftp.ebi.ac.uk/pub/databases/embl/release/)
2. Download the NCBI taxonomy (available from: ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz)
3. Format them into the ecoPCR format (set :doc:`obiconvert <scripts/obiconvert>` for how you can produce
ecoPCR compatible files)
4. Use ecoPCR to simulate amplification and build a reference database on the basis of putatively
1. Download the whole set of EMBL sequences (available from:
ftp://ftp.ebi.ac.uk/pub/databases/embl/release/)
2. Download the NCBI taxonomy (available from:
ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz)
3. Format them into the ecoPCR format (see :doc:`obiconvert <scripts/obiconvert>` for how
you can produce ecoPCR compatible files)
4. Use ecoPCR to simulate amplification and build a reference database based on putatively
amplified barcodes together with their recorded taxonomic information
As step 1 and step 3 can be really time consuming (about one day) we provide the reference database
produce by the following commands so that you can skip the following steps. Note that as both the EMBL database
and the taxonomic data can evolve dayly, if you run the following commands you may end up with quite different
results.
As step 1 and step 3 can be really time-consuming (about one day), we alredy provide the
reference database produced by the following commands so that you can skip its
construction. Note that as the EMBL database and taxonomic data can evolve daily, if you
run the following commands you may end up with quite different results.
Note that any utility that allows downloading of files from a ftp site can be used. In the following commands,
we use the commonly used ``wget`` *Unix* command.
Any utility allowing file downloading from a ftp site can be used. In the following
commands, we use the commonly used ``wget`` *Unix* command.
Download the sequences
......................
......@@ -419,8 +454,8 @@ Format the data
> obiconvert --embl -t ./TAXO --ecopcrDB-output=embl_last ./EMBL/*.dat
Retrieve the sequences
......................
Use ecoPCR to simulate an in silico` PCR
........................................
.. code-block:: bash
......@@ -428,20 +463,20 @@ Retrieve the sequences
TTAGATACCCCACTATGC TAGAACAGGCTCCTCTAG > v05.ecopcr
Note that the primers must be in the same order both
in ``wolf_diet_ngsfilter.txt`` and in the :doc:`ecoPCR <scripts/ecoPCR>` command.
Note that the primers must be in the same order both in ``wolf_diet_ngsfilter.txt`` and in
the :doc:`ecoPCR <scripts/ecoPCR>` command.
Clean the database
..................
1. filter the sequences so that they have a good taxonomic description at the species,
1. filter sequences so that they have a good taxonomic description at the species,
genus, and family levels (:doc:`obigrep <scripts/obigrep>` command below).
2. remove redundant sequences (:doc:`obiuniq <scripts/obiuniq>` command below).
3. ensure that the dereplicated sequences have a taxid at the family level
(:doc:`obigrep <scripts/obigrep>` command below).
4. ensure that sequences each have a uniq identification (:doc:`obiannotate <scripts/obiannotate>`
command below)
4. ensure that sequences each have a unique identification
(:doc:`obiannotate <scripts/obiannotate>` command below)
.. code-block:: bash
......@@ -458,15 +493,16 @@ Clean the database
.. warning::
From now, for the sake of clarity, the following commands will use the filenames of the provided data.
If you decided to run the last steps and use the files you have produced, you'll have to use
``db_v05.ecopcr`` instead of ``db_v05_r117.ecopcr`` and ``embl_last`` instead of ``embl_r117``
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 ``db_v05.fasta`` instead of
``db_v05_r117.fasta`` and ``embl_last`` instead of ``embl_r117``
Assign each sequence to a taxon
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Once the reference database is built, the taxonomic assignment can be carried out using
Once the reference database is built, taxonomic assignment can be carried out using
the :doc:`ecotag <scripts/ecotag>` command.
.. code-block:: bash
......@@ -475,11 +511,15 @@ the :doc:`ecotag <scripts/ecotag>` command.
wolf.ali.assigned.uniq.c10.l80.clean.tag.fasta
The :doc:`ecotag <scripts/ecotag>` adds several `key=value` attributes, among them are :
The :doc:`ecotag <scripts/ecotag>` adds several `key=value` attributes in the sequence
record header, among them:
- best_match=ACCESSION where ACCESSION is the id of one the sequence in the reference database that best align to the query sequence;
- best_identity=FLOAT where FLOAT*100 is the percentage identity between the best match sequence and the query sequence;
- taxid=TAXID where TAXID is the final assignation of the sequence by :doc:`ecotag <scripts/ecotag>`
- best_match=ACCESSION where ACCESSION is the id of hte sequence in the reference database
that best aligns to the query sequence;
- best_identity=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
:doc:`ecotag <scripts/ecotag>`
- scientific_name=NAME where NAME is the scientific name of the assigned taxid.
The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.fasta`` is:
......@@ -511,20 +551,21 @@ The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.fasta``
Generate the final result table
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Some unuseful attributes can be removed.
Some unuseful attributes can be removed at this stage.
.. code-block:: bash
> obiannotate --delete-tag=scientific_name_by_db --delete-tag=obiclean_samplecount \
--delete-tag=obiclean_count --delete-tag=obiclean_singletoncount \
--delete-tag=obiclean_cluster --delete-tag=obiclean_internalcount \
--delete-tag=obiclean_head --delete-tag=taxid_by_db --delete-tag=obiclean_headcount \
--delete-tag=id_status --delete-tag=rank_by_db --delete-tag=order_name \
--delete-tag=order wolf.ali.assigned.uniq.c10.l80.clean.tag.fasta > \
--delete-tag=obiclean_count --delete-tag=obiclean_singletoncount \
--delete-tag=obiclean_cluster --delete-tag=obiclean_internalcount \
--delete-tag=obiclean_head --delete-tag=taxid_by_db --delete-tag=obiclean_headcount \
--delete-tag=id_status --delete-tag=rank_by_db --delete-tag=order_name \
--delete-tag=order wolf.ali.assigned.uniq.c10.l80.clean.tag.fasta > \
wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.fasta
The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.fasta`` is:
The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.fasta`` is
then:
.. code-block:: bash
......@@ -539,15 +580,15 @@ The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.fast
ttagccctaaacacaagtaattaatataacaaaattattcgccagagtactaccggcaat
agcttaaaactcaaaggacttggcggtgctttataccctt
The sequences can be sorted in decreasing order of `count`.
The sequences can be sorted by decreasing order of `count`.
.. code-block:: bash
> obisort -k count -r wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.fasta > \
wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.sort.fasta
The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.sort.fasta`` is:
The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.sort.fasta`` is then:
.. code-block:: bash
......@@ -562,7 +603,7 @@ The first sequence record of ``wolf.ali.assigned.uniq.c10.l80.clean.tag.ann.sort
ttagccctaaacacaagtaattaatataacaaaattattcgccagagtactaccggcaat
agcttaaaactcaaaggacttggcggtgctttataccctt
Finally, a tab delimited file that can be open by excel or R is generated.
Finally, a tab-delimited file that can be open by excel or R is generated.
.. code-block:: bash
......@@ -575,7 +616,9 @@ This file contains 26 sequences. You can deduce the diet of each sample:
- 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 despite that a wolf blocking primer was used.
Note that we also obtained a few wolf sequences although a wolf-blocking oligonucleotide
was used.
References
......@@ -590,12 +633,12 @@ References
Acids Research, 39, e145.
- Seguritan V, Rohwer F. (2001) FastGroup: a program to dereplicate libraries of
16S rDNA sequences. BMC Bioinformatics. 2001;2:9. Epub 2001 Oct 16.
Contact
-------
For any suggestion and improvement of this tutorial, please contact :
For any suggestion or improvement, please contact :
- eric.coissac@metabarcoding.org
- frederic.boyer@metabarcoding.org
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment