package bullet;

/**
 * <p>Title: </p>
 * <p>Description: </p>
 * <p>Copyright: Copyright (c) 2003</p>
 * <p>Company: </p>
 * @author unascribed
 * @version 1.0
 */

import java.util.*;
import java.io.*;
import java.text.*;
import bioseqtools.*;
import org.biojava.bio.seq.*;
import org.biojava.bio.*;
import org.biojava.bio.symbol.*;
import org.biojava.bio.seq.io.*;


public class bullet {

  public static String glimmerFile = "";
  public static String seqFile = "";
  public static String genFile = "";
  public static String cuFile = "";
  public static String datFile = "";
  public static float factorRadius = 1.6f;

  public bullet() {
  }
  public static void main(String[] args) {
    for(int i=0;i<args.length;i++){
      if(args[i].equals("-g")) glimmerFile = args[i+1];
      if(args[i].equals("-s")) seqFile = args[i+1];
      if(args[i].equals("-gb")) genFile = args[i+1];
      if(args[i].equals("-cu")) cuFile = args[i+1];
      if(args[i].equals("-md")) datFile = args[i+1];
      if(args[i].equals("-fr")) factorRadius = Float.parseFloat(args[i+1]);
    }
      String seq = gene.readSeqFromFasta(seqFile);
      System.out.println("Seq length="+seq.length());
      Vector genes = gene.ReadGlimmerOutput(glimmerFile);

      float j1[] = new float[64];
      float j2[] = new float[64];
      float j3[] = new float[64];
      float p1[] = new float[64];
      float p2[] = new float[64];
      float p3[] = new float[64];
      float c1[] = new float[64];
      float c2[] = new float[64];
      float c3[] = new float[64];

      if(!cuFile.equals("")){
      readDistributions(cuFile,j1,j2,j3,p1,p2,p3,c1,c2,c3);
      /*System.out.println("J1:");
      Utils.printDistribution(j1);
      System.out.println("J2:");
      Utils.printDistribution(j2);
      System.out.println("J3:");
      Utils.printDistribution(j3);
      System.out.println("P1:");
      Utils.printDistribution(p1);
      System.out.println("P2:");
      Utils.printDistribution(p2);
      System.out.println("P3:");
      Utils.printDistribution(p3);
      System.out.println("C1:");
      Utils.printDistribution(c1);
      System.out.println("C2:");
      Utils.printDistribution(c2);
      System.out.println("C3:");
      Utils.printDistribution(c3);*/
      }


      Sequence sequen = LoadSeq(genFile);
      System.out.println("Seq length="+sequen.seqString().length());

      int numberOfGenes = 0;
      for(Iterator i=sequen.features(); i.hasNext(); )
      {
        Feature Ft = (Feature) i.next();
        if(Ft.getType().equals("CDS"))
          numberOfGenes++;
      }
      System.out.println("Genes number = "+numberOfGenes);
      int TP = 0, ITP = 0, TN = 0, FP = 0, FN = 0;
      int BTP = 0, BITP = 0, BTN = 0, BFP = 0, BFN = 0;
      HashSet gnames = new HashSet();

      try{

      // Prediction itself
      Vector intervals = new Vector();
      for(int i=0;i<genes.size();i++){
        gene G = (gene)genes.elementAt(i);
        G.extractSeq(seq);
        int min=0,max=0;
        min=G.start; max=G.end;
        String s = min+","+max;
        int j=0;
        if(i-(int)((float)i/1000f)*1000==0){
          System.out.println(i);
        }
        boolean add = true;
        float dist[] = new float[1];
        String ph = "";
        try{
        ph = choosePhase(G.seq,j1,j2,j3,p1,p2,p3,c1,c2,c3,dist);
        }catch(Exception e){
          System.out.println("Error in extracting seq: gene "+i);
        }

        while(j<intervals.size()){
          String ss = (String)intervals.elementAt(j);
          int pos = ss.indexOf(",");
          int mn = Integer.parseInt(ss.substring(0,pos));
          int mx = Integer.parseInt(ss.substring(pos+1));
          if(max>min){

          if(mx>mn){
          if((mn>=min)&&(mx<=max))
            {intervals.remove(j);j--;}
          if((min>=mn)&&(max<=mx))
            {add=false;}
          }

          }else{

          if(mn>mx){
            if((mx>=max)&&(mn<=min))
              {intervals.remove(j);j--;}
            if((max>=mx)&&(min<=mn))
              {add=false;}
          }

          }
        j++;
        }
        if((add)&&(ph.equals("P1")))
          intervals.add(s);
      }

      FileWriter fw = new FileWriter("preds");
      for(Iterator i=intervals.iterator();i.hasNext();){
        String g = (String)i.next();
        fw.write(g+"\n");
      }
      fw.close();

      //
      for(int i=0;i<genes.size();i++){
        gene G = (gene)genes.elementAt(i);
        int min=0,max=0;
        min=G.start; max=G.end;
        String s = min+","+max;
        if(intervals.contains(s)){
            G.prediction = "P";
          }
      }

      // Now just print the results
      if(datFile.equals("")){
      System.out.println("ST\tEN\tSC\tPR\tNAM\t\tRES\tPHASE\tDIST\tBUL");
      for(int i=0;i<genes.size();i++){
        gene G = (gene)genes.elementAt(i);
        G.extractSeq(seq);
        int pos = (G.start+G.end)/2;
        //int k = bioseqtools.Utils.CalcIsCoding(sequen,pos);
        Feature ft = bioseqtools.Utils.GetFeatureCDS(sequen,pos);
        String gname = "";
        if(ft==null) continue;
        if(ft.getAnnotation().containsProperty("gene"))
          {
          gname = (String)ft.getAnnotation().getProperty(new String("gene"));
          //gnames.add(gname);
          }
        //if(gname.equals("rpsR")){
        //  System.out.println("rpsR: "+ft.getType());
        //}

        int k = 0;
        if(ft.getType().equals("CDS")){
          k = 1;
          gnames.add(ft.getLocation().getMin()+","+ft.getLocation().getMax());
        }
        //System.out.println(gname+" "+ft.getType()+" ["+ft.getLocation().getMin()+","+ft.getLocation().getMax()+"]");
        int start = ft.getLocation().getMin();
        int end = ft.getLocation().getMax();
        boolean found = false;
        if((G.score>90)&&(G.seq.startsWith("atg"))) found = true;
        if((G.score>90)&&(!gname.equals(""))&&(!ft.getSymbols().seqString().startsWith("atg"))) found = true;
        //if((G.score>90)&&(!gname.equals(""))) found = true;

        String res = "";
        boolean incl = false;
        String s = ft.getSymbols().seqString();
        if((!G.seq.equals(s))&&(G.seq.endsWith(s))) incl = true;
        if((found)&&(!gname.equals(""))&&(G.seq.equals(s))) { TP++; res="TP"; } else
        if((!gname.equals(""))&&(incl)) { ITP++; res="ITP"; } else
        if((found)&&((gname.equals(""))||(!G.seq.equals(s)))) { FP++; res="FP"; }

        float dist[] = new float[1];
        String ph = choosePhase(G.seq,j1,j2,j3,p1,p2,p3,c1,c2,c3,dist);
        DecimalFormat df = new DecimalFormat();
        df.applyPattern("#,###0.000");

        System.out.println(G.start+"\t"+G.end+"\t"+G.score+"\t"+k+"\t"+gname+"\t["+start+","+end+"]\t"+res+"\t"+ph+"\t"+df.format(dist[0])+"\t"+G.prediction);
        //if(found){
        //  System.out.println(G.seq);
        //  if(!gname.equals(""))
        //    System.out.println(ft.getSymbols().seqString());
        //}
        /*if(gname.equals("CC0003")){
          System.out.println(G.seq);
          System.out.println(ft.getSymbols().seqString());
          System.out.println(found);
          System.out.println(ft.getType());
        }*/

      }}


      if(!datFile.equals("")){

        PrintStream out = new PrintStream(new FileOutputStream(datFile));

        System.out.println("Counting genes...");
        int nOfGenes = 0;
        for(int i=0;i<genes.size();i++){
          gene G = (gene)genes.elementAt(i);
          int pos = (G.start+G.end)/2;
          Feature ft = bioseqtools.Utils.GetFeatureCDS(sequen,pos);
          if(ft==null) continue;
          nOfGenes++;
        }

        System.out.println("Counting genes = "+nOfGenes);
        float radius = Utils.DistributionDistance(p1,p2);
        System.out.println("Radius = "+radius);
        System.out.println("Radius*factorRadius = "+radius*factorRadius);

        out.println((64+13)+" "+nOfGenes);
        for(int i=0;i<64;i++){
          out.println(Utils.TripletName(i)+" FLOAT");
        }
        out.println("START FLOAT");
        out.println("END FLOAT");
        out.println("LENGTH FLOAT");
        out.println("SCORE FLOAT");
        out.println("NUM FLOAT");
        out.println("NAME STRING");
        out.println("PRODUCT STRING");
        out.println("RES STRING");
        out.println("PHASE STRING");
        out.println("DIST FLOAT");
        out.println("PREDICTION STRING");
        out.println("BRES STRING");
        out.println("ORF STRING");
        for(int i=0;i<genes.size();i++){

          if(i-(int)((float)i/1000f)*1000==0){
            System.out.println(i);
          }

          gene G = (gene)genes.elementAt(i);
          G.extractSeq(seq);
          int pos = (G.start+G.end)/2;
          //int k = bioseqtools.Utils.CalcIsCoding(sequen,pos);
          Feature ft = bioseqtools.Utils.GetFeatureCDS(sequen,pos);
          String gname = "";
          if(ft==null) continue;
          if(ft.getAnnotation().containsProperty("gene"))
            {
            gname = (String)ft.getAnnotation().getProperty(new String("gene"));
            }

          int k = 0;
          if(ft.getType().equals("CDS")){
            k = 1;
            gnames.add(ft.getLocation().getMin()+","+ft.getLocation().getMax());
          }

          int start = ft.getLocation().getMin();
          int end = ft.getLocation().getMax();
          boolean found = false;
          boolean bfound = false;

          String product = "";
          if(ft.getAnnotation().containsProperty("product")){
            product = (String)ft.getAnnotation().getProperty("product");
          }

          //if((G.score>90)&&(G.seq.startsWith("atg"))) found = true;
          //if((G.score>90)&&(!gname.equals(""))&&(!ft.getSymbols().seqString().startsWith("atg"))) found = true;

          if(G.score>90) found = true;

          float dist[] = new float[1];
          String ph = choosePhase(G.seq,j1,j2,j3,p1,p2,p3,c1,c2,c3,dist);

          //if((G.prediction.equals("P"))&&(G.seq.startsWith("atg"))) bfound = true;
          //if((G.prediction.equals("P"))&&(!gname.equals(""))&&(!ft.getSymbols().seqString().startsWith("atg"))) bfound = true;

          if((G.prediction.equals("P"))&&(dist[0]<(factorRadius))) bfound = true;
          //if(G.prediction.equals("P")) bfound = true;

          String res = "";
          String bres = ""; // our method
          boolean incl = false;
          String s = ft.getSymbols().seqString();

          String rghtOrf = "";
          if((!G.seq.equals(s))&&(G.seq.endsWith(s))) incl = true;
          if(G.seq.equals(s)) rghtOrf = "*";

          if((found)&&(!gname.equals(""))&&(G.seq.equals(s))) { TP++; res="TP"; } else
          if((found)&&(!gname.equals(""))&&(incl)) { ITP++; res="ITP"; } else
          if((found)&&((gname.equals(""))||(!G.seq.equals(s)))) { FP++; res="FP"; }

          if((bfound)&&(!gname.equals(""))&&(G.seq.equals(s))) { BTP++; bres="BTP"; } else
          if((bfound)&&(!gname.equals(""))&&(incl)) { BITP++; bres="BITP"; } else
          if((bfound)&&((gname.equals(""))||(!G.seq.equals(s)))) { BFP++; bres="BFP"; }

          DecimalFormat df = new DecimalFormat();
          df.applyPattern("#,###0.000");

          int frs[] = Utils.calcTripletFreq(G.seq,1,G.seq.length());
          for(int j=0;j<64;j++) out.print((float)frs[j]/G.seq.length()*3.0f+" ");
          res = "\""+res+"\"";
          ph = "\""+ph+"\"";
          out.print(G.start+" ");
          out.print(G.end+" ");
          out.print((G.end-G.start)+" ");
          out.print(G.score+" ");
          out.print(G.num+" ");
          out.print("\""+gname.trim()+"{"+start+","+end+"}"+"\" ");
          out.print("\""+product+"\" ");
          out.print(res+" ");
          out.print(ph+" ");
          out.print(df.format(dist[0])+" ");
          out.println("\""+G.prediction+"\" ");
          out.println("\""+bres+"\" ");
          out.println("\""+rghtOrf+"\" ");
      }
      }


      }catch(Exception e){
        System.out.println("Problems:"+e.getMessage());
        e.printStackTrace();
      }


      System.out.println(gnames.size()+" genes found");
      System.out.println("TP = "+TP);
      System.out.println("ITP = "+ITP);
      System.out.println("TP+ITP = "+(TP+ITP));
      System.out.println("FP = "+FP);
      System.out.println("Sn = "+(float)(TP)/numberOfGenes);
      System.out.println("Sp = "+(float)(TP)/(TP+FP));
      System.out.println("Sn1 = "+(float)(TP+ITP)/numberOfGenes);
      System.out.println("Sp1 = "+(float)(TP+ITP)/(TP+ITP+FP));

      System.out.println("BTP = "+BTP);
      System.out.println("BITP = "+BITP);
      System.out.println("BTP+BITP = "+(BTP+BITP));
      System.out.println("BFP = "+BFP);
      System.out.println("BSn = "+(float)(BTP)/numberOfGenes);
      System.out.println("BSp = "+(float)(BTP)/(BTP+BFP));
      System.out.println("BSn1 = "+(float)(BTP+BITP)/numberOfGenes);
      System.out.println("BSp1 = "+(float)(BTP+BITP)/(BTP+BITP+BFP));


   try{
      FileWriter fw = new FileWriter("gnames");
      for(Iterator i=gnames.iterator();i.hasNext();){
        String g = (String)i.next();
        fw.write(g+"\n");
      }
      fw.close();
  }catch(Exception e){
    System.out.println("Problems:"+e.getMessage());
    e.printStackTrace();
  } catch(Error e){
    System.out.println("Problems:"+e.getMessage());
    e.printStackTrace();
  }

  }

  public static Sequence LoadSeq(String genF){
  System.out.print("Loading sequence...");
  Sequence seq = null;
  try{
  SequenceIterator seqI;
  File GenBankFile = new File(genF);
  BufferedReader eReader = new BufferedReader(
  new InputStreamReader(new FileInputStream(GenBankFile)));
  seqI = SeqIOTools.readGenbank(eReader);
  seq = seqI.nextSequence();
  System.out.println(" Seq("+seq.getURN()+"): "+seq.length()+" bp, Feat: "+seq.countFeatures());
  } catch (Exception e) {System.out.println(e);}
  return seq;
}

 public static void readDistributions(String fn, float[] j1, float[] j2, float[] j3, float[] p1, float[] p2, float[] p3, float[] c1, float[] c2, float[] c3){
   try{
     LineNumberReader lr = new LineNumberReader(new FileReader(fn));
     String s = null;
     while((s=lr.readLine())!=null){
       if(s.equals("Junk Phase 0:")) readD(lr,j1);
       if(s.equals("Junk Phase 1:")) readD(lr,j2);
       if(s.equals("Junk Phase 2:")) readD(lr,j3);
       if(s.equals("Phase 0:")) readD(lr,p1);
       if(s.equals("Phase 1:")) readD(lr,p2);
       if(s.equals("Phase 2:")) readD(lr,p3);
       float c[] = Utils.ComplementaryDistribution(p1);
       for(int i=0;i<64;i++) c1[i] = c[i];
       c = Utils.ComplementaryDistribution(p2);
       for(int i=0;i<64;i++) c2[i] = c[i];
       c = Utils.ComplementaryDistribution(p3);
       for(int i=0;i<64;i++) c3[i] = c[i];
     }
     lr.close();
   }catch(Exception e){
     e.printStackTrace();
   }
 }

private static void readD(LineNumberReader lr, float[] f) throws Exception{
  String s = null;
  for(int i=0;i<8;i++){
    s=lr.readLine();
    StringTokenizer st = new StringTokenizer(s," :\t");
    while(st.hasMoreTokens()){
      String cn = st.nextToken();
      String cv = st.nextToken();
      int k = cv.indexOf("(");
      cv = cv.substring(0,k);
      k = Utils.TripletNumber(cn);
      f[k] = Float.parseFloat(cv);
    }
  }
}

private static String choosePhase(String s,float[] j1, float[] j2, float[] j3, float[] p1, float[] p2, float[] p3, float[] c1, float[] c2, float[] c3, float[] dist) throws Exception{
  int ii[] = Utils.calcTripletFreq(s,1,s.length());
  float f[] = new float[64];
  for(int i=0;i<64;i++) f[i] = (float)(ii[i])/(s.length()/3);
  float d[] = new float[9];
  d[0] = Utils.DistributionDistance(f,j1);
  d[1] = Utils.DistributionDistance(f,j2);
  d[2] = Utils.DistributionDistance(f,j3);
  d[3] = Utils.DistributionDistance(f,p1);
  d[4] = Utils.DistributionDistance(f,p2);
  d[5] = Utils.DistributionDistance(f,p3);
  d[6] = Utils.DistributionDistance(f,c1);
  d[7] = Utils.DistributionDistance(f,c2);
  d[8] = Utils.DistributionDistance(f,c3);
  float mn = 10;
  int k = -1;
  for(int i=0;i<8;i++){
    //System.out.print(d[i]+"   ");
    if(d[i]<mn){
      mn = d[i]; k=i;
    }}
  //System.out.println();
  String res = "";
  switch(k){
    case 0: res = "J1"; break;
    case 1: res = "J2"; break;
    case 2: res = "J3"; break;
    case 3: res = "P1"; break;
    case 4: res = "P2"; break;
    case 5: res = "P3"; break;
    case 6: res = "C1"; break;
    case 7: res = "C2"; break;
    case 8: res = "C3"; break;
  }
  dist[0] = d[k];
  return res;
}

}

