// cjt
// v3.1.5
// 10-12-16


import java.util.*;
import java.io.*;
import java.util.regex.*;
import java.security.*;
import java.math.*;

import jargs.CmdLineParser;

public class generateDatabase
{
    private static CmdLineParser parser;

    //-mandatory
    private static String opt_outpath;
    private static String opt_dbname;
    private static String opt_inputsequence;

    //-optional
    private static double opt_evalue;
    private static int opt_wvalue;
    private static int opt_fragmentlength;
    private static int opt_fragmentoverlap;
    private static int opt_minHitLength;
    private static int opt_maxHitLength;




    private static void printInfo()
    {
        System.err.println();
        System.out.println("SplamiR - Prediction of spliced miRNAs in plants.");
        System.out.println("generateDatabase, v3.1.5");
        System.out.println("-------------------------------------------------");
        System.out.println("");
    }

    private static void printUsage()
    {
        System.err.println();
        System.err.println(
            "Usage: java -Xms512M -Xmx1024M generateDatabase [-d,--output-dir] [-i,--input-sequence] [{-e,--e-value}] [{-w,--w-value}] [{-l,--fragment-length}] [{-o,--fragment-overlap}] [{-n,--minHitLength}] [{-m,--maxHitLength}]");
        System.err.println();
        System.err.println("Mandatory parameters:");
        System.err.println("  output-dir: directory to store the database");
        System.err.println("  input-sequence: input FASTA file");
        System.err.println();
        System.err.println("Optional parameters:");
        System.err.println("  e-value: BLAST threshold (default: 10.0)");
        System.err.println("  w-value: BLAST initial word length (default: 9)");
        System.err.println("  fragment-length: size of the fragments, limits pri-u-miRNA length (default: 30,000)");
        System.err.println("  fragment-overlap: share of successive fragments (default: 10,000)");
        System.err.println("  minHitLength: min length of fragments to be stored in database (default: 16)");
        System.err.println("  maxHitLength: max length of fragments to be stored in database (default: 500)");
    }





    public static void main (String[] args)
    {
        printInfo();

        parser = new CmdLineParser();

        integrityCheck();
        readParams(args);

        System.out.print("Preprocess FASTA input file");
        fastaInputFileToGenomicSequence();
        System.out.println("...done.");

        System.out.print("Generate fragments");
        genomicSequenceFileToFragments();
        System.out.println("...done.");

        System.out.print("BLAST fragments");
        blastFragments();
        System.out.println("...done.");

        System.out.print("Processing results");
        processBlastOutput();
        System.out.println("...done.");

        System.out.print("Formatting database");
        generateDatabase();
        System.out.println("...done.");

        System.out.println("\nFinished.");
    }





    /*
     *  check for blast, etc.
     */
    private static void integrityCheck()
    {
        //check if all tools are present
        File checkFile=new File("RNAblast/formatdbrna");
        if (!checkFile.isFile())
        {
            System.out.println("ERROR: File RNAblast/formatdbrna is missing. Wrong directory?");
            System.exit(3);
        }
        checkFile=new File("RNAblast/bl2seqrna");
        if (!checkFile.isFile())
        {
            System.out.println("ERROR: File RNAblast/bl2seqrna is missing. Wrong directory?");
            System.exit(3);
        }
    }




    /*
     *  Helper procedure to load instructions from cmd line
     */
    private static void readParams(String[] args)
    {

        //define cmd line options
        //-mandatory
        CmdLineParser.Option opt_outpath1 = parser.addStringOption('d', "output-dir");
        CmdLineParser.Option opt_inputsequence1 = parser.addStringOption('i', "input-sequence");

        //-optional
        CmdLineParser.Option opt_evalue1 = parser.addDoubleOption('e', "e-value");
        CmdLineParser.Option opt_wvalue1 = parser.addIntegerOption('w', "w-value");
        CmdLineParser.Option opt_fragmentlength1 = parser.addIntegerOption('l', "fragment-length");
        CmdLineParser.Option opt_fragmentoverlap1 = parser.addIntegerOption('o', "fragment-overlap");
        CmdLineParser.Option opt_minHitLength1 = parser.addIntegerOption('n', "minHitLength");
        CmdLineParser.Option opt_maxHitLength1 = parser.addIntegerOption('m', "maxHitLength");

        //read params
        try
        {
            parser.parse(args);
        }
        catch ( CmdLineParser.OptionException e )
        {
            System.err.println("ERROR parsing parameters:");
            System.err.println(e.getMessage());
            printUsage();
            System.exit(2);
        }


        //-mandatory
        opt_outpath=(String)parser.getOptionValue(opt_outpath1);
        opt_inputsequence=(String)parser.getOptionValue(opt_inputsequence1);

        //-optional
        opt_evalue=(Double)parser.getOptionValue(opt_evalue1, new Double(10.0));
        opt_wvalue=(Integer)parser.getOptionValue(opt_wvalue1, new Integer(9));
        opt_fragmentlength=(Integer)parser.getOptionValue(opt_fragmentlength1, new Integer(30000));
        opt_fragmentoverlap=(Integer)parser.getOptionValue(opt_fragmentoverlap1, new Integer(10000));
        opt_minHitLength=(Integer)parser.getOptionValue(opt_minHitLength1, new Integer(16));
        opt_maxHitLength=(Integer)parser.getOptionValue(opt_maxHitLength1, new Integer(500));




        //something missing?
        if(opt_outpath==null)
        {
            System.err.println("ERROR: No output directory specified.");
            printUsage();
            System.exit(2);
        }
        if(opt_inputsequence==null)
        {
            System.err.println("ERROR: No input sequence specified.");
            printUsage();
            System.exit(2);
        }


        if (opt_outpath.substring(opt_outpath.length()-1)!="/")
        {
            opt_outpath+="/";
        }

        File f = new File(opt_inputsequence);
        opt_dbname = f.getName().substring(Math.max(f.getName().lastIndexOf("/"),0), f.getName().lastIndexOf("."));


        //avoid overwriting files
        File path=new File(opt_outpath);
        if (path.isDirectory())
        {
            File[] listOfFiles = path.listFiles();
            if(listOfFiles.length>=1)
            {
                System.err.println("ERROR: Output path ('"+opt_outpath+"') is not empty.");
                printUsage();
                System.exit(5);
            }
        }
        else
        {
            //create folder
            path.mkdirs();
        }

        File fastafile=new File(opt_inputsequence);
        if (!fastafile.isFile())
        {
            System.err.println("ERROR: Input sequence ('"+fastafile+"') is missing.");
            printUsage();
            System.exit(5);
        }


        //Output the list to the console
        System.out.println("Configuration:");
        System.out.println("\tdb-name= "+opt_dbname);
        System.out.println("\toutput-dir= "+opt_outpath);
        System.out.println("\tinput-sequence= "+opt_inputsequence);

        System.out.println("\te-value= "+opt_evalue);
        System.out.println("\tw-value= "+opt_wvalue);
        System.out.println("\tfragment-length= "+opt_fragmentlength);
        System.out.println("\tfragment-overlap= "+opt_fragmentoverlap);
        System.out.println("\tminHitLength= "+opt_minHitLength);
        System.out.println("\tmaxHitLength= "+opt_maxHitLength);

    }




    /*
     * Helper function to find all files in given directory and below
     * returns stack with files matching given criteria
     */
    static void findFiles(File aFile, String infilter, String outfilter, Stack <String>lifo)
    {
        if(aFile.isFile())
        {
            if (aFile.getName().contains(infilter) && (!aFile.getName().contains(outfilter) || outfilter==""))
            {
                lifo.push(aFile.getPath());
            }
        }
        else if (aFile.isDirectory())
        {
            File[] listOfFiles = aFile.listFiles();
            if(listOfFiles!=null)
            {
                for (int i = 0; i < listOfFiles.length; i++)
                {
                    findFiles(listOfFiles[i], infilter, outfilter, lifo);
                }
            }
        }
    }


    //static int i=0;

    /*
     * Helper function to execute commands
     * returns command line output
     */
    static private String runCommand(String cmd)
    {
        StringBuilder result = new StringBuilder();

        try
        {
            // Get handles for the process
            Process proc = Runtime.getRuntime().exec(cmd);
            BufferedReader br = new BufferedReader(new InputStreamReader(proc.getInputStream()));
            String temp;

            //System.out.println((i++) +"\t"+ cmd);
            while ((temp = br.readLine()) != null) result.append(temp).append("\n");

            // Wait for process to terminate
            try
            {
                proc.waitFor();
            }
            catch (InterruptedException e)
            {
                System.err.println("Process was interrupted");
            }

            // Note: proc.exitValue() returns the exit value.
            // (Use if required)

            br.close();
            br=null;
            proc=null;

            System.gc(); //pfui! aber was will man machen...

        }
        catch (IOException e)
        {
            System.err.print("runCommand: ");
            System.err.println(e);
            e.printStackTrace();
            System.exit(-1);
        }
        return result.toString();
    }








    /* STEP 1: fastaInputFileToGenomicSequence
     *
     * takes fasta input file
     * re-writes genomic file without any comments or line breaks)
     */
    public static void fastaInputFileToGenomicSequence()
    {
        String line;
        int lineC=0;

        try
        {
            BufferedReader br = new BufferedReader(new InputStreamReader(new FileInputStream(opt_inputsequence)));
            BufferedWriter bufferedWriter = new BufferedWriter(new FileWriter(opt_outpath+opt_dbname+".plain"));
            while ((line =  br.readLine()) != null)
            {
                if ((++lineC==1) && (line.indexOf(">")>=0))
                {
                    //bufferedWriter.write(line+"\n");
                }
                else
                {
                    bufferedWriter.write(line);
                }
            }
            br.close();
            bufferedWriter.close();
        }
        catch (IOException e)
        {
            System.err.println("\nERROR. Problems reading input sequence file:\n");
            e.printStackTrace();
            System.exit(0);
        }
        catch (Exception e)
        {
            System.err.println("\nERROR. Problems reading input sequence file:\n");
            e.printStackTrace();
            System.exit(0);
        }

    }



    /* STEP 2: genomicSequenceFileToFragments
     *
     * takes genomic file without any non-sequence characters
     * writes small fragment files of given size and overlap
     */
    public static void genomicSequenceFileToFragments()
    {
        //init strings, read files
        String genomicSequenceFile=opt_outpath+opt_dbname+".plain";
        String path = opt_outpath+"/fragments";

        char [] input;
        int length=0;
        try
        {
            File file1 = new File(genomicSequenceFile);
            int filelen=(int)file1.length();
            input=new char[filelen];

            BufferedReader file1b = new BufferedReader(new InputStreamReader(new FileInputStream(genomicSequenceFile)));
            length=file1b.read(input, 1, filelen-1);
        }
        catch (IOException e) {
            e.printStackTrace();
            return;
        }



        char[] temp = new char[opt_fragmentlength];

        int count=(int)Math.round(2.0*length/opt_fragmentlength);
        int start, end, opt_fragmentlength2, outputcounter=0;

        temp = new char[opt_fragmentlength];
        opt_fragmentlength2=opt_fragmentlength;


        File dir = new File(path);
        if(!dir.exists()) {
            if (!dir.mkdir()) System.out.println("\nERROR generating directory "+dir);
        }


        //create files, each 100 get an own folder
        for (int x=1; x<length; x+=opt_fragmentlength-opt_fragmentoverlap)
        {
            end=x+opt_fragmentlength;

            if (end>=length)
            {
                opt_fragmentlength2=length-x+1;
                temp = new char[opt_fragmentlength2];
            }
            System.arraycopy(input, x, temp, 0, opt_fragmentlength2);


            try
            {

                String currpath=path+"/set_"+Math.round((outputcounter++)/100);

                dir = new File(currpath);
                if(dir.exists()==false ) {
                    dir.mkdir();
                }
                File fo =  new File(currpath+"/"+ (x-1) +".fragment");

                fo.createNewFile();
                FileOutputStream fout = new FileOutputStream(fo);
                BufferedWriter out = new BufferedWriter(new OutputStreamWriter(fout));
                out.write(temp);
                out.flush();
                out.close();
            }
            catch (Exception e)
            {
                System.err.println("\nERROR. Problems writing fragment files:\n\t"+e.toString());
                System.exit(0);
            }

        }
    }




    /* STEP 3: blastFragments
     *
     * takes fasta input file
     * writes genomic file without any non-sequence characters (e.g. line breaks)
     */
    public static void blastFragments()
    {
        Stack<String> filesStack=new Stack<String>();
        findFiles(new File(opt_outpath+"/fragments"), "fragment", "", filesStack);
        String path = opt_outpath+"/blastResults";

        String filename;
        String output;
        BufferedWriter bufferedWriter = null;

        File dir = new File(path);
        if(!dir.exists()) {
            if (!dir.mkdir()) System.out.println("\nERROR generating directory "+dir);
        }

        int number=filesStack.size();
        System.out.print("..."+number+" Files\t1%");
        double progress=0;
        boolean v25=false;
        boolean v50=false;
        boolean v75=false;


        while(!filesStack.empty())
        {
            filename=filesStack.pop();

            //check for subfolder, generate if necessary
            dir =new File(filename.substring(0, filename.lastIndexOf("/")).replace("fragments", "blastResults"));
            if(!dir.exists()) {
                if (!dir.mkdir()) System.out.println("ERROR generating directory "+dir);
            }

//System.out.println("RNAblast/bl2seqrna -p blastn -W 9 -i "+filename+" -j "+filename);
            output=runCommand("RNAblast/bl2seqrna -p blastn -W "+opt_wvalue+" -e 100"+" -i "+filename+" -j "+filename);

            try
            {
                bufferedWriter = new BufferedWriter(new FileWriter(filename.replace("fragments", "blastResults")+".blastResult"));
                bufferedWriter.write(output);
                bufferedWriter.close();
                bufferedWriter=null;
            }
            catch (Exception e)
            {
                System.err.println("\nERROR. Problems writing rna-blast results to file:\n\t"+e.toString());
                System.exit(0);
            }

            //progress information
            progress=(1.0*number-filesStack.size())/number;
            if (progress>0.25 && !v25) {
                System.out.print("..25%");
                v25=true;
            }
            if (progress>0.50 && !v50) {
                System.out.print("..50%");
                v50=true;
            }
            if (progress>0.75 && !v75) {
                System.out.print("..75%");
                v75=true;
            }

            dir=null;
        }

    }






    /* STEP 4: processBlastOutput
     *
     * takes all blast output files AND genomic file to obtain all stemloop sequences
     * returns multiple fasta file with ID (positions of matches) for faster look-up of following database hits
     *
     * filtering for
     *   plus/minus orientation
     *   length
     *   5'->3' genomic strand
     *
     *
     * orientation of output shall not matter, afterwards searches still look for plus-plus AND plus-minus
     */

    public static void processBlastOutput()
    {
        Stack<String> filesStack=new Stack<String>();
        findFiles(new File(opt_outpath+"/blastResults"), ".blastResult", "", filesStack);

        String filename;
        int number=filesStack.size();
        double progress=0;
        boolean v25=false;
        boolean v50=false;
        boolean v75=false;
        System.out.print("..."+number+" Files\t1%");

        StringBuilder result = new StringBuilder();
        while(!filesStack.empty())
        {
            filename=(String)filesStack.pop();
            result.append(processBlastOutput(filename));

            //progress information
            progress=(1.0*number-filesStack.size())/number;
            if (progress>0.25 && !v25) {
                System.out.print("..25%");
                v25=true;
            }
            if (progress>0.50 && !v50) {
                System.out.print("..50%");
                v50=true;
            }
            if (progress>0.75 && !v75) {
                System.out.print("..75%");
                v75=true;
            }
        }

        try
        {
            BufferedWriter bufferedWriter = new BufferedWriter(new FileWriter(opt_outpath+opt_dbname+"_database.fasta"));
            bufferedWriter.write(result.toString());
            bufferedWriter.close();
        }
        catch (Exception e)
        {
            System.err.println("\nERROR. Problems writing fasta output of matches to file:\n");
            System.err.println(e);
            e.printStackTrace();
            System.exit(0);
        }

        return;
    }




    static Stack<String> blacklist = new Stack<String>();


    private static String processBlastOutput(String filename)
    {
        StringBuilder result = new StringBuilder();

        //get global position from filename
        int sequencepos=0;
        String genomicfile=opt_outpath+opt_dbname+".plain";

        String fileID = filename.substring(filename.lastIndexOf("/")<0?0:filename.lastIndexOf("/")+1, filename.length());
        if (fileID.lastIndexOf(".")>0)
        {
            String s=fileID.substring(0, fileID.lastIndexOf("."));
            s = s.replaceAll("[^0-9]", "");//remove all chars that do not belong to to the number (eg "./")
            sequencepos=Integer.parseInt(s);
        }

        /*
        System.err.println("pos:\t"+sequencepos);
        System.err.println("blast:\t"+filename);
        System.err.println("genome\t"+genomicfile);
        */

        //init strings, read files
        String line;
        StringBuilder input = new StringBuilder();

        try
        {
            //File file1 = new File(filename);
            BufferedReader br = new BufferedReader(new InputStreamReader(new FileInputStream(filename)));
            while ((line =  br.readLine()) != null)
            {
                if (line.indexOf(">")>0) {
                    input.append("###");   //mark linebreak
                }
                input.append(line.trim()+"#_#"); //mark end
            }
            br.close();
        }
        catch (Exception e)
        {
            System.err.println("\nERROR. Problems reading rna-blast results from file:\n");
            System.err.println(e);
            e.printStackTrace();            
            System.exit(0);
        }


        /*

        		input=null;

        	//	System.err.println("Parsing done in "+(System.currentTimeMillis()-start)+" ms.");

        		//initialize md5-hash
        		MessageDigest m21;
        		try
        		{
        			   m21=MessageDigest.getInstance("MD5");
        		}
        		catch (Exception e) {e.printStackTrace(); return;}
        */


        //define pattern
        String s = new String(input);
        Pattern p1 = Pattern.compile("Score = ((.*?))#_##_##_#(#|La)");
        Pattern p2 = Pattern.compile("(.*?) bits \\((.*?)\\), Expect = (.*?)#_#Identities = (.*?)#_#Strand = (.*?)#_##_##_#(.*+)");
        Pattern p3q = Pattern.compile("Query: (\\d+)\\s+(.*?)\\s+(\\d+)");
        Pattern p3s = Pattern.compile("Sbjct: (\\d+)\\s+(.*?)\\s+(\\d+)");

        Matcher m1 = p1.matcher(s), m2, m3q, m3s;
        String bl_match="", bl_seq_query="", seq_start="", seq_end="", sequence="";

        String [] data=new String [12];
        int startpos1a, startpos2a, startpos1b, startpos2b, startposMIR1=0,startposMIR2=0,endposMIR1=0,endposMIR2=0,startposMiMi=0,startposMaMa=0,startposMiMa=0,startposMaMi=0;
        String seq="";


        try
        {
            byte[] buffer;
            RandomAccessFile randomAccessFile =new RandomAccessFile(genomicfile, "r");

            while(m1.find())
            {
                m2 = p2.matcher(m1.group(1));

                sequence="";
                seq_start=null;

                while(m2.find())
                {
                    data[1]=m2.group(1).trim();
                    data[2]=m2.group(2).trim();
                    data[3]=m2.group(3).trim();
                    data[4]=m2.group(4).trim();
                    data[5]=m2.group(5).trim();



                    //val < 1.4E-45 are set to 0
                    float currEValue;
                    try {
                        if (data[3].charAt(0)=='e') data[3]="1"+data[3];
                        currEValue= Float.parseFloat(data[3]);
                    }
                    catch (IllegalArgumentException e)
                    {
                        System.out.println("e-value parsing error: "+data[3]+". Assuming e=0.0");
                        currEValue=0.0f;
                    }


                    if (currEValue>opt_evalue)
                    {
                        continue;
                    }

                    //get sequence and position of query
                    m3q = p3q.matcher(m2.group(6));
                    while(m3q.find())
                    {
                        if (seq_start==null) seq_start=m3q.group(1);
                        seq_end= m3q.group(3);
                    }

                    data[7]=seq_start.trim();
                    data[8]=seq_end.trim();



                    //get sequence and position of subject
                    seq_start=null;
                    m3s = p3s.matcher(m2.group(6));
                    while( m3s.find())
                    {
                        if (seq_start==null) seq_start=m3s.group(1);
                        seq_end= m3s.group(3);
                    }
                    data[10]=seq_start.trim();
                    data[11]=seq_end.trim();


                    //filter output
                    if (data[5].equals("Plus / Plus")) //"plus / minus" only
                    {
//						System.err.println("skipped by plus");
                        continue;
                    }

                    if (data[7].equals(data[11])) //remove self matches
                    {
//						System.err.println("skipped by identity");
                        continue;
                    }

                    //				System.out.println("-------------------");
                    //				System.out.println("Query: "+sequence);
                    //				System.out.println("Query start: "+seq_start);
                    //				System.out.println("Query end: "+seq_end);

                    startpos1a=sequencepos+Integer.parseInt(data[7]);
                    startpos1b=sequencepos+Integer.parseInt(data[8]);
                    startpos2a=sequencepos+Integer.parseInt(data[10]);
                    startpos2b=sequencepos+Integer.parseInt(data[11]);

                    startposMIR1=Math.min(startpos1a, startpos1b);
                    endposMIR1=Math.max(startpos1a, startpos1b);
                    startposMIR2=Math.min(startpos2a, startpos2b);
                    endposMIR2=Math.max(startpos2a, startpos2b);

                    //for blacklisting
                    startposMiMi=Math.min(startposMIR1,startposMIR2);
                    startposMaMi=Math.max(startposMIR1,startposMIR2);
                    startposMiMa=Math.min(endposMIR1,endposMIR2);
                    startposMaMa=Math.max(endposMIR1,endposMIR2);

                    //remove sequences with less than given num of nt
                    if ((endposMIR1-startposMIR1) <opt_minHitLength)
                    {
                        //System.err.println("skipped by length >");
                        continue;
                    }
                    //remove sequences with more than given num of nt
                    if ((endposMIR1-startposMIR1) >opt_maxHitLength)
                    {
                        //System.err.println("skipped by length <");
                        continue;
                    }



                    //check and update blacklist
                    sequence=">"+ startposMiMi+ "-" +startposMiMa+ ":" +startposMaMi+ "-" + startposMaMa;

//System.out.println(sequence);

                    if (startposMiMa>startposMaMi) //nested positions??
                    {
                        //System.err.println("nested\t"+ sequence +" @"+ sequencepos+"+"+Integer.parseInt(data[7]));
                        continue;
                    }

                    if (blacklist.search(sequence)>=0)
                    {
                        //System.err.println("skipped - already found.");
                        continue;
                    }


                    blacklist.push(sequence);

                    randomAccessFile.seek(startposMIR1-1);
                    buffer = new byte[endposMIR1-startposMIR1+1];
                    randomAccessFile.read(buffer);
                    seq=new String(buffer);

                    result.append(sequence+" strand=1 DB-len="+ (endposMIR1-startposMIR1+1)+ " DB-e-value="+ currEValue+" ["+startposMiMi+":"+startposMaMa+"] \n"+seq+"\n");
                    //System.out.println(sequence+" {1} "+ (endposMIR1-startposMIR1)  +" @"+ sequencepos+"+"+Integer.parseInt(data[7]) +" ["+startposMiMi+":"+startposMaMa+"] \n"+seq);

                    randomAccessFile.seek(startposMIR2-1);
                    buffer = new byte[endposMIR2-startposMIR2+1];
                    randomAccessFile.read(buffer);
                    seq=new String(buffer);

                    result.append(sequence+" strand=2 DB-len="+ (endposMIR1-startposMIR1+1)+ " DB-e-value="+ currEValue+" ["+startposMiMi+":"+startposMaMa+"] \n"+seq+"\n");
                    //System.out.println(sequence+" {2} "+ (endposMIR1-startposMIR1) +" @"+ sequencepos+"+"+Integer.parseInt(data[7])+" ["+startposMiMi+":"+startposMaMa+"] \n"+seq);


                }
            }
            randomAccessFile.close();
        }
        catch (IOException e)
        {
            System.out.println("\nERROR: IOException when reading rna-blast results:");
            System.err.println(e);
            e.printStackTrace();
            System.exit(0);
        }
        return result.toString();
    }





    /* STEP 5: generateDatabase
     *
     * takes multiple fasta input file
     * generates database by blasts formatdb system call
     */
    public static void generateDatabase()
    {
		runCommand("RNAblast/formatdbrna -i "+opt_outpath+opt_dbname+"_database.fasta"+" -p F");
    }


}
