不知道有什么好办法可以让primer允许漏配，现在仅仅是允许错配，还是有一些没有配上，454数据有些primer漏配了一些，下一步解决这个问题
1 #include <cstdio> 2 #include <cstdlib> 3 #include <string> 4 #include <iostream> 5 #include <fstream> 6 #include <iomanip> 7 #include <getopt.h> 8 using namespace std; 9 10 int max_mis = 1; 11 int N = 0; 12 string primer = ""; 13 string* barcode = NULL; 14 string* barcodeName = NULL; 15 void usage () 16 { 17 cout << " Usage: splitByBarcode [options] <read.fasta> <barcode.txt> <output prefix> " 18 << " -m <int> allow the max mis num, default " << max_mis << " " 19 << " -p <string> primer, required "<<endl 20 << " -n <int> barcode num, required"<<endl; 21 exit(1); 22 } 23 struct COMPARE{ 24 int ifCompared; 25 int posStart; 26 int posEnd; 27 }; 28 COMPARE compare(string reads, string barcode, int maxmis,int b,int e); 29 int main(int argc, char *argv[]) 30 { 31 if (argc < 4) 32 { 33 usage(); 34 } 35 int c; 36 while ((c=getopt(argc,argv,"m:n:p:")) != -1) 37 { 38 switch (c) 39 { 40 case 'm' : max_mis = atoi(optarg);break; 41 case 'n' : N = atoi(optarg); break; 42 case 'p' : primer = optarg;break; 43 default : cout<<"error:"<<(char)c<<endl;usage(); 44 } 45 } 46 if(N<=0){ 47 cerr<<"please input the required barcode num "; 48 exit(1); 49 } 50 if(primer == ""){ 51 cerr<<"please input the required primer "; 52 exit(1); 53 } 54 string seq_file = argv[optind++]; 55 string barcode_file = argv[optind++]; 56 string outPrefix = argv[optind++]; 57 string out_file = outPrefix + ".list"; 58 barcode = new string[N]; 59 barcodeName = new string[N]; 60 61 for(int i=0;i<N;i++){ 62 barcodeName[i] = ""; 63 barcode[i] = ""; 64 } 65 66 ifstream reads; 67 reads.open(seq_file.c_str()); 68 if(!reads){ 69 cerr<<"fail to open input file"<<seq_file<<endl; 70 } 71 ifstream barcode_in; 72 barcode_in.open(barcode_file.c_str()); 73 if(!barcode_in){ 74 cerr<<"fail to open barcode file" <<barcode_file<<endl; 75 } 76 77 string text; 78 int k=0; 79 while( getline( barcode_in, text, ' ') ) 80 { 81 int startPos = 0, endPos; 82 endPos = text.find(" "); 83 barcodeName[k] = text.substr(startPos, endPos-startPos+1); 84 startPos = endPos + 1; 85 endPos = text.find(" ",startPos); 86 barcode[k] = text.substr(startPos, endPos-startPos+1); 87 //cout << barcodeName[k] <<" "<< barcodeForward[k] <<" "<< barcodeReverse[k] << endl; 88 k++; 89 } 90 91 ofstream splitlist; 92 splitlist.open(out_file.c_str()); 93 if(!splitlist){ 94 cerr<<"fail to creat output file" << out_file <<endl; 95 } 96 int line_num = 1; 97 splitlist<<"barcodeName id id_line "; 98 while( getline( reads, text, ' ') ) 99 { 100 string seq = "", id = ""; 101 id = text; 102 getline( reads, text, ' '); 103 seq = text; 104 105 COMPARE p1 ,p2; 106 p1.ifCompared = p2.ifCompared = 0; 107 for(int k=0;k<N;k++){ 108 p1 = compare(seq,primer,max_mis,0,seq.length()-primer.length()+1); 109 if( p1.ifCompared ) 110 { 111 p2 = compare(seq,barcode[k],0,0,p1.posStart); 112 if ( p2.ifCompared ){ 113 //cout<<seq1<<" "<<barcodeForward[k]<<" "<<seq2<<" "<<barcodeReverse[k]<<endl; 114 splitlist<<barcodeName[k]<<" "<<id<<" "<<line_num<<endl; 115 break; 116 } 117 } 118 } 119 line_num = line_num + 2; 120 } 121 delete[] barcode; 122 delete[] barcodeName; 123 splitlist.close(); 124 barcode_in.close(); 125 reads.close(); 126 } 127 COMPARE compare(string reads, string barcode, int maxmis, int b, int e) 128 { 129 COMPARE p; 130 //cout<<b<<" "<<e<<endl; 131 // cout<<barcode<<" "<<barcode.length()<<endl; 132 for(int i=b;i<e;i++) 133 { 134 int mismatch = 0; 135 int pos = 0; 136 while(mismatch <= maxmis ){ 137 if(reads[i+pos] != barcode[pos]) mismatch++; 138 pos++; 139 if(pos == barcode.length()){ 140 p.posStart = i; 141 p.posEnd = i+pos; 142 // cout<<p.posStart<<" "<<p.posEnd<<endl; 143 p.ifCompared = 1; 144 return p; 145 } 146 } 147 } 148 p.ifCompared = 0; 149 return p; 150 }