diff --git a/base/DNA_string.h b/base/DNA_string.h
new file mode 100644
index 0000000..8f24d97
--- /dev/null
+++ b/base/DNA_string.h
@@ -0,0 +1,395 @@
+#pragma once
+
+#include "global.h"
+#include "nthash.hpp"
+
+namespace DNA{
+	//Map bases to 2 bit representation: A|a:00; C|c:01; G|g:10; T|t:11; N|n:arbitrary.
+	//Every 8 bits store 4 bases.
+	const uint8_t base_mask[4]={0xC0, 0x30, 0x0C, 0x03};
+	const uint8_t base2bits[]={0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x0F
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x1F
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x2F
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x3F
+													0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x4F
+													0x00, 0x00, 0x00, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x5F
+													0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x6F
+													0x00, 0x00, 0x00, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x7F
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x8F
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0x9F
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0xAF
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0xBF
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0xCF
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0xDF
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,//0xEF
+													0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};//0xFF
+	const char bits2base[]={'A', 'C', 'G', 'T', 'C', ' ', ' ', ' ', 'G', ' ', ' ', ' ', 'T', ' ', ' ', ' ',//0x0F
+													'C', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x1F
+													'G', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x2F
+													'T', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x3F
+													'C', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x4F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x5F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x6F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x7F
+													'G', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x8F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x9F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xAF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xBF
+													'T', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xCF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xDF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xEF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' '};//0xFF
+	const char bits2RCbase[]={'T', 'G', 'C', 'A', 'G', ' ', ' ', ' ', 'C', ' ', ' ', ' ', 'A', ' ', ' ', ' ',//0x0F
+													'G', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x1F
+													'C', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x2F
+													'A', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x3F
+													'G', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x4F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x5F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x6F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x7F
+													'C', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x8F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0x9F
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xAF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xBF
+													'A', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xCF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xDF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ',//0xEF
+													' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' '};//0xFF											
+	//const char bits2base[]={'A','C','G','T'};
+}
+
+class DNAString{
+public:
+	DNAString():_dna_base_num(0), _offset(0){}
+
+	DNAString(vector<uint8_t> d, size_t dna_base_num):_data(d.begin(), d.end()), _dna_base_num(dna_base_num){
+		this->_offset = (d.size()*4 - dna_base_num - 1)*2;
+	}
+
+	DNAString(const DNAString& dnaStr){
+		this->_data.assign(dnaStr._data.begin(), dnaStr._data.end());
+		this->_dna_base_num = dnaStr._dna_base_num;
+		this->_offset = dnaStr._offset;
+	}
+
+	DNAString& operator= (const DNAString& dnaStr){
+		this->_data.assign(dnaStr._data.begin(), dnaStr._data.end());
+		this->_dna_base_num = dnaStr._dna_base_num;
+		this->_offset = dnaStr._offset;
+		return *this;
+	}
+
+	DNAString& operator= (const string& str){
+		this->clear();
+		this->_dna_base_num = str.size();
+		size_t anchor, buf_anchor, buf_size;
+		if(_dna_base_num % 4){
+			buf_size = str.size()/4+1;
+		}else{
+			buf_size = str.size()/4;
+		} 
+		this->_offset = (buf_size*4-this->_dna_base_num-1)*2;
+		_data.resize(buf_size, 0x00);
+		anchor=buf_anchor=0;
+		while(anchor < str.size()){
+			uint8_t c=0x00;
+			for(int idx = 0; idx < 4 && anchor < str.size(); idx++, anchor++){
+				//cout<<str[anchor]<<":"<<DNA::base2bits[anchor]<<(6-idx*2)<<endl;
+				c |= (DNA::base2bits[str[anchor]]<<(6-idx*2));
+			}
+			_data[buf_anchor++] =c;
+		}
+		return *this;
+	}
+
+	const uint8_t* data() const{
+		return this->_data.data();
+	}
+
+	char operator[](size_t n){
+		assert(n < this->_dna_base_num);
+		return DNA::bits2base[this->_data[n/4] & DNA::base_mask[n%4]];
+	}
+
+	DNAString& clear(){
+		this->_data.clear();
+		this->_dna_base_num = 0;
+		this->_offset = 0;
+		return *this;
+	}
+
+	DNAString(const string& str):_dna_base_num(str.size()){
+		size_t anchor, buf_anchor, buf_size;
+		if(_dna_base_num % 4){
+			buf_size = str.size()/4+1;
+		}else{
+			buf_size = str.size()/4;
+		} 
+		this->_offset = (buf_size*4-this->_dna_base_num-1)*2;
+		_data.resize(buf_size, 0x00);
+		anchor=buf_anchor=0;
+		while(anchor < str.size()){
+			uint8_t c=0x00;
+			for(int idx = 0; idx < 4 && anchor < str.size(); idx++, anchor++){
+				//cout<<str[anchor]<<":"<<DNA::base2bits[anchor]<<(6-idx*2)<<endl;
+				c |= (DNA::base2bits[str[anchor]]<<(6-idx*2));
+			}
+			_data[buf_anchor++] =c;
+		}
+	}
+
+	DNAString make_DNAString(const string& str){
+		DNAString tmp;
+		tmp._dna_base_num = str.size();
+		
+		size_t anchor, buf_anchor, buf_size;
+		if(tmp._dna_base_num % 4){
+			buf_size = str.size()/4+1;
+		}else{
+			buf_size = str.size()/4;
+		} 
+		tmp._offset = (buf_size*4 - tmp._dna_base_num-1)*2;
+		tmp._data.resize(buf_size, 0x00);
+		anchor=buf_anchor=0;
+		while(anchor < str.size()){
+			uint8_t c=0x00;
+			for(int idx = 0; idx < 4 && anchor < str.size(); idx++, anchor++){
+				//cout<<str[anchor]<<":"<<DNA::base2bits[anchor]<<(6-idx*2)<<endl;
+				c |= (DNA::base2bits[str[anchor]]<<(6-idx*2));
+			}
+			tmp._data[buf_anchor++] =c;
+		}
+		return tmp;
+	}
+
+	DNAString& append(const char& base){
+		if(this->_offset == -2){
+			_data.push_back(DNA::base2bits[base]<<6);
+			this->_offset = 4;
+		}else{
+			_data.back() = (_data.back() | (DNA::base2bits[base]<<this->_offset));
+			this->_offset -= 2;
+		}
+		_dna_base_num++;
+		return *this;
+	}
+
+	DNAString& operator +=(const char& base){
+		this->append(base);
+		return *this;
+	}
+
+	size_t dna_base_num() const{
+		return _dna_base_num;
+	}
+
+	size_t length() const{
+		return _dna_base_num;
+	}
+
+	size_t data_size() const{
+		return this->_data.size();
+	}
+
+	DNAString& RC(){
+		vector<uint8_t> _rc_data(_data.size(), 0x00);
+		size_t _rc_slot_idx; _rc_slot_idx = 0;
+		size_t base_2_process = this->_dna_base_num;
+		size_t backward_slot_idx = (base_2_process-1)/4;
+		int backward_offset = (base_2_process % 4 + 3)%4;
+
+		while(base_2_process > 0){
+			uint8_t tmp=0x00;
+			for(int x = 0; x< 4 && base_2_process > 0; x++, base_2_process--){
+				if(backward_offset == -1){
+					backward_slot_idx--;
+					backward_offset = 3;
+				}
+				uint8_t extract = ((~(this->_data[backward_slot_idx])) & DNA::base_mask[backward_offset]);
+				
+				if(backward_offset == x){
+					tmp |= extract;
+				}else if(backward_offset < x){
+					tmp |= (extract >> ((x - backward_offset)*2));
+				}else{
+					tmp |= (extract << ((backward_offset - x)*2));
+				}
+				backward_offset--;
+			}
+			_rc_data[_rc_slot_idx++] = tmp;
+		}
+		this->_data.assign(_rc_data.begin(), _rc_data.end());
+
+		return *this;
+	}
+
+	DNAString get_RC(){
+		DNAString tmp(*this);
+		tmp.RC();
+		return tmp;
+	}
+
+	//pos is 0-based coordinate
+	void replace(const size_t pos, const char c){
+		assert(pos < this->_dna_base_num);
+		uint8_t d = this->_data[pos/4];
+		d &= ~(DNA::base_mask[pos%4]);
+		d |= (DNA::base2bits[c]<<((3 - pos%4)*2));
+		this->_data[pos/4] = d;
+	}
+
+	//start_pos is 0-based coordinate
+	DNAString substr(size_t start_pos, size_t len=string::npos) const{
+		assert(start_pos < this->_dna_base_num);
+
+		if(len == string::npos){
+			len = this->_dna_base_num - start_pos;
+		}else if(start_pos + len > this->_dna_base_num){
+			len = this->_dna_base_num - start_pos;
+		}
+
+		size_t byte_num = len/4;
+		if(len%4){
+			byte_num++;
+		} 
+
+		vector<uint8_t> data(byte_num, 0x00);
+		size_t data_slot_idx = 0;
+		size_t _data_slot_idx = start_pos/4;
+		size_t _data_slot_offset = start_pos%4;
+
+		if(start_pos==0){
+			while(data_slot_idx <= (len-1)/4 ){
+				data[data_slot_idx] = this->_data[data_slot_idx];
+				data_slot_idx++;
+			}
+		}else{
+			while(_data_slot_idx <= (start_pos+len-1)/4){
+				uint8_t slot = (this->_data[_data_slot_idx]<<(_data_slot_offset*2));
+				if((_data_slot_idx+1) <= (start_pos+len-1)/4 && _data_slot_idx < this->_data.size()){
+					slot |= (this->_data[_data_slot_idx+1] >> (8 - _data_slot_offset*2));
+				}
+				data[data_slot_idx++] = slot;
+				_data_slot_idx++;
+			}
+		}
+		uint8_t mask = 0x00;
+		for(int x = 0; x<len%4; x++){
+			mask |= DNA::base_mask[x];
+		}
+		data[data_slot_idx-1] &= mask;//reset the rest of the bits
+
+		return DNAString(data, len);
+	}
+
+	string to_str() const{
+		string str(this->_dna_base_num, ' ');
+		size_t base_idx=0;
+		for(size_t slot=0; slot < this->_data.size(); slot++){
+			for(size_t x=0; x < 4 && base_idx < _dna_base_num ; x++, base_idx++){
+				str[base_idx]=DNA::bits2base[(this->_data[slot]) & DNA::base_mask[x]];
+			}
+		}
+		return str;
+	}
+
+	bool is_palindrome() const{
+		if(this->_dna_base_num % 2){
+			return false;
+		}
+		return is_hairpin(this->_dna_base_num/2);
+	}
+
+	bool is_hairpin(size_t len=0) const{
+		if(len==0){
+			len = this->_dna_base_num / 2;
+		}else{
+			assert(len <= this->_dna_base_num / 2);
+		}
+		size_t forward_block_idx, backward_block_idx;
+		size_t forward_offset, backward_offset;//index(0-3) of slot to compare inside each data block.
+		forward_block_idx = 0; backward_block_idx = this->_data.size()-1;
+		forward_offset = 0; backward_offset = (this->_dna_base_num-1) % 4;
+		for(int x = 0; x < len; x++){
+			// cout<<x<<":"<<endl;
+			// cout<<forward_block_idx<<":"<<forward_offset<<":"<<DNA::bits2base[DNA::base_mask[forward_offset] & this->_data[forward_block_idx]]<<endl;
+			// cout<<backward_block_idx<<":"<<backward_offset<<":"<<DNA::bits2RCbase[DNA::base_mask[backward_offset] & this->_data[backward_block_idx]]<<endl;
+			if(DNA::bits2base[DNA::base_mask[forward_offset] & this->_data[forward_block_idx]] != DNA::bits2RCbase[DNA::base_mask[backward_offset] & this->_data[backward_block_idx]]){
+				return false;
+			}
+			if(forward_offset == 3){
+				forward_block_idx++;
+				forward_offset=0;	
+			}else{
+				forward_offset++;
+			}
+			if(backward_offset == 0){
+				backward_block_idx--;
+				backward_offset=3;
+			}else{
+				backward_offset--;
+			}
+		}
+		return true;
+	}
+
+	friend bool operator== (const DNAString& lhs, const DNAString& rhs);
+	friend bool operator< (const DNAString& lhs, const DNAString& rhs);
+	friend ostream& operator<<( ostream &output, const DNAString& dnaStr);
+
+private:
+	vector<uint8_t> _data;
+	size_t _dna_base_num;
+	int8_t _offset;//Number of bits shift to the left in order to store the next base. 
+};
+
+bool operator== (const DNAString& lhs, const DNAString& rhs){
+	//cout<<"called"<<endl;
+	if(lhs._data == rhs._data && lhs._dna_base_num == rhs._dna_base_num){
+		return true;
+	}
+	return false;
+}
+
+bool operator< (const DNAString& lhs, const DNAString& rhs){
+	if(lhs._data < rhs._data){
+		return true;
+	}
+	return false;
+}
+
+std::ostream& operator<<(std::ostream &output, const DNAString& dnaStr){
+	// output<<dnaStr._data.size()<<endl;
+	// for(const auto ele:dnaStr._data){
+	// 	output<<ele<<endl;
+	// }
+	output<<dnaStr.to_str();
+	return output;
+}
+
+DNAString operator + (const DNAString& str, const char& base){
+	DNAString tmp(str);
+	tmp.append(base);
+	return tmp;
+}
+DNAString operator + (const char& base, const DNAString& str){
+	DNAString tmp(string(1, base) + str.to_str());
+	return tmp;
+}
+
+//canonical nthash
+inline uint64_t NTPC64(const DNAString& dnaStr, const unsigned k, uint64_t& fhVal, uint64_t& rhVal) {
+	fhVal=0, rhVal=0;
+	const uint8_t* data = dnaStr.data();
+	size_t slot = 0;
+  for(unsigned i=0; i<dnaStr.length(); ) {
+  	for(int x = 0; x<4 && i < dnaStr.length(); x++, i++){
+  		char c = DNA::bits2base[data[slot]&DNA::base_mask[x]];
+    	fhVal ^= msTab[(unsigned char)c][(k-1-i)%64];
+    	rhVal ^= msTab[(unsigned char)c & cpOff][i%64];
+  	}
+  	slot++;
+  }
+  return (rhVal<fhVal)? rhVal : fhVal;
+}
+
diff --git a/base/Hash.h b/base/Hash.h
index 90e1883..c80bc2f 100644
--- a/base/Hash.h
+++ b/base/Hash.h
@@ -17,7 +17,7 @@
 // 2. It will not produce the same results on little-endian and big-endian
 //    machines.    */
 
-uint32_t MurmurHash2 ( const void * key, int len, uint32_t seed )
+uint32_t MurmurHash2 ( const void * key, int len, uint32_t seed=0 )
 {
   /* 'm' and 'r' are mixing constants generated offline.
      They're not really 'magic', they just happen to work well.  */
diff --git a/base/Utility.h b/base/Utility.h
index 41268d5..a99802f 100644
--- a/base/Utility.h
+++ b/base/Utility.h
@@ -5,6 +5,7 @@
 
 #include "global.h"
 #include "Hash.h"
+#include "DNA_string.h"
 
 using namespace std;
 
@@ -14,6 +15,21 @@ typedef vector<double> vec_double;
 typedef vector<vector<int>> matrix_int;
 typedef vector<int> vec_int;
 
+
+struct Contig
+{
+  DNAString seq;
+  double median_abundance;//median abundances of k-mers in
+  Contig(string s="", double a=0){
+    seq = s;
+    median_abundance = a;
+  }
+  Contig(DNAString s, double a=0){
+    seq = s;
+    median_abundance = a;
+  }
+};
+
 /*Compute median from a vector of int
 */
 double median(vector<int>& nums){
diff --git a/base/global.h b/base/global.h
index e60218f..1903319 100644
--- a/base/global.h
+++ b/base/global.h
@@ -47,6 +47,7 @@ using std::cin;
 using std::cout;
 using std::endl;
 using std::string;
+using std::ostream;
 using std::ifstream;
 using std::ofstream;
 using std::fstream;
@@ -85,15 +86,6 @@ const char DNA_bases[4] = {'A', 'C', 'G','T'};
 
 const std::string currentDateTime();
 
-struct Contig
-{
-  string seq;
-  double median_abundance;//median abundances of k-mers in
-  Contig(string s="", double a=0){
-    seq = s;
-    median_abundance = a;
-  }
-};
 
 #if 0
 std::ostream& operator<<( std::ostream& dest, __int128_t value ){
diff --git a/src/compare_contigs.cpp b/src/compare_contigs.cpp
index 044ed14..c384311 100644
--- a/src/compare_contigs.cpp
+++ b/src/compare_contigs.cpp
@@ -58,9 +58,12 @@ int main(int argc, char* argv[]){
     }
     seqs2.push_back(seq);
   }
-
   fin1.close();
   fin2.close();
+
+  //compare seqs in two sets by considering only exact match 
+
+
   fout.close();
 
   return 0;
diff --git a/src/contig_assembly.cpp b/src/contig_assembly.cpp
index ac32245..ebdd248 100644
--- a/src/contig_assembly.cpp
+++ b/src/contig_assembly.cpp
@@ -6,11 +6,11 @@
 #include "base/vector_mt.h"
 #include "base/Params.h"
 #include "base/Hash.h"
+#include "base/DNA_string.h"
 //#include "base/nthash.h"
 #include <boost/iostreams/filter/gzip.hpp>
 #include <boost/iostreams/copy.hpp>
 #include <boost/iostreams/filtering_streambuf.hpp>
-
 #include <tbb/concurrent_unordered_set.h>
 #include <tbb/concurrent_vector.h>
 #include <tbb/concurrent_queue.h>
@@ -19,23 +19,56 @@
 using tbb::concurrent_vector;
 using tbb::concurrent_queue;
 
-struct str_hasher
-{
-    size_t operator()(const string& val)const
-    {
-        return MurmurHash2(val);
-    }
-};
-
-struct str_equalto 
-{
-    bool operator()(const string& one, const string& two) const
+// struct DNAString_hasher
+// {
+//     size_t operator()(const DNAString& val)const //A is 00
+//     {
+//         return MurmurHash2(val.data(), val.length());
+//     }
+// };
+
+namespace std{
+  template<>
+  struct hash<DNAString>
+  {
+     size_t operator () (const DNAString& x) const
+     {
+        //cout<<"hash being called"<<endl;
+        return MurmurHash2(x.data(), x.data_size());
+     }
+  };
+
+  template<>
+  struct equal_to<DNAString>
+  {
+    bool operator()(const DNAString &lhs, const DNAString &rhs) const 
     {
-        return (one == two);
+        //cout<<"equal being called"<<endl;
+        return lhs == rhs;
     }
-};
+  };
+}
 
-typedef tbb::concurrent_unordered_set<string, str_hasher, str_equalto> unordered_set_mt;
+namespace tbb{
+  template<>
+  struct tbb_hash<DNAString>
+  { 
+    size_t operator () (const DNAString& x) const
+     {
+        return MurmurHash2(x.data(), x.data_size());
+     }
+  };
+}
+// struct DNAString_equalto 
+// {
+//     bool operator()(const DNAString& one, const DNAString& two) const
+//     {
+//         return (one == two);
+//     }
+// };
+
+typedef tbb::concurrent_unordered_set<DNAString> unordered_set_mt;
+//typedef tbb::concurrent_unordered_set<DNAString, DNAString_hasher, DNAString_equalto> unordered_set_mt;
 //using namespace boost::program_options;
 //namespace po = boost::program_options;
 //using boost::program_options::variables_map;
@@ -49,7 +82,7 @@ Params get_opts(int argc, char* argv[]){
     ("input,i", po::value<string>()->required(), "a file containing list of input file name(s), should be absolute address or file names when in the running directory.")
     ("format,f", po::value<char>()->default_value('f'), "format of the input: g(gzip); b(bzip2); f(plain fastq)")
     ("cqf,c", po::value<string>()->required(), "the counting quotient filter built with the same 'k'")
-    ("abundance_min,s", po::value<int>()->default_value(2), "minimum coverage of k-mers used to extend the assembly") 
+    ("abundance_min,s", po::value<int>()->default_value(1), "minimum coverage of k-mers used to extend the assembly") 
     ("solid_abundance_min,x", po::value<int>()->default_value(2), "minimum coverage of a solid k-mer to start the assembly")
     ("solid_abundance_max,X", po::value<int>()->default_value(100), "maximum coverage of a solid k-mer to start the assembly")
     (",t", po::value<int>()->default_value(16), "number of threads")
@@ -117,6 +150,51 @@ void get_unitig_forward(CQF_mt& cqf, const Params& options, concurrent_vector<Co
 //void get_unitig_backward(CQF_mt& cqf, const Params& options, vector_mt<Contig>& contigs, unordered_map_mt<string, int>& startKmer2unitig, WorkQueue* work_queue, int contig_id);
 
 int main(int argc, char* argv[]){
+  // /*test of DNA string 
+  // */
+  // string tmp_s("ATGCAGGATGCAT");
+  // DNAString dna_seq = tmp_s;
+  // dna_seq.append('T');
+  // cout<<dna_seq<<endl;
+  // cout<<dna_seq.substr(0, 5)<<endl;
+  // cout<<dna_seq.substr(7,5)<<endl;
+  // if(dna_seq.substr(0,5) == dna_seq.substr(7,5)){
+  //   cout<<"Equal operator works"<<endl;
+  // }
+  // dna_seq.RC();
+  // cout<<dna_seq<<endl;
+  // //cout<<dna_seq<<endl;
+  // //cout<<dna_seq.dna_base_num()<<endl;
+  
+  // DNAString dna_seq = string("AAATT");
+  // cout<<dna_seq.is_palindrome()<<endl;
+  // cout<<dna_seq.is_hairpin()<<endl;
+  
+  // unordered_map<DNAString, int> kmerSet;
+  
+  // DNAString dna_seq = string("AAATTT");
+  
+  // char DNA_bases[4]={'A','C','G','T'};
+  // DNAString fix = dna_seq.substr(0, 5);
+  // DNAString dna_seq1;
+  // for(int x = 0; x < 4; x++){
+  //   dna_seq1 = dna_seq.substr(0,5).append(DNA_bases[x]);
+  //   kmerSet[dna_seq1]=0;
+  //   cout<<"Insert "<<dna_seq1<<endl;
+  // }
+
+  // for(int x = 0; x < 4; x++){
+  //   dna_seq1 = fix;
+  //   dna_seq1.append(DNA_bases[x]);
+  //   auto it = kmerSet.find(dna_seq1);
+  //   if(it != kmerSet.end()){
+  //     cout<<dna_seq1<<" found."<<endl;
+  //     cout<<it->first<<": "<<it->second<<endl;
+  //   }
+  // }
+  
+  // return 0;
+
   Params options = get_opts(argc, argv);
   
   vector<string> seqFileNames;
@@ -132,7 +210,7 @@ int main(int argc, char* argv[]){
   seqFile_batch seqFiles(seqFileNames, ftype, options.fmode);   
 
   DisplayCurrentDateTime(); 
-  cout<<"[CQF] loading cqf from disk"<<endl;
+  cout<<"[CQF] load cqf from disk"<<endl;
   CQF_mt cqf_mt;
   cqf_mt.load(options.cqfFile);
   cout<<"[CQF] cqf loaded!"<<endl;
@@ -151,25 +229,29 @@ int main(int argc, char* argv[]){
   // return 0;
 
   DisplayCurrentDateTime();
-  cout<<"[Unitig] constructing unitigs"<<endl;
+  cout<<"[Unitig] find unitigs"<<endl<<std::flush;
   concurrent_vector<Contig> contigs;
   contigs.resize(1);
   find_unitigs_mt_master(cqf_mt, seqFiles, options, contigs);
   
   //bulid the graph
+  cout<<"[Unitig] build compacted DBG graph with unitigs as nodes."<<endl<<std::flush;
   uint64_t totalUnitigs_len, duplicateUnitigs_len; totalUnitigs_len = duplicateUnitigs_len=0;
   int duplicateUnitig_num=0;
   auto contig_num = contigs.size();
-  unordered_map<string, int> startKmer2unitig(contig_num*2);
+  unordered_map<DNAString, int, std::hash<DNAString>> startKmer2unitig;
   vector<bool> isDuplicate(contig_num, false);
+  vector<bool> isHairpin(contig_num, false);//whether the non-duplicate unitigs are hairpin(Palindrome excluded).
   vector<int> id2id_afterRemoveDuplicate(contig_num, -1);
 
-  string first_kmer, last_kmer_RC;
+  DNAString first_kmer, last_kmer_RC;
   int noduplicate_contig_num=0;
+  int palindrome_contig_num, hairpin_contig_num;
+  palindrome_contig_num = hairpin_contig_num = 0;
   for(int contig_id = 1; contig_id<contigs.size(); contig_id++){
     totalUnitigs_len += contigs[contig_id].seq.length();
     first_kmer = contigs[contig_id].seq.substr(0,options.K);
-    last_kmer_RC = RC_DNA(contigs[contig_id].seq.substr(contigs[contig_id].seq.length()-options.K));
+    last_kmer_RC = contigs[contig_id].seq.substr(contigs[contig_id].seq.length()-options.K).RC();
     auto it = startKmer2unitig.find(first_kmer);
     if(it != startKmer2unitig.end()){
       // auto tmp = abs(it->second);
@@ -192,66 +274,118 @@ int main(int argc, char* argv[]){
       noduplicate_contig_num++;
       id2id_afterRemoveDuplicate[contig_id]=noduplicate_contig_num;
       if(first_kmer == last_kmer_RC){
+        if(contigs[contig_id].seq.is_palindrome()){
+          //cout<<"[Warning] contig "<<contig_id<<"("<<contigs[contig_id].seq.length()<<" bp) is a palindrome seq (+ and - strand): "<<contigs[contig_id].seq<<endl;
+          palindrome_contig_num++;
+          startKmer2unitig[first_kmer] = noduplicate_contig_num;
+        }else{//Hairpin is the seq where >= K bases at both ends are complementary, but no palindrome.
+          //cout<<"[Warning] contig "<<contig_id<<"("<<contigs[contig_id].seq.length()<<" bp) is a hairpin seq (+ and - strand): "<<contigs[contig_id].seq<<endl;
+          isHairpin[noduplicate_contig_num] = true;
+          hairpin_contig_num++;
+          startKmer2unitig[first_kmer] = noduplicate_contig_num;
+        }
         //dont't know how to do  
-        cout<<"[Warning] contig "<<contig_id<<" seems like a palindrome seq (+ and - strand): "<<contigs[contig_id].seq<<endl;
+        //cout<<"[Warning] contig "<<contig_id<<" seems like a palindrome seq (+ and - strand): "<<contigs[contig_id].seq<<endl;
+      }else{
+        startKmer2unitig[first_kmer] = noduplicate_contig_num;
+        startKmer2unitig[last_kmer_RC] = -noduplicate_contig_num;
       }
-      startKmer2unitig[first_kmer] = noduplicate_contig_num;
-      startKmer2unitig[last_kmer_RC] = -noduplicate_contig_num;
     }
+
+    // if(startKmer2unitig.find(first_kmer) == startKmer2unitig.end()){
+    //   cerr<<"[error] why not found!"<<endl;
+    // }
+    // if(startKmer2unitig.find(last_kmer_RC) == startKmer2unitig.end()){
+    //   cerr<<"[error] why not found!"<<endl;
+    // }
   }
   cout<<"[Unitig] "<<contigs.size()<<" unitigs reported of length "<<totalUnitigs_len<<" bp in total"<<endl;
   cout<<"[Unitig] Among them, "<<duplicateUnitig_num<<" duplicate unitigs found of length "<<duplicateUnitigs_len<<" bp in total."<<endl;
   cout<<"[Unitig] After removing duplicates, we have "<<noduplicate_contig_num<<" unitigs of length "<<totalUnitigs_len- duplicateUnitigs_len<<" bp in total."<<endl;
+  cout<<"[Unitig] Among all non-duplicate unitigs, there are "<<palindrome_contig_num<<" palindromes among "<<hairpin_contig_num+palindrome_contig_num<<" hairpins."<<endl;
+  
   //check whether can reconstruct the full sequences from the graph.
-  // string refFile = "/public/hshi/tools/minia/test2/read50x_ref10K_e001.contigs.fa";
-  // ifstream refFin;
-  // refFin.open(refFile, ios::in);
-  // string refSeq="";
-  // getline(refFin, line);
-  // while(getline(refFin, line)){
-  //   refSeq += line;
+  // if(true){
+  //   string refFile = "/public/hshi/tools/SH-assembly/test/case1/genome10K.fasta";
+  //   ifstream refFin;
+  //   refFin.open(refFile, ios::in);
+  //   string refSeq="";
+  //   getline(refFin, line);
+  //   while(getline(refFin, line)){
+  //     refSeq += line;
+  //   }
+  //   to_upper_DNA(refSeq);
+  //   cout<<"[Test] length of reference sequence: "<<refSeq.length()<<" bp."<<endl;
+  //   int constructed_len=0;
+  //   DNAString kmer=refSeq.substr(0,options.K);
+  //   // while(true){
+  //   //   auto it = startKmer2unitig.find(kmer);
+  //   //   if(it == startKmer2unitig.end()){
+  //   //     cout<<"[Test] kmer not found: "<<kmer<<endl;  
+  //   //     if(constructed_len != 0){
+  //   //       cout<<"[Test] last kmer in the traversed sequence: "<<refSeq.substr(constructed_len-1, options.K)<<endl;
+  //   //       cout<<"[Test] traversal of reference seq broken at "<<constructed_len+options.K-1<<endl;
+  //   //     }
+  //   //     break;
+  //   //   }
+  //   //   //get the real ID
+  //   //   auto id = abs(it->second);
+  //   //   while(id2id_afterRemoveDuplicate[id] < abs(it->second)){
+  //   //     id++;
+  //   //   }
+  //   //   constructed_len += (contigs[id].seq.length() - options.K +1);
+  //   //   cout<<"[Test] find kmer: "<<kmer<<" in contig "<<it->second<<" of length "<<contigs[id].seq.length()<<" bp."<<endl;
+  //   //   cout<<"[Test] contig "<< it->second<<": ";
+  //   //   if(it->second < 0)
+  //   //     cout<<contigs[id].seq.get_RC();
+  //   //   else
+  //   //     cout<<contigs[id].seq;
+  //   //   cout<<" median_abundance: "<<contigs[id].median_abundance<<endl;
+  //   //   kmer = refSeq.substr(constructed_len, options.K);
+  //   //   if(constructed_len+options.K-1 >= refSeq.length()){
+  //   //     cout<<"[Test] reference seq is encoded in the DBG of unitigs"<<endl;
+  //   //     break;
+  //   //   }
+  //   // }
+
+  //   //check whether k-mer are in the CQF with proper coverage
+  //   kmer = refSeq.substr(constructed_len, options.K);
+  //   uint64_t kmer_hash, kmer_RC_hash, kmer_count;
+  //   kmer_hash = NTPC64(kmer, options.K, kmer_hash, kmer_RC_hash);
+  //   kmer_count = cqf_mt.count(kmer_hash%cqf_mt.qf->metadata->range);
+  //   cout<<"[Test] count: "<<kmer_count<<" of kmer: "<<kmer<<endl;
+  // }
+  
+  //output all kmers
+  // for(auto ele : startKmer2unitig){
+  //   cout<<ele.first<<":"<<ele.second<<endl;
   // }
-  // to_upper_DNA(refSeq);
-  // cout<<"[Test] length of reference sequence: "<<refSeq.length()<<" bp."<<endl;
-  // int constructed_len=0;
-  // string kmer=refSeq.substr(0,options.K);
-  // while(true){
-  //   auto it = startKmer2unitig.find(kmer);
-  //   if(it == startKmer2unitig.end()){
-  //     cout<<"[Test] kmer not found: "<<kmer<<endl;
-  //     cout<<"[Test] last kmer in the traversed sequence: "<<refSeq.substr(constructed_len-1, options.K)<<endl;
-  //     cout<<"[Test] traversal of reference seq broken at "<<constructed_len+options.K-1<<endl;
-  //     break;
+  // for(int contig_id = 1; contig_id < 2; contig_id++){
+  //   cout<<"Contig "<<contig_id<<":"<<contigs[contig_id].seq<<endl;
+  //   if(isDuplicate[contig_id]){
+  //     continue;
   //   }
-  //   //get the real ID
-  //   auto id = abs(it->second);
-  //   while(id2id_afterRemoveDuplicate[id] < abs(it->second)){
-  //     id++;
+  //   DNAString first_kmer = contigs[contig_id].seq.substr(0,options.K);
+  //   DNAString last_kmer_RC = contigs[contig_id].seq.substr(contigs[contig_id].seq.length()-options.K).RC();
+  //   if(startKmer2unitig.find(first_kmer) == startKmer2unitig.end()){
+  //     cerr<<"[error] why not found: "<<first_kmer<<endl;
+  //     for(auto ele : startKmer2unitig){
+  //       if(ele.first == first_kmer){
+  //         cerr<<"[error] Hmmm..."<<endl;
+  //       }
+  //     }
+  //   }else{
+  //     cerr<<"[error] "<<startKmer2unitig.find(first_kmer)->first<<endl;
   //   }
-  //   constructed_len += (contigs[id].seq.length() - options.K +1);
-  //   cout<<"[Test] find kmer: "<<kmer<<" in contig "<<it->second<<" of length "<<contigs[id].seq.length()<<" bp."<<endl;
-  //   cout<<"[Test] contig "<< it->second<<": ";
-  //   if(it->second < 0)
-  //     cout<<RC_DNA(contigs[id].seq);
-  //   else
-  //     cout<<contigs[id].seq;
-  //   cout<<" median_abundance: "<<contigs[id].median_abundance<<endl;
-  //   kmer = refSeq.substr(constructed_len, options.K);
-  //   if(constructed_len+options.K-1 >= refSeq.length()){
-  //     cout<<"[Test] reference seq is encoded in the DBG of unitigs"<<endl;
-  //     break;
+  //   if(startKmer2unitig.find(last_kmer_RC) == startKmer2unitig.end()){
+  //     cerr<<"[error] why not found!"<<endl;
   //   }
   // }
-  // //check whether k-mer are in the CQF with proper coverage
-  // kmer = refSeq.substr(constructed_len, options.K);
-  // uint64_t kmer_hash, kmer_RC_hash, kmer_count;
-  // kmer_hash = NTPC64(kmer.c_str(), options.K, kmer_hash, kmer_RC_hash);
-  // kmer_count = cqf_mt.count(kmer_hash%cqf_mt.qf->metadata->range);
-  // cout<<"[Test] count: "<<kmer_count<<" of kmer: "<<kmer<<endl;
 
   //output graph
   ofstream fout;
-  string header, seq, kmer_fix, kmer;
+  string header, seq;
+  DNAString kmer_fix, kmer;
   char bases[4]={'A','C','G','T'};
   fout.open(options.output, ios::out);
   for(int contig_id = 1; contig_id<contigs.size(); contig_id++){
@@ -259,14 +393,18 @@ int main(int argc, char* argv[]){
       continue;
     }
     fout<<">"<<(id2id_afterRemoveDuplicate[contig_id]-1)<<" LN:i:"<<contigs[contig_id].seq.length()<<" KC:i:"<<contigs[contig_id].median_abundance*(contigs[contig_id].seq.length()-options.K+1)<<" km:f:"<<contigs[contig_id].median_abundance;
-    kmer_fix=contigs[contig_id].seq.substr(contigs[contig_id].seq.length()-options.K+1);
+    kmer_fix=contigs[contig_id].seq.substr(contigs[contig_id].seq.length()-options.K+1, options.K-1);
+
     for(int x=0; x<4; x++){
-      kmer = kmer_fix+bases[x];
+      kmer = kmer_fix; kmer.append(bases[x]);
       auto it = startKmer2unitig.find(kmer);
       if(it != startKmer2unitig.end()){
         auto tmp = it->second;
         if(tmp>0){
           fout<<" L:+:"<<tmp-1<<":+";
+          if(isHairpin[tmp]){
+            fout<<" L:+:"<<tmp-1<<":-";
+          }
         }else{
           fout<<" L:+:"<<-tmp-1<<":-";
         }
@@ -274,14 +412,17 @@ int main(int argc, char* argv[]){
         continue;
       }
     }
-    kmer_fix=RC_DNA(contigs[contig_id].seq.substr(0,options.K-1));
+    kmer_fix=contigs[contig_id].seq.substr(0, options.K-1).RC();
     for(int x=0; x<4; x++){
-      kmer = kmer_fix+bases[x];
+      kmer = kmer_fix; kmer.append(bases[x]);
       auto it = startKmer2unitig.find(kmer);
       if(it != startKmer2unitig.end()){
         auto tmp = it->second;
         if(tmp>0){
           fout<<" L:-:"<<tmp-1<<":+";
+          if(isHairpin[tmp]){
+            fout<<" L:-:"<<tmp-1<<":-";
+          }
         }else{
           fout<<" L:-:"<<-tmp-1<<":-";
         }
@@ -604,7 +745,7 @@ void processDataChunk(CQF_mt& cqf, const Params& options, concurrent_vector<Cont
       startKmers.insert(kmer);//fake id 
 
       get_unitig_forward(cqf, options, contigs, startKmers, work_queue, iter);
-      iter->seq = RC_DNA(iter->seq);
+      iter->seq.RC();
       get_unitig_forward(cqf, options, contigs, startKmers, work_queue, iter);
     }
     //skip two lines
@@ -671,7 +812,7 @@ void find_unitigs_mt_master(CQF_mt& cqf, seqFile_batch& seqFiles, const Params&
         //work_queue->add_skip_work(1);
         auto iter = contigs.push_back(Contig(kmer, kmer_count));
         //startKmer2unitig.insert_mt(kmer, contig_id); //may not be the start k-mer
-        startKmers.insert(kmer);//fake id 
+        startKmers.insert(DNAString(kmer));//fake id 
         //startKmer2unitig.insert_mt(kmer_RC, -contig_id); //may not be the end k-mer
         // Contig master_contig(kmer, kmer_count);
         // get_unitig_forward(cqf, options, contigs, startKmer2unitig, work_queue, master_contig);
@@ -682,7 +823,7 @@ void find_unitigs_mt_master(CQF_mt& cqf, seqFile_batch& seqFiles, const Params&
         //break;
         get_unitig_forward(cqf, options, contigs, startKmers, work_queue, iter);
         //contigs.set(contig_id, Contig(RC_DNA(contigs[contig_id].seq), contigs[contig_id].median_abundance));
-        iter->seq = RC_DNA(iter->seq);
+        iter->seq.RC();
         get_unitig_forward(cqf, options, contigs, startKmers, work_queue, iter);
         //get_unitig_backward(cqf, options, contigs, startKmer2unitig, work_queue, contig_id);
 
@@ -718,6 +859,7 @@ void find_unitigs_mt_master(CQF_mt& cqf, seqFile_batch& seqFiles, const Params&
   prod_threads.join_all();
 
   //contigs.insert(contigs.end(), master_contigs.begin(), master_contigs.end());
+  startKmers.clear();
 }
 
 /*
@@ -1392,18 +1534,20 @@ void get_unitig_forward(CQF_mt& cqf, const Params& options, concurrent_vector<Co
   int candidates_before_num, candidates_after_num;
   int nodes_before_num, nodes_after_num;
   uint64_t kmer_hash, kmer_RC_hash, current_kmer_hash, current_kmer_RC_hash;
-  string kmer, kmer_RC, current_kmer, current_kmer_RC, current_kmer_fix;
+  DNAString kmer, kmer_RC, current_kmer, current_kmer_RC, current_kmer_fix;
   uint64_t kmer_count, kmer_RC_count;
   int idx;
   
-  string contig_seq = contigIter->seq;
+  DNAString& contig_seq = contigIter->seq;
+  //string contig_seq = contigIter->seq.to_str();
   current_kmer = contig_seq.substr(contig_seq.length()-K);
-  current_kmer_RC = RC_DNA(current_kmer);
+  current_kmer_RC = current_kmer.get_RC();
 
   std::vector<int> abundances;
   abundances.push_back(int(contigIter->median_abundance));
 
-  NTPC64(current_kmer.c_str(), K, current_kmer_hash, current_kmer_RC_hash);
+  //NTPC64(current_kmer.c_str(), K, current_kmer_hash, current_kmer_RC_hash);
+  NTPC64(current_kmer, K, current_kmer_hash, current_kmer_RC_hash);
   int node_after_x, node_before_x;//useful only when there is only one node after and without candidates after.
   while(true){
     //NTPC64(current_kmer.c_str(), K, current_kmer_hash, current_kmer_RC_hash);
@@ -1451,7 +1595,7 @@ void get_unitig_forward(CQF_mt& cqf, const Params& options, concurrent_vector<Co
       //kmer_hash = current_kmer_hash; kmer_RC_hash = current_kmer_RC_hash;
       //kmer_hash = NTPC64(current_kmer_RC[0], DNA_bases[x], K, kmer_hash, kmer_RC_hash);
       //kmer = current_kmer_RC;
-      kmer[K-1] = DNA_bases[x];
+      kmer.replace(K-1, DNA_bases[x]);
       kmer_hash = current_kmer_hash; kmer_RC_hash = current_kmer_RC_hash;
       kmer_hash = NTPC64('T', DNA_bases[x], options.K, kmer_RC_hash, kmer_hash);
       //kmer_hash = NTPC64(kmer.c_str(), K, kmer_hash, kmer_RC_hash);
@@ -1482,7 +1626,7 @@ void get_unitig_forward(CQF_mt& cqf, const Params& options, concurrent_vector<Co
       // }
       startKmers.insert(current_kmer_RC);
       //contigs.set_mt(contig_id, Contig(contig_seq, median(abundances)));
-      contigIter->seq = contig_seq;
+      //contigIter->seq = contig_seq;
       contigIter->median_abundance = median(abundances);
 
       for(int x= 0; x < 4; x++){
@@ -1498,7 +1642,7 @@ void get_unitig_forward(CQF_mt& cqf, const Params& options, concurrent_vector<Co
       kmer = current_kmer_RC;
       for(int x = 0; x < 4; x++){
         if(candidates_before[x]){
-          kmer[K-1] = DNA_bases[x];
+          kmer.replace(K-1, DNA_bases[x]);
           if(startKmers.find(kmer) == startKmers.end()){
             startKmers.insert(kmer);
             auto it = contigs.push_back(Contig(kmer, kmer_abundance_befores[x]));
@@ -1530,7 +1674,7 @@ void get_unitig_forward(CQF_mt& cqf, const Params& options, concurrent_vector<Co
       NTPC64('T', DNA_bases[node_after_x], options.K, current_kmer_hash, current_kmer_RC_hash);
     }else{ //stop
       startKmers.insert(current_kmer_RC);
-      contigIter->seq = contig_seq;
+      //contigIter->seq = contig_seq;
       contigIter->median_abundance = median(abundances);
       //contigs.set_mt(contig_id, Contig(contig_seq, median(abundances)));
       break;