#!/usr/bin/env sh

help="[Usage] ./prepare_data.sh infasta ingtf species outfileprefix  species: hs38 mm10"
infa=$1
ingtf=$2
species=$3
outfile=$4

if [[ $# -eq 0 ]];then
    echo ${help}
    exit 2
fi

## should judge where species is legal

ingtf2="ref_Hs.gtf"
infa2="ref_filled.fa"
# gff2gtf
# Convert NC... to chrxxx; remove non-chr records; remove records located in more than one chromosome; get short transcript name.
# Keep transcript both start with NM_ and BestRefSeq level.
# Random(with out replace) pick 5000 transcript and get it fasta and gtf file. 
gff2gtf ${ingtf} -o ${ingtf2}
cat ${ingtf2} | grep "transcript_id" > ${ingtf2}.tmp
../7.Misc/ncbi_chr.py -s ${species} --strict ${ingtf2}.tmp > ${ingtf2}
cat ${ingtf2} | grep "BestRefSeq" | cut -d ";" -f 2| cut -d "\"" -f 2| uniq|sort|uniq > ref_gff.txt
cat ${infa} | grep ">" | cut -d "|" -f4 | sort | uniq >ref_NM.txt
../../base/filter.py ref_NM.txt ref_gff.txt -n1 -n2 > common_NM.txt

../7.Misc/illegal_filter.py ${infa} > ${infa2}
../7.Misc/subset_fasta.py -k common_NM.txt ${infa2} -l 200 -n 5000 --row_len 0 -o ${outfile}.fa
cat ${outfile}.fa | awk -F "|" '{if($0~/>/)print(">"$4);else print$0}'> tmp
mv tmp ${outfile}.fa
transeq -frame F ${outfile}.fa ${outfile}.p.fa
cat ${outfile}.fa | grep ">" | cut -d '>' -f2 | cut -d '|' -f1 >${outfile}.txt
cat ${ingtf2} |cgat gtf2gtf -m filter --filter-method=transcript --map-tsv-file=${outfile}.txt |grep -v "#" > ${outfile}.gtf
../7.Misc/del_multi_chromosome_transcripts.py ${outfile}.gtf > tmp
mv tmp ${outfile}.gtf
if [[ ${species} == "human" ]];then
    ../7.Misc/gtf38to19.sh ${outfile}.gtf ${outfile}_hg19.gtf
fi

rm ref_NM.txt
rm common_NM.txt
rm ref_gff.txt
rm ${ingtf2}
rm ${ingtf2}.tmp
rm ${infa2}