gclib/fqtrim/fqtrim.cpp

#include "GArgs.h"
#include "GStr.h"
#include "GHash.hh"
#include "GList.hh"
#include <ctype.h>
#include "GAlnExtend.h"

#define USAGE "Usage:\n\
fqtrim [{-5 <5adapter> -3 <3adapter>|-f <adapters_file>}] [-a <min_matchlen>]\\\n\
   [-R] [-q <minq> [-t <trim_max_len>]] [-p {64|33}] [-o <outsuffix>]\\\n\
   [-l <minlen>] [-C] [-D] [-Q] [-n <rename_prefix>] [-r <discarded.lst>]\\\n\
    <input.fq>[,<input_mates.fq>\n\
 \n\
 Trim low quality bases at the 3' end and can trim adapter sequence(s), filter\n\
 for low complexity and collapse duplicate reads.\n\
 If read pairs should be trimmed and kept together (i.e. without discarding\n\
 one read in a pair), the two file names should be given delimited by a comma\n\
 or a colon character\n\
\n\
Options:\n\
-n  rename all the reads using the <prefix> followed by a read counter;\n\
    if -C option was given, the suffix \"_x<N>\" is appended, with <N> being\n\
    the read duplication count\n\
-o  unless this parameter is '-', write the trimmed/filtered reads to \n\
    file(s) named <input>.<outsuffix> which will be created in the \n\
    current (working) directory; (writes to stdout if -o- is given);\n\
    a suffix ending with .gz, .gzip or .bz2 will enforce compression\n\
-f  file with adapter sequences to trim, each line having this format:\n\
    <5'-adapter-sequence> <3'-adapter-sequence>\n\
-5  trim the given adapter or primer sequence at the 5' end of each read\n\
    (e.g. -5 CGACAGGTTCAGAGTTCTACAGTCCGACGATC)\n\
-3  trim the given adapter sequence at the 3' end of each read\n\
    (e.g. -3 TCGTATGCCGTCTTCTGCTTG)\n\
-A  disable polyA/T trimming (enabled by default)\n\
-y  minimum length of exact match to adaptor sequence at the proper end (6)\n\
-q  trim bases with quality value lower than <minq> (starting at the 3' end)\n\
-t  for -q option, maximum trimming at the 3' end is limited to <trim_max_len>\n\
-m  maximum percentage of Ns allowed in a read after trimming (default 7)\n\
-l  minimum \"clean\" length after trimming that a read must have\n\
    in order to pass the filter (default: 16)\n\
-r  write a simple \"trash report\" file listing the discarded reads with a\n\
    one letter code for the reason why they were trashed\n\
-D  apply a low-complexity (dust) filter and discard any read that has over\n\
    50% of its length detected as low complexity\n\
-C  collapse duplicate reads and add a prefix with their count to the read\n\
    name (e.g. for microRNAs)\n\
-p  input is phred64/phred33 (use -p64 or -p33)\n\
-Q  convert quality values to the other Phred qv type\n\
-V  verbose processing\n\
 "

//-z  for -o option, the output stream(s) will be first piped into the given\n
//   <zcmd> command, which must output to stdout (e.g. -z 'bzip2 -9 -c')\n


// example 3' adapter for miRNAs: TCGTATGCCGTCTTCTGCTTG

//For paired reads sequencing:
//3' : ACACTCTTTCCCTACACGACGCTCTTCCGATCT
//5' : GATCGGAAGAGCGGTTCAGCAGGAATGCCGAG
//FILE* f_out=NULL; //stdout if not provided
//FILE* f_out2=NULL; //for paired reads
//FILE* f_in=NULL; //input fastq (stdin if not provided)
//FILE* f_in2=NULL; //for paired reads

FILE* freport=NULL;

bool debug=false;
bool verbose=false;
bool doCollapse=false;
bool doDust=false;
bool doPolyTrim=true;
bool fastaOutput=false;
bool trashReport=false;
bool revCompl=false; //also reverse complement adaptor sequences
int min_read_len=16;
double max_perc_N=7.0;
int dust_cutoff=16;
bool isfasta=false;
bool convert_phred=false;
GStr outsuffix; // -o
GStr prefix;
GStr zcmd;
int num_trimmed5=0;
int num_trimmed3=0;
int num_trimmedN=0;
int num_trimmedQ=0;
int min_trimmed5=INT_MAX;
int min_trimmed3=INT_MAX;

int qvtrim_qmin=0;
int qvtrim_max=0;  //for -q, do not trim at 3'-end more than this number of bases
int qv_phredtype=0; // could be 64 or 33 (0 means undetermined yet)
int qv_cvtadd=0; //could be -31 or +31

// adaptor matching metrics -- for X-drop ungapped extension
//const int match_reward=2;
//const int mismatch_penalty=3;
const int match_reward=2;
const int mismatch_penalty=4;
const int Xdrop=10;

const int poly_m_score=2; //match score
const int poly_mis_score=-3; //mismatch
const int poly_dropoff_score=7;
int poly_minScore=12; //i.e. an exact match of 6 bases at the proper ends WILL be trimmed

const char *polyA_seed="AAAA";
const char *polyT_seed="TTTT";

struct CASeqData {
   //positional data for every possible hexamer in an adapter
   GVec<uint16>* pz[4096]; //0-based coordinates of all possible hexamers in the adapter sequence
   GVec<uint16>* pzr[4096]; //0-based coordinates of all possible hexamers for the reverse complement of the adapter sequence
   GStr seq; //actual adapter sequence data
   GStr seqr; //reverse complement sequence
   int amlen; //fraction of adapter length matching that's
              //enough to consider the alignment
   bool use_reverse;
   CASeqData(bool rev=false):seq(),seqr(),
             amlen(0), use_reverse(rev) {
     for (int i=0;i<4096;i++) {
        pz[i]=NULL;
        pzr[i]=NULL;
        }
     }

   void update(const char* s) {
         seq=s;
         table6mers(seq.chars(), seq.length(), pz);
         amlen=iround(double(seq.length())*0.8);
         if (amlen<12)
                amlen=12;
         if (!use_reverse) return;
         //reverse complement
         seqr=s;
         int slen=seq.length();
         for (int i=0;i<slen;i++)
           seqr[i]=ntComplement(seq[slen-i-1]);
         table6mers(seqr.chars(), seqr.length(), pzr);
   }

   ~CASeqData() {
     for (int i=0;i<4096;i++) {
       delete pz[i];
       delete pzr[i];
     }
   }
};

GVec<CASeqData> adapters5;
GVec<CASeqData> adapters3;

CGreedyAlignData* gxmem_l=NULL;
CGreedyAlignData* gxmem_r=NULL;

// element in dhash:
class FqDupRec {
 public:
   int count; //how many of these reads are the same
   int len; //length of qv
   char* firstname; //optional, only if we want to keep the original read names
   char* qv;
   FqDupRec(GStr* qstr=NULL, const char* rname=NULL) {
     len=0;
     qv=NULL;
     firstname=NULL;
     count=0;
     if (qstr!=NULL) {
       qv=Gstrdup(qstr->chars());
       len=qstr->length();
       count++;
       }
     if (rname!=NULL) firstname=Gstrdup(rname);
     }
   ~FqDupRec() {
     GFREE(qv);
     GFREE(firstname);
     }
   void add(GStr& d) { //collapse another record into this one
     if (d.length()!=len)
       GError("Error at FqDupRec::add(): cannot collapse reads with different length!\n");
     count++;
     for (int i=0;i<len;i++)
       qv[i]=(qv[i]+d[i])/2;
     }
 };

void openfw(FILE* &f, GArgs& args, char opt) {
  GStr s=args.getOpt(opt);
  if (!s.is_empty()) {
      if (s=='-') f=stdout;
      else {
       f=fopen(s.chars(),"w");
       if (f==NULL) GError("Error creating file: %s\n", s.chars());
       }
     }
}

#define FWCLOSE(fh) if (fh!=NULL && fh!=stdout) fclose(fh)
#define FRCLOSE(fh) if (fh!=NULL && fh!=stdin) fclose(fh)

GHash<FqDupRec> dhash; //hash to keep track of duplicates

void addAdapter(GVec<CASeqData>& adapters, GStr& seq);
int loadAdapters(const char* fname);

void setupFiles(FILE*& f_in, FILE*& f_in2, FILE*& f_out, FILE*& f_out2,
                       GStr& s, GStr& infname, GStr& infname2);
// uses outsuffix to generate output file names and open file handles as needed

void writeRead(FILE* f_out, GStr& rname, GStr& rinfo, GStr& rseq, GStr& rqv, int& outcounter);
void trash_report(char trashcode, GStr& rname, FILE* freport);

bool getFastxRec(GLineReader& fq, GStr& rseq, GStr& rqv,
          GStr& rname, GStr& rinfo, GStr& infname);

char process_read(GStr& rname, GStr& rseq, GStr& rqv, int &l5, int &l3);
//returns 0 if the read was untouched, 1 if it was trimmed and a trash code if it was trashed

bool ntrim(GStr& rseq, int &l5, int &l3); //returns true if any trimming occured
bool qtrim(GStr& qvs, int &l5, int &l3); //return true if any trimming occured
int dust(GStr& seq);
bool trim_poly5(GStr &seq, int &l5, int &l3, const char* poly_seed); //returns true if any trimming occured
bool trim_poly3(GStr &seq, int &l5, int &l3, const char* poly_seed);
bool trim_adapter5(GStr& seq, int &l5, int &l3); //returns true if any trimming occured
bool trim_adapter3(GStr& seq, int &l5, int &l3);

void convertPhred(char* q, int len);
void convertPhred(GStr& q);

int main(int argc, char * const argv[]) {
  GArgs args(argc, argv, "YQDCRVAl:d:3:5:m:n:r:p:q:f:t:o:z:a:");
  int e;
  if ((e=args.isError())>0) {
      GMessage("%s\nInvalid argument: %s\n", USAGE, argv[e]);
      exit(224);
      }
  debug=(args.getOpt('Y')!=NULL);
  verbose=(args.getOpt('V')!=NULL);
  convert_phred=(args.getOpt('Q')!=NULL);
  doCollapse=(args.getOpt('C')!=NULL);
  doDust=(args.getOpt('D')!=NULL);
  revCompl=(args.getOpt('R')!=NULL);
  if (args.getOpt('A')) doPolyTrim=false;
  /*
  rawFormat=(args.getOpt('R')!=NULL);
  if (rawFormat) {
    GError("Sorry, raw qseq format parsing is not implemented yet!\n");
    }
  */
  prefix=args.getOpt('n');
  GStr s=args.getOpt('l');
  if (!s.is_empty())
     min_read_len=s.asInt();
  s=args.getOpt('m');
  if (!s.is_empty())
     max_perc_N=s.asDouble();
  s=args.getOpt('d');
  if (!s.is_empty()) {
     dust_cutoff=s.asInt();
     doDust=true;
     }
  s=args.getOpt('q');
  if (!s.is_empty()) {
     qvtrim_qmin=s.asInt();
     }
  s=args.getOpt('t');
  if (!s.is_empty()) {
     qvtrim_max=s.asInt();
     }
  s=args.getOpt('p');
  if (!s.is_empty()) {
     int v=s.asInt();
     if (v==33) {
        qv_phredtype=33;
        qv_cvtadd=31;
        }
      else if (v==64) {
        qv_phredtype=64;
        qv_cvtadd=-31;
        }
       else
         GMessage("%s\nInvalid value for -p option (can only be 64 or 33)!\n",USAGE);
     }
  s=args.getOpt('f');
  if (!s.is_empty()) {
   loadAdapters(s.chars());
   }
  bool fileAdapters=adapters5.Count()+adapters3.Count();
  s=args.getOpt('5');
  if (!s.is_empty()) {
    if (fileAdapters)
      GError("Error: options -5 and -f cannot be used together!\n");
    s.upper();
    addAdapter(adapters5, s);
    }
  s=args.getOpt('3');
  if (!s.is_empty()) {
    if (fileAdapters)
      GError("Error: options -3 and -f cannot be used together!\n");
      s.upper();
      addAdapter(adapters3, s);
    }
  s=args.getOpt('y');
  if (!s.is_empty()) {
     int minmatch=s.asInt();
     poly_minScore=minmatch*poly_m_score;
     }

  if (args.getOpt('o')!=NULL) outsuffix=args.getOpt('o');
                         else outsuffix="-";
  trashReport=  (args.getOpt('r')!=NULL);
  int fcount=args.startNonOpt();
  if (fcount==0) {
    GMessage(USAGE);
    exit(224);
    }
   if (fcount>1 && doCollapse) {
    GError("%s Sorry, the -C option only works with a single input.\n", USAGE);
    }
  //openfw(f_out, args, 'o');
  //if (f_out==NULL) f_out=stdout;
  if (trashReport)
    openfw(freport, args, 'r');
  char* infile=NULL;

  if (adapters5.Count()>0)
    //gxmem_l=new CGreedyAlignData(match_reward, mismatch_penalty, Xdrop-2);
        gxmem_l=new CGreedyAlignData(match_reward, mismatch_penalty, Xdrop);
  if (adapters3.Count()>0)
    gxmem_r=new CGreedyAlignData(match_reward, mismatch_penalty, Xdrop);

  while ((infile=args.nextNonOpt())!=NULL) {
    //for each input file
    int incounter=0; //counter for input reads
    int outcounter=0; //counter for output reads
    int trash_s=0; //too short from the get go
    int trash_Q=0;
    int trash_N=0;
    int trash_D=0;
    int trash_poly=0;
    int trash_A3=0;
    int trash_A5=0;
    s=infile;
    GStr infname;
    GStr infname2;
    FILE* f_in=NULL;
    FILE* f_in2=NULL;
    FILE* f_out=NULL;
    FILE* f_out2=NULL;
    bool paired_reads=false;
    setupFiles(f_in, f_in2, f_out, f_out2, s, infname, infname2);
    GLineReader fq(f_in);
    GLineReader* fq2=NULL;
    if (f_in2!=NULL) {
       fq2=new GLineReader(f_in2);
       paired_reads=true;
       }
    GStr rseq, rqv, seqid, seqinfo;
    GStr rseq2, rqv2, seqid2, seqinfo2;
    while (getFastxRec(fq, rseq, rqv, seqid, seqinfo, infname)) {
       incounter++;
       int a5=0, a3=0, b5=0, b3=0;
       char tcode=0, tcode2=0;
       tcode=process_read(seqid, rseq, rqv, a5, a3);
       if (fq2!=NULL) {
            getFastxRec(*fq2, rseq2, rqv2, seqid2, seqinfo2, infname2);
            if (seqid.substr(0,seqid.length()-1)!=seqid2.substr(0,seqid2.length()-1)) {
               GError("Error: no paired match for read %s vs %s (%s,%s)\n",
                  seqid.chars(), seqid2.chars(), infname.chars(), infname2.chars());
               }
            tcode2=process_read(seqid2, rseq2, rqv2, b5, b3);
            //decide what to do with this pair and print rseq2 if the pair makes it
            if (tcode>1 && tcode2<=1) {
               //"untrash" rseq
               if (a3-a5+1<min_read_len) {
                   a5=1;
                   if (a3<min_read_len) { a3= GMIN(rseq.length()-1, min_read_len+1); }
                   }
               tcode=1;
               }
             else if (tcode<=1 && tcode2>1) {
               //"untrash" rseq2
               if (b3-b5+1<min_read_len) {
                   b5=1;
                   if (b3<min_read_len) { b3= GMIN((rseq2.length()-1),(min_read_len+1)); }
                   }
               tcode2=1;
               }
            if (tcode<=1) { //trimmed or left intact -- write it!
               if (tcode>0) {
                 rseq2=rseq2.substr(b5,b3-b5+1);
                 if (!rqv2.is_empty()) rqv2=rqv2.substr(b5,b3-b5+1);
                 }
               int nocounter=0;
               writeRead(f_out2, seqid2, seqinfo2, rseq2, rqv2, nocounter);
               }
            } //pair read
       if (tcode>1) { //trashed
         #ifdef GDEBUG
         GMessage(" !!!!TRASH code = %c\n",tcode);
         #endif
          if (tcode=='s') trash_s++;
          else if (tcode=='A' || tcode=='T') trash_poly++;
            else if (tcode=='Q') trash_Q++;
              else if (tcode=='N') trash_N++;
               else if (tcode=='D') trash_D++;
                else if (tcode=='3') trash_A3++;
                 else if (tcode=='5') trash_A5++;
          if (trashReport) trash_report(tcode, seqid, freport);
          }
         else if (!doCollapse) { //write it
          if (tcode>0) {
            rseq=rseq.substr(a5,a3-a5+1);
            if (!rqv.is_empty()) rqv=rqv.substr(a5,a3-a5+1);
            }
         #ifdef GDEBUG
            GMessage("  After trimming:\n");
            GMessage("%s\n",rseq.chars());
         #endif
          writeRead(f_out, seqid, seqinfo, rseq, rqv, outcounter);
          }
       } //for each fastq record
    delete fq2;
    FRCLOSE(f_in);
    FRCLOSE(f_in2);
    if (doCollapse) {
       outcounter=0;
       int maxdup_count=1;
       char* maxdup_seq=NULL;
       dhash.startIterate();
       FqDupRec* qd=NULL;
       char* seq=NULL;
       while ((qd=dhash.NextData(seq))!=NULL) {
         GStr rseq(seq);
         //do the dusting here
         if (doDust) {
            int dustbases=dust(rseq);
            if (dustbases>(rseq.length()>>1)) {
               if (trashReport && qd->firstname!=NULL) {
                 fprintf(freport, "%s_x%d\tD\n",qd->firstname, qd->count);
                 }
               trash_D+=qd->count;
               continue;
               }
            }
         outcounter++;
         if (qd->count>maxdup_count) {
            maxdup_count=qd->count;
            maxdup_seq=seq;
            }
         if (isfasta) {
           if (prefix.is_empty()) {
             fprintf(f_out, ">%s_x%d\n%s\n", qd->firstname, qd->count,
                           rseq.chars());
             }
           else { //use custom read name
             fprintf(f_out, ">%s%08d_x%d\n%s\n", prefix.chars(), outcounter,
                        qd->count, rseq.chars());
             }
           }
         else { //fastq format
          if (convert_phred) convertPhred(qd->qv, qd->len);
          if (prefix.is_empty()) {
            fprintf(f_out, "@%s_x%d\n%s\n+\n%s\n", qd->firstname, qd->count,
                           rseq.chars(), qd->qv);
            }
          else { //use custom read name
            fprintf(f_out, "@%s%08d_x%d\n%s\n+\n%s\n", prefix.chars(), outcounter,
                        qd->count, rseq.chars(), qd->qv);
            }
           }
         }//for each element of dhash
       if (maxdup_count>1) {
         GMessage("Maximum read multiplicity: x %d (read: %s)\n",maxdup_count, maxdup_seq);
         }
       } //collapse entries
    if (verbose) {
       if (paired_reads) {
           GMessage(">Input files : %s , %s\n", infname.chars(), infname2.chars());
           GMessage("Number of input pairs :%9d\n", incounter);
           GMessage("         Output pairs :%9d\n", outcounter);
           }
         else {
           GMessage(">Input file : %s\n", infname.chars());
           GMessage("Number of input reads :%9d\n", incounter);
           GMessage("         Output reads :%9d\n", outcounter);
           }
       GMessage("------------------------------------\n");
       if (num_trimmed5)
          GMessage("           5' trimmed :%9d  (min. trim: %d)\n", num_trimmed5, min_trimmed5);
       if (num_trimmed3)
          GMessage("           3' trimmed :%9d  (min. trim: %d)\n", num_trimmed3, min_trimmed3);
       GMessage("------------------------------------\n");
       if (trash_s>0)
         GMessage("     Trashed by length:%9d\n", trash_s);
       if (trash_N>0)
         GMessage("         Trashed by N%%:%9d\n", trash_N);
       if (trash_Q>0)
         GMessage("Trashed by low quality:%9d\n", trash_Q);
       if (trash_poly>0)
         GMessage("   Trashed by poly-A/T:%9d\n", trash_poly);
       if (trash_A5>0)
         GMessage(" Trashed by 5' adapter:%9d\n", trash_A5);
       if (trash_A3>0)
         GMessage(" Trashed by 3' adapter:%9d\n", trash_A3);
       }
    if (trashReport) {
          FWCLOSE(freport);
          }
    FWCLOSE(f_out);
    FWCLOSE(f_out2);
   } //while each input file
 delete gxmem_l;
 delete gxmem_r;
 //getc(stdin);
}

class NData {
 public:
   int NPos[1024]; //there should be no reads longer than 1K ?
   int NCount;
   int end5;
   int end3;
   int seqlen;
   double perc_N; //percentage of Ns in end5..end3 range only!
   const char* seq;
   bool valid;
   NData() {
    seqlen=0;
    NCount=0;
    end5=0;
    end3=0;
    seq=NULL;
    perc_N=0;
    valid=true;
    }
   void init(GStr& rseq) {
     NCount=0;
     perc_N=0;
     seqlen=rseq.length();
     seq=rseq.chars();
     end5=1;
     end3=seqlen;
     for (int i=0;i<seqlen;i++)
        if (seq[i]=='N') {// if (!ichrInStr(rseq[i], "ACGT")
           NPos[NCount]=i;
           NCount++;
           }
     perc_N=(NCount*100.0)/seqlen;
     valid=true;
     }
  void N_calc() { //only in the region after trimming
     int n=0;
     for (int i=end3-1;i<end5;i++) {
       if (seq[i]=='N') n++;
       }
     perc_N=(n*100.0)/(end5-end3+1);
     }
 };

static NData feat;
int perc_lenN=12; // incremental distance from ends, in percentage of
          // sequence length, where N-trimming is done (default:12 %) (autolimited to 20)

void N_analyze(int l5, int l3, int p5, int p3) {
/* assumes feat was filled properly */
 int old_dif, t5,t3,v;
 if (l3<l5+2 || p5>p3 ) {
   feat.end5=l5+1;
   feat.end3=l3+1;
   return;
   }

 t5=feat.NPos[p5]-l5;
 t3=l3-feat.NPos[p3];
 old_dif=p3-p5;
 v=(int)((((double)(l3-l5))*perc_lenN)/100);
 if (v>20) v=20; /* enforce N-search limit for very long reads */
 if (t5 < v ) {
   l5=feat.NPos[p5]+1;
   p5++;
   }
 if (t3 < v) {
   l3=feat.NPos[p3]-1;
   p3--;
   }
 /* restNs=p3-p5; number of Ns in the new CLR */
 if (p3-p5==old_dif) { /* no change, return */
           /* don't trim if this may shorten the read too much */
           //if (l5-l3<min_read_len) return;
           feat.end5=l5+1;
           feat.end3=l3+1;
           return;
           }
    else
      N_analyze(l5,l3, p5,p3);
}


bool qtrim(GStr& qvs, int &l5, int &l3) {
if (qvtrim_qmin==0 || qvs.is_empty()) return false;
if (qv_phredtype==0) {
  //try to guess the Phred type
  int vmin=256, vmax=0;
  for (int i=0;i<qvs.length();i++) {
     if (vmin>qvs[i]) vmin=qvs[i];
     if (vmax<qvs[i]) vmax=qvs[i];
     }
  if (vmin<64) { qv_phredtype=33; qv_cvtadd=31; }
  if (vmax>95) { qv_phredtype=64; qv_cvtadd=-31; }
  if (qv_phredtype==0) {
    GError("Error: couldn't determine Phred type, please use the -p33 or -p64 !\n");
    }
  if (verbose)
    GMessage("Input reads have Phred-%d quality values.\n", (qv_phredtype==33 ? 33 : 64));
  } //guessing Phred type
for (l3=qvs.length()-1;l3>2;l3--) {
  if (qvs[l3]-qv_phredtype>=qvtrim_qmin && qvs[l3-1]-qv_phredtype>=qvtrim_qmin) break;
  }
//just in case, check also the 5' the end (?)
for (l5=0;l5<qvs.length()-3;l5++) {
  if (qvs[l5]-qv_phredtype>=qvtrim_qmin && qvs[l5+1]-qv_phredtype>=qvtrim_qmin) break;
  }
if (qvtrim_max>0) {
  if (qvs.length()-1-l3>qvtrim_max) l3=qvs.length()-1-qvtrim_max;
  if (l5>qvtrim_max) l5=qvtrim_max;
  }
return (l5>0 || l3<qvs.length()-1);
}

bool ntrim(GStr& rseq, int &l5, int &l3) {
 //count Ns in the sequence, trim N-rich ends
 feat.init(rseq);
 l5=feat.end5-1;
 l3=feat.end3-1;
 N_analyze(feat.end5-1, feat.end3-1, 0, feat.NCount-1);
 if (l5==feat.end5-1 && l3==feat.end3-1) {
    if (feat.perc_N>max_perc_N) {
           feat.valid=false;
           l3=l5+1;
           return true;
           }
      else {
       return false; //nothing changed
       }
    }
 l5=feat.end5-1;
 l3=feat.end3-1;
 if (l3-l5+1<min_read_len) {
   feat.valid=false;
   return true;
   }
 feat.N_calc();

 if (feat.perc_N>max_perc_N) {
      feat.valid=false;
      l3=l5+1;
      return true;
      }
 return true;
 }

//--------------- dust functions ----------------
class DNADuster {
 public:
  int dustword;
  int dustwindow;
  int dustwindow2;
  int dustcutoff;
  int mv, iv, jv;
  int counts[32*32*32];
  int iis[32*32*32];
  DNADuster(int cutoff=16, int winsize=32, int wordsize=3) {
    dustword=wordsize;
    dustwindow=winsize;
    dustwindow2 = (winsize>>1);
    dustcutoff=cutoff;
    mv=0;
    iv=0;
    jv=0;
    }
  void setWindowSize(int value) {
    dustwindow = value;
    dustwindow2 = (dustwindow >> 1);
    }
  void setWordSize(int value) {
    dustword=value;
    }
void wo1(int len, const char* s, int ivv) {
  int i, ii, j, v, t, n, n1, sum;
  int js, nis;
  n = 32 * 32 * 32;
  n1 = n - 1;
  nis = 0;
  i = 0;
  ii = 0;
  sum = 0;
  v = 0;
  for (j=0; j < len; j++, s++) {
        ii <<= 5;
        if (*s<=32) {
           i=0;
           continue;
           }
        ii |= *s - 'A'; //assume uppercase!
        ii &= n1;
        i++;
        if (i >= dustword) {
              for (js=0; js < nis && iis[js] != ii; js++) ;
              if (js == nis) {
                    iis[nis] = ii;
                    counts[ii] = 0;
                    nis++;
              }
              if ((t = counts[ii]) > 0) {
                    sum += t;
                    v = 10 * sum / j;
                    if (mv < v) {
                          mv = v;
                          iv = ivv;
                          jv = j;
                    }
              }
              counts[ii]++;
        }
  }
}

int wo(int len, const char* s, int* beg, int* end) {
      int i, l1;
      l1 = len - dustword + 1;
      if (l1 < 0) {
            *beg = 0;
            *end = len - 1;
            return 0;
            }
      mv = 0;
      iv = 0;
      jv = 0;
      for (i=0; i < l1; i++) {
            wo1(len-i, s+i, i);
            }
      *beg = iv;
      *end = iv + jv;
      return mv;
 }

void dust(const char* seq, char* seqmsk, int seqlen, int cutoff=0) { //, maskFunc maskfn) {
  int i, j, l, a, b, v;
  if (cutoff==0) cutoff=dustcutoff;
  a=0;b=0;
  //GMessage("Dust cutoff=%d\n", cutoff);
  for (i=0; i < seqlen; i += dustwindow2) {
        l = (seqlen > i+dustwindow) ? dustwindow : seqlen-i;
        v = wo(l, seq+i, &a, &b);
        if (v > cutoff) {
           //for (j = a; j <= b && j < dustwindow2; j++) {
           for (j = a; j <= b; j++) {
                    seqmsk[i+j]='N';//could be made lowercase instead
                    }
           }
         }
//return first;
 }
};

static DNADuster duster;

int dust(GStr& rseq) {
 char* seq=Gstrdup(rseq.chars());
 duster.dust(rseq.chars(), seq, rseq.length(), dust_cutoff);
 //check the number of Ns:
 int ncount=0;
 for (int i=0;i<rseq.length();i++) {
   if (seq[i]=='N') ncount++;
   }
 GFREE(seq);
 return ncount;
 }

struct SLocScore {
  int pos;
  int score;
  SLocScore(int p=0,int s=0) {
    pos=p;
    score=s;
    }
  void set(int p, int s) {
    pos=p;
    score=s;
    }
  void add(int p, int add) {
    pos=p;
    score+=add;
    }
};

bool trim_poly3(GStr &seq, int &l5, int &l3, const char* poly_seed) {
 if (!doPolyTrim) return false;
 int rlen=seq.length();
 l5=0;
 l3=rlen-1;
 int32 seedVal=*(int32*)poly_seed;
 char polyChar=poly_seed[0];
 //assumes N trimming was already done
 //so a poly match should be very close to the end of the read
 // -- find the initial match (seed)
 int lmin=GMAX((rlen-16), 0);
 int li;
 for (li=rlen-4;li>lmin;li--) {
   if (seedVal==*(int*)&(seq[li])) {
      break;
      }
   }
 if (li<=lmin) return false;
 //seed found, try to extend it both ways
 //extend right
 int ri=li+3;
 SLocScore loc(ri, poly_m_score<<2);
 SLocScore maxloc(loc);
 //extend right
 while (ri<rlen-1) {
   ri++;
   if (seq[ri]==polyChar) {
                loc.add(ri,poly_m_score);
                }
   else if (seq[ri]=='N') {
                loc.add(ri,0);
                }
   else { //mismatch
        loc.add(ri,poly_mis_score);
        if (maxloc.score-loc.score>poly_dropoff_score) break;
        }
   if (maxloc.score<=loc.score) {
      maxloc=loc;
      }
   }
 ri=maxloc.pos;
 if (ri<rlen-6) return false; //no trimming wanted, too far from 3' end
 //ri = right boundary for the poly match
 //extend left
 loc.set(li, maxloc.score);
 maxloc.pos=li;
 while (li>0) {
    li--;
    if (seq[li]==polyChar) {
                 loc.add(li,poly_m_score);
                 }
    else if (seq[li]=='N') {
                 loc.add(li,0);
                 }
    else { //mismatch
         loc.add(li,poly_mis_score);
         if (maxloc.score-loc.score>poly_dropoff_score) break;
         }
    if (maxloc.score<=loc.score) {
       maxloc=loc;
       }
    }
li=maxloc.pos;
if ((maxloc.score==poly_minScore && ri==rlen-1) ||
    (maxloc.score>poly_minScore && ri>=rlen-3) ||
    (maxloc.score>(poly_minScore*3) && ri>=rlen-8)) {
  //trimming this li-ri match at 3' end
    l3=li-1;
    if (l3<0) l3=0;
    return true;
    }
return false;
}

bool trim_poly5(GStr &seq, int &l5, int &l3, const char* poly_seed) {
 if (!doPolyTrim) return false;
 int rlen=seq.length();
 l5=0;
 l3=rlen-1;
 int32 seedVal=*(int32*)poly_seed;
 char polyChar=poly_seed[0];
 //assumes N trimming was already done
 //so a poly match should be very close to the end of the read
 // -- find the initial match (seed)
 int lmax=GMIN(12, rlen-4);//how far from 5' end to look for 4-mer seeds
 int li;
 for (li=0;li<=lmax;li++) {
   if (seedVal==*(int*)&(seq[li])) {
      break;
      }
   }
 if (li>lmax) return false;
 //seed found, try to extend it both ways
 //extend left
 int ri=li+3; //save rightmost base of the seed
 SLocScore loc(li, poly_m_score<<2);
 SLocScore maxloc(loc);
 while (li>0) {
    li--;
    if (seq[li]==polyChar) {
                 loc.add(li,poly_m_score);
                 }
    else if (seq[li]=='N') {
                 loc.add(li,0);
                 }
    else { //mismatch
         loc.add(li,poly_mis_score);
         if (maxloc.score-loc.score>poly_dropoff_score) break;
         }
    if (maxloc.score<=loc.score) {
       maxloc=loc;
       }
    }
 li=maxloc.pos;
 if (li>5) return false; //no trimming wanted, too far from 5' end
 //li = right boundary for the poly match

 //extend right
 loc.set(ri, maxloc.score);
 maxloc.pos=ri;
 while (ri<rlen-1) {
   ri++;
   if (seq[ri]==polyChar) {
                loc.add(ri,poly_m_score);
                }
   else if (seq[ri]=='N') {
                loc.add(ri,0);
                }
   else { //mismatch
        loc.add(ri,poly_mis_score);
        if (maxloc.score-loc.score>poly_dropoff_score) break;
        }
   if (maxloc.score<=loc.score) {
      maxloc=loc;
      }
   }
ri=maxloc.pos;
if ((maxloc.score==poly_minScore && li==0) ||
     (maxloc.score>poly_minScore && li<2)
     || (maxloc.score>(poly_minScore*3) && li<8)) {
    //adjust l5 to reflect this trimming of 5' end
    l5=ri+1;
    if (l5>rlen-1) l5=rlen-1;
    return true;
    }
return false;
}

bool trim_adapter3(GStr& seq, int&l5, int &l3) {
 if (adapters3.Count()==0) return false;
 int rlen=seq.length();
 l5=0;
 l3=rlen-1;
 bool trimmed=false;
 GStr wseq(seq);
 int wlen=rlen;
 GXSeqData seqdata;
 for (int ai=0;ai<adapters3.Count();ai++) {
  seqdata.update(adapters3[ai].seq.chars(), adapters3[ai].seq.length(),
       adapters3[ai].pz, wseq.chars(), wlen, adapters3[ai].amlen);
  GXAlnInfo* bestaln=match_RightEnd(seqdata, gxmem_r, 86);
  if (bestaln) {
     trimmed=true;
     //keep unmatched region on the left OR right (the longer one)
     if (bestaln->sl > wlen-bestaln->sr) {
         //keep left side
         l3-=(wlen-bestaln->sl+1);
         if (l3<0) l3=0;
         }
     else { //keep right side
         l5+=bestaln->sr;
         if (l5>=rlen) l5=rlen-1;
         }
     delete bestaln;
     if (l3-l5+1<min_read_len) return true;
     wseq=seq.substr(l5,l3-l5+1);
     wlen=wseq.length();
     }
  }//for each 5' adapter
  return trimmed;
 }

bool trim_adapter5(GStr& seq, int&l5, int &l3) {
 if (adapters5.Count()==0) return false;
 int rlen=seq.length();
 l5=0;
 l3=rlen-1;
 bool trimmed=false;
 GStr wseq(seq);
 int wlen=rlen;
 GXSeqData seqdata;
 for (int ai=0;ai<adapters5.Count();ai++) {
   seqdata.update(adapters5[ai].seq.chars(), adapters5[ai].seq.length(),
       adapters5[ai].pz, wseq.chars(), wlen, adapters5[ai].amlen);
  GXAlnInfo* bestaln=match_LeftEnd(seqdata, gxmem_l, 90);
  if (bestaln) {
     trimmed=true;
     if (bestaln->sl-1 > wlen-bestaln->sr) {
         //keep left side
         l3-=(wlen-bestaln->sl+1);
         if (l3<0) l3=0;
         }
     else { //keep right side
         l5+=bestaln->sr;
         if (l5>=rlen) l5=rlen-1;
         }
     delete bestaln;
     if (l3-l5+1<min_read_len) return true;
     wseq=seq.substr(l5,l3-l5+1);
     wlen=wseq.length();
     }
  }//for each 5' adapter
  return trimmed;
}

//convert qvs to/from phred64 from/to phread33
void convertPhred(GStr& q) {
 for (int i=0;i<q.length();i++) q[i]+=qv_cvtadd;
}

void convertPhred(char* q, int len) {
 for (int i=0;i<len;i++) q[i]+=qv_cvtadd;
}

bool getFastxRec(GLineReader& fq, GStr& rseq, GStr& rqv,
          GStr& rname, GStr& rinfo, GStr& infname) {
 rseq="";
 rqv="";
 rname="";
 rinfo="";
 if (fq.eof()) return false;
 char* l=fq.getLine();
 while (l!=NULL && (l[0]==0 || isspace(l[0]))) l=fq.getLine(); //ignore empty lines
 if (l==NULL) return false;
 /* if (rawFormat) {
      //TODO: implement raw qseq parsing here
      //if (raw type=N) then continue; //skip invalid/bad records
      } //raw qseq format
 else { // FASTQ or FASTA */
 isfasta=(l[0]=='>');
 if (!isfasta && l[0]!='@') GError("Error: fasta/fastq record marker not found(%s)\n%s\n",
      infname.chars(), l);
 GStr s(l);
 rname=&(l[1]);
 for (int i=0;i<rname.length();i++)
    if (rname[i]<=' ') {
       if (i<rname.length()-2) rinfo=rname.substr(i+1);
       rname.cut(i);
       break;
       }
  //now get the sequence
 if ((l=fq.getLine())==NULL)
      GError("Error: unexpected EOF after header for read %s (%s)\n",
                   rname.chars(), infname.chars());
 rseq=l; //this must be the DNA line
 while ((l=fq.getLine())!=NULL) {
      //seq can span multiple lines
      if (l[0]=='>' || l[0]=='+') {
           fq.pushBack();
           break; //
           }
      rseq+=l;
      } //check for multi-line seq
 if (!isfasta) { //reading fastq quality values, which can also be multi-line
    if ((l=fq.getLine())==NULL)
        GError("Error: unexpected EOF after sequence for %s\n", rname.chars());
    if (l[0]!='+') GError("Error: fastq qv header marker not detected!\n");
    if ((l=fq.getLine())==NULL)
        GError("Error: unexpected EOF after qv header for %s\n", rname.chars());
    rqv=l;
    //if (rqv.length()!=rseq.length())
    //  GError("Error: qv len != seq len for %s\n", rname.chars());
    while (rqv.length()<rseq.length() && ((l=fq.getLine())!=NULL)) {
      rqv+=l; //append to qv string
      }
    }// fastq
 // } //<-- FASTA or FASTQ
 rseq.upper();
 return true;
}

#ifdef GDEBUG
void showTrim(GStr& s, int l5, int l3) {
  if (l5>0 || l3==0) {
    color_bg(c_red);
    }
  for (int i=0;i<s.length()-1;i++) {
    if (i && i==l5) color_resetbg();
    fprintf(stderr, "%c", s[i]);
    if (i && i==l3) color_bg(c_red);
   }
  fprintf(stderr, "%c", s[s.length()-1]);
  color_reset();
  fprintf(stderr, "\n");
}
#endif

char process_read(GStr& rname, GStr& rseq, GStr& rqv, int &l5, int &l3) {
 //returns 0 if the read was untouched, 1 if it was just trimmed
 // and a trash code if it was trashed
 l5=0;
 l3=rseq.length()-1;
 #ifdef GDEBUG
   //rseq.reverse();
   GMessage(">%s\n", rname.chars());
   GMessage("%s\n",rseq.chars());
 #endif
 if (l3-l5+1<min_read_len) {
   return 's';
   }
GStr wseq(rseq.chars());
GStr wqv(rqv.chars());
int w5=l5;
int w3=l3;
//first do the q-based trimming
if (qvtrim_qmin!=0 && !wqv.is_empty() && qtrim(wqv, w5, w3)) { // qv-threshold trimming
   if (w3-w5+1<min_read_len) {
     return 'Q'; //invalid read
     }
    //-- keep only the w5..w3 range
   l5=w5;
   l3=w3;
   wseq=wseq.substr(w5, w3-w5+1);
   if (!wqv.is_empty())
      wqv=wqv.substr(w5, w3-w5+1);
   } //qv trimming
// N-trimming on the remaining read seq
if (ntrim(wseq, w5, w3)) {
   //GMessage("before: %s\n",wseq.chars());
   //GMessage("after : %s (%d)\n",wseq.substr(w5,w3-w5+1).chars(),w3-w5+1);
   l5+=w5;
   l3-=(wseq.length()-1-w3);
   if (w3-w5+1<min_read_len) {
     return 'N'; //to be trashed
     }
    //-- keep only the w5..w3 range
   wseq=wseq.substr(w5, w3-w5+1);
   if (!wqv.is_empty())
      wqv=wqv.substr(w5, w3-w5+1);
   w5=0;
   w3=wseq.length()-1;
   }
char trim_code;
do {
  trim_code=0;
  if (trim_poly5(wseq, w5, w3, polyA_seed)) {
      trim_code='A';
      }
  else if (trim_poly5(wseq, w5, w3, polyT_seed)) {
      trim_code='T';
      }
  else if (trim_adapter5(wseq, w5, w3)) {
      trim_code='5';
      }
  if (trim_code) {
     #ifdef GDEBUG
      GMessage("#### TRIM by '%c' code ( w5-w3 = %d-%d ):\n",trim_code, w5,w3);
      showTrim(wseq, w5, w3);
     #endif
     int trimlen=wseq.length()-(w3-w5+1);
     num_trimmed5++;
     if (trimlen<min_trimmed5)
         min_trimmed5=trimlen;
     l5+=w5;
     l3-=(wseq.length()-1-w3);
     if (w3-w5+1<min_read_len) {
         return trim_code;
         }
      //-- keep only the w5..w3 range
      wseq=wseq.substr(w5, w3-w5+1);
      if (!wqv.is_empty())
         wqv=wqv.substr(w5, w3-w5+1);
      }// trimmed at 5' end
} while (trim_code);

do {
  trim_code=0;
  if (trim_poly3(wseq, w5, w3, polyA_seed)) {
      trim_code='A';
      }
  else if (trim_poly3(wseq, w5, w3, polyT_seed)) {
      trim_code='T';
      }
  else if (trim_adapter3(wseq, w5, w3)) {
      trim_code='3';
      }
  if (trim_code) {
     #ifdef GDEBUG
     GMessage("#### TRIM by '%c' code ( w5-w3 = %d-%d ):\n",trim_code, w5,w3);
     showTrim(wseq, w5, w3);
     #endif
     int trimlen=wseq.length()-(w3-w5+1);
     num_trimmed3++;
     if (trimlen<min_trimmed3)
         min_trimmed3=trimlen;
     l5+=w5;
     l3-=(wseq.length()-1-w3);
     if (w3-w5+1<min_read_len) {
         return trim_code;
         }
      //-- keep only the w5..w3 range
      wseq=wseq.substr(w5, w3-w5+1);
      if (!wqv.is_empty())
         wqv=wqv.substr(w5, w3-w5+1);
      }//trimming at 3' end
} while (trim_code);


if (doCollapse) {
   //keep read for later
   FqDupRec* dr=dhash.Find(wseq.chars());
   if (dr==NULL) { //new entry
          //if (prefix.is_empty())
             dhash.Add(wseq.chars(),
                  new FqDupRec(&wqv, rname.chars()));
          //else dhash.Add(wseq.chars(), new FqDupRec(wqv.chars(),wqv.length()));
         }
      else
         dr->add(wqv);
   } //collapsing duplicates
 else { //not collapsing duplicates
   //apply the dust filter now
   if (doDust) {
     int dustbases=dust(wseq);
     if (dustbases>(wseq.length()>>1)) {
        return 'D';
        }
     }
   } //not collapsing duplicates
return (l5>0 || l3<rseq.length()-1) ? 1 : 0;
}

void printHeader(FILE* f_out, char recmarker, GStr& rname, GStr& rinfo) {
 //GMessage("printing Header..%c%s\n",recmarker, rname.chars());
 if (rinfo.is_empty()) fprintf(f_out, "%c%s\n",recmarker,rname.chars());
  else fprintf(f_out, "%c%s %s\n",recmarker, rname.chars(), rinfo.chars());
 }

void writeRead(FILE* f_out, GStr& rname, GStr& rinfo, GStr& rseq, GStr& rqv, int& outcounter) {
   outcounter++;
   bool asFasta=(rqv.is_empty() || fastaOutput);
   if (asFasta) {
    if (prefix.is_empty()) {
       printHeader(f_out, '>',rname,rinfo);
       fprintf(f_out, "%s\n", rseq.chars()); //plain one-line fasta for now
       }
      else {
       fprintf(f_out, ">%s%08d\n%s\n", prefix.chars(), outcounter,
                          rseq.chars());
       }
     }
   else {  //fastq
    if (convert_phred) convertPhred(rqv);
    if (prefix.is_empty()) {
       printHeader(f_out, '@', rname, rinfo);
       fprintf(f_out, "%s\n+\n%s\n", rseq.chars(), rqv.chars());
       }
      else
       fprintf(f_out, "@%s_%08d\n%s\n+\n%s\n", prefix.chars(), outcounter,
                          rseq.chars(),rqv.chars() );
     }
}

void trash_report(char trashcode, GStr& rname, FILE* freport) {
 if (freport==NULL || trashcode<=' ') return;
 if (trashcode=='3' || trashcode=='5') {
   fprintf(freport, "%s\ta%c\n",rname.chars(),trashcode);
   }
 else {
   fprintf(freport, "%s\t%c\n",rname.chars(),trashcode);
   }
 //tcounter++;
}

GStr getFext(GStr& s, int* xpos=NULL) {
 //s must be a filename without a path
 GStr r("");
 if (xpos!=NULL) *xpos=0;
 if (s.is_empty() || s=="-") return r;
 int p=s.rindex('.');
 int d=s.rindex('/');
 if (p<=0 || p>s.length()-2 || p<s.length()-7 || p<d) return r;
 r=s.substr(p+1);
 if (xpos!=NULL) *xpos=p+1;
 r.lower();
 return r;
 }

void baseFileName(GStr& fname) {
 //remove all known extensions, like .txt,fq,fastq,fasta,fa)(.gz .gzip .bz2 .bzip2) .
 int xpos=0;
 GStr fext=getFext(fname, &xpos);
 if (xpos<=1) return;
 bool extdel=false;
 GStr f2;
 if (fext=="z" || fext=="zip") {
   extdel=true;
   }
  else if (fext.length()>=2) {
   f2=fext.substr(0,2);
   extdel=(f2=="gz" || f2=="bz");
   }
 if (extdel) {
   fname.cut(xpos-1);
   fext=getFext(fname, &xpos);
   if (xpos<=1) return;
   }
 extdel=false;
 if (fext=="f" || fext=="fq" || fext=="txt" || fext=="seq" || fext=="sequence") {
   extdel=true;
   }
  else if (fext.length()>=2) {
   extdel=(fext.substr(0,2)=="fa");
   }
 if (extdel) fname.cut(xpos-1);
 GStr fncp(fname);
 fncp.lower();
 fncp.chomp("seq");
 fncp.chomp("sequence");
 fncp.trimR("_.");
 if (fncp.length()<fname.length()) fname.cut(fncp.length());
}

FILE* prepOutFile(GStr& infname, GStr& pocmd) {
  FILE* f_out=NULL;
  GStr fname(getFileName(infname.chars()));
  //eliminate known extensions
  baseFileName(fname);
  if (outsuffix.is_empty() || outsuffix=="-") { return stdout; }
    else if (pocmd.is_empty()) {
               GStr oname(fname);
               oname.append('.');
               oname.append(outsuffix);
               f_out=fopen(oname.chars(),"w");
               if (f_out==NULL) GError("Error: cannot create '%s'\n",oname.chars());
               }
            else {
              GStr oname(">");
              oname.append(fname);
              oname.append('.');
              oname.append(outsuffix);
              pocmd.append(oname);
              f_out=popen(pocmd.chars(), "w");
              if (f_out==NULL) GError("Error: cannot popen '%s'\n",pocmd.chars());
              }
 return f_out;
}

void guess_unzip(GStr& fname, GStr& picmd) {
 GStr fext=getFext(fname);
 if (fext=="gz" || fext=="gzip" || fext=="z") {
    picmd="gzip -cd ";
    }
   else if (fext=="bz2" || fext=="bzip2" || fext=="bz" || fext=="bzip") {
    picmd="bzip2 -cd ";
    }
}

void addAdapter(GVec<CASeqData>& adapters, GStr& seq) {
//TODO: prepare CASeqData here, and collect hexamers as well
 CASeqData adata(revCompl);
 int idx=adapters.Add(adata);
 if (idx<0) GError("Error: failed to add adaptor!\n");
 adapters[idx].update(seq.chars());
}


int loadAdapters(const char* fname) {
  GLineReader lr(fname);
  char* l;
  while ((l=lr.nextLine())!=NULL) {
   if (lr.length()<=3 || l[0]=='#') continue;
   if ( l[0]==' ' || l[0]=='\t' || l[0]==',' ||
        l[0]==';'|| l[0]==':' ) {
      int i=1;
      while (l[i]!=0 && isspace(l[i])) {
        i++;
        }
      if (l[i]!=0) {
        GStr s(&(l[i]));
      #ifdef GDEBUG
          //s.reverse();
      #endif
        addAdapter(adapters3, s);
        continue;
        }
      }
    else {
      GStr s(l);
      s.startTokenize("\t ;,:");
      GStr a5,a3;
      if (s.nextToken(a5))
         s.nextToken(a3);
      a5.upper();
      a3.upper();
     #ifdef GDEBUG
     //   a5.reverse();
     //   a3.reverse();
     #endif
      addAdapter(adapters5, a5);
      addAdapter(adapters3, a3);
      }
   }
   return adapters5.Count()+adapters3.Count();
}

void setupFiles(FILE*& f_in, FILE*& f_in2, FILE*& f_out, FILE*& f_out2,
                       GStr& s, GStr& infname, GStr& infname2) {
// uses outsuffix to generate output file names and open file handles as needed
 infname="";
 infname2="";
 f_in=NULL;
 f_in2=NULL;
 f_out=NULL;
 f_out2=NULL;
 //analyze outsuffix intent
 GStr pocmd;
 if (outsuffix=="-") {
    f_out=stdout;
    }
   else {
    GStr ox=getFext(outsuffix);
    if (ox.length()>2) ox=ox.substr(0,2);
    if (ox=="gz") pocmd="gzip -9 -c ";
        else if (ox=="bz") pocmd="bzip2 -9 -c ";
    }
 if (s=="-") {
    f_in=stdin;
    infname="stdin";
    f_out=prepOutFile(infname, pocmd);
    return;
    } // streaming from stdin
 s.startTokenize(",:");
 s.nextToken(infname);
 bool paired=s.nextToken(infname2);
 if (fileExists(infname.chars())==0)
    GError("Error: cannot find file %s!\n",infname.chars());
 GStr fname(getFileName(infname.chars()));
 GStr picmd;
 guess_unzip(fname, picmd);
 if (picmd.is_empty()) {
   f_in=fopen(infname.chars(), "r");
   if (f_in==NULL) GError("Error opening file '%s'!\n",infname.chars());
   }
  else {
   picmd.append(infname);
   f_in=popen(picmd.chars(), "r");
   if (f_in==NULL) GError("Error at popen %s!\n", picmd.chars());
   }
 if (f_out==stdout) {
   if (paired) GError("Error: output suffix required for paired reads\n");
   return;
   }
 f_out=prepOutFile(infname, pocmd);
 if (!paired) return;
 if (doCollapse) GError("Error: sorry, -C option cannot be used with paired reads!\n");
 // ---- paired reads:-------------
 if (fileExists(infname2.chars())==0)
     GError("Error: cannot find file %s!\n",infname2.chars());
 picmd="";
 GStr fname2(getFileName(infname2.chars()));
 guess_unzip(fname2, picmd);
 if (picmd.is_empty()) {
   f_in2=fopen(infname2.chars(), "r");
   if (f_in2==NULL) GError("Error opening file '%s'!\n",infname2.chars());
   }
  else {
   picmd.append(infname2);
   f_in2=popen(picmd.chars(), "r");
   if (f_in2==NULL) GError("Error at popen %s!\n", picmd.chars());
   }
 f_out2=prepOutFile(infname2, pocmd);
}
Revision:	177
Committed:	Thu Feb 16 07:10:32 2012 UTC (7 years, 7 months ago) by gpertea
File size:	43339 byte(s)
Log Message:	collectSeeds unified