#!/usr/bin/env bash

# This is the adapted version of panEDTA
#       Each genome will be reannotated by the panEDTA library to get the final TE annotation.
#       The panEDTA library is generated by the panEDTA.pl script in the previous step.
#       The panEDTA library is a non-redundant library of TEs from all input genomes.

## Test files explained
# Col.test.fa   Genome 1
# Ler.test.fa   Genome 2
# genome.cds.fa CDS sequence. Could be absent, or from all input genomes, or from one of the genomes, or from a closely related species.
#               If users provided only one CDS file, then use it for all genomes. Each genome could use a different CDS.
#               This file is optional, but having it will improve the annotation by purging gene sequences from TEs.
# athrep.ref.formatted  Existing TE library from this species. Optional. Providing a curated library will significantly improve
#               the annotation and also harvest existing knowledge.


# Help info
helpFunction()
{
   echo ""
   echo "Usage: $0 -genomes genome_list.txt -cds cds.fasta -threads 10"
   echo -e "\t-g        A list of genome files with paths accessible from the working directory.
                        Required: You can provide only a list of genomes in this file (one column, one genome each row).
                        Optional: You can also provide both genomes and CDS files in this file (two columns, one genome and one CDS each row).
                                  Missing of CDS files (eg, for some or all genomes) is allowed."
   echo -e "\t-c                Optional. Coding sequence files in fasta format.
                        The CDS file provided via this parameter will fill in the missing CDS files in the genome list.
                        If no CDS files are provided in the genome list, then this CDS file will be used on all genomes."
   echo -e "\t-l        Optional. A manually curated, non-redundant library following the RepeatMasker naming format."
   echo -e "\t-f        Minimum number of full-length TE copies in individual genomes to be kept as candidate TEs for the pangenome.
                        Lower is more inclusive, and will ↑ library size, ↑ sensitivity, and ↑ inconsistency.
                        Higher is more stringent, and will ↓ library size, ↓ sensitivity, and ↓ inconsistency.
                        Default: 3."
   echo -e "\t-t        Number of CPUs to run panEDTA. Default: 10."
   echo -e "\t-s        Start line of files"
   echo -e "\t-e        End line of files"
   echo ""
   exit 1 # Exit script after printing help
}

# Preset parameters
genome_list=''
cds=''
curatedlib=''
fl_copy=3
threads=10
sta=1
end=1

# Read user inputs
while getopts "g:c::l::f::t::s::e::" opt
do
   case "$opt" in
      g ) genome_list="$OPTARG" ;;
      c ) cds="$OPTARG" ;;
      l ) curatedlib="$OPTARG" ;;
      f ) fl_copy="$OPTARG" ;;
      t ) threads="$OPTARG" ;;
      ? ) helpFunction ;; # Print helpFunction in case parameter is non-existent
      s ) str="$OPTARG" ;;
      e ) end="$OPTARG" ;; 
  esac
done

# Check parameters
if [ ! -s "$genome_list" ]; then
   echo "ERROR: The genomes $genomes_list file is not found or is empty";
   helpFunction
fi

if [ "$cds" != '' ] & [ ! -s "$cds" ]; then
   echo "ERROR: The cds $cds file is not found or is empty"
   helpFunction
fi

if [ "$curatedlib" != '' ] & [ ! -s "$curatedlib" ]; then
   echo "ERROR: The curated library $curatedlib file is not found or is empty"
   helpFunction
fi

# Set paths
path=$(dirname "$0") #program path
dir=$(pwd) #current work directory


### Begin panEDTA and print all parameters
echo `date`
echo "Pan-genome Extensive de-novo TE Annotator (panEDTA)"
echo -e "\tGenome files: $genome_list"
echo -e "\tCoding sequences: $cds"
echo -e "\tCurated library: $curatedlib"
echo -e "\tCopy cutoff: $fl_copy"
echo -e "\tCPUs: $threads"

## Step 3, re-annotate all genomes with the panEDTA library, consider to submit each RepeatMasker and EDTA job to different nodes.
for i in `sed -n "$sta,$end p" $genome_list`; do
        genome=`echo $i|awk '{print $1}'`
        accession=`echo $i|awk '{print $3}'`
        cds_ind=`echo $i|awk '{print $2}'`

        mkdir ~/pan_genome_of_rice/assembly/EDTA/$accession
        cd ~/pan_genome_of_rice/assembly/EDTA/$accession
        echo "Reannotate genome $genome with the panEDTA library - homology"
        RepeatMasker -pa $threads -q -div 40 -lib ~/pan_genome_of_rice/assembly/EDTA/panEDTA.TElib.fa -cutoff 225 -gff $genome.mod >/dev/null
        perl -i -nle 's/\s+DNA\s+/\tDNA\/unknown\t/; print $_'  $genome.mod.out 
        echo "Reannotate genome $genome with the panEDTA library - structural"
        perl $path/EDTA.pl --genome $genome -t $threads --step final --anno 1 --curatedlib ~/pan_genome_of_rice/assembly/EDTA/panEDTA.TElib.fa --cds $cds_ind --rmout $genome.mod.out

done