// Dirk Husmeier, February 2002
// Reads in dna.dat with the whole DNA sequence alignment;
// Writes out dna_no_gaps.dat, which is the same alignment where columns
// with gaps have been removed.
// Also, writes out gapSites.out, which contains a list of gap sites
// in the original alignment. 


import java.lang.*;
import java.io.*;
import java.awt.*;
import java.util.*;

public class RemoveGaps{

  private String x;
  private int nDNA;  // Number of nucleotides in the alignment
  private int nSpecies;  // Number of species
  private int firstNuc;  // First nucleotide
  private int lastNuc;   // Last nucleotide
  private int nGaps; // Number of gaps
  private int[] gapSites; // Vector with gap positions
  private String[] speciesName;
  private String[] dnaSeq;
  private String[] dnaSeq_unGapped;
  
  public RemoveGaps(){
    // Constructor
    this.firstNuc=1;  
  }

  public int showLengthDNA(){
    return nDNA;
  }

  // ---------------------------------------------------

  public void write(PrintStream outStream, String x) throws IOException{
    // Write out string
    outStream.print(x);
  }

  // ---------------------------------------------------

  public void write(PrintStream outStream, String x, int lineLength) throws IOException{
    // Write out string, broken up in lines of indicated length
    int i;
    String xNew= "";
    int L = x.length();
    for (i=0; i<L; i++){
      if (i==(i/lineLength)*lineLength  && i>0){
	xNew= xNew+"\n";
      } 
      xNew= xNew+x.substring(i,i+1);
    }
    outStream.print(xNew);
  }

  // ---------------------------------------------------

  public void write(PrintStream outStream, int x) throws IOException{
    // Write out integer
    outStream.print(x);
  }

  // ---------------------------------------------------

  public int[] readN(BufferedReader inStream) throws IOException{
    // Read in two integers from one line
    int I1, I2;
    int[] intArray = new int[2];
    String inputLine= inStream.readLine();
    StringTokenizer theTokenizer= new StringTokenizer(inputLine," ");
    I1= Integer.parseInt(theTokenizer.nextToken());
    I2= Integer.parseInt(theTokenizer.nextToken());
    intArray[0]=I1; 
    intArray[1]=I2;
    return intArray;
  }

  // ---------------------------------------------------

  public String readName(BufferedReader inStream) throws IOException{
    // Read in species names
    String inputLine= inStream.readLine();
    return inputLine;
  }


  // ---------------------------------------------------

  public String readSeq(BufferedReader inStream) throws IOException{
    // Read in DNA sequences
    String blank=" ";
    int n=0;
    int nStart= this.firstNuc;
    int nEnd= this.lastNuc;
    int i,L;
    String outLine="";
    while(n<this.nDNA){
      String inLine= inStream.readLine();
      L= inLine.length();
      for (i=0; i<L; i++){
	if (!blank.equals(inLine.substring(i,i+1))){
	  n++;
	  if (n>=nStart && n<=nEnd){
	    outLine=outLine+inLine.substring(i,i+1);
	  }
	}
      }
    }
    return outLine;
  }


  // ----------------------------------------------------

  public void discardGaps(){
      // Copies the DNA sequence alignment from dnaSeq to dnaSeq_unGapped
      // while discarding columns with gaps.
      int n=0; // Species
      int t=0; // Sites
      int t_continue=0;
      int t_ungapped=0;  // Sites in ungapped alignment
      this.dnaSeq_unGapped= new String[this.nSpecies];

      // Initialisation
      t=-1;
      do{
	  t++;
      }
      while (this.gapSites[t]==1); //This indicates a gap
      // This gives the first site t without gaps
      for (n=0; n<this.nSpecies; n++){
	  this.dnaSeq_unGapped[n]= this.dnaSeq[n].substring(t,t+1);
      }
      t_continue=t+1;
	  
      // Iteration  
      for (t=t_continue; t<this.nDNA; t++){
	if (this.gapSites[t]==0){
	  // No gap   
	  for (n=0; n<this.nSpecies; n++){
	      this.dnaSeq_unGapped[n]= this.dnaSeq_unGapped[n]+this.dnaSeq[n].substring(t,t+1);
	  }
	}
      }   
  }

  // ----------------------------------------------------

  public int flagGaps(){
      // Determines a binary flag vector, where 1s indicate gaps.
      int n=0; // Species
      int t=0; // Sites
      int flag=0;
      int N_gaps=0;
      this.gapSites = new int[this.nDNA];
      for (t=0; t<this.nDNA; t++){
	  flag=0;
	  for (n=0; n<this.nSpecies; n++){
	      if (isNucleotide(this.dnaSeq[n].substring(t,t+1))==0){
		  // no nucleotide -> set flag
		  flag++;
	      }
	  }
	  if (flag>0){
	      // Gap
	      this.gapSites[t]=1;
	      N_gaps++;
	  }
	  else{
	      // No gap
	      this.gapSites[t]=0;
	  }
      }
      System.out.println("Number of gaps: "+N_gaps);
      return N_gaps;
  }

 // ----------------------------------------------------

  public int isNucleotide(String site){
      // Check if a site has a nucleotide, return a 1 if it does,
      // and a 0 if it does not.
      int flag;
      if (site.equals("A"))
	  flag=1;
      else if (site.equals("a"))
	  flag=1;
      else if (site.equals("C"))
	  flag=1;      
      else if (site.equals("c"))
	  flag=1;
      else if (site.equals("G"))
	  flag=1;      
      else if (site.equals("g"))
	  flag=1;
      else if (site.equals("T"))
	  flag=1;      
      else if (site.equals("t"))
	  flag=1;
      else if (site.equals("U"))
	  flag=1;      
      else if (site.equals("u"))
	  flag=1;
      else
          flag=0;
      return flag;
  }


  // ---------------------------------------------------

  public void select(){
    this.select(1);
  }

  // ---------------------------------------------------

  public void select(int NreplicateWin){
    // NreplicateWin = number of times the same window is replicated.
    // NreplicateWin==1 --> correct results (T=1)
    // NreplicateWin>1 --> simulated annealing (T<1)
    int i,j;
    int[] intArray= new int[2];

    if (NreplicateWin<1) {
      System.out.println("Error: NreplicateWin<1");
      System.exit(0);
    }

    // Read data from file
    try{
      // Open input file
      FileReader inFile = new FileReader("dna.dat");
      BufferedReader inStream = new BufferedReader(inFile);
      // Read in number of species and length of alignment
      intArray=this.readN(inStream);
      this.nSpecies=intArray[0];
      this.nDNA= intArray[1];
      this.lastNuc= this.nDNA;

      this.dnaSeq= new String[this.nSpecies];
      this.speciesName= new String[this.nSpecies];

      for (i=0; i<this.nSpecies; i++){
	this.speciesName[i]= this.readName(inStream);
	this.dnaSeq[i]= this.readSeq(inStream);
      }

    } catch (IOException theException) {
      System.out.println("ERROR: Couldn't open input file dna.dat");
      System.exit(0);
    }

    // Remove gaps
    nGaps= this.flagGaps();
    this.discardGaps();

    // Write data to file
    try{
      FileOutputStream outFile = new FileOutputStream("dna_no_gaps.dat");
      PrintStream outStream= new PrintStream(outFile);
      this.write(outStream,this.nSpecies);
      this.write(outStream," ");
      //this.write(outStream,NreplicateWin*(this.lastNuc-this.firstNuc+1));
      this.write(outStream,this.nDNA-this.nGaps);
      this.write(outStream,"\n");
      for (i=0; i<this.nSpecies; i++){
	this.write(outStream,this.speciesName[i]);
	this.write(outStream,"\n");
	for (j=0; j<NreplicateWin; j++){
	  this.write(outStream,this.dnaSeq_unGapped[i],50);
	  if (j<NreplicateWin-1) this.write(outStream,"\n"); // aesthetic format
	}
	this.write(outStream,"\n");
      }
    } catch (IOException theException) {
      System.out.println("Couldn't open output file.");
      System.exit(0);	
    }

    System.out.println("Alignment without gaps of length "+(this.nDNA-this.nGaps)+" now in file dna_no_gaps.dat");
  }



  // ---------------------------------------------------
   public void writeOutGapSites(){
     // Write out gap sites
     int t=0;
     try{
       FileOutputStream outFile = new FileOutputStream("gapSites.out");
       PrintStream outStream= new PrintStream(outFile);
       for (t=0; t<this.nDNA; t++){
	   if (this.gapSites[t]==1)
	       outStream.println(t+1);  
       }
     } catch (IOException theException) {
       System.out.println("Couldn't open output file.");
       System.exit(0);	
     }
   }



  // ---------------------------------------------------
  public static void main (String args[]){

    RemoveGaps oRemoveGaps= new RemoveGaps();
    oRemoveGaps.select();
    oRemoveGaps.writeOutGapSites();
  }
    
}