Added fasta program

Author: pmelsted

Makefile

@@ -7,7 +7,7 @@ CXXFLAGS+= -O3 -DSEQAN_ENABLE_TESTING=0 -DSEQAN_ENABLE_DEBUG=0 -DSEQAN_HAS_ZLIB=
 LDLIBS=-lz -lssl -lcrypto 
 
 
-TARGET = bamhash_checksum_bam bamhash_checksum_fastq
+TARGET = bamhash_checksum_bam bamhash_checksum_fastq bamhash_checksum_fasta
 all: $(TARGET)
 
 bamhash_checksum_bam: bamhash_checksum_common.o bamhash_checksum_bam.o
@@ -16,5 +16,9 @@ bamhash_checksum_bam: bamhash_checksum_common.o bamhash_checksum_bam.o
 bamhash_checksum_fastq: bamhash_checksum_common.o bamhash_checksum_fastq.o
 	 $(CXX) $(LDFLAGS) -o $@ $^ $(LDLIBS)
 
+bamhash_checksum_fasta: bamhash_checksum_common.o bamhash_checksum_fasta.o
+	 $(CXX) $(LDFLAGS) -o $@ $^ $(LDLIBS)
+
+
 clean:
 	$(RM) *.o *~ $(TARGET)  

bamhash_checksum_fasta.cpp

@@ -0,0 +1,148 @@
+#include <iostream>
+#include <seqan/stream.h>
+#include <seqan/seq_io.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <cmath>
+#include <cstdlib>
+#include <stdint.h>
+#include <vector>
+#include <seqan/arg_parse.h>
+
+#include "bamhash_checksum_common.h"
+
+struct Fastainfo {
+  std::vector<std::string> fastafiles;
+  bool debug;
+  bool noReadNames;
+
+  Fastainfo() : debug(false), noReadNames(false) {}
+
+};
+
+seqan::ArgumentParser::ParseResult
+parseCommandLine(Fastainfo& options, int argc, char const **argv) {
+  // Setup ArgumentParser.
+  seqan::ArgumentParser parser("bamhash_checksum_fasta");
+
+  setShortDescription(parser, "Checksum of a set of fasta files");
+  setVersion(parser, BAMHASH_VERSION);
+  setDate(parser, "May 2015");
+
+  addUsageLine(parser, "[\\fIOPTIONS\\fP] \\fI<in1.fasta>\\fP [\\fIin2.fasta ... \\fP]");
+  addDescription(parser, "Program for checksum of sequence reads. ");
+
+  addArgument(parser, seqan::ArgParseArgument(seqan::ArgParseArgument::INPUTFILE,"fastafiles", true));
+
+  setValidValues(parser, 0,"fa fa.gz fasta fasta.gz");
+
+  addSection(parser, "Options");
+  //add debug option:
+  addOption(parser, seqan::ArgParseOption("d", "debug", "Debug mode. Prints full hex for each read to stdout"));
+  addOption(parser, seqan::ArgParseOption("R", "no-readnames", "Do not use read names as part of checksum"));
+
+  // Parse command line.
+  seqan::ArgumentParser::ParseResult res = seqan::parse(parser, argc, argv);
+  if (res != seqan::ArgumentParser::PARSE_OK) {
+    return res;
+  }
+
+  options.debug = seqan::isSet(parser, "debug");
+  options.noReadNames = seqan::isSet(parser, "no-readnames");
+
+
+  options.fastafiles = getArgumentValues(parser, 0);
+
+  
+  return seqan::ArgumentParser::PARSE_OK;
+}
+
+
+int main(int argc, char const **argv) {
+  Fastainfo info; // Define structure variable
+  seqan::ArgumentParser::ParseResult res = parseCommandLine(info, argc, argv); // Parse the command line.
+
+  if (res != seqan::ArgumentParser::PARSE_OK) {
+    return res == seqan::ArgumentParser::PARSE_ERROR;
+  }
+
+  // Define:
+  uint64_t sum = 0;
+  unsigned count = 0;
+  seqan::StringSet<seqan::CharString> idSub;
+  seqan::CharString string2hash;
+
+  seqan::CharString id;
+  seqan::CharString seq;
+  hash_t hex;
+
+
+  // Open GZStream
+  seqan::Stream<seqan::GZFile> gzStream;
+
+
+
+  for (int i = 0; i < info.fastafiles.size(); i++) {
+    const char* fasta = info.fastafiles[i].c_str();
+    
+    if (!open(gzStream, fasta, "r")) {
+      std::cerr << "ERROR: Could not open the file: " << fasta << " for reading.\n";
+      return 1;
+    }
+
+    //Setup RecordReader for reading FASTA file from gzip-compressed file
+    seqan::RecordReader<seqan::Stream<seqan::GZFile>, seqan::SinglePass<> > reader(gzStream);
+
+
+
+    // Read record
+    while (!atEnd(reader)) {
+      if (readRecord(id, seq, reader, seqan::Fasta()) != 0) {
+        if (atEnd(reader)) {
+          std::cerr << "WARNING: Could not continue reading " << fasta <<  " at line: " << count+1 << ".\n";
+          return 1;
+        }
+        std::cerr << "ERROR: Could not read from " << fasta << "\n";
+        return 1;
+      }
+
+      count +=1;
+
+
+      // cut away after first space
+      seqan::strSplit(idSub, id, ' ', false, 1);
+
+
+      if (!info.noReadNames) {
+        seqan::append(string2hash, idSub[0]);
+        seqan::append(string2hash, "/1"); // to be consistent with BAM and FASTQ
+      }
+      seqan::append(string2hash, seq);
+
+      // Get MD5 hash
+      hex = str2md5(toCString(string2hash), length(string2hash));
+
+      if (info.debug) {
+        std::cout << string2hash << " " <<   std::hex << hex.p.low << "\n";
+      } else {
+        hexSum(hex, sum);
+      }
+
+      seqan::clear(string2hash);
+      seqan::clear(idSub);
+
+
+    }
+
+  }
+
+  if (!info.debug) {
+    std::cout << std::hex << sum << "\t";
+    std::cout << std::dec << count << "\n";
+  }
+
+    
+  return 0;
+}
+