#!/bin/sh

# utilisation de script / librairies de FROGS: 
#   solution 1 :utiliser un environnement FROGS (exemple ci-dessous de l'activation de l'environnement sur le cluster Genologin
module load system/Miniconda3-4.7.10
module load bioinfo/FROGS-v3.2.3
conda_lib_dir=`echo $(dirname $(dirname $(which preprocess.py)))/libexec`
export PATH=$conda_lib_dir:$PATH
#   solution 2 utiliser directement les scripts envoyés

# il sera également nécessaire d'avoir RDP Classifier et NCBI BLast+ (éventuellement FROGS pour vérifier que la base fonctionne)

###########
# Téléchargement de la base
# URL : https://data.inrae.fr/dataset.xhtml?persistentId=doi:10.15454/TOMBYZ
# commentaire : certaines cellule du fichier Excel contiennent des retours à la ligne. Ajout au préalable en premier colonne du tag new_line, parser le fichier et réécrire les lignes complètes en supprimant les retours à la ligne


###########
dos2unix 2021-06-25-Diat.barcode-release-version_10.1.txt

### suppression des retours à la ligne
awk '{if($1=="new_line"){if(s!=""){print s} s=$0}else{s=s""$0}}' 2021-06-25-Diat.barcode-release-version_10.1.txt > 2021-06-25-Diat.barcode-release-version_10.1_remove_new_line.txt 
# check nombre de lignes /colonnes
awk -F "\t" '{print NF}' 2021-06-25-Diat.barcode-release-version_10.1_remove_new_line.txt | sort | uniq -c
#    8136 121

###########
# description base
cut -f 5 2021-06-25-Diat.barcode-release-version_10.1_remove_new_line.txt |sort | uniq -c
    #~ 135 18s
   #~ 2426 18S
      #~ 3 18s-V4
     #~ 82 18S-V4
    #~ 315 28S
      #~ 1 Amplified region
    #~ 293 cox1
     #~ 83 ITS2
   #~ 4545 rbcl
    #~ 253 Rbcl-312bp


# Les bases préformatées mothur / DADA2 contiennent 1969 séquences rbCL 263.
####    ==>  EXPLICATION DE FREDERIC RIMET
#   - Sur les nombreuses séquences qui sont présentent dans Diat.barcode, certaines sont identiques (mais provenant de cultures différentes).
#   - D’autre part, les bases Mothur et DADA2 sont trimées sur le barcode court utilisé pour le métabarcoding (pour avoir des temps de calcul plus courts). Donc certaines séquences différentes avec 1500 bp deviennent identiques avec 312 bp. Il faut donc les supprimer.
#   - C’est de cette manière qu’on supprimer pas mal de séquences.

# on part de toutes les séquences RBCL + RBCL 312 qui sont des séquences de 312 bp qui proviennent pas de cultures (comme les 4543 autres séquences) mais de données de métabarcoding, elles sont importantes car elles couvrent un part de la diversité non cultivée en eau douce

# Il y a 2 taxonomies disponibles avec un nombre de genre différent, on conserve la plus ancienne, qui est plus dichotomique et qui ne respecte pas les résultats phylogénétiques des publis récentes, mais qui « fonctionne » mieux. La plus récente est dérivée de l’évolution paraphylétique (en escalier), et du coup beaucoup de rang supérieurs sont très pauvres en taxons, donc assez peu informative.

FASTA=diat_barcode_v10_all_rbcl_frogs.fasta
TAX=diat_barcode_v10_all_rbcl_frogs_tax.txt

awk -F "\t" -v F=$FASTA -v T=$TAX '{
    if ($5=="rbcl" || $5=="Rbcl-312bp"){
        print ">"$2 >> F
        print $4 >> F
        print $2"\t"$41";"$40";"$39";"$38";"$37";"$36";"$35";"$34";"$33 >> T
    }
}' 2021-06-25-Diat.barcode-release-version_10.1_remove_new_line.txt

# Nombre de séquences
grep -c '>' $FASTA
# 4798
wc -l $TAX
# 4798

# check de la présence de doublon
cut -f 1 $TAX | sort | uniq -c | awk '$1>1' | wc -l

# check du même nombre de niveau pour toutes les séquences
cut -f 2 $TAX |awk -F ";" '{print NF}' | sort -u 
# 9

# vérifier que chaque taxon est non vide
sed 's/\t/;/' $TAX | awk -F ";" '{for(i=2;i<=NF;i++){if($i==""){print $1 }}}'
# ==> ok

##################################   PREPARATION DES FICHIERS POUR RDP ##############################
DB_NAME=Diat.barcode_rbcL
DB_VERSION=10.1

mkdir ${DB_NAME}_${DB_VERSION}

# formatage de la base pour RDP
# l'option -r précise le nom des rangs taxonomiques, ici 9 rangs différents

    # si utilisation d'un environnement FROGS
fasta2RDP.py -d $FASTA -t $TAX -r E K SK P C O F G S --rdp-taxonomy ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.tax --rdp-fasta ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.fasta
    # sinon avec les scripts dans le répertoire courant
#~ python3 fasta2RDP.py -d $FASTA -t $TAX -r E K SK P C O F G S --rdp-taxonomy ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.tax --rdp-fasta ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.fasta


##################################   OPTIONNEL VERIFICATION DE L'UNICITE D'UTILISATION DES TAXON ##############################
# recherche de taxons ambigus retrouvés plusieurs fois :
#   - soit à des rangs différents 
#   - soit au même rangs mais avec une parenté différents

cut -f 2 -d '*' ${DB_NAME}_${DB_VERSION}/diat_barcode_v10_all_rbcl_frogs.tax |  cut -f 1 -d ' ' | sort | uniq -c | awk '$1!=1 && $2 !="unclassified"' > problematic_taxon
if [[ -e problematic_taxon_and_rank_number ]]
then
rm problematic_taxon_and_rank_number
fi
for i in `awk '{print $2}' problematic_taxon `
do 
exp=`echo \"\*$i \[id\"`
grep "\*"$i" \[id" ${DB_NAME}_${DB_VERSION}/diat_barcode_v10_all_rbcl_frogs.tax| sed 's/ \[id/*/'| awk -F "*" '{print $2,$5}' >> problematic_taxon_and_rank_number
done

# nombre de taxon présent au même rang mais  avec une parenté différente
cat problematic_taxon_and_rank_number | sort | uniq -c | awk '$1>1' | wc -l
# 1

cat problematic_taxon_and_rank_number | sort | uniq -c | awk '$1>1'
#   2 Madinithidium 8

# nombre de taxon présent à des rang différents
sort -u problematic_taxon_and_rank_number | cut -f 1 -d " "| uniq -c | awk '$1>1' | wc -l
# 0 

##################################  ENTRAINEMENT RDP, INDEX BLAST ET TEST ##############################
### RDP (besoins du chemin du fichier jar de RDP classifier)
sbatch -J ${DB_NAME}_RDP -o %x.out -e %x.err -p workq --mem=80G --wrap="java -Xmx60g -jar /usr/local/bioinfo/src/RDPTools/RDPTools-v2.0.2/classifier.jar train -o ${DB_NAME}_${DB_VERSION} -s ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.fasta -t  ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.tax"
# copie d'un fichier "properties", on peut prendre n'importe lequel et le renommer FASTA.properties
cp /save/frogs/galaxy_databanks/SILVA/16S/silva_128/silva_128_16S.fasta.properties ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.fasta.properties

### BLAST
# binaire Blast disponible dans le PATH (sur genologin /usr/local/bioinfo/src/NCBI_Blast+/ncbi-blast-2.7.1+/bin/)
sbatch -J ${DB_NAME}_Blast -o %x.out -e %x.err --wrap="makeblastdb -in ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.fasta -dbtype nucl"


### TEST D'UTILISATION DE LA BASE en précisant les 9 noms noms dans rang taxonomique
sbatch -J test_${DB_NAME}  -o %x.out -e %x.err --mem=100G --cpus-per-task=30 --wrap="affiliation_OTU.py --rdp --reference ${DB_NAME}_${DB_VERSION}/diat_barcode_v10.1_all_rbcl_frogs.fasta --input-biom swarm.biom --input-fasta swarm.fasta --output-biom affiliation.biom --summary affiliation.summary --nb-cpus 30 --java-mem 15 --taxonomy-ranks Empire Kingdom Subkingdom Phylum Class Order Family Genus Species"

##################################  FINALISATION ##############################

# ajout d'un readme (copier coller du précédent avec éventuellement quelques mises à jour)
touch ${DB_NAME}_${DB_VERSION}/readme.txt

# faire une archive
tar -cvzf ${DB_NAME}_${DB_VERSION}.tar.gz ${DB_NAME}_${DB_VERSION}
