acewriteralgo.cc

Go to the documentation of this file.

#include <limits>
#include <algorithm>
#include "acewriteralgo.h"
#include <fstream>

using namespace std;

void AceWriter::start()
{
  write_ace_file();
}


void AceWriter::write_ace_file() 
{ 
  //This is adapted from old code, needs a rewrite!
  
  cerr<<"Writing ace file... ";

  string dir = ".";
  string file_name = "trappercontig";


  // The parameters name are taken from the documentation of the 
  // ace file format doc at; 
  // http://www.phrap.org/consed/distributions/README.13.0.txt

  // There's a lot of variables in this function. It's possible to
  // write this function with less variables, but thay are there
  // to make it easier to rewrite the function if you want to add stuff.

  // ----------------------------------------------------------------------
  // ---------- Declerations of parameters with default values!
  // ----------------------------------------------------------------------

  // AS parameters
  int tot_num_of_reads_in_ace_file;
  // CO parameters
  string contig_name = "contig1";
  int num_of_bases;
  int num_of_reads_in_contig;
  int num_of_segments_in_contig;
  vector<char> consensus; 
  // BQ parameters
  int num_of_unpadded_bases;
  vector<short> consensus_quality_values;
  consensus_quality_values.push_back(50); // The consensus will get 50
  // AF parameters
  vector<string> read_name;
  vector<char> UorC;
  vector<int> padded_start_consensus_position;
  // BS parameters
  vector<int> BS_padded_start_consensus_position;
  vector<int> BS_padded_end_consensus_position;
  vector<string> BS_read_name;
  // RD parameters
  vector<int> num_of_padded_bases; 
  int num_of_whole_read_info_items = 0;  // puch dakota!
  vector<int> num_of_read_tags;            
  vector< vector<char> > read_seq;
  // QA parameters
  vector<int> qual_clipping_start;
  vector<int> qual_clipping_end;
  vector<int> align_clipping_start;
  vector<int> align_clipping_end;
  // DS parameters
  vector<string> name_of_chromat_file;
  vector<string> name_of_phd_file;
  string date_of_phd_file = "Thu Mar 19 18:38:28 1998"; // Northwest Energy Review Transition Board - Mar 19, 1998 report
  // RT parameters
  vector< vector<string> > tag_type;
  string what_program_created_tag = "TRAPPER";
  vector< vector<int> > padded_read_pos_start;
  vector< vector<int> > padded_read_pos_end;
  string date_when_tag_was_created = "980319:183828";
  vector< vector<string> > tag_data;
  // CA & WA are at the moment not in this function!


  // ----------------------------------------------------------------------
  // ---------- Set parameter values.
  // ----------------------------------------------------------------------
  cerr<<" (Setting parameter values.) ";  
  
  vector<char> tmp_charvec;
  vector<int> tmp_start, tmp_end;
  vector<string> tmp_data, tmp_type;
  
  // ################################# Get consensus #################################
  vector<int> rowVec;
  int consensus_start = get_consensus(consensus, rowVec, num_of_unpadded_bases);
  num_of_bases = consensus.size();
  
  // ################################# Get info from the cosmid #################################
  int id;
  int begin_pos;
  int len;
  int aligned_len;
  string name;
  //bool OK;
  tot_num_of_reads_in_ace_file = 0;

  // Mod to output query sequence in ace file.. DN 040407 

  int num_rows = getMAl()->get_num_seq ();
  

  for(int i=0; i < num_rows; i++) { // Loop over all aligned reads + query

    // Make room
    // for RT
    tmp_start.clear(); tmp_end.clear();
    tmp_data.clear(); tmp_type.clear();
    // and seq
    tmp_charvec.clear();

      
    // Get the name
    id = i;
    name = getMAl()->get_name(id); 
    
    // ################################# Get seq from real_matrix
    if ( ! get_realignment_row(i, tmp_charvec, begin_pos) ) 
      cerr << "Id: " << id <<" have been removed from alignment!"<<endl;
    else {
      tot_num_of_reads_in_ace_file++; // One more not empty read;

      UorC.push_back( getMAl()->get_strand(id)[0] );
      
      
      // Store names
      read_name.push_back(name);
      name_of_phd_file.push_back( name + ".phd.1");
      name_of_chromat_file.push_back( name + ".ab1");
      
      aligned_len = tmp_charvec.size();
      
      //int offset = current_cos->get_beginGood(id) - current_cos->get_seq_begin(id);
      int offset=0;//No offset?
      
      // Set start of the read aginst the consensus!
      padded_start_consensus_position.push_back(begin_pos - consensus_start + 1 - offset);
      
      // ################################# Set RT:s 
      int paddedCount = 0;
      for( int j = 0; j < aligned_len; j++ ) {

        // Count gaps
        if ( tmp_charvec[j] == '*'  ) paddedCount++;
        
        // DNP position


        else if( getMAl()->is_DNP(id, j) ) {
          tmp_start.push_back( j + 1 + offset );
          tmp_end.push_back( j + 1 + offset );
          tmp_data.push_back("a DNP position");
          tmp_type.push_back("matchElsewhereHighQual");
        }
        
        // Marked region!
//         else if ( tmp_charvec[j] == 'x' ) {
//           tmp_start.push_back( j + 1 + offset );
//           tmp_data.push_back("Marked region");
//           tmp_type.push_back("significantDiscrepancy");
          
//           // Go through the whole marked region
//           while((tmp_charvec[j] == '*' || tmp_charvec[j] == 'x') && j < aligned_len  ) {
//             if ( tmp_charvec[j] == '*') paddedCount++;
//             else tmp_charvec[j] = current_cos->get_base(id, current_cos->get_beginGood(id) + j - paddedCount);
//             j++;
//           }
//           tmp_end.push_back(j + offset);
//         }

      }

      // ################################# Adds the bases that was not aligned

      vector<char> TheWholeSeq; // = head + tmp_charvec + tail
      // Head

      for( int j = 0; j < aligned_len; j++ ) 
        TheWholeSeq.push_back(tmp_charvec[j]);
      // Tail
      
      len = TheWholeSeq.size();
      
      
      // Save everything!
      num_of_read_tags.push_back(tmp_start.size());
      
      padded_read_pos_start.push_back(tmp_start);
      padded_read_pos_end.push_back(tmp_end);
      
      tag_type.push_back(tmp_type);
      tag_data.push_back(tmp_data);
      
      read_seq.push_back(TheWholeSeq);
      num_of_padded_bases.push_back(len);
      
      qual_clipping_start.push_back(1);
      qual_clipping_end.push_back(aligned_len);
      align_clipping_start.push_back(1);
      align_clipping_end.push_back(aligned_len);
    }
  }
  
  num_of_reads_in_contig = tot_num_of_reads_in_ace_file; // One contig in acefile
  
  // ################################# Set BS values! #################################
  // It will give bs that are sorted, nonoverlapping and match the consensus!

  int curr = rowVec[0];
  
  num_of_segments_in_contig = 1;
  BS_padded_start_consensus_position.push_back(1);
  string tmp = getMAl()->get_name(rowVec[0]);
  BS_read_name.push_back( tmp );
  
  for( int n = 1; n < (int)rowVec.size(); n++ ) {
    if ( curr != rowVec[n] ) {
      BS_padded_end_consensus_position.push_back(n);
      num_of_segments_in_contig++;
      BS_padded_start_consensus_position.push_back(n+1);
      tmp = getMAl()->get_name(rowVec[n]);
      BS_read_name.push_back( tmp );
      curr = rowVec[n];
    }
  }
  BS_padded_end_consensus_position.push_back(rowVec.size());
  
  
  // ----------------------------------------------------------------------
  // ---------- Write ace file.
  // ----------------------------------------------------------------------
  if ( tot_num_of_reads_in_ace_file == 0 )
    return;
  else {
    // Disk full?
    if(access(dir.data(), W_OK) == -1) { cerr << "Cannot write ace-files to " << dir << endl; exit(-1); }
    
    ofstream outFile((dir + "/" + file_name + ".ace").data());
    
    cerr<<" (Printing to file.) ";
    
    // write AS
    int num_contigs = 1;
    outFile<<"AS "<<num_contigs<<" "<<tot_num_of_reads_in_ace_file<<endl;
    outFile<<endl;
    
    // write CO
    outFile<<"CO "<<contig_name<<" "<<num_of_bases
           <<" "<<num_of_reads_in_contig<<" "<<num_of_segments_in_contig
           <<" U"<<endl;
    for( int i = 0; i < num_of_bases; i++ ) {
      if ( i != 0 && i % 50 == 0 )
        outFile<<endl;
      outFile<<consensus[i];
    }
    outFile<<endl;
    outFile<<endl;
    
    // write BQ
    outFile<<"BQ "<<endl;
    for( int i = 0; i < num_of_unpadded_bases; i++ ) {
      if ( i != 0 && i % 50 == 0 )
        outFile<<endl;
      outFile<<' '<<consensus_quality_values[0];
    }
    outFile<<endl;
    outFile<<endl;
  
    // write AF
    for( int i = 0; i < num_of_reads_in_contig; i++ ) {
      // Change parameters value!
      outFile<<"AF "<<read_name[i]<<" "<<UorC[i]
             <<" "<<padded_start_consensus_position[i]<<endl;
    }
    
    // write BS
    for( int i = 0; i < num_of_segments_in_contig; i++ ) {
      // Change parameters value!
      outFile<<"BS "<<BS_padded_start_consensus_position[i]
             <<" "<<BS_padded_end_consensus_position[i]
             <<" "<<BS_read_name[i]<<endl;
    }
    outFile<<endl;
    
    // write RD, QA & DS
    for( int i = 0; i < num_of_reads_in_contig; i++ ) {
      // Change parameters value!
      
      // write RD
      outFile<<"RD "<<read_name[i]<<" "<<num_of_padded_bases[i]
             <<" "<<num_of_whole_read_info_items
             <<" "<<num_of_read_tags[i]<<endl;
      for( int j = 0; j < num_of_padded_bases[i]; j++ ) {
        if ( j != 0 && j % 50 == 0 )
          outFile<<endl;
        outFile<<read_seq[i][j];
      }
      outFile<<endl;
      outFile<<endl;
      
      // write QA
      outFile<<"QA "<<qual_clipping_start[i]
             <<" "<<qual_clipping_end[i]
             <<" "<<align_clipping_start[i]
             <<" "<<align_clipping_end[i]<<endl;
      outFile<<endl;
      
      // write DS
      outFile<<"DS CHROMAT_FILE: "<<name_of_chromat_file[i]
             <<" PHD_FILE: "<<name_of_phd_file[i]
             <<" TIME: "<<date_of_phd_file<<endl;
      outFile<<endl;
      
      // write RT
      for( int j = 0; j < (int)padded_read_pos_start[i].size(); j++ ) {
        outFile<<"RT{"<<endl;
        outFile<<read_name[i]<<' '<<tag_type[i][j]<<' '
               <<what_program_created_tag<<' '
               <<padded_read_pos_start[i][j]<<' '
               <<padded_read_pos_end[i][j]<<' '
               <<date_when_tag_was_created<<endl;
        outFile<<tag_data[i][j]<<endl;
        outFile<<'}'<<endl;
        outFile<<endl;
      }
    }

    cerr<<" Done. "<<endl;
    outFile.close();
  }
  
}


int AceWriter::get_consensus(vector <char> &cons, vector<int> &rowVec, int &num) {

  cons.clear();
  rowVec.clear();
  num = 0;
  
  int min_row_value = numeric_limits<int>::max();
  int num_cols = 0;
  
  for(int i=0; i < getMAl()->get_num_seq(); i++){
    if( num_cols < getMAl()->get_seq_end_global(i))
      num_cols = getMAl()->get_seq_end_global(i);
  }
  
  for( int c = 0; c < num_cols; c++ ) {   
    short Acount(0), Tcount(0), Gcount(0), Ccount(0), GAPcount(0);
    int Arow(1), Trow(1), Grow(1), Crow(1), GAProw(1);
    
    for( int r = 0; r < getMAl()->get_num_seq(); r++ ) {
      
      if ( min_row_value > getMAl()->get_seq_begin_global(r) )
        min_row_value = getMAl()->get_seq_begin_global(r);
      
      if( c >= getMAl()->get_seq_begin_global(r) && c < getMAl()->get_seq_end_global(r)  ){
        char base = getMAl()->get_base_global(r, c);
        
        if( base == 'a' )     { Acount++; if ( Acount == 1 ) Arow = r; }
        else if( base == 't' ){ Tcount++; if ( Tcount == 1 ) Trow = r; }
        else if( base == 'g' ){ Gcount++; if ( Gcount == 1 ) Grow = r; }
        else if( base == 'c' ){ Ccount++; if ( Ccount == 1 ) Crow = r; }
        else if( base == '*' || base == 'n'){ GAPcount++; if ( GAPcount == 1 ) GAProw = r; }
        
      }
    }
    if( Acount != 0 && Acount >= Tcount && Acount >= Gcount && Acount >= Ccount && Acount >= GAPcount )
      { cons.push_back( 'a' ); rowVec.push_back( Arow ); }//
    else if( Tcount != 0 && Tcount >= Gcount && Tcount >= Ccount && Tcount >= GAPcount )
      { cons.push_back( 't' ); rowVec.push_back( Trow ); }//
    else if( Gcount != 0 && Gcount >= Ccount && Gcount >= GAPcount )
      { cons.push_back( 'g' ); rowVec.push_back( Grow ); }//
    else if( Ccount != 0 && Ccount >= GAPcount )
      { cons.push_back( 'c' ); rowVec.push_back( Crow ); }//
    else if( GAPcount != 0 ) {
      num++;
      cons.push_back( '*' ); 
      rowVec.push_back( GAProw );
    }
  }
  assert( cons.size() == rowVec.size());
  num = cons.size() - num; // Number of unpadded bases!
  return min_row_value;
  
}

int AceWriter::get_realignment_row(const int row, vector<char> &aligned_read, int &row_start) {

  string tmp = getMAl()->get_seq(row);
  copy(tmp.begin(), tmp.end(), back_inserter(aligned_read));
  row_start = getMAl()->get_seq_begin_global(row);
  return getMAl()->get_seq_end_global(row);
}

---