gclib/tophat_cpp/reads.cpp

/*
 *  reads.cpp
 *  TopHat
 *
 *  Created by Cole Trapnell on 9/2/48.
 *  Copyright 2448 Cole Trapnell. All rights reserved.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <cassert>
#include <string>
#include <algorithm>
#include <cstring>
#include <cstdlib>

#include <seqan/find.h>
#include <seqan/file.h>
#include <seqan/modifier.h>

#include "reads.h"
#include "bwt_map.h"
#include "tokenize.h"

using namespace std;

char* FLineReader::nextLine() {
   if(!file) return NULL;
   if (pushed) { pushed=false; return buf; }
   //reads a char at a time until \n and/or \r are encountered
   len=0;
   int c=0;
   while ((c=getc(file))!=EOF) {
     if (len>=allocated-1) {
        allocated+=512;
        buf=(char*)realloc(buf,allocated);
     }
     if (c=='\n' || c=='\r') {
       buf[len]='\0';
       if (c=='\r') { //DOS file: double-char line terminator, skip the second one
          if ((c=getc(file))!='\n')
              ungetc(c,file);
          }
       lcount++;
       return buf;
       }
     buf[len]=(char)c;
     len++;
     }
   if (c==EOF) {
     isEOF=true;
     if (len==0) return NULL;
     }
   buf[len]='\0';
   lcount++;
   return buf;
}

void skip_lines(FLineReader& fr)
{
  if (fr.fhandle() == NULL) return;
  char* buf = NULL;
  while ((buf = fr.nextLine()) != NULL) {
    if (buf[0] == '\0') continue;
    if (buf[0] == '>' || buf[0] == '@')
      {
        fr.pushBack();
        break;
      }
  }
}

bool next_fasta_record(FLineReader& fr,
                       string& defline, 
                       string& seq,
                       ReadFormat reads_format)

{
  seq.clear();
  defline.clear();
  char* buf=NULL;
  while ((buf=fr.nextLine())!=NULL) {
    if (buf[0]==0) continue; //skip empty lines
    if ((reads_format == FASTA && buf[0] == '>') || (reads_format == FASTQ && (buf[0] == '+' || buf[0] == '@'))) { //next record
        if (seq.length()>0) { //current record ending
           fr.pushBack();
           return true;
           }
        defline=buf+1;
        string::size_type space_pos = defline.find_first_of(" \t");
        if (space_pos != string::npos) {
            defline.resize(space_pos);
            }
        continue;
        } //defline
    // sequence line
    seq.append(buf);
    } //line reading loop

  replace(seq.begin(), seq.end(), '.', color ? '4' : 'N'); //shouldn't really be needed for FASTA files
  return !(seq.empty());
}

bool next_fastq_record(FLineReader& fr,
                       const string& seq,
                       string& alt_name,
                       string& qual,
                       ReadFormat reads_format)
{
  alt_name.clear();
  qual.clear();
  char* fline=fr.nextLine();
  if (fline==NULL) return false;
  while (fline[0]==0) { //skip empty lines
    fline=fr.nextLine();
    if (fline==NULL) return false;
    }
  //must be on '+' line here
  if (fline==NULL || (reads_format == FASTQ && fline[0] != '+') ||
      (reads_format == FASTA && quals && fline[0] != '>')) {
     err_exit("Error: '+' not found for fastq record %s\n",fline);
     return false;
     }
  alt_name=fline+1;
  string::size_type space_pos = alt_name.find_first_of(" \t");
  if (space_pos != string::npos) alt_name.resize(space_pos);
   //read qv line(s) now:
  while ((fline=fr.nextLine())!=NULL) {
    if (integer_quals)
      {
      vector<string> integer_qual_values;
      tokenize(string(fline), " ", integer_qual_values);

      string temp_qual;
      for (size_t i = 0; i < integer_qual_values.size(); ++i)
        {
          int qual_value = atoi(integer_qual_values[i].c_str());
          if (qual_value < 0) qual_value = 0;
          temp_qual.push_back((char)(qual_value + 33));
        }

      qual.append(temp_qual);
      }
    else
      qual.append(fline);
      if (qual.length()>=seq.length()-1) break;
     }
  // final check
  if (color) {
     if (seq.length()==qual.length()) {
        //discard first qv
        qual=qual.substr(1);
        }
     if (seq.length()!=qual.length()+1) {
        err_exit("Error: length of quality string does not match seq length (%d) for color read %s!\n",
           seq.length(), alt_name.c_str());
        }
     }
  else {
    if (seq.length()!=qual.length()) {
           err_exit("Error: qual string length (%d) differs from seq length (%d) for read %s!\n",
               qual.length(), seq.length(), alt_name.c_str());
           //return false;
           }
    }
  //
  return !(qual.empty());
}

bool next_fastx_read(FLineReader& fr, Read& read, ReadFormat reads_format, 
                        FLineReader* frq) {
  read.clear();
  char* buf=NULL;
  while ((buf=fr.nextLine())!=NULL) {
    if (buf[0]==0) continue; //skip empty lines
    if ((reads_format == FASTA && buf[0] == '>') ||
          (reads_format == FASTQ && (buf[0] == '+' || buf[0] == '@'))) { //next record
        if (read.seq.length()>0) { //current record ending
           fr.pushBack();
           break;
           }
        read.name=buf+1;
        string::size_type space_pos = read.name.find_first_of(" \t");
        if (space_pos != string::npos) {
            read.name.resize(space_pos);
            }
        continue;
        } //defline
    // sequence line
    read.seq.append(buf);
    } //line reading loop

  replace(read.seq.begin(), read.seq.end(), '.', color ? '4' : 'N'); //shouldn't really be needed for FASTA files
  if (reads_format != FASTQ && frq==NULL)
      return (!read.seq.empty());
  if (frq==NULL) frq=&fr; //FASTQ
  //FASTQ or quals in a separate file -- now read quality values
  buf=frq->nextLine();
  if (buf==NULL) return false;
  while (buf[0]==0) { //skip empty lines
    buf=frq->nextLine();
    if (buf==NULL) return false;
    }
  //must be on '+' line here
  if (buf==NULL || (reads_format == FASTQ && buf[0] != '+') ||
           (reads_format == FASTA && buf[0] != '>')) {
     err_exit("Error: beginning of quality values record not found! (%s)\n",buf);
     return false;
     }
  read.alt_name=buf+1;
  string::size_type space_pos = read.alt_name.find_first_of(" \t");
  if (space_pos != string::npos) read.alt_name.resize(space_pos);
   //read qv line(s) now:
  while ((buf=frq->nextLine())!=NULL) {
    if (integer_quals)
      {
        vector<string> integer_qual_values;
        tokenize(string(buf), " ", integer_qual_values);
        string temp_qual;
        for (size_t i = 0; i < integer_qual_values.size(); ++i)
          {
            int qual_value = atoi(integer_qual_values[i].c_str());
            if (qual_value < 0) qual_value = 0;
            temp_qual.push_back((char)(qual_value + 33));
          }
        read.qual.append(temp_qual);
      }
    else {
        read.qual.append(buf);
      }
    if (read.qual.length()>=read.seq.length()-1)
          break;
    } //while qv lines
  
  // final check
  if (color) {
     if (read.seq.length()==read.qual.length()) {
        //discard first qv
        read.qual=read.qual.substr(1);
        }
     if (read.seq.length()!=read.qual.length()+1) {
        err_exit("Error: length of quality string does not match sequence length (%d) for color read %s!\n",
            read.seq.length(), read.alt_name.c_str());
        }
     }
  else {
   if (read.seq.length()!=read.qual.length()) {
           err_exit("Error: qual length (%d) differs from seq length (%d) for fastq record %s!\n",
               read.qual.length(), read.seq.length(), read.alt_name.c_str());
           return false;
           }
    }
  return !(read.seq.empty());
}

// This could be faster.
void reverse_complement(string& seq)
{
        //fprintf(stderr,"fwd: %s\n", seq.c_str());
        for (string::size_type i = 0; i < seq.length(); ++i)
        {
                switch(seq[i])
                {
                        case 'A' : seq[i] = 'T'; break;
                        case 'T' : seq[i] = 'A'; break;
                        case 'C' : seq[i] = 'G'; break;
                        case 'G' : seq[i] = 'C'; break;
                        default: seq[i]   = 'N'; break;
                }
        }
        reverse(seq.begin(), seq.end());
        //fprintf(stderr, "rev: %s\n", seq.c_str());
}

string convert_color_to_bp(const string& color)
{
  if (color.length() <= 0)
    return "";

  char base = color[0];
  string bp;
  for (string::size_type i = 1; i < color.length(); ++i)
    {
      char next = color[i];
      switch(base)
        {
          // 'A0':'A', 'A1':'C', 'A2':'G', 'A3':'T', 'A4':'N', 'A.':'N',
        case 'A':
          {
            switch(next)
              {
              case '0': next = 'A'; break;
              case '1': next = 'C'; break;
              case '2': next = 'G'; break;
              case '3': next = 'T'; break;
              default: next = 'N'; break;
              }
          }
          break;
        case 'C':
          {
            // 'C0':'C', 'C1':'A', 'C2':'T', 'C3':'G', 'C4':'N', 'C.':'N',
             switch(next)
              {
              case '0': next = 'C'; break;
              case '1': next = 'A'; break;
              case '2': next = 'T'; break;
              case '3': next = 'G'; break;
              default: next = 'N'; break;
              }
          }
          break;
        case 'G':
          {
            // 'G0':'G', 'G1':'T', 'G2':'A', 'G3':'C', 'G4':'N', 'G.':'N',
             switch(next)
              {
              case '0': next = 'G'; break;
              case '1': next = 'T'; break;
              case '2': next = 'A'; break;
              case '3': next = 'C'; break;
              default: next = 'N'; break;
              }
          }
          break;
        case 'T':
          {
            // 'T0':'T', 'T1':'G', 'T2':'C', 'T3':'A', 'T4':'N', 'T.':'N',
             switch(next)
              {
              case '0': next = 'T'; break;
              case '1': next = 'G'; break;
              case '2': next = 'C'; break;
              case '3': next = 'A'; break;
              default: next = 'N'; break;
              }
          }
          break;
        default: next = 'N'; break;
        }

      bp.push_back(next);
      base = next;
    }

  return bp;
}

// daehwan - reduce code redundancy!
seqan::String<char> convert_color_to_bp(char base, const seqan::String<char>& color)
{
  if (seqan::length(color) <= 0)
    return "";

  string bp;
  for (string::size_type i = 0; i < seqan::length(color); ++i)
    {
      char next = color[i];
      switch(base)
        {
          // 'A0':'A', 'A1':'C', 'A2':'G', 'A3':'T', 'A4':'N', 'A.':'N',
        case 'A':
          {
            switch(next)
              {
              case '0': next = 'A'; break;
              case '1': next = 'C'; break;
              case '2': next = 'G'; break;
              case '3': next = 'T'; break;
              default: next = 'N'; break;
              }
          }
          break;
        case 'C':
          {
            // 'C0':'C', 'C1':'A', 'C2':'T', 'C3':'G', 'C4':'N', 'C.':'N',
             switch(next)
              {
              case '0': next = 'C'; break;
              case '1': next = 'A'; break;
              case '2': next = 'T'; break;
              case '3': next = 'G'; break;
              default: next = 'N'; break;
              }
          }
          break;
        case 'G':
          {
            // 'G0':'G', 'G1':'T', 'G2':'A', 'G3':'C', 'G4':'N', 'G.':'N',
             switch(next)
              {
              case '0': next = 'G'; break;
              case '1': next = 'T'; break;
              case '2': next = 'A'; break;
              case '3': next = 'C'; break;
              default: next = 'N'; break;
              }
          }
          break;
        case 'T':
          {
            // 'T0':'T', 'T1':'G', 'T2':'C', 'T3':'A', 'T4':'N', 'T.':'N',
             switch(next)
              {
              case '0': next = 'T'; break;
              case '1': next = 'G'; break;
              case '2': next = 'C'; break;
              case '3': next = 'A'; break;
              default: next = 'N'; break;
              }
          }
          break;
        default: next = 'N'; break;
        }

      bp.push_back(next);
      base = next;
    }

  return bp;
}


#define check_color(b1, b2, c1, c2) ((b1 == c1 && b2 == c2) || (b1 == c2 && b2 == c1))

#define two_bps_to_color(b1, b2, c)                                   \
  if (((b1) == 'A' || (b1) == 'G' || (b1) == 'C' || (b1) == 'T') && (b1) == (b2)) \
  c = '0'; \
  else if (check_color((b1), (b2), 'A', 'C') || check_color((b1), (b2), 'G', 'T')) \
  c = '1'; \
  else if (check_color((b1), (b2), 'A', 'G') || check_color((b1), (b2), 'C', 'T')) \
  c = '2'; \
  else if (check_color((b1), (b2), 'A', 'T') || check_color((b1), (b2), 'C', 'G')) \
  c = '3'; \
  else \
  c = '4';
 

string convert_bp_to_color(const string& bp, bool remove_primer)
{
  if (bp.length() <= 1)
    return "";

  char base = toupper(bp[0]);
  string color;
  if (!remove_primer)
    color.push_back(base);
  
  for (string::size_type i = 1; i < bp.length(); ++i)
    {
      char next = toupper(bp[i]);

      char c = '0';
      two_bps_to_color(base, next, c);
      color.push_back(c);

      base = next;
    }

  return color;
}

// daehwan - check this - 
seqan::String<char> convert_bp_to_color(const seqan::String<char>& bp, bool remove_primer)
{
  if (seqan::length(bp) <= 1)
    return "";

  char base = toupper(bp[0]);
  string color;
  if (!remove_primer)
    color.push_back(base);
  
  for (string::size_type i = 1; i < seqan::length(bp); ++i)
    {
      char next = toupper(bp[i]);

      char c = '0';
      two_bps_to_color(base, next, c);
      color.push_back(c);

      base = next;
    }

  return color;
}

/*
 */
void BWA_decode(const string& color, const string& qual, const string& ref, string& decode)
{
  assert(color.length() == ref.length() - 1);
  
  const size_t max_length = 256;
  const unsigned int max_value = max_length * 0xff;
  size_t length = color.length();
  if (length < 1 || length + 1 > max_length)
    {
      return;
    }

  unsigned int f[max_length * 4];
  char ptr[max_length * 4];

  unsigned int q_prev = 0;
  for (unsigned int i = 0; i < length + 1; ++i)
    {
      unsigned int q = (unsigned int) (qual.length() <= i ? 'I' : qual[i]) - 33;
      for (unsigned int j = 0; j < 4; ++j)
        {
          size_t i_j = i * 4 + j;
          if (i == 0)
            {
              f[i_j] = "ACGT"[j] == ref[i] ? 0 : q;
              ptr[i_j] = 4;
              continue;
            }

          f[i_j] = max_value;
          char base = "ACGT"[j];
          for (unsigned int k = 0; k < 4; ++k)
            {
              char base_prev = "ACGT"[k];
              char ref_color;
              two_bps_to_color(base_prev, base, ref_color);

              char base_prev_prev = "ACGTN"[ptr[(i-1)*4 + k]];
              char ref_color_prev;
              two_bps_to_color(base_prev_prev, base_prev, ref_color_prev);

              char color_curr = color[i-1];
              char color_prev = i >= 2 ? color[i-2] : '4';            

              int q_hat = 0;
              if (color_prev == ref_color_prev && color_prev != '4')
                {
                  if (color_curr == ref_color)
                    q_hat = q + q_prev;
                  else
                    q_hat = q_prev - q;
                }
              else if (color_curr == ref_color)
                {
                  q_hat = q - q_prev;
                }

              unsigned int f_k = f[(i-1) * 4 + k] +
                (base == ref[i] ? 0 : q_hat) +
                (color_curr == ref_color ? 0 : q);

              if (f_k < f[i_j])
                {
                  f[i_j] = f_k;
                  ptr[i_j] = k;
                }
            }
        }

      q_prev = q;
    }

  unsigned int min_index = 0;
  unsigned int min_f = f[length * 4];
  for (unsigned int i = 1; i < 4; ++i)
    {
      unsigned int temp_f = f[length * 4 + i];
      if (temp_f < min_f)
        {
          min_f = temp_f;
          min_index = i;
        }
    }

  decode.resize(length + 1);
  decode[length] = "ACGT"[min_index];
  for (unsigned int i = length; i > 0; --i)
    {
      min_index = ptr[i * 4 + min_index];
      decode[i-1] = "ACGT"[min_index];
    }
}


bool ReadStream::next_read(Read& r, ReadFormat read_format) {
  FLineReader fr(fstream.file);
  while (read_pq.size()<100000 && !r_eof) {
    //keep the queue topped off
    Read rf;
    if (!next_fastx_read(fr, rf, read_format)) {
        r_eof=true;
        break;
        }
    //Read read=Read(rf);
    uint64_t id = (uint64_t)atol(rf.name.c_str());
    read_pq.push(make_pair(id, rf));
    }
  if (read_pq.size()==0)
     return false;
  const pair<uint64_t, Read>& t = read_pq.top();
  r=t.second; //copy strings
  //free(t.second);
  read_pq.pop();
  return true;
}

// reads must ALWAYS requested in increasing order of their ID
bool ReadStream::getRead(uint64_t r_id,
            Read& read,
            ReadFormat read_format,
            bool strip_slash,
            FILE* um_out, //unmapped reads output
            bool um_write_found) {
  if (!fstream.file)
       err_die("Error: calling ReadStream::getRead() with no file handle!");
  if (r_id<last_id)
      err_die("Error: ReadStream::getRead() called with out-of-order id#!");
  last_id=r_id;
  bool found=false;
  while (!found) {
    read.clear();
      // Get the next read from the file
    if (!next_read(read, read_format))
        break;
    if (strip_slash) {
       string::size_type slash = read.name.rfind("/");
       if (slash != string::npos)
          read.name.resize(slash);
       }
    if ((uint64_t)atoi(read.name.c_str()) == r_id) {
       found=true;
       }
    if (um_out && (um_write_found || !found)) {
     //write unmapped reads
      fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(),
                              read.seq.c_str(), read.qual.c_str());
      }
    //rt.get_id(read.name, ref_str);
    } //while reads
  return found;
}


bool get_read_from_stream(uint64_t insert_id,
        FILE* reads_file,
        ReadFormat reads_format,
        bool strip_slash,
        Read& read,
        FILE* um_out,
        bool um_write_found
        ) {
  FLineReader fr(reads_file);
  bool found=false;
  while(!found && !fr.isEof())
    {
    read.clear();
      // Get the next read from the file
    if (!next_fastx_read(fr, read, reads_format))
        break;
    if (strip_slash)
      {
        string::size_type slash = read.name.rfind("/");
        if (slash != string::npos)
          read.name.resize(slash);
      }
    if ((uint64_t)atoi(read.name.c_str()) == insert_id)
      {
      //return true;
      found=true;
      }
    if (um_out && (um_write_found || !found)) {
     //write unmapped reads
      fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(),
                              read.seq.c_str(), read.qual.c_str());
      }
    //rt.get_id(read.name, ref_str);
    } //while reads
  return found;
}


bool get_read_from_stream(uint64_t insert_id,
                          FILE* reads_file,
                          ReadFormat reads_format,
                          bool strip_slash,
                          char read_name [], 
                          char read_seq  [],
                          char read_alt_name [], 
                          char read_qual [],
                          FILE* um_out)
{
  Read read;
  FLineReader fr(reads_file);
  while(!fr.isEof())
    {
    read.clear();

      // Get the next read from the file
    if (!next_fastx_read(fr, read, reads_format))
        break;

    if (strip_slash)
      {
        string::size_type slash = read.name.rfind("/");
        if (slash != string::npos)
          read.name.resize(slash);
      }

    if ((uint64_t)atoi(read.name.c_str()) == insert_id)
      {
        if (read_name) strcpy(read_name, read.name.c_str());
        if (read_seq) strcpy(read_seq, read.seq.c_str());
        if (read_alt_name) strcpy(read_alt_name, read.alt_name.c_str());
        if (read_qual) strcpy(read_qual, read.qual.c_str());
        return true;
      }
    else if (um_out!=NULL) {
     //write unmapped reads
      fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(), read.seq.c_str(),
                          read.qual.c_str());
      }
      
        //rt.get_id(read.name, ref_str);
    } //while reads
  
  return false;
}
Revision:	135
Committed:	Mon Dec 12 22:28:38 2011 UTC (12 years, 8 months ago) by gpertea
File size:	18715 byte(s)
Log Message:	wip - SplicedSAMHitFactory() still not implemented
Line	File contents
1	/*
2	* reads.cpp
3	* TopHat
4	*
5	* Created by Cole Trapnell on 9/2/48.
6	* Copyright 2448 Cole Trapnell. All rights reserved.
7	*
8	*/
9
10	#ifdef HAVE_CONFIG_H
11	#include <config.h>
12	#endif
13
14	#include <iostream>
15	#include <cassert>
16	#include <string>
17	#include <algorithm>
18	#include <cstring>
19	#include <cstdlib>
20
21	#include <seqan/find.h>
22	#include <seqan/file.h>
23	#include <seqan/modifier.h>
24
25	#include "reads.h"
26	#include "bwt_map.h"
27	#include "tokenize.h"
28
29	using namespace std;
30
31	char* FLineReader::nextLine() {
32	if(!file) return NULL;
33	if (pushed) { pushed=false; return buf; }
34	//reads a char at a time until \n and/or \r are encountered
35	len=0;
36	int c=0;
37	while ((c=getc(file))!=EOF) {
38	if (len>=allocated-1) {
39	allocated+=512;
40	buf=(char*)realloc(buf,allocated);
41	}
42	if (c=='\n' \|\| c=='\r') {
43	buf[len]='\0';
44	if (c=='\r') { //DOS file: double-char line terminator, skip the second one
45	if ((c=getc(file))!='\n')
46	ungetc(c,file);
47	}
48	lcount++;
49	return buf;
50	}
51	buf[len]=(char)c;
52	len++;
53	}
54	if (c==EOF) {
55	isEOF=true;
56	if (len==0) return NULL;
57	}
58	buf[len]='\0';
59	lcount++;
60	return buf;
61	}
62
63	void skip_lines(FLineReader& fr)
64	{
65	if (fr.fhandle() == NULL) return;
66	char* buf = NULL;
67	while ((buf = fr.nextLine()) != NULL) {
68	if (buf[0] == '\0') continue;
69	if (buf[0] == '>' \|\| buf[0] == '@')
70	{
71	fr.pushBack();
72	break;
73	}
74	}
75	}
76
77	bool next_fasta_record(FLineReader& fr,
78	string& defline,
79	string& seq,
80	ReadFormat reads_format)
81
82	{
83	seq.clear();
84	defline.clear();
85	char* buf=NULL;
86	while ((buf=fr.nextLine())!=NULL) {
87	if (buf[0]==0) continue; //skip empty lines
88	if ((reads_format == FASTA && buf[0] == '>') \|\| (reads_format == FASTQ && (buf[0] == '+' \|\| buf[0] == '@'))) { //next record
89	if (seq.length()>0) { //current record ending
90	fr.pushBack();
91	return true;
92	}
93	defline=buf+1;
94	string::size_type space_pos = defline.find_first_of(" \t");
95	if (space_pos != string::npos) {
96	defline.resize(space_pos);
97	}
98	continue;
99	} //defline
100	// sequence line
101	seq.append(buf);
102	} //line reading loop
103
104	replace(seq.begin(), seq.end(), '.', color ? '4' : 'N'); //shouldn't really be needed for FASTA files
105	return !(seq.empty());
106	}
107
108	bool next_fastq_record(FLineReader& fr,
109	const string& seq,
110	string& alt_name,
111	string& qual,
112	ReadFormat reads_format)
113	{
114	alt_name.clear();
115	qual.clear();
116	char* fline=fr.nextLine();
117	if (fline==NULL) return false;
118	while (fline[0]==0) { //skip empty lines
119	fline=fr.nextLine();
120	if (fline==NULL) return false;
121	}
122	//must be on '+' line here
123	if (fline==NULL \|\| (reads_format == FASTQ && fline[0] != '+') \|\|
124	(reads_format == FASTA && quals && fline[0] != '>')) {
125	err_exit("Error: '+' not found for fastq record %s\n",fline);
126	return false;
127	}
128	alt_name=fline+1;
129	string::size_type space_pos = alt_name.find_first_of(" \t");
130	if (space_pos != string::npos) alt_name.resize(space_pos);
131	//read qv line(s) now:
132	while ((fline=fr.nextLine())!=NULL) {
133	if (integer_quals)
134	{
135	vector<string> integer_qual_values;
136	tokenize(string(fline), " ", integer_qual_values);
137
138	string temp_qual;
139	for (size_t i = 0; i < integer_qual_values.size(); ++i)
140	{
141	int qual_value = atoi(integer_qual_values[i].c_str());
142	if (qual_value < 0) qual_value = 0;
143	temp_qual.push_back((char)(qual_value + 33));
144	}
145
146	qual.append(temp_qual);
147	}
148	else
149	qual.append(fline);
150	if (qual.length()>=seq.length()-1) break;
151	}
152	// final check
153	if (color) {
154	if (seq.length()==qual.length()) {
155	//discard first qv
156	qual=qual.substr(1);
157	}
158	if (seq.length()!=qual.length()+1) {
159	err_exit("Error: length of quality string does not match seq length (%d) for color read %s!\n",
160	seq.length(), alt_name.c_str());
161	}
162	}
163	else {
164	if (seq.length()!=qual.length()) {
165	err_exit("Error: qual string length (%d) differs from seq length (%d) for read %s!\n",
166	qual.length(), seq.length(), alt_name.c_str());
167	//return false;
168	}
169	}
170	//
171	return !(qual.empty());
172	}
173
174	bool next_fastx_read(FLineReader& fr, Read& read, ReadFormat reads_format,
175	FLineReader* frq) {
176	read.clear();
177	char* buf=NULL;
178	while ((buf=fr.nextLine())!=NULL) {
179	if (buf[0]==0) continue; //skip empty lines
180	if ((reads_format == FASTA && buf[0] == '>') \|\|
181	(reads_format == FASTQ && (buf[0] == '+' \|\| buf[0] == '@'))) { //next record
182	if (read.seq.length()>0) { //current record ending
183	fr.pushBack();
184	break;
185	}
186	read.name=buf+1;
187	string::size_type space_pos = read.name.find_first_of(" \t");
188	if (space_pos != string::npos) {
189	read.name.resize(space_pos);
190	}
191	continue;
192	} //defline
193	// sequence line
194	read.seq.append(buf);
195	} //line reading loop
196
197	replace(read.seq.begin(), read.seq.end(), '.', color ? '4' : 'N'); //shouldn't really be needed for FASTA files
198	if (reads_format != FASTQ && frq==NULL)
199	return (!read.seq.empty());
200	if (frq==NULL) frq=&fr; //FASTQ
201	//FASTQ or quals in a separate file -- now read quality values
202	buf=frq->nextLine();
203	if (buf==NULL) return false;
204	while (buf[0]==0) { //skip empty lines
205	buf=frq->nextLine();
206	if (buf==NULL) return false;
207	}
208	//must be on '+' line here
209	if (buf==NULL \|\| (reads_format == FASTQ && buf[0] != '+') \|\|
210	(reads_format == FASTA && buf[0] != '>')) {
211	err_exit("Error: beginning of quality values record not found! (%s)\n",buf);
212	return false;
213	}
214	read.alt_name=buf+1;
215	string::size_type space_pos = read.alt_name.find_first_of(" \t");
216	if (space_pos != string::npos) read.alt_name.resize(space_pos);
217	//read qv line(s) now:
218	while ((buf=frq->nextLine())!=NULL) {
219	if (integer_quals)
220	{
221	vector<string> integer_qual_values;
222	tokenize(string(buf), " ", integer_qual_values);
223	string temp_qual;
224	for (size_t i = 0; i < integer_qual_values.size(); ++i)
225	{
226	int qual_value = atoi(integer_qual_values[i].c_str());
227	if (qual_value < 0) qual_value = 0;
228	temp_qual.push_back((char)(qual_value + 33));
229	}
230	read.qual.append(temp_qual);
231	}
232	else {
233	read.qual.append(buf);
234	}
235	if (read.qual.length()>=read.seq.length()-1)
236	break;
237	} //while qv lines
238
239	// final check
240	if (color) {
241	if (read.seq.length()==read.qual.length()) {
242	//discard first qv
243	read.qual=read.qual.substr(1);
244	}
245	if (read.seq.length()!=read.qual.length()+1) {
246	err_exit("Error: length of quality string does not match sequence length (%d) for color read %s!\n",
247	read.seq.length(), read.alt_name.c_str());
248	}
249	}
250	else {
251	if (read.seq.length()!=read.qual.length()) {
252	err_exit("Error: qual length (%d) differs from seq length (%d) for fastq record %s!\n",
253	read.qual.length(), read.seq.length(), read.alt_name.c_str());
254	return false;
255	}
256	}
257	return !(read.seq.empty());
258	}
259
260	// This could be faster.
261	void reverse_complement(string& seq)
262	{
263	//fprintf(stderr,"fwd: %s\n", seq.c_str());
264	for (string::size_type i = 0; i < seq.length(); ++i)
265	{
266	switch(seq[i])
267	{
268	case 'A' : seq[i] = 'T'; break;
269	case 'T' : seq[i] = 'A'; break;
270	case 'C' : seq[i] = 'G'; break;
271	case 'G' : seq[i] = 'C'; break;
272	default: seq[i] = 'N'; break;
273	}
274	}
275	reverse(seq.begin(), seq.end());
276	//fprintf(stderr, "rev: %s\n", seq.c_str());
277	}
278
279	string convert_color_to_bp(const string& color)
280	{
281	if (color.length() <= 0)
282	return "";
283
284	char base = color[0];
285	string bp;
286	for (string::size_type i = 1; i < color.length(); ++i)
287	{
288	char next = color[i];
289	switch(base)
290	{
291	// 'A0':'A', 'A1':'C', 'A2':'G', 'A3':'T', 'A4':'N', 'A.':'N',
292	case 'A':
293	{
294	switch(next)
295	{
296	case '0': next = 'A'; break;
297	case '1': next = 'C'; break;
298	case '2': next = 'G'; break;
299	case '3': next = 'T'; break;
300	default: next = 'N'; break;
301	}
302	}
303	break;
304	case 'C':
305	{
306	// 'C0':'C', 'C1':'A', 'C2':'T', 'C3':'G', 'C4':'N', 'C.':'N',
307	switch(next)
308	{
309	case '0': next = 'C'; break;
310	case '1': next = 'A'; break;
311	case '2': next = 'T'; break;
312	case '3': next = 'G'; break;
313	default: next = 'N'; break;
314	}
315	}
316	break;
317	case 'G':
318	{
319	// 'G0':'G', 'G1':'T', 'G2':'A', 'G3':'C', 'G4':'N', 'G.':'N',
320	switch(next)
321	{
322	case '0': next = 'G'; break;
323	case '1': next = 'T'; break;
324	case '2': next = 'A'; break;
325	case '3': next = 'C'; break;
326	default: next = 'N'; break;
327	}
328	}
329	break;
330	case 'T':
331	{
332	// 'T0':'T', 'T1':'G', 'T2':'C', 'T3':'A', 'T4':'N', 'T.':'N',
333	switch(next)
334	{
335	case '0': next = 'T'; break;
336	case '1': next = 'G'; break;
337	case '2': next = 'C'; break;
338	case '3': next = 'A'; break;
339	default: next = 'N'; break;
340	}
341	}
342	break;
343	default: next = 'N'; break;
344	}
345
346	bp.push_back(next);
347	base = next;
348	}
349
350	return bp;
351	}
352
353	// daehwan - reduce code redundancy!
354	seqan::String<char> convert_color_to_bp(char base, const seqan::String<char>& color)
355	{
356	if (seqan::length(color) <= 0)
357	return "";
358
359	string bp;
360	for (string::size_type i = 0; i < seqan::length(color); ++i)
361	{
362	char next = color[i];
363	switch(base)
364	{
365	// 'A0':'A', 'A1':'C', 'A2':'G', 'A3':'T', 'A4':'N', 'A.':'N',
366	case 'A':
367	{
368	switch(next)
369	{
370	case '0': next = 'A'; break;
371	case '1': next = 'C'; break;
372	case '2': next = 'G'; break;
373	case '3': next = 'T'; break;
374	default: next = 'N'; break;
375	}
376	}
377	break;
378	case 'C':
379	{
380	// 'C0':'C', 'C1':'A', 'C2':'T', 'C3':'G', 'C4':'N', 'C.':'N',
381	switch(next)
382	{
383	case '0': next = 'C'; break;
384	case '1': next = 'A'; break;
385	case '2': next = 'T'; break;
386	case '3': next = 'G'; break;
387	default: next = 'N'; break;
388	}
389	}
390	break;
391	case 'G':
392	{
393	// 'G0':'G', 'G1':'T', 'G2':'A', 'G3':'C', 'G4':'N', 'G.':'N',
394	switch(next)
395	{
396	case '0': next = 'G'; break;
397	case '1': next = 'T'; break;
398	case '2': next = 'A'; break;
399	case '3': next = 'C'; break;
400	default: next = 'N'; break;
401	}
402	}
403	break;
404	case 'T':
405	{
406	// 'T0':'T', 'T1':'G', 'T2':'C', 'T3':'A', 'T4':'N', 'T.':'N',
407	switch(next)
408	{
409	case '0': next = 'T'; break;
410	case '1': next = 'G'; break;
411	case '2': next = 'C'; break;
412	case '3': next = 'A'; break;
413	default: next = 'N'; break;
414	}
415	}
416	break;
417	default: next = 'N'; break;
418	}
419
420	bp.push_back(next);
421	base = next;
422	}
423
424	return bp;
425	}
426
427
428	#define check_color(b1, b2, c1, c2) ((b1 == c1 && b2 == c2) \|\| (b1 == c2 && b2 == c1))
429
430	#define two_bps_to_color(b1, b2, c) \
431	if (((b1) == 'A' \|\| (b1) == 'G' \|\| (b1) == 'C' \|\| (b1) == 'T') && (b1) == (b2)) \
432	c = '0'; \
433	else if (check_color((b1), (b2), 'A', 'C') \|\| check_color((b1), (b2), 'G', 'T')) \
434	c = '1'; \
435	else if (check_color((b1), (b2), 'A', 'G') \|\| check_color((b1), (b2), 'C', 'T')) \
436	c = '2'; \
437	else if (check_color((b1), (b2), 'A', 'T') \|\| check_color((b1), (b2), 'C', 'G')) \
438	c = '3'; \
439	else \
440	c = '4';
441
442
443	string convert_bp_to_color(const string& bp, bool remove_primer)
444	{
445	if (bp.length() <= 1)
446	return "";
447
448	char base = toupper(bp[0]);
449	string color;
450	if (!remove_primer)
451	color.push_back(base);
452
453	for (string::size_type i = 1; i < bp.length(); ++i)
454	{
455	char next = toupper(bp[i]);
456
457	char c = '0';
458	two_bps_to_color(base, next, c);
459	color.push_back(c);
460
461	base = next;
462	}
463
464	return color;
465	}
466
467	// daehwan - check this -
468	seqan::String<char> convert_bp_to_color(const seqan::String<char>& bp, bool remove_primer)
469	{
470	if (seqan::length(bp) <= 1)
471	return "";
472
473	char base = toupper(bp[0]);
474	string color;
475	if (!remove_primer)
476	color.push_back(base);
477
478	for (string::size_type i = 1; i < seqan::length(bp); ++i)
479	{
480	char next = toupper(bp[i]);
481
482	char c = '0';
483	two_bps_to_color(base, next, c);
484	color.push_back(c);
485
486	base = next;
487	}
488
489	return color;
490	}
491
492	/*
493	*/
494	void BWA_decode(const string& color, const string& qual, const string& ref, string& decode)
495	{
496	assert(color.length() == ref.length() - 1);
497
498	const size_t max_length = 256;
499	const unsigned int max_value = max_length * 0xff;
500	size_t length = color.length();
501	if (length < 1 \|\| length + 1 > max_length)
502	{
503	return;
504	}
505
506	unsigned int f[max_length * 4];
507	char ptr[max_length * 4];
508
509	unsigned int q_prev = 0;
510	for (unsigned int i = 0; i < length + 1; ++i)
511	{
512	unsigned int q = (unsigned int) (qual.length() <= i ? 'I' : qual[i]) - 33;
513	for (unsigned int j = 0; j < 4; ++j)
514	{
515	size_t i_j = i * 4 + j;
516	if (i == 0)
517	{
518	f[i_j] = "ACGT"[j] == ref[i] ? 0 : q;
519	ptr[i_j] = 4;
520	continue;
521	}
522
523	f[i_j] = max_value;
524	char base = "ACGT"[j];
525	for (unsigned int k = 0; k < 4; ++k)
526	{
527	char base_prev = "ACGT"[k];
528	char ref_color;
529	two_bps_to_color(base_prev, base, ref_color);
530
531	char base_prev_prev = "ACGTN"[ptr[(i-1)*4 + k]];
532	char ref_color_prev;
533	two_bps_to_color(base_prev_prev, base_prev, ref_color_prev);
534
535	char color_curr = color[i-1];
536	char color_prev = i >= 2 ? color[i-2] : '4';
537
538	int q_hat = 0;
539	if (color_prev == ref_color_prev && color_prev != '4')
540	{
541	if (color_curr == ref_color)
542	q_hat = q + q_prev;
543	else
544	q_hat = q_prev - q;
545	}
546	else if (color_curr == ref_color)
547	{
548	q_hat = q - q_prev;
549	}
550
551	unsigned int f_k = f[(i-1) * 4 + k] +
552	(base == ref[i] ? 0 : q_hat) +
553	(color_curr == ref_color ? 0 : q);
554
555	if (f_k < f[i_j])
556	{
557	f[i_j] = f_k;
558	ptr[i_j] = k;
559	}
560	}
561	}
562
563	q_prev = q;
564	}
565
566	unsigned int min_index = 0;
567	unsigned int min_f = f[length * 4];
568	for (unsigned int i = 1; i < 4; ++i)
569	{
570	unsigned int temp_f = f[length * 4 + i];
571	if (temp_f < min_f)
572	{
573	min_f = temp_f;
574	min_index = i;
575	}
576	}
577
578	decode.resize(length + 1);
579	decode[length] = "ACGT"[min_index];
580	for (unsigned int i = length; i > 0; --i)
581	{
582	min_index = ptr[i * 4 + min_index];
583	decode[i-1] = "ACGT"[min_index];
584	}
585	}
586
587
588	bool ReadStream::next_read(Read& r, ReadFormat read_format) {
589	FLineReader fr(fstream.file);
590	while (read_pq.size()<100000 && !r_eof) {
591	//keep the queue topped off
592	Read rf;
593	if (!next_fastx_read(fr, rf, read_format)) {
594	r_eof=true;
595	break;
596	}
597	//Read read=Read(rf);
598	uint64_t id = (uint64_t)atol(rf.name.c_str());
599	read_pq.push(make_pair(id, rf));
600	}
601	if (read_pq.size()==0)
602	return false;
603	const pair<uint64_t, Read>& t = read_pq.top();
604	r=t.second; //copy strings
605	//free(t.second);
606	read_pq.pop();
607	return true;
608	}
609
610	// reads must ALWAYS requested in increasing order of their ID
611	bool ReadStream::getRead(uint64_t r_id,
612	Read& read,
613	ReadFormat read_format,
614	bool strip_slash,
615	FILE* um_out, //unmapped reads output
616	bool um_write_found) {
617	if (!fstream.file)
618	err_die("Error: calling ReadStream::getRead() with no file handle!");
619	if (r_id<last_id)
620	err_die("Error: ReadStream::getRead() called with out-of-order id#!");
621	last_id=r_id;
622	bool found=false;
623	while (!found) {
624	read.clear();
625	// Get the next read from the file
626	if (!next_read(read, read_format))
627	break;
628	if (strip_slash) {
629	string::size_type slash = read.name.rfind("/");
630	if (slash != string::npos)
631	read.name.resize(slash);
632	}
633	if ((uint64_t)atoi(read.name.c_str()) == r_id) {
634	found=true;
635	}
636	if (um_out && (um_write_found \|\| !found)) {
637	//write unmapped reads
638	fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(),
639	read.seq.c_str(), read.qual.c_str());
640	}
641	//rt.get_id(read.name, ref_str);
642	} //while reads
643	return found;
644	}
645
646
647	bool get_read_from_stream(uint64_t insert_id,
648	FILE* reads_file,
649	ReadFormat reads_format,
650	bool strip_slash,
651	Read& read,
652	FILE* um_out,
653	bool um_write_found
654	) {
655	FLineReader fr(reads_file);
656	bool found=false;
657	while(!found && !fr.isEof())
658	{
659	read.clear();
660	// Get the next read from the file
661	if (!next_fastx_read(fr, read, reads_format))
662	break;
663	if (strip_slash)
664	{
665	string::size_type slash = read.name.rfind("/");
666	if (slash != string::npos)
667	read.name.resize(slash);
668	}
669	if ((uint64_t)atoi(read.name.c_str()) == insert_id)
670	{
671	//return true;
672	found=true;
673	}
674	if (um_out && (um_write_found \|\| !found)) {
675	//write unmapped reads
676	fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(),
677	read.seq.c_str(), read.qual.c_str());
678	}
679	//rt.get_id(read.name, ref_str);
680	} //while reads
681	return found;
682	}
683
684
685	bool get_read_from_stream(uint64_t insert_id,
686	FILE* reads_file,
687	ReadFormat reads_format,
688	bool strip_slash,
689	char read_name [],
690	char read_seq [],
691	char read_alt_name [],
692	char read_qual [],
693	FILE* um_out)
694	{
695	Read read;
696	FLineReader fr(reads_file);
697	while(!fr.isEof())
698	{
699	read.clear();
700
701	// Get the next read from the file
702	if (!next_fastx_read(fr, read, reads_format))
703	break;
704
705	if (strip_slash)
706	{
707	string::size_type slash = read.name.rfind("/");
708	if (slash != string::npos)
709	read.name.resize(slash);
710	}
711
712	if ((uint64_t)atoi(read.name.c_str()) == insert_id)
713	{
714	if (read_name) strcpy(read_name, read.name.c_str());
715	if (read_seq) strcpy(read_seq, read.seq.c_str());
716	if (read_alt_name) strcpy(read_alt_name, read.alt_name.c_str());
717	if (read_qual) strcpy(read_qual, read.qual.c_str());
718	return true;
719	}
720	else if (um_out!=NULL) {
721	//write unmapped reads
722	fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(), read.seq.c_str(),
723	read.qual.c_str());
724	}
725
726	//rt.get_id(read.name, ref_str);
727	} //while reads
728
729	return false;
730	}