diff -rwu pIRS_111_original/src/pirs/src/global.h pIRS_111/src/pirs/src/global.h
--- pIRS_111_original/src/pirs/src/global.h	2013-09-26 14:52:51.000000000 +0200
+++ pIRS_111/src/pirs/src/global.h	2015-09-30 15:39:16.000000000 +0200
@@ -4,6 +4,8 @@
 using namespace std;
 using namespace boost; 
 
+#include <vector>
+
 
 typedef struct{
   int Read_length;
@@ -17,7 +19,7 @@
   int Q_shift;
   int Mask_quality_mode;
   int Output_type;
-  double Coverage;
+  int NPairs;
   double Error_rate;
 	string Input_ref1;
   string Input_ref2;
@@ -25,6 +27,10 @@
   string GC_depth_profile;
   string InDel_error_profile;
   string Output_prefix;
+  vector<string> adapter_1;
+  vector<string> adapter_2;
+  vector<string> barcode_1;
+  vector<string> barcode_2;
 }PARAMETER;
 
 #endif
diff -rwu pIRS_111_original/src/pirs/src/simulate_Illumina_reads.cpp pIRS_111/src/pirs/src/simulate_Illumina_reads.cpp
--- pIRS_111_original/src/pirs/src/simulate_Illumina_reads.cpp	2013-09-26 14:52:52.000000000 +0200
+++ pIRS_111/src/pirs/src/simulate_Illumina_reads.cpp	2015-09-30 15:39:16.000000000 +0200
@@ -15,6 +15,30 @@
 
 using namespace std;
 
+
+//                                                               Barcode:
+//                                                               NNNNNN
+const std::string ADAPTER_1 = "AGATCGGAAGAGCACACGTCTGAACTCCAGTCACCACCTAATCTCGTATGCCGTCTTCTGCTTG";
+const std::string ADAPTER_2 = "AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT";
+
+
+void add_adapter(const std::string& adapter, int length, std::string& dst)
+{
+    if (dst.size() < length) {
+        dst.reserve(length);
+
+        for (int pos = 0; dst.size() < length; ++pos) {
+            if (pos < adapter.size()) {
+                dst.push_back(adapter.at(pos));
+            } else {
+                dst.push_back('A');
+            }
+        }
+    }
+}
+
+
+
 /*parameter variable:
 
   int Read_length;
@@ -28,7 +52,7 @@
   int Q_shift;
   int Mask_quality_mode;
   int Output_type;
-  double Coverage;
+  int NPairs;
   double Error_rate;
 	string Input_ref1;
   string Input_ref2;
@@ -37,7 +61,7 @@
   string InDel_error_profile;
   string Output_prefix;
  */
-PARAMETER InputParameter ={100,500,-1,1,0,1,1,1,64,0,1,5,-1,"","","","","","Illumina"};
+PARAMETER InputParameter ={100,500,-1,1,0,1,1,1,64,0,1,10000,-1,"","","","","","Illumina"};
 
 int Ref_Base_num = 0;  //ATCG: 4
 int Statistical_Cycle_num = 0; //the cycle number in Base-calling profile
@@ -94,13 +118,19 @@
 	cout<<"you can get another diploid genome sequence by the command \"pirs diploid\", but remember that heterozygosis SNP rate and heterozygosis Indel rate only exist in diploid. \n";
 	cout<<endl<<"Program: pirs simulate"<<endl;
 	cout<<endl<<"Usage:\t./pirs simulate [options]"<<endl;
+
+    cout<<"\t-1  <string>  Adapter sequence appended to mate 1 reads,default:"<<ADAPTER_1<<endl;
+    cout<<"\t-2  <string>  Adapter sequence appended to mate 2 reads,default:"<<ADAPTER_2<<endl;
+    cout<<"\t-8  <string>  Barcode sequence appended to mate 1 reads"<<endl;
+    cout<<"\t-9  <string>  Barcode sequence appended to mate 2 reads"<<endl;
+
 	cout<<"\t-i  <string>  input_ref1, input reference genome sequence *.fa/*.fa.gz, no default vaule"<<endl;
 	cout<<"\t-I  <string>  input_ref2, for diploid genome, input another reference genome sequence which was generated by command \"pirs diploid\""<<endl;
 	cout<<"\t-s  <string>  Base-calling profile, input Base-calling profile for simulating substitution-error and quality score,default: (exe_path)"<<BASE_CALLING_PROFILE<<endl;
 	cout<<"\t-d  <string>  GC content-coverage profile, input GC content-coverage file for simulating GC bias, the default profile are determined based on the twice of read length"<<endl;
 	cout<<"\t-b  <string>  InDel-error profile, input InDel-error profile for simulating InDel-error of reads, default:(exe_path)"<<INDEL_ERROR_PROFILE<<endl;
 	cout<<"\t-l  <int>     read_len, set length of read, read1 and read2 have the same length,default:"<<InputParameter.Read_length<<endl;
-	cout<<"\t-x  <double>  coverage, set the sequencing coverage(sometimes called depth),default:"<<InputParameter.Coverage<<endl;
+    cout<<"\t-x  <int>     number of reads-pairs to generate,default:"<<InputParameter.NPairs<<endl;
 	cout<<"\t-m  <int>     insertsize_mean, set the average value of insert size,default:"<<InputParameter.Insertsize_mean<<endl;
 	cout<<"\t-v  <int>     insertsize_sd, set the standard deviation of insert sizes, default:insertsize_mean/20"<<endl;
 	cout<<"\t-e  <double>  substitution-error rate, set the average substitution-error rate(0 or 0.0001~0.63) over all cycles, default=average substitution-error rate of Base-calling profile"<<endl;
@@ -117,8 +147,8 @@
 	cout<<endl<<"Example:"<<endl;
 	cout<<"\t1. ./pirs simulate -i ref_sequence.fa"<<endl;
 	cout<<"\t  Every parameter use the default value."<<endl;
-	cout<<"\t2. ./pirs simulate -i ref_sequence.fa -l 100 -x 20 -o human_500_100"<<endl;
-	cout<<"\t  Just set read length and coverage you needed."<<endl;
+    cout<<"\t2. ./pirs simulate -i ref_sequence.fa -l 100 -x 2000000 -o human_500_100"<<endl;
+    cout<<"\t  Just set read length and #reads you needed."<<endl;
 	cout<<"\t3. ./pirs simulate -i ref_sequence.fa -o human -m 600 -v 30"<<endl;
 	cout<<"\t  Set insertsize distribution."<<endl;
 	cout<<"\t4. ./pirs simulate -i ref_sequence.fa -I ref_seq.snp.indel.inversion.fa.gz -o human "<<endl;
@@ -137,16 +167,20 @@
 
 void SimReads_Getopt(int argc,char *argv[]){
 	int c;
-	while ((c=getopt(argc,argv,"i:I:s:d:b:l:x:m:v:e:a:f:g:q:M:Q:E:c:o:h"))!=-1)
+    while ((c=getopt(argc,argv,"1:2:8:9:i:I:s:d:b:l:x:m:v:e:a:f:g:q:M:Q:E:c:o:h"))!=-1)
 	{
 		switch(c){
+            case '1': InputParameter.adapter_1.push_back(optarg);break;
+            case '2': InputParameter.adapter_2.push_back(optarg);break;
+            case '8': InputParameter.barcode_1.push_back(optarg);break;
+            case '9': InputParameter.barcode_2.push_back(optarg);break;
 			case 'i': InputParameter.Input_ref1=optarg;break;
 			case 'I': InputParameter.Input_ref2=optarg;break;
 			case 's': InputParameter.BaseCalling_profile=optarg;break;
 			case 'd': InputParameter.GC_depth_profile=optarg;break;
 			case 'b': InputParameter.InDel_error_profile=optarg;break;
 			case 'l': InputParameter.Read_length=atoi(optarg);break;
-			case 'x': InputParameter.Coverage=atof(optarg);break;
+            case 'x': InputParameter.NPairs=atoi(optarg);break;
 			case 'm': InputParameter.Insertsize_mean=atoi(optarg);break;
 			case 'v': InputParameter.Insertsize_sd=atoi(optarg);break;
 			case 'e': InputParameter.Error_rate=atof(optarg);break;
@@ -185,8 +219,7 @@
 	//check parameter
 	if(InputParameter.Input_ref1 == ""){cerr<<"Error: there is not default value with option -i, please input reference sequence!"<<endl;exit(-1);}
 	if(InputParameter.Read_length <= 0){cerr<<"Error: read length should be set bigger than 0, please check option -l !"<<endl;exit(-1);}
-	if(InputParameter.Coverage <= 0){cerr<<"Error: coverage should be set bigger than 0, please check option -x !"<<endl;exit(-1);}
-	if(InputParameter.Insertsize_mean < InputParameter.Read_length){cerr<<"Error: insertize mean should be set bigger than read_length, please check option -m !"<<endl;exit(-1);}
+    if(InputParameter.NPairs <= 0){cerr<<"Error: NReads should be set bigger than 0, please check option -x !"<<endl;exit(-1);}
 	if(InputParameter.Insertsize_sd < 0){cerr<<"Error: insertsize_sd should be set bigger than 0, please check option -v !"<<endl;exit(-1);}
 	if(InputParameter.Error_rate != -1 && InputParameter.Error_rate != 0 && InputParameter.Error_rate < 0.0001 || InputParameter.Error_rate > 0.63 ){cerr<<"Error: error_rate should be set 0 or between 0.0001 and 0.63, you can also set -1 to simulate default error rate according with error profile, please check option -e !"<<endl;exit(-1);}
 	if(InputParameter.Is_cyclization != 0 && InputParameter.Is_cyclization != 1){cerr<<"Error: Is_cyclization should be set 0 or 1, please check option -f !"<<endl;exit(-1);}
@@ -197,6 +230,22 @@
 	if(InputParameter.Mask_quality_mode != 0 && InputParameter.Mask_quality_mode != 1 && InputParameter.Mask_quality_mode != 2){cerr<<"Error: Mask_quality_mode should be set 0 , 1 or 2, please check option -E !"<<endl;exit(-1);}
 	if(InputParameter.Output_type != 0 && InputParameter.Output_type != 1){cerr<<"Error: output_type should be set 0 or 1, please check option -c !"<<endl;exit(-1);}
 
+    if (InputParameter.adapter_1.empty() && InputParameter.adapter_2.empty()) {
+        InputParameter.adapter_1.push_back(ADAPTER_1);
+        InputParameter.adapter_2.push_back(ADAPTER_2);
+    }
+
+    if (InputParameter.barcode_1.empty() && InputParameter.barcode_2.empty()) {
+        InputParameter.barcode_1.push_back("");
+        InputParameter.barcode_2.push_back("");
+    }
+
+    if (InputParameter.adapter_1.size() != InputParameter.adapter_2.size()) {
+        cerr<<"Unequal number of adapter -1 and adapter -2 sequences"<<endl;exit(-1);
+    } else if (InputParameter.barcode_1.size() != InputParameter.barcode_2.size()) {
+        cerr<<"Unequal number of barcode -8 and barcode -9 sequences"<<endl;exit(-1);
+    }
+
 	//set the simulate cycle number
 	Simulate_Cycle_num = InputParameter.Read_length*2;
 	//initialize error position distribution table
@@ -423,7 +472,7 @@
 		<<"#simulate InDel-error in reads: "<<Is_Indel<<endl
 		<<"#InDel-error profile: "<<indel_error_profile<<endl
 		<<"#read length: "<<InputParameter.Read_length<<endl
-		<<"#data coverage: "<<InputParameter.Coverage<<endl
+        <<"#read pairs: "<<InputParameter.NPairs<<endl
 		<<"#substitution-error rate setting by user: "<<input_error_rate<<" (you can find the real substitution-error rate of output data in file *.error_rate.distr)"<<endl
 		<<"#mean of insertsize: "<<InputParameter.Insertsize_mean<<endl
 		<<"#standard deviation of insert sizes: "<<InputParameter.Insertsize_sd<<endl
@@ -439,7 +488,7 @@
 			<<"#mode of mask quality: "<<mode_of_Mask<<endl;
 	}
 
-	Infor_outfile<<endl<<"#readId\t-i/-I\tchr\tposition\t+/-\tinsertSize\tmaskEndLen\tsubstitution\tinsertion\tdeletion"<<endl;
+    Infor_outfile<<endl<<"#readId\t-i/-I\tchr\tposition\t+/-\tinsertSize\tmaskEndLen\tsubstitution\tinsertion\tdeletion\tadapter\tbarcode_1\tbarcode_2"<<endl;
 	
 	///////////////////////////get genome seq and start to simulate reads/////////////////////////////
 	
@@ -649,17 +698,9 @@
 	to_upper(sequence);
 
 	uint64_t sequence_length=sequence.size();
-	uint64_t reads_pair_num=0;
+    Total_read_pair = Total_read_pair + InputParameter.NPairs;
 	
-	if (InputParameter.Input_ref2 != "")
-	{
-		reads_pair_num=(uint64_t)(sequence_length*InputParameter.Coverage/InputParameter.Read_length/2/2);
-	}else{
-		reads_pair_num=(uint64_t)(sequence_length*InputParameter.Coverage/InputParameter.Read_length/2);
-	}
-	Total_read_pair = Total_read_pair + (uint64_t)(sequence_length*InputParameter.Coverage/InputParameter.Read_length/2);
-	
-	readonchr = simulate_reads(sequence,sequence_length,reads_pair_num, id, 1, read_genome);
+    readonchr = simulate_reads(sequence,sequence_length,InputParameter.NPairs, id, 1, read_genome);
 
 	//simulate reads of another diploid genome 
 	if (InputParameter.Input_ref2 != "") 
@@ -669,11 +710,10 @@
 		{
 			return 0;
 		}
-		reads_pair_num=(uint64_t)(sequence_length*InputParameter.Coverage/InputParameter.Read_length/2/2);
 
 		uint64_t readonchr2=read_genome+readonchr;
 		
-		readonchr += simulate_reads(sequence2,sequence_length,reads_pair_num, id, 2, readonchr2);
+        readonchr += simulate_reads(sequence2,sequence_length,InputParameter.NPairs, id, 2, readonchr2);
 	}
 	return readonchr;
 }
@@ -698,15 +738,11 @@
 	{
 		//simulate insertsize
 		int insertsize = simulate_insertsize(InputParameter.Insertsize_mean,InputParameter.Insertsize_sd);
-		if (insertsize<InputParameter.Read_length){pair_count--;continue;}
-		if (seqlen<insertsize){return reads_count;}
+        if (insertsize < 0){pair_count--;continue;}
 		uint64_t N = 1000000000000;
 		uint64_t pos = (uint64_t)(((uint64_t)((double)rand() / double(RAND_MAX) * N)) % seqlen);
 		if (pos+insertsize>seqlen){pair_count--;continue;}
 			
-		//get insert seq
-		string sub_str=seq.substr(pos,insertsize);
-
 		map<int,string,less<int> > indel1;
 		map<int,string,less<int> > indel2;
 		int r1_slen = 0;
@@ -714,27 +750,36 @@
 		if(InputParameter.Is_simulate_InDel)
 		{
 			get_reads_indel(InputParameter.Read_length, indel1, indel2, r1_slen, r2_slen, InDel_max_len, InDel_error_matrix, InDel_num);
-  		//fixed in v1.1.1
-  		if(InputParameter.Read_length-r1_slen > sub_str.size() || InputParameter.Read_length-r2_slen > sub_str.size())
-  		{
-  			pair_count--;continue;
-  		}
 		}
 		
 		string ref_read1, ref_read2;
 		int selection=int(rand()%2); //0 or 1, for selecting output file randomly and deciding read +/-
 		
 		int read1_pos, read2_pos;
+
+        const int barcode_idx = rand() / (RAND_MAX / InputParameter.barcode_1.size() + 1);
+        const std::string barcode_1 = InputParameter.barcode_1.at(barcode_idx);
+        const std::string barcode_2 = InputParameter.barcode_2.at(barcode_idx);
+
+        //get insert seq
 		if(selection == 0)
 		{
+            const string sub_str = barcode_1 + seq.substr(pos,insertsize) + reversecomplementary(barcode_2);
+            const string sub_str_rev(sub_str.rbegin(), sub_str.rend());
+
 			ref_read1 = sub_str.substr(0, InputParameter.Read_length-r1_slen);
-			ref_read2 = sub_str.substr(insertsize-InputParameter.Read_length+r2_slen, InputParameter.Read_length-r2_slen);
+            ref_read2 = sub_str_rev.substr(0, InputParameter.Read_length-r2_slen);
+            ref_read2 = std::string(ref_read2.rbegin(), ref_read2.rend());
 			read1_pos = pos+1;
-			read2_pos = pos+insertsize-InputParameter.Read_length-r2_slen+1;
+            read2_pos = pos+ref_read2.size()+1;
 		}else{
-			ref_read1 = sub_str.substr(insertsize-InputParameter.Read_length+r1_slen, InputParameter.Read_length-r1_slen);
+            const string sub_str = barcode_2 + seq.substr(pos,insertsize) + reversecomplementary(barcode_1);
+            const string sub_str_rev(sub_str.rbegin(), sub_str.rend());
+
+            ref_read1 = sub_str_rev.substr(0, InputParameter.Read_length-r1_slen);
+            ref_read1 = std::string(ref_read1.rbegin(), ref_read1.rend());
 			ref_read2 = sub_str.substr(0, InputParameter.Read_length-r2_slen);
-			read1_pos = pos+insertsize-InputParameter.Read_length-r1_slen+1;		
+            read1_pos = pos+ref_read1.size()+1;
 			read2_pos = pos+1;		
 		}
 		
@@ -747,8 +792,10 @@
 		//simulate GC bias
 		if(InputParameter.Is_simulate_GC_bias){
 			string check_seq = ref_read1+ref_read2;
+            if (!check_seq.empty()) {
 			if(simulate_GC_bias(check_seq,GC_bias_abundance)){pair_count--;continue;}
 		}
+		}
 		
 		//insertsize statistics
 		if(InsertSize_distr[insertsize]>0)
@@ -790,6 +837,10 @@
 			exit(-1);
 		}
 		
+        int adapter_idx = rand() / (RAND_MAX / InputParameter.adapter_1.size() + 1);
+        add_adapter(InputParameter.adapter_1.at(adapter_idx), InputParameter.Read_length - r1_slen, ref_read1);
+        add_adapter(InputParameter.adapter_2.at(adapter_idx), InputParameter.Read_length - r2_slen, ref_read2);
+
 		bool* Is_insertion_pos1;
 		bool* Is_insertion_pos2;
 		Is_insertion_pos1 = new bool[InputParameter.Read_length];
@@ -1234,7 +1285,9 @@
     	}
   	}
   	
-  	Infor_outfile<<endl;
+   Infor_outfile<<"\t"<<adapter_idx<<"\t"
+                << InputParameter.barcode_1.at(barcode_idx) << "\t"
+                << InputParameter.barcode_2.at(barcode_idx) << endl;
   	
   	//output read2 information
 		Infor_outfile<<id_header<<"read_"<<InputParameter.Insertsize_mean<<"_"<<reads_count+reads_all<<"/"<<2<<"\t"<<i_or_I<<"\t"<<id_seq<<"\t"<<read2_pos<<"\t"<<read2_order<<"\t"<<insertsize<<"\t"<<mask_num2<<"\t";
@@ -1299,7 +1352,9 @@
     	}
   	}
   	
-  	Infor_outfile<<endl;
+   Infor_outfile<<"\t"<<adapter_idx<<"\t"
+                << InputParameter.barcode_1.at(barcode_idx) << "\t"
+                << InputParameter.barcode_2.at(barcode_idx) << endl;
   	
   	//output reads
 		if(!InputParameter.Output_type){