#!/bin/sh
# zagor

# see https://www.cyberciti.biz/howto/question/general/compress-file-unix-linux-cheat-sheet.php in the case you need to uncompress some data


####  ####  ####  ####  ####
####  all steps in comments:
####  ####  ####  ####  ####

# https://github.com/aaschaffer/vecscreen_plus_taxonomy

# find -maxdepth 10 -iname "*vecscreen_output*"
# find -maxdepth 10 -iname "*vecscreen_to_summary.pl"

####  replace appropriately
# export VECPLUSDIR="/home/administrator/Software/vecscreen_plus_taxonomy"
# export PATH="$VECPLUSDIR/scripts:$PATH"
# export PERL5LIB="$VECPLUSDIR"
# export BLASTDB=/DATA/majaz/myDB/

# echo $VECPLUSDIR
# echo $PERL5LIB
# echo $PATH
# echo $BLASTDB

####  source ~/.bashrc

# cd $VECPLUSDIR/scripts
# gunzip vecscreen.gz
# gunzip srcchk.gz
# chmod +x vecscreen
# chmod +x srcchk


# makeblastdb -in UniVec -input_type 'fasta' -dbtype 'nucl'
# wget ftp://ftp.ncbi.nih.gov/pub/UniVec/UniVec_Core
# makeblastdb -in UniVec_Core -input_type fasta -dbtype nucl -parse_seqids -out UniVecDB


# $VECPLUSDIR/scripts/vecscreen -version
## vecscreen: 2.0.0


####  get VECPLUSDIR DB
# cd /DATA/majaz/myDB
# wget ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdmp.zip
# unzip taxdmp.zip
# cut -f1,3,5 nodes.dmp > taxonomy_tree.txt
# $VECPLUSDIR/scripts/assign_levels_to_taxonomy.pl --input_taxa taxonomy_tree.txt --outfile taxonomy_tree_wlevels.txt

mkdir ../intermediate
whereTo=../intermediate
myDB=/DATA/majaz/myDB




$VECPLUSDIR/scripts/from_vecscreen_to_summary.pl \
--input_fasta ../input/Rywal.tr.fasta \
--input_taxa $myDB/taxonomy_tree_wlevels.txt \
--output_root $whereTo/Rywal \
--verbose \
--combine_output


$VECPLUSDIR/scripts/from_vecscreen_to_summary.pl \
--input_fasta ../input/PW363.tr.fasta \
--input_taxa $myDB/taxonomy_tree_wlevels.txt \
--output_root $whereTo/PW363 \
--verbose \
--combine_output


$VECPLUSDIR/scripts/from_vecscreen_to_summary.pl \
--input_fasta ../input/Desiree.tr.fasta \
--input_taxa $myDB/taxonomy_tree_wlevels.txt \
--output_root $whereTo/Desiree \
--verbose \
--combine_output

cut -f1 $whereTo/Rywal.output_combined_wtaxonomy.txt | uniq | sort > $whereTo/Rywal_contaminated_IDs.txt;
cut -f1 $whereTo/PW363.output_combined_wtaxonomy.txt | uniq | sort > $whereTo/PW363_contaminated_IDs.txt;
cut -f1 $whereTo/Desiree.output_combined_wtaxonomy.txt | uniq | sort > $whereTo/Desiree_contaminated_IDs.txt;

xargs faidx -d ' ' ../input/Desiree.tr.fasta < $whereTo/Desiree_contaminated_IDs.txt > $whereTo/Desiree_contaminated_IDs.fasta
xargs faidx -d ' ' ../input/PW363.tr.fasta < $whereTo/PW363_contaminated_IDs.txt > $whereTo/PW363_contaminated_IDs.fasta
xargs faidx -d ' ' ../input/Rywal.tr.fasta < $whereTo/Rywal_contaminated_IDs.txt > $whereTo/Rywal_contaminated_IDs.fasta


export BLASTDB="/NFSi-bioDBs/blast/701-NCBI/" # or wherever it is for you

blastn -query $whereTo/Desiree_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/Desiree_contaminated_nt.out -num_threads 24
./ENCH_sumablastplus.pl -empties $whereTo/Desiree_contaminated_nt.out $whereTo/Desiree_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/Desiree_contaminated_nt.out $whereTo/Desiree_contaminated_nt_top1.out


blastn -query $whereTo/PW363_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/PW363_contaminated_nt.out -num_threads 24
./ENCH_sumablastplus.pl -empties $whereTo/PW363_contaminated_nt.out $whereTo/PW363_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/PW363_contaminated_nt.out $whereTo/PW363_contaminated_nt_top1.out


blastn -query $whereTo/Rywal_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/Rywal_contaminated_nt.out -num_threads 24
./ENCH_sumablastplus.pl -empties $whereTo/Rywal_contaminated_nt.out $whereTo/Rywal_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/Rywal_contaminated_nt.out $whereTo/Rywal_contaminated_nt_top1.out

rm $whereTo/*_contaminated_nt.out



####  calculate seq lengths
seqkit fx2tab --length --name --header-line $whereTo/Desiree_contaminated_IDs.fasta > $whereTo/Desiree_contaminated_seq_len
seqkit fx2tab --length --name --header-line $whereTo/PW363_contaminated_IDs.fasta > $whereTo/PW363_contaminated_seq_len
seqkit fx2tab --length --name --header-line $whereTo/Rywal_contaminated_IDs.fasta > $whereTo/Rywal_contaminated_seq_len


rm $whereTo/*fasta
rm $whereTo/*_IDs.txt
rm $whereTo/*fai


############
## utrorf ##
############

mkdir ../intermediate/utrorf
whereTo=../intermediate/utrorf

# Rywal
grep "utrorf "  ../../_A_01_evigene/intermediate_tr2aacds_Rywal/okayset/*.aa > $whereTo/Rywal.all.aa.utrorf.txt;
awk -F" " '{if ($1) print $1}' $whereTo/Rywal.all.aa.utrorf.txt | cut -d '>' -f2 > $whereTo/retrieveIDs_Rywal.txt;
sed 's/utrorf//' $whereTo/retrieveIDs_Rywal.txt | sponge $whereTo/retrieveIDs_Rywal.txt;
sort $whereTo/retrieveIDs_Rywal.txt | sponge $whereTo/retrieveIDs_Rywal.txt;
# PW363
grep "utrorf "  .../../_A_01_evigene/intermediate_tr2aacds_PW363/okayset/*.aa > $whereTo/PW363.all.aa.utrorf.txt;
awk -F" " '{if ($1) print $1}' $whereTo/PW363.all.aa.utrorf.txt | cut -d '>' -f2 > $whereTo/retrieveIDs_PW363.txt;
sed 's/utrorf//' $whereTo/retrieveIDs_PW363.txt | sponge $whereTo/retrieveIDs_PW363.txt;
sort $whereTo/retrieveIDs_PW363.txt | sponge $whereTo/retrieveIDs_PW363.txt;
# Desiree
grep "utrorf "  ../../_A_01_evigene/intermediate_tr2aacds_Desiree/okayset/*.aa > $whereTo/Desiree.all.aa.utrorf.txt;
awk -F" " '{if ($1) print $1}' $whereTo/Desiree.all.aa.utrorf.txt | cut -d '>' -f2 > $whereTo/retrieveIDs_Desiree.txt;
sed 's/utrorf//' $whereTo/retrieveIDs_Desiree.txt | sponge $whereTo/retrieveIDs_Desiree.txt;
sort $whereTo/retrieveIDs_Desiree.txt | sponge $whereTo/retrieveIDs_Desiree.txt;

xargs faidx -d ' ' ../../_A_01_evigene/output/Rywal.tr.fasta < $whereTo/retrieveIDs_Rywal.txt > $whereTo/ignore.fasta 2>  $whereTo/Rywal_retrieve.tr.err.txt;
xargs faidx -d ' ' ../../_A_01_evigene/output/PW363.tr.fasta < $whereTo/retrieveIDs_PW363.txt > $whereTo/ignore.fasta 2>  $whereTo/PW363_retrieve.tr.err.txt;
xargs faidx -d ' ' ../../_A_01_evigene/output/Desiree.tr.fasta < $whereTo/retrieveIDs_Desiree.txt > $whereTo/ignore.fasta 2>  $whereTo/Desiree_retrieve.tr.err.txt;

rm $whereTo/ignore.fasta
rm $whereTo/*all.aa.utrorf.txt
rm $whereTo/retrieveIDs_*.txt

cut -f2 -d " " $whereTo/Rywal_retrieve.tr.err.txt > $whereTo/retrieveIDs_Rywal.txt
cut -f2 -d " " $whereTo/PW363_retrieve.tr.err.txt > $whereTo/retrieveIDs_PW363.txt
cut -f2 -d " " $whereTo/Desiree_retrieve.tr.err.txt > $whereTo/retrieveIDs_Desiree.txt

rm $whereTo/*tr.err.txt

# Rywal
xargs faidx -d ' ' ../../_A_01_evigene/intermediate_tr2aacds_Rywal/all.tr < $whereTo/retrieveIDs_Rywal.txt > $whereTo/Rywal.all.aa.utrorf_rawIDs_tr.fasta 2>  $whereTo/Rywal.all.aa.utrorf_rawIDs.fasta.tr.err.txt;
fasta_formatter -i  $whereTo/Rywal.all.aa.utrorf_rawIDs_tr.fasta -o $whereTo/Rywal.linear.utrorf_tr -w 0;
rm $whereTo/Rywal.all.aa.utrorf_rawIDs_tr.fasta;
# PW363
xargs faidx -d ' ' .../../_A_01_evigene/intermediate_tr2aacds_PW363/all.tr < $whereTo/retrieveIDs_PW363.txt > $whereTo/PW363.all.aa.utrorf_rawIDs_tr.fasta 2>  $whereTo/PW363.all.aa.utrorf_rawIDs.fasta.tr.err.txt;
fasta_formatter -i  $whereTo/PW363.all.aa.utrorf_rawIDs_tr.fasta -o $whereTo/PW363.linear.utrorf_tr -w 0;
rm $whereTo/PW363.all.aa.utrorf_rawIDs_tr.fasta;
# Desiree
xargs faidx -d ' ' ../../_A_01_evigene/intermediate_tr2aacds_Desiree/all.tr < $whereTo/retrieveIDs_Desiree.txt > $whereTo/Desiree.all.aa.utrorf_rawIDs_tr.fasta 2>  $whereTo/Desiree.all.aa.utrorf_rawIDs.fasta.tr.err.txt;
fasta_formatter -i  $whereTo/Desiree.all.aa.utrorf_rawIDs_tr.fasta -o $whereTo/Desiree.linear.utrorf_tr -w 0;
rm $whereTo/Desiree.all.aa.utrorf_rawIDs_tr.fasta;

rm $whereTo/retrieveIDs_*.txt



$VECPLUSDIR/scripts/from_vecscreen_to_summary.pl \
--input_fasta $whereTo/Rywal.linear.utrorf_tr \
--input_taxa $myDB/taxonomy_tree_wlevels.txt \
--output_root $whereTo/Rywal.utrorf \
--verbose \
--combine_output


$VECPLUSDIR/scripts/from_vecscreen_to_summary.pl \
--input_fasta $whereTo/PW363.linear.utrorf_tr \
--input_taxa $myDB/taxonomy_tree_wlevels.txt \
--output_root $whereTo/PW363.utrorf \
--verbose \
--combine_output


$VECPLUSDIR/scripts/from_vecscreen_to_summary.pl \
--input_fasta $whereTo/Desiree.linear.utrorf_tr \
--input_taxa $myDB/taxonomy_tree_wlevels.txt \
--output_root $whereTo/Desiree.utrorf \
--verbose \
--combine_output

cut -f1 $whereTo/Rywal.utrorf.output_combined_wtaxonomy.txt | uniq | sort > $whereTo/Rywal.utrorf_contaminated_IDs.txt;
cut -f1 $whereTo/PW363.utrorf.output_combined_wtaxonomy.txt | uniq | sort > $whereTo/PW363.utrorf_contaminated_IDs.txt;
cut -f1 $whereTo/Desiree.utrorf.output_combined_wtaxonomy.txt | uniq | sort > $whereTo/Desiree.utrorf_contaminated_IDs.txt;

xargs faidx -d ' ' $whereTo/Desiree.linear.utrorf_tr < $whereTo/Desiree.utrorf_contaminated_IDs.txt > $whereTo/Desiree.utrorf_contaminated_IDs.fasta
xargs faidx -d ' ' $whereTo/PW363.linear.utrorf_tr < $whereTo/PW363.utrorf_contaminated_IDs.txt > $whereTo/PW363.utrorf_contaminated_IDs.fasta
xargs faidx -d ' ' $whereTo/Rywal.linear.utrorf_tr < $whereTo/Rywal.utrorf_contaminated_IDs.txt > $whereTo/Rywal.utrorf_contaminated_IDs.fasta

rm $whereTo/*_IDs.txt
rm $whereTo/*fai
cp $whereTo/Rywal.linear.utrorf_tr  ../../_A_02.8_DIAMOND/input/Rywal.linear.utrorf_tr
cp $whereTo/PW363.linear.utrorf_tr ../../_A_02.8_DIAMOND/input/PW363.linear.utrorf_tr
cp $whereTo/Desiree.linear.utrorf_tr ../../_A_02.8_DIAMOND/input/Desiree.linear.utrorf_tr
rm $whereTo/*linear.utrorf_tr

####  calculate seq lengths
seqkit fx2tab --length --name --header-line $whereTo/Desiree.utrorf_contaminated_IDs.fasta > $whereTo/Desiree.utrorf_contaminated_seq_len
seqkit fx2tab --length --name --header-line $whereTo/PW363.utrorf_contaminated_IDs.fasta > $whereTo/PW363.utrorf_contaminated_seq_len
seqkit fx2tab --length --name --header-line $whereTo/Rywal.utrorf_contaminated_IDs.fasta > $whereTo/Rywal.utrorf_contaminated_seq_len


blastn -query $whereTo/Rywal.utrorf_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/Rywal.utrorf_contaminated_nt.out -num_threads 24
./ENCH_sumablastplus.pl -empties $whereTo/Rywal.utrorf_contaminated_nt.out $whereTo/Rywal.utrorf_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/Rywal.utrorf_contaminated_nt.out $whereTo/Rywal.utrorf_contaminated_nt_top1.out

blastn -query $whereTo/PW363.utrorf_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/PW363.utrorf_contaminated_nt.out -num_threads 24
# ./ENCH_sumablastplus.pl -empties $whereTo/PW363.utrorf_contaminated_nt.out $whereTo/PW363.utrorf_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/PW363.utrorf_contaminated_nt.out $whereTo/PW363.utrorf_contaminated_nt_top1.out

blastn -query $whereTo/Rywal.utrorf_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/Rywal.utrorf_contaminated_nt.out -num_threads 24
# ./ENCH_sumablastplus.pl -empties $whereTo/Rywal.utrorf_contaminated_nt.out $whereTo/Rywal.utrorf_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/Rywal.utrorf_contaminated_nt.out $whereTo/Rywal.utrorf_contaminated_nt_top1.out

blastn -query $whereTo/Desiree.utrorf_contaminated_IDs.fasta -db $BLASTDB/nt \
-out $whereTo/Desiree.utrorf_contaminated_nt.out -num_threads 24
# ./ENCH_sumablastplus.pl -empties $whereTo/Desiree.utrorf_contaminated_nt.out $whereTo/Desiree.utrorf_contaminated_nt_tsv.out
./ENCH_sumablastplus.pl -maxhits 1 -empties $whereTo/Desiree.utrorf_contaminated_nt.out $whereTo/Desiree.utrorf_contaminated_nt_top1.out

rm $whereTo/*fasta
rm $whereTo/*_nt.out


rm ../input/*fasta
rm ../input/*fai



####################################
## quick and dirty - manual labor ##
# merge files to output
# add colnames
# VecScreen column names:
# Descriptions of the 11 columns:
# Column  1: Accession of query
# Column  2: Genus of query if known, or 1 otherwise
# Column  3: Species of query if known, or 1 otherwise
# Column  4: Lower end of the alignment in the vector
# Column  5: Upper end of the alignment in the vector
# Column  6: Matching vector, starting with uv|
# Column  7: One end of the alignment in the vector
# Column  8: The other end of the alignment in the vector
# Column  9: The strength of this vecscreen match
# Column 10: The strength of the strongest vecscreen match for this query
# Column 11: Whether there is any dangling part (called "Suspect" by vecscreen) at either end of the query
# ENCH_sumablastplus.pl column names
# Query_ID
# Query_length
# Target_ID
# Target_length
# Target_description
# E-value
# Identity%
# Aligned_seq_length
# Score
# Positives%
# Gaps%


tar -zcvf ../intermediate.tar.gz ../intermediate/
rm -rf ../intermediate/

