#!/bin/bash -e

#SBATCH --job-name       MR                  # The job name
#SBATCH --account        ####                # The account code
#SBATCH --time           24:00:00            # The walltime
#SBATCH --mem            16G                 #
#SBATCH -D ./                                # The initial directory
#SBATCH -o %x_out.txt                        # The output file
#SBATCH -e %x_err.txt                        # The error file

## Load required programmes/modules
module load cutadapt/3.5-gimkl-2020a-Python-3.8.2
module load VSEARCH/2.21.1-GCC-11.3.0
module load BLAST/2.13.0-GCC-11.3.0
module load BLASTDB/2023-07

## Forward and (reverse-complemented) reverse primer sequences for cutadapt:
## 16S F
F_16S=CCTACGGGNGGCWGCAG # 341F
## 16S R reverse complement
RC_16S=ATTAGAWACCCBNGTAGTCC # 806R
## ITS F
F_ITS=GTGARTCATCGAATCTTTG # fITS7
## ITS R reverse complement
RC_ITS=GCATATCAATAAGCGGAGGA # ITS4
## AMF F
F_AMF=TTGGAGGGCAAGTCTGGTGCC # NS31
## AMF R reverse complement 
RC_AMF=GGAAACCAAAGTGTTTGGGTTC # AML2
## COI F
F_COI=GGWACWGGWTGAACWGTWTAYCCYCC # mlCOIintF
## COI R reverse complement
RC_COI=TGATTTTTTGGTCACCCTGAAGTTTA # HCO2198
## 18S
F_18S=GYGGTGCATGGCCGTTSKTRGTT # #3F
## 18S R reverse complement
RC_18S=GTCCCTGVCCTTTGTRCACAC # #5RC

###############################################################################
## Test effect of truncating reads at different quality thresholds for several different samples
## For datasets with poor merging success rates - 16S, ITS, AMF
## Q7 best for 16S, Q8 for ITS; but no benefit for AMF  
# tqs=(5 6 7 8 9 10)
# g="COI"
# mkdir -p processing_test 
# for i in ${!tqs[@]}; do
  # tq=${tqs[$i]}
  # for r1 in $( ls MR_${g}/MRustS-${g}__25_1.fastq.gz ); do 
    # echo $r1
    # r2=$(sed -e "s/_1.fastq/_2.fastq/" <<< "$r1") # Get reverse filename
    # echo $r2
    # s=$( basename $r1 | cut -d_ -f3 ) # Get sample ID
    # echo Merging $g $s
    # echo truncating at Q $tq
    # # Trim low quality regions before merging
    # vsearch -fastx_filter $r1 -reverse $r2 -fastq_truncqual ${tq} -fastqout processing_test/${s}_r1_trunc.fastq -fastqout_rev processing_test/${s}_r2_trunc.fastq 
    # vsearch -fastq_mergepairs processing_test/${s}_r1_trunc.fastq -reverse processing_test/${s}_r2_trunc.fastq -fastqout processing_test/merged_temp_${s}.fastq -fastq_allowmergestagger -relabel "S$s|"
  # done
# done

###############################################################################
# Process full dats with optimal quality truncation
# groups=("18S" "COI")
# tqs=(8 9) 
groups=("16S" "18S" "ITS" "COI")
tqs=(7 8 8 9) 
for i in ${!groups[@]}; do
  g=${groups[$i]}
  mkdir -p processed_${g}
  tq=${tqs[$i]}
  for r1 in $( ls MR_${g}/*_1.fastq.gz ); do
    r2=$(sed -e "s/_1.fastq/_2.fastq/" <<< "$r1") # Get reverse filename
    echo $r1
    echo $r2
    s=$( basename $r1 | cut -d_ -f3 ) # Get sample ID
    echo Merging $g $s after truncating at Q $tq
    # Trim low quality regions before merging
    vsearch -fastx_filter $r1 -reverse $r2 -fastq_truncqual ${tq} -fastqout processed_${g}/${s}_r1_trunc.fastq -fastqout_rev processed_${g}/${s}_r2_trunc.fastq 
    vsearch -fastq_mergepairs processed_${g}/${s}_r1_trunc.fastq -reverse processed_${g}/${s}_r2_trunc.fastq -fastqout processed_${g}/merged_temp_${s}.fastq -fastq_allowmergestagger -relabel "S$s|"
  done
  cat processed_${g}/merged_temp*.fastq > processed_${g}/all_merged_${g}.fastq
#for i in ${!groups[@]}; do
#  g=${groups[$i]}
  primer_F=F_${g}
  primer_R=RC_${g}
  cutadapt -a ^${!primer_F}...${!primer_R}$ processed_${g}/all_merged_${g}.fastq -o processed_${g}/trimmed_${g}.fastq --discard-untrimmed --cores 4
  vsearch -fastq_stats processed_${g}/trimmed_${g}.fastq -log processed_${g}/trimmed_fq_stats_${g}.txt 
done

######################################################################################
# # AMF doesn't overlap by much and has very low merging rates. So, truncate low-quality ends and join reads instead?
# Need to retain all sequences during truncation to get consistent numbers for fastq_join.
groups2=("AMF")
for i in ${!groups2[@]}; do
  g=${groups2[$i]}
  mkdir processed_joined_"${g}"
  for r1 in $( ls MR_${g}/*_1.fastq.gz ); do
    echo $r1
    r2=$(sed -e "s/_1.fastq/_2.fastq/" <<< "$r1") # Get reverse filename
    echo $r2
    s=$( basename $r1 | cut -d_ -f3 ) # Get sample ID
    echo $s
    echo Truncating and joining ${g} $s
    vsearch -fastx_filter $r1 -fastq_trunclen_keep 260 -fastqout processed_joined_"${g}"/R1_trunc_"$s".fastq
    vsearch -fastx_filter $r2 -fastq_trunclen_keep 160 -fastqout processed_joined_"${g}"/R2_trunc_"$s".fastq
    vsearch -fastq_join processed_joined_"${g}"/R1_trunc_"$s".fastq -reverse processed_joined_"${g}"/R2_trunc_"$s".fastq -fastqout processed_joined_"${g}"/joined_temp_"$s".fastq -relabel "$s|"
  done
  cat processed_joined_"${g}"/joined_temp_*.fastq > processed_joined_"${g}"/all_joined_"${g}".fastq
  primer_F=F_"${g}"
  primer_R=RC_"${g}"
  cutadapt -a ^${!primer_F}...${!primer_R}$ processed_joined_"${g}"/all_joined_"${g}".fastq -o processed_joined_"${g}"/joined_trimmed_"${g}".fastq --discard-untrimmed --cores 4
  vsearch -fastq_stats processed_joined_"${g}"/joined_trimmed_"${g}".fastq -log processed_joined_"${g}"/joined_trimmed_fq_stats_"${g}".txt 
done

###############################################################################
## Process merged sequences into ASVs and OTUs
## (run as separate jobs per amplicon)

for i in ${!groups[@]}; do
  g=${groups[$i]}
  echo processing $g merged
  ## Set min and max lengths for filtering step
  if [ $g == "16S" ]; then 
    #minlen=260
    minlen=380
    maxlen=450
  elif [ $g == "18S" ]; then # Strangely variable lengths?
    minlen=200
    maxlen=450
  elif [ $g == "ITS" ]; then # Variable but expected
    minlen=180
    maxlen=470
  elif [ $g == "AMF" ]; then 
    minlen=390
    maxlen=550
  elif [ $g == "COI" ]; then 
    minlen=300
    maxlen=330
  fi
  ## Filter for errors and length, and dereplicate
  vsearch -fastx_filter processed_${g}/trimmed_${g}.fastq -fastq_maxee 1.0 -fastq_maxns 0 -fastq_minlen $minlen -fastq_maxlen $maxlen -fastaout processed_${g}/filtered_${g}.fasta
  vsearch -derep_fulllength processed_${g}/filtered_${g}.fasta -relabel Uniq -sizeout -output processed_${g}/uniques_${g}.fasta
  vsearch -sortbysize processed_${g}/uniques_${g}.fasta -minsize 2 -sizeout -output processed_${g}/uniques_min2_${g}.fasta
  ## Convert trimmed fastq to fasta for read mapping?
  vsearch -fastx_filter processed_${g}/trimmed_${g}.fastq -fastaout processed_${g}/trimmed_${g}.fasta
  ## Denoising
  vsearch -cluster_unoise processed_${g}/uniques_${g}.fasta -centroids processed_${g}/ASVs_temp_${g}.fasta -relabel ASV
  vsearch -uchime3_denovo processed_${g}/ASVs_temp_${g}.fasta -nonchimeras processed_${g}/ASVs_${g}.fasta
  vsearch -usearch_global processed_${g}/trimmed_${g}.fasta -db processed_${g}/ASVs_${g}.fasta -id 0.97 -otutabout processed_${g}/ASVtab_${g}.txt
  ## OTU clustering
  #vsearch -cluster_fast processed_${g}/uniques_min2_${g}.fasta -id 0.97 -centroids processed_${g}/OTUs_temp_${g}.fasta -relabel otu
  #vsearch -uchime_denovo processed_${g}/OTUs_temp_${g}.fasta -nonchimeras processed_${g}/OTUs_${g}.fasta
  #vsearch -usearch_global processed_${g}/trimmed_${g}.fasta -db processed_${g}/OTUs_${g}.fasta -id 0.97 -otutabout processed_${g}/OTUtab_${g}.txt
done

# Identify the ASVs and OTUs
for i in ${!groups[@]}; do
  g=${groups[$i]}
  echo identifying ${g}
  if [ $g == "16S" ]; then 
    java -Xmx8g -jar bioinformatics/rdp_classifier_2.13/dist/classifier.jar -o processed_16S/ASVs_16S_rdp.txt processed_16S/ASVs_16S.fasta
    java -Xmx8g -jar bioinformatics/rdp_classifier_2.13/dist/classifier.jar -o processed_16S/OTUs_16S_rdp.txt processed_16S/OTUs_16S.fasta
  elif [ $g == "ITS" ]; then
# Identification of ITS via RDP and UNITE classifier v2.0 (based on UNITE + INSD v8.3; https://github.com/terrimporter/UNITE_ITSClassifier)
    java -Xmx16g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/UNITE_v2_RDP/mydata_trained/rRNAClassifier.properties -o processed_ITS/ASVs_ITS_rdp_UNITEv2.txt processed_ITS/ASVs_ITS.fasta
    #java -Xmx16g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/UNITE_v2_RDP/mydata_trained/rRNAClassifier.properties -o processed_ITS/OTUs_ITS_rdp_UNITEv2.txt processed_ITS/OTUs_ITS.fasta
  elif [ $g == "18S" ]; then 
# Identification of 18S via RDP and 18S classifier v4.1 (based on SILVA 138 SSU Ref Nr99; https://github.com/terrimporter/18SClassifier)
    java -Xmx8g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/18Sv4.1_mydata_trained/rRNAClassifier.properties -o processed_18S/ASVs_18S_rdp_SILVA138.txt processed_18S/ASVs_18S.fasta
    #java -Xmx8g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/18Sv4.1_mydata_trained/rRNAClassifier.properties -o processed_18S/OTUs_18S_rdp_SILVA138.txt processed_18S/OTUs_18S.fasta
  elif [ $g == "COI" ]; then 
# Identification of COI via RDP and COI classifier v4 (based on BOLD and GenBank > 500 bp with binomial species names - thus excluding most LBI barcodes)
    java -Xmx8g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/RDP_COIv5.1.0/rRNAClassifier.properties -o processed_COI/ASVs_COI_rdp_BOLDGBv5.txt processed_COI/ASVs_COI.fasta
    #java -Xmx8g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/RDP_COIv5.1.0/rRNAClassifier.properties -o processed_COI/OTUs_COI_rdp_BOLDGBv5.txt processed_COI/OTUs_COI.fasta
# Also try MIDORI COI database server?
  elif [ $g == "AMF" ]; then 
# RDP identification of AMF via RDP and SILVA database? (or Marjaam?)
  java -Xmx8g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/18Sv4.1_mydata_trained/rRNAClassifier.properties -o processed_joined_AMF/ASVs_AMF_rdp_SILVA138.txt processed_joined_AMF/ASVs_AMF.fasta
  #java -Xmx8g -jar rdp_classifier_2.13/dist/classifier.jar -t databases/18Sv4.1_mydata_trained/rRNAClassifier.properties -o processed_joined_AMF/OTUs_AMF_rdp_SILVA138.txt processed_joined_AMF/OTUs_AMF.fasta
  fi
done

# BLAST ID for COI
BLASTOPTS="-evalue 0.001 -max_target_seqs 1 -perc_identity 70"
blastn $BLASTOPTS -db nt -query processed_COI/ASVs_COI.fasta -outfmt '6 qseqid qstart qend sgi sacc sstart send staxids sscinames stitle length pident evalue bitscore' -out processed_COI/ASVs_COI_nt_blastnx1_70.txt -num_threads 4
blastn $BLASTOPTS -db nt -query processed_COI/OTUs_COI.fasta -outfmt '6 qseqid qstart qend sgi sacc sstart send staxids sscinames stitle length pident evalue bitscore' -out processed_COI/OTUs_COI_nt_blastnx1_70.txt -num_threads 4

# Generate match list for lulu curation
for i in ${!groups[@]}; do
  g=${groups[$i]}
  vsearch --usearch_global processed_${g}/ASVs_${g}.fasta --db processed_${g}/ASVs_${g}.fasta --self --id .84 --iddef 1 --userout processed_${g}/match_list.txt -userfields query+target+id --maxaccepts 0 --query_cov .9 --maxhits 10
done
vsearch --usearch_global processed_joined_AMF/ASVs_AMF.fasta --db processed_joined_AMF/ASVs_AMF.fasta --self --id .84 --iddef 1 --userout processed_joined_AMF/match_list_AMF.txt -userfields query+target+id --maxaccepts 0 --query_cov .9 --maxhits 10