build_bdr.sh

##########################################################################
## Codes for scientific papers related to metabarcoding studies
##
## AUTHORS
## =======
## * Pierre-Edouard Guerin   | pierre-edouard.guerin@cefe.cnrs.fr
## * Virginie Marques        | virginie.marques@etu.umontpellier.fr
## * CNRS/CEFE, CNRS/MARBEC  | Montpellier, France
## * 2018-2020
##
##
## Inspired by the Eric Coissac et al. obitools
## Molecular Ecology Resources 2015
##
## DESCRIPTION
## ===========
##
## Build a reference database for "miseq" sequences data.
## We recommand to launch each command line by line into linux terminal
## to check each step results
##
##
##
## USAGE
## =====
## bash buil_bdr.sh
##
##
##########################################################################

## load an environment with obitools
SINGULARITY_SIMG="/media/superdisk/utils/conteneurs/obitools.simg"
singularity shell --bind /media/superdisk:/media/superdisk $SINGULARITY_SIMG


## configure arguments value
source ./config.sh
## download the sequences
mkdir EMBL
cd EMBL
wget ftp://ftp.ebi.ac.uk/pub/databases/ena/sequence/release/std/*
gzip -d *
cd ..
## download taxonomy
mkdir TAXO
cd TAXO
wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
tar -zxvf taxdump.tar.gz
cd ..

## add mitofish sequences (by default we skip this step)
if [ $MITOFISH == 'y' ]
then
 echo "adding sequences from mitofish..."
 bash scripts/add_sequences_from_mitofish.sh
 obiconvert --skip-on-error --fasta -t ./TAXO --ecopcrdb-output=mitofish/"${rd_prefix}" mitofish/mitogene_12S.fasta
else
 echo "skip adding sequences from mitofish"
fi

## format the data
obiconvert --skip-on-error --embl -t ./TAXO --ecopcrdb-output="${rd_prefix}" EMBL/rel_std_*.dat
## ecoPCR to simulate an in silico PCR
#### 50 :: change to 20 : lost of lamproie
ecoPCR -d "${rd_prefix}" -e "${ecoPCR_e}" -l "${ecoPCR_l}" -L "${ecoPCR_L}" "${primer5}" "${primer3}" > v_"${rd_prefix}".ecopcr
## clean the database
## filter sequences so that they have a good taxonomic description at the species genus and family levels
obigrep -d "${rd_prefix}" --require-rank=species --require-rank=genus --require-rank=family v_"${rd_prefix}".ecopcr > v_"${rd_prefix}"_clean.fasta
## remove redundant sequences
obiuniq -d "${rd_prefix}" v_"${rd_prefix}"_clean.fasta > v_"${rd_prefix}"_clean_uniq.fasta
## ensure that the dereplicated sequences have a taxid at the family level
obigrep -d "${rd_prefix}" --require-rank=family v_"${rd_prefix}"_clean_uniq.fasta > v_"${rd_prefix}"_clean_uniq_clean.fasta
## ensure that sequences each have a unique identification
obiannotate --uniq-id v_"${rd_prefix}"_clean_uniq_clean.fasta > db_"${rd_prefix}".fasta
## your reference database is built !



#obiconvert --skip-on-error --fasta -t ./TAXO --ecopcrdb-output=mitofish/"${rd_prefix}" mitofish/mitogene_12S.fasta
#obiconvert --skip-on-error --embl -t ./TAXO --ecopcrdb-output="${rd_prefix}" EMBL/rel_std_*.dat