C++ EList::push_back方法代码示例

/**
 * Bowtie main function.  It is placed in a separate source file to
 * make it slightly easier to compile Bowtie as a library.
 *
 * If the user specifies -A <file> as the first two arguments, main
 * will interpret that file as having one set of command-line arguments
 * per line, and will dispatch each batch of arguments one at a time to
 * bowtie.
 */
int main(int argc, const char **argv) {
	if(argc > 2 && strcmp(argv[1], "-A") == 0) {
		const char *file = argv[2];
		ifstream in;
		in.open(file);
		char buf[4096];
		int lastret = -1;
		while(in.getline(buf, 4095)) {
			EList<string> args;
			args.push_back(string(argv[0]));
			tokenize(buf, " \t", args);
			const char **myargs = (const char**)malloc(sizeof(char*)*args.size());
			for(size_t i = 0; i < args.size(); i++) {
				myargs[i] = args[i].c_str();
			}
			if(args.size() == 1) continue;
			lastret = bowtie((int)args.size(), myargs);
			free(myargs);
		}
		if(lastret == -1) {
			cerr << "Warning: No arg strings parsed from " << file << endl;
			return 0;
		}
		return lastret;
	} else {
		return bowtie(argc, argv);
	}
}

EList ListGraph::getAdj(NodeID u) const{
	 EList lst;
	 node* temp = ary[u].next;
	 while(temp != NULL){
		 lst.push_back(temp->p);
		 temp = temp->next;
	 };
	 return lst;
 };

   // @note The general idea of this method comes from https://github.com/ScottDVincent/HW05_vincensd_v2/
	std::list<NWPair> MatrixGraph::getAdj(NodeID u) const {
		if (0 <= u < M.size()) {
			EList list; 
			for(int i = 0; i < M.at(u).size(); i++) {
				if (M.at(u).at(i) != 0 ) {
					NWPair pair(i, M.at(u).at(i));
					if (pair.second != 0) 
						list.push_back(pair);
				}
			}
			return list;
		}
	}

std::list<NWPair> ListGraph::getAdj(NodeID u) const
{
    EList temp;
    EList::const_iterator it;
    for(it = edgeList[u].begin(); it != edgeList[u].end(); it++)
    {
        NWPair theEdge = *it;
        if(theEdge.first != NULL)
            temp.push_back(NWPair(theEdge.first, theEdge.second));
    }
    return temp;

}

 EList MatrixGraph::getAdj(NodeID u) const{

	 EList lst;
	 //for the number of nodes
	 for(int i = 0; i < num_nodes; i++){
		 // if they are atached
		 if(ary[u][i] != 0){
			 //add to the vector
			lst.push_back(NWPair(i, ary[u][i]));
		 };
	 };
	 return lst;
 };

/**
 * Calculate a vector containing the sizes of all of the patterns in
 * all of the given input files, in order.  Returns the total size of
 * all references combined.  Rewinds each istream before returning.
 */
std::pair<size_t, size_t>
fastaRefReadSizes(
	EList<FileBuf*>& in,
	EList<RefRecord>& recs,
	const RefReadInParams& rparms,
	BitpairOutFileBuf* bpout,
	int& numSeqs)
{
	uint32_t unambigTot = 0;
	uint32_t bothTot = 0;
	RefReadInParams rpcp = rparms;
	assert_gt(in.size(), 0);
	// For each input istream
	for(size_t i = 0; i < in.size(); i++) {
		bool first = true;
		assert(!in[i]->eof());
		// For each pattern in this istream
		while(!in[i]->eof()) {
			RefRecord rec = fastaRefReadSize(*in[i], rparms, first, bpout);
			if((unambigTot + rec.len) < unambigTot) {
				cerr << "Error: Reference sequence has more than 2^32-1 characters!  Please divide the" << endl
				     << "reference into batches or chunks of about 3.6 billion characters or less each" << endl
				     << "and index each independently." << endl;
				throw 1;
			}
			// Add the length of this record.
			if(rec.first) numSeqs++;
			unambigTot += rec.len;
			bothTot += rec.len;
			bothTot += rec.off;
			first = false;
			if(rec.len == 0 && rec.off == 0 && !rec.first) continue;
			recs.push_back(rec);
		}
		// Reset the input stream
		in[i]->reset();
		assert(!in[i]->eof());
#ifndef NDEBUG
		// Check that it's really reset
		int c = in[i]->get();
		assert_eq('>', c);
		in[i]->reset();
		assert(!in[i]->eof());
#endif
	}
	assert_geq(bothTot, 0);
	assert_geq(unambigTot, 0);
	return make_pair(
		unambigTot, // total number of unambiguous DNA characters read
		bothTot); // total number of DNA characters read, incl. ambiguous ones
}

/**
 * Merge second argument into the first.  Assume both are sorted to
 * begin with.
 */
void Edit::merge(EList<Edit>& dst, const EList<Edit>& src) {
	size_t di = 0, si = 0;
	while(di < dst.size()) {
		if(src[si].pos < dst[di].pos) {
			dst.insert(src[si], di);
			si++; di++;
		} else if(src[si].pos == dst[di].pos) {
			// There can be two inserts at a given position, but we
			// can't merge them because there's no way to know their
			// order
			assert(src[si].isReadGap() != dst[di].isReadGap());
			if(src[si].isReadGap()) {
				dst.insert(src[si], di);
				si++; di++;
			} else if(dst[di].isReadGap()) {
				di++;
			}
		}
	}
	while(si < src.size()) dst.push_back(src[si++]);
}

/**
 * Given the values for all of the various arguments used to specify
 * the read and quality input, create a list of pattern sources to
 * dispense them.
 */
PairedPatternSource* PairedPatternSource::setupPatternSources(
	const EList<string>& si,   // singles, from argv
	const EList<string>& m1,   // mate1's, from -1 arg
	const EList<string>& m2,   // mate2's, from -2 arg
	const EList<string>& m12,  // both mates on each line, from --12 arg
#ifdef USE_SRA
    const EList<string>& sra_accs,
#endif
	const EList<string>& q,    // qualities associated with singles
	const EList<string>& q1,   // qualities associated with m1
	const EList<string>& q2,   // qualities associated with m2
	const PatternParams& p,    // read-in parameters
    size_t nthreads,
	bool verbose)              // be talkative?
{
	EList<PatternSource*>* a  = new EList<PatternSource*>();
	EList<PatternSource*>* b  = new EList<PatternSource*>();
	EList<PatternSource*>* ab = new EList<PatternSource*>();
	// Create list of pattern sources for paired reads appearing
	// interleaved in a single file
	for(size_t i = 0; i < m12.size(); i++) {
		const EList<string>* qs = &m12;
		EList<string> tmp;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmp;
			tmp.push_back(m12[i]);
			assert_eq(1, tmp.size());
		}
		ab->push_back(PatternSource::patsrcFromStrings(p, *qs, nthreads));
		if(!p.fileParallel) {
			break;
		}
	}
    
#ifdef USE_SRA
    for(size_t i = 0; i < sra_accs.size(); i++) {
        const EList<string>* qs = &sra_accs;
        EList<string> tmp;
        if(p.fileParallel) {
            // Feed query files one to each PatternSource
            qs = &tmp;
            tmp.push_back(sra_accs[i]);
            assert_eq(1, tmp.size());
        }
        ab->push_back(PatternSource::patsrcFromStrings(p, *qs, nthreads));
        if(!p.fileParallel) {
            break;
        }
    }
#endif

	// Create list of pattern sources for paired reads
	for(size_t i = 0; i < m1.size(); i++) {
		const EList<string>* qs = &m1;
		EList<string> tmpSeq;
		EList<string> tmpQual;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmpSeq;
			tmpSeq.push_back(m1[i]);
			assert_eq(1, tmpSeq.size());
		}
		a->push_back(PatternSource::patsrcFromStrings(p, *qs, nthreads));
		if(!p.fileParallel) {
			break;
		}
	}

	// Create list of pattern sources for paired reads
	for(size_t i = 0; i < m2.size(); i++) {
		const EList<string>* qs = &m2;
		EList<string> tmpSeq;
		EList<string> tmpQual;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmpSeq;
			tmpSeq.push_back(m2[i]);
			assert_eq(1, tmpSeq.size());
		}
		b->push_back(PatternSource::patsrcFromStrings(p, *qs, nthreads));
		if(!p.fileParallel) {
			break;
		}
	}
	// All mates/mate files must be paired
	assert_eq(a->size(), b->size());

	// Create list of pattern sources for the unpaired reads
	for(size_t i = 0; i < si.size(); i++) {
		const EList<string>* qs = &si;
		PatternSource* patsrc = NULL;
		EList<string> tmpSeq;
		EList<string> tmpQual;
		if(p.fileParallel) {
//.........这里部分代码省略.........

static void driver(
                   const string& infile,
                   EList<string>& infiles,
                   const string& snpfile,
                   const string& htfile,
                   const string& ssfile,
                   const string& exonfile,
                   const string& svfile,
                   const string& outfile,
                   bool packed,
                   int reverse)
{
    initializeCntLut();
    initializeCntBit();
	EList<FileBuf*> is(MISC_CAT);
	bool bisulfite = false;
	RefReadInParams refparams(false, reverse, nsToAs, bisulfite);
	assert_gt(infiles.size(), 0);
	if(format == CMDLINE) {
		// Adapt sequence strings to stringstreams open for input
		stringstream *ss = new stringstream();
		for(size_t i = 0; i < infiles.size(); i++) {
			(*ss) << ">" << i << endl << infiles[i].c_str() << endl;
		}
		FileBuf *fb = new FileBuf(ss);
		assert(fb != NULL);
		assert(!fb->eof());
		assert(fb->get() == '>');
		ASSERT_ONLY(fb->reset());
		assert(!fb->eof());
		is.push_back(fb);
	} else {
		// Adapt sequence files to ifstreams
		for(size_t i = 0; i < infiles.size(); i++) {
			FILE *f = fopen(infiles[i].c_str(), "r");
			if (f == NULL) {
				cerr << "Error: could not open "<< infiles[i].c_str() << endl;
				throw 1;
			}
			FileBuf *fb = new FileBuf(f);
			assert(fb != NULL);
			if(fb->peek() == -1 || fb->eof()) {
				cerr << "Warning: Empty fasta file: '" << infile.c_str() << "'" << endl;
				continue;
			}
			assert(!fb->eof());
			assert(fb->get() == '>');
			ASSERT_ONLY(fb->reset());
			assert(!fb->eof());
			is.push_back(fb);
		}
	}
	if(is.empty()) {
		cerr << "Warning: All fasta inputs were empty" << endl;
		throw 1;
	}
    filesWritten.push_back(outfile + ".1." + gfm_ext);
    filesWritten.push_back(outfile + ".2." + gfm_ext);
	// Vector for the ordered list of "records" comprising the input
	// sequences.  A record represents a stretch of unambiguous
	// characters in one of the input sequences.
	EList<RefRecord> szs(MISC_CAT);
	std::pair<size_t, size_t> sztot;
	{
		if(verbose) cerr << "Reading reference sizes" << endl;
		Timer _t(cerr, "  Time reading reference sizes: ", verbose);
		if(!reverse && (writeRef || justRef)) {
			filesWritten.push_back(outfile + ".3." + gfm_ext);
			filesWritten.push_back(outfile + ".4." + gfm_ext);
			sztot = BitPairReference::szsFromFasta(is, outfile, bigEndian, refparams, szs, sanityCheck);
		} else {
			sztot = BitPairReference::szsFromFasta(is, string(), bigEndian, refparams, szs, sanityCheck);
		}
	}
	if(justRef) return;
	assert_gt(sztot.first, 0);
	assert_gt(sztot.second, 0);
	assert_gt(szs.size(), 0);
    
	// Construct index from input strings and parameters	
    filesWritten.push_back(outfile + ".5." + gfm_ext);
    filesWritten.push_back(outfile + ".6." + gfm_ext);
    filesWritten.push_back(outfile + ".7." + gfm_ext);
    filesWritten.push_back(outfile + ".8." + gfm_ext);
	TStr s;
	HGFM<TIndexOffU> hGFM(
                          s,
                          packed,
                          1,  // TODO: maybe not?
                          lineRate,
                          offRate,      // suffix-array sampling rate
                          ftabChars,    // number of chars in initial arrow-pair calc
                          localOffRate,
                          localFtabChars,
                          nthreads,
                          snpfile,
                          htfile,
                          ssfile,
                          exonfile,
                          svfile,
//.........这里部分代码省略.........

/**
 * A way of feeding simply tests to the seed alignment infrastructure.
 */
int main(int argc, char **argv) {

    EList<string> strs;
    //                            GCTATATAGCGCGCTCGCATCATTTTGTGT
    strs.push_back(string("CATGTCAGCTATATAGCGCGCTCGCATCATTTTGTGTGTAAACCA"
                          "NNNNNNNNNN"
                          "CATGTCAGCTATATAGCGCGCTCGCATCATTTTGTGTGTAAACCA"));
    //                            GCTATATAGCGCGCTTGCATCATTTTGTGT
    //                                           ^
    bool packed = false;
    int color = 0;
	pair<GFM*, GFM*> gfms = GFM::fromStrings<SString<char> >(
		strs,
		packed,
		REF_READ_REVERSE,
		Ebwt::default_bigEndian,
		Ebwt::default_lineRate,
		Ebwt::default_offRate,
		Ebwt::default_ftabChars,
		".aligner_seed2.cpp.tmp",
		Ebwt::default_useBlockwise,
		Ebwt::default_bmax,
		Ebwt::default_bmaxMultSqrt,
		Ebwt::default_bmaxDivN,
		Ebwt::default_dcv,
		Ebwt::default_seed,
		false,  // verbose
		false,  // autoMem
		false); // sanity
    
    gfms.first->loadIntoMemory (-1, true, true, true, true, false);
    gfms.second->loadIntoMemory(1, true, true, true, true, false);
	
	int testnum = 0;

	// Query is longer than ftab and matches exactly twice
    for(int rc = 0; rc < 2; rc++) {
		for(int i = 0; i < 2; i++) {
			cerr << "Test " << (++testnum) << endl;
			cerr << "  Query with length greater than ftab" << endl;
			DescentMetrics mets;
			PerReadMetrics prm;
			DescentDriver dr;
			
			// Set up the read
			BTDnaString seq ("GCTATATAGCGCGCTCGCATCATTTTGTGT", true);
			BTString    qual("ABCDEFGHIabcdefghiABCDEFGHIabc");
			if(rc) {
				seq.reverseComp();
				qual.reverse();
			}
			dr.initRead(Read("test", seq.toZBuf(), qual.toZBuf()), -30, 30);

			// Set up the DescentConfig
			DescentConfig conf;
			conf.cons.init(GFM::default_ftabChars, 1.0);
			conf.expol = DESC_EX_NONE;
			
			// Set up the search roots
			dr.addRoot(
				conf,   // DescentConfig
				(i == 0) ? 0 : (seq.length() - 1), // 5' offset into read of root
				(i == 0) ? true : false,           // left-to-right?
				rc == 0,   // forward?
				0.0f);   // root priority
			
			// Do the search
			Scoring sc = Scoring::base1();
			dr.go(sc, *gfms.first, *gfms.second, mets, prm);
			
			// Confirm that an exact-matching alignment was found
			assert_eq(1, dr.sink().nrange());
			assert_eq(2, dr.sink().nelt());
		}
	}
	
	// Query has length euqal to ftab and matches exactly twice
    for(int i = 0; i < 2; i++) {
		cerr << "Test " << (++testnum) << endl;
		cerr << "  Query with length equal to ftab" << endl;
        DescentMetrics mets;
		PerReadMetrics prm;
        DescentDriver dr;
        
        // Set up the read
        BTDnaString seq ("GCTATATAGC", true);
        BTString    qual("ABCDEFGHIa");
		dr.initRead(Read("test", seq.toZBuf(), qual.toZBuf()), -30, 30);
        
        // Set up the DescentConfig
        DescentConfig conf;
        conf.cons.init(GFM::default_ftabChars, 1.0);
        conf.expol = DESC_EX_NONE;
        
        // Set up the search roots
        dr.addRoot(
            conf,   // DescentConfig
//.........这里部分代码省略.........

/**
 * Merman main driver function.  Does the following:
 *
 * 1. Parses command-line options
 */
int merman(int argc, char **argv) {
	reset();
	try {
		parseCommandLine(argc, argv);
		Timer tov(cerr, "Overall time: ", timing);
		EList<string> refstrs;
		ReferenceSet refs;
		EList<string> refnames;
		EList<size_t> reflens;
		string refstr = argv[optind++];
		tokenize(refstr, ",", refstrs);
		auto_ptr<MerIndex> ind(
			new MerIndex(ap, rp, readLen, seedWidth, nk.first, nk.second,
			             specificity, begin, naiveCheck, nthreads));
		{
			Timer t(cerr, "... ", timing);
			if(timing) cerr << "Reading reference sequences..." << endl;
			for(size_t i = 0; i < refstrs.size(); i++) {
				if(timing) {
					cerr << "  Sequence " << (i+1) << " of " << refstrs.size() << endl;
				}
				if(refIsStr) {
					refs.addOrigReferenceString(refstrs[i].c_str(), rp);
				} else {
					refs.addOrigReferenceFasta(refstrs[i].c_str(), rp);
				}
			}
			for(size_t i = 0; i < refs.numRefs(); i++) {
				refnames.push_back(string(refs[i].name.toZBuf()));
				reflens.push_back(refs[i].seq.length(color));
			}
			if(refs.numRefs() == 0) {
				cerr << "Warning: No references were found" << endl;
			}
			if(rp.genCrick) {
				if(timing) {
					cerr << "  Crickizing" << endl;
				}
				// Add the crick strand.  If there were bisulfite
				// transformations to the Watson strand, they are
				// removed from the Watson strand before the Crick copy
				// is made.  Transformations are then applied to the
				// new Crick strand.  This has the effect of correctly
				// producing either Watson / Crick in the non-bisulfite
				// case, or BS Watson / BS Crick in the bisulfite case.
				refs.addReferenceRevComps(rp, false, 1, 0);
			}
			if(rp.genRevcomps) {
				if(timing) {
					cerr << "  Adding reverse comps" << endl;
				}
				// Add reverse complements of all existing references
				// (after the transformations have already been
				// applied).
				refs.addReferenceRevComps(rp, true, -1, 1);
			}
			assert(refs.repOk());
		}

		pair<size_t, size_t> mers = make_pair(0, 0);
		EList<MerIndexThread> threads;
		{
			Timer t(cerr, "... ", timing);
			if(timing) cerr << "Preparing to extract sub-sequences..." << endl;
			// Instantiate and run index threads
			assert_gt(nthreads, 0);
			threads.resize(nthreads);
			for(int i = 0; i < nthreads; i++) {
				threads[i].runCount(&refs, ind.get(), i, nthreads, color);
			}
			for(int i = 0; i < nthreads; i++) {
				pair<size_t, size_t> mrs = threads[i].join();
				mers.first += mrs.first;
				mers.second += mrs.second;
			}
			ind->allocateMers();
		}
		if(timing || verbose || justBlowup) {
			cerr << "Expecting index footprint of ";
			printBytes(mers.first * sizeof(mer_ent), cerr);
			cerr << endl;
			if(mers.first > mers.second) {
				cerr.setf(ios::fixed);
				cerr << "  base footprint is ";
				printBytes(mers.second * sizeof(mer_ent), cerr);
				cerr << endl
				     << "  blowup factor: " << setprecision(2) << ((double)mers.first / (double)mers.second) << endl;
			}
			if(justBlowup) throw 0;
		}
		{
			Timer t(cerr, "... ", timing);
			if(timing) cerr << "Extracting index sub-sequences..." << endl;
			// Instantiate and run index threads
			for(int i = 0; i < nthreads; i++) {
//.........这里部分代码省略.........

static void driver(
	const char * type,
	const string& bt2indexBase,
	const string& cf_out)
{
	if(gVerbose || startVerbose)  {
		cerr << "Entered driver(): "; logTime(cerr, true);
	}

    //initializeCntLut();  // FB: test commenting

	// Vector of the reference sequences; used for sanity-checking
	EList<SString<char> > names, os;
	EList<size_t> nameLens, seqLens;

	// Initialize Ebwt object and read in header
	if(gVerbose || startVerbose) {
		cerr << "About to initialize fw Ebwt: "; logTime(cerr, true);
	}
	adjIdxBase = adjustEbwtBase(argv0, bt2indexBase, gVerbose);
	Ebwt<index_t> ebwt(
		adjIdxBase,
	    0,        // index is colorspace
		-1,       // fw index
	    true,     // index is for the forward direction
	    /* overriding: */ offRate,
		0, // amount to add to index offrate or <= 0 to do nothing
	    useMm,    // whether to use memory-mapped files
	    useShmem, // whether to use shared memory
	    mmSweep,  // sweep memory-mapped files
	    !noRefNames, // load names?
		true,        // load SA sample?
		true,        // load ftab?
		true,        // load rstarts?
	    gVerbose, // whether to be talkative
	    startVerbose, // talkative during initialization
	    false /*passMemExc*/,
	    sanityCheck);
	//Ebwt<index_t>* ebwtBw = NULL;


	EList<size_t> reflens;
	EList<string> refnames;
	readEbwtRefnames<index_t>(adjIdxBase, refnames);
	map<uint32_t,pair<string,uint64_t> > speciesID_to_name_len;
	for(size_t i = 0; i < ebwt.nPat(); i++) {
		// cerr << "Push back to reflens: "<<  refnames[i] << " is so long: " << ebwt.plen()[i] << endl;
		reflens.push_back(ebwt.plen()[i]);

		// extract numeric id from refName
		const string& refName = refnames[i];
		uint64_t id = extractIDFromRefName(refName);
		uint32_t speciesID = (uint32_t)(id >> 32);

		// extract name from refName
		const string& name_part = refName.substr(refName.find_first_of(' '));

		//uint32_t genusID = (uint32_t)(id & 0xffffffff);
		speciesID_to_name_len[speciesID] = pair<string,uint64_t>(name_part,ebwt.plen()[i]);

	}
//	EList<string> refnames;
//	readEbwtRefnames<index_t>(adjIdxBase, refnames);

	// Read Centrifuge output file
	ifstream infile(cf_out.c_str());

	string line;
	map<uint32_t,uint32_t> species_to_score;

	while (getline(infile,line)) {
		string rd_name;
		uint32_t genusID;
		uint32_t speciesID;
		uint32_t score;
		uint32_t secbest_score;

		istringstream iss(line);
		iss >> rd_name >> genusID >> speciesID >> score >> secbest_score;
		// cerr << rd_name << " -> " << genusID << " -> " << speciesID << " -> " << score << " -> " << secbest_score << "\n";
		species_to_score[speciesID] += score;
	}

	// Sort the species by their score
	vector<pair<uint32_t,uint32_t> > species_to_score_v(species_to_score.begin(), species_to_score.end());

	sort(species_to_score_v.begin(),species_to_score_v.end(),Pair2ndComparator<uint32_t>());

	cout << "Name\tTaxonID\tLength\tSummed Score\tNormalized Score\n";
	// Output the summed species scores
	for (vector<pair<uint32_t,uint32_t> >::iterator species_score = species_to_score_v.begin();
			species_score != species_to_score_v.end();
			++species_score) {
		uint32_t speciesID = species_score->first;
		pair<string,uint64_t> name_len = speciesID_to_name_len[speciesID];
		uint64_t slength = name_len.second;
		uint64_t sumscore = species_score->second;

		cout << name_len.first << "\t" <<
				speciesID << "\t" <<
//.........这里部分代码省略.........

/**
 * Given the values for all of the various arguments used to specify
 * the read and quality input, create a list of pattern sources to
 * dispense them.
 */
PairedPatternSource* PairedPatternSource::setupPatternSources(
	const EList<string>& si,   // singles, from argv
	const EList<string>& m1,   // mate1's, from -1 arg
	const EList<string>& m2,   // mate2's, from -2 arg
	const EList<string>& m12,  // both mates on each line, from --12 arg
	const EList<string>& q,    // qualities associated with singles
	const EList<string>& q1,   // qualities associated with m1
	const EList<string>& q2,   // qualities associated with m2
	const PatternParams& p,    // read-in parameters
	bool verbose)              // be talkative?
{
	//std::cout << "setupPatternSources\n";
	EList<PatternSource*>* a  = new EList<PatternSource*>();
	EList<PatternSource*>* b  = new EList<PatternSource*>();
	EList<PatternSource*>* ab = new EList<PatternSource*>();
	// Create list of pattern sources for paired reads appearing
	// interleaved in a single file
	for(size_t i = 0; i < m12.size(); i++) {
		const EList<string>* qs = &m12;
		EList<string> tmp;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmp;
			tmp.push_back(m12[i]);
			assert_eq(1, tmp.size());
		}
		ab->push_back(PatternSource::patsrcFromStrings(p, *qs));
		if(!p.fileParallel) {
			break;
		}
	}

	// Create list of pattern sources for paired reads
	for(size_t i = 0; i < m1.size(); i++) {
		const EList<string>* qs = &m1;
		EList<string> tmpSeq;
		EList<string> tmpQual;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmpSeq;
			tmpSeq.push_back(m1[i]);
			assert_eq(1, tmpSeq.size());
		}
		a->push_back(PatternSource::patsrcFromStrings(p, *qs));
		if(!p.fileParallel) {
			break;
		}
	}

	// Create list of pattern sources for paired reads
	for(size_t i = 0; i < m2.size(); i++) {
		const EList<string>* qs = &m2;
		EList<string> tmpSeq;
		EList<string> tmpQual;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmpSeq;
			tmpSeq.push_back(m2[i]);
			assert_eq(1, tmpSeq.size());
		}
		b->push_back(PatternSource::patsrcFromStrings(p, *qs));
		if(!p.fileParallel) {
			break;
		}
	}
	// All mates/mate files must be paired
	assert_eq(a->size(), b->size());

	// Create list of pattern sources for the unpaired reads
	for(size_t i = 0; i < si.size(); i++) {
		const EList<string>* qs = &si;
		PatternSource* patsrc = NULL;
		EList<string> tmpSeq;
		EList<string> tmpQual;
		if(p.fileParallel) {
			// Feed query files one to each PatternSource
			qs = &tmpSeq;
			tmpSeq.push_back(si[i]);
			assert_eq(1, tmpSeq.size());
		}
		patsrc = PatternSource::patsrcFromStrings(p, *qs);
		assert(patsrc != NULL);
		a->push_back(patsrc);
		b->push_back(NULL);
		if(!p.fileParallel) {
			break;
		}
	}

	PairedPatternSource *patsrc = NULL;
	if(m12.size() > 0) {
		patsrc = new PairedSoloPatternSource(ab, p);
		for(size_t i = 0; i < a->size(); i++) delete (*a)[i];
		for(size_t i = 0; i < b->size(); i++) delete (*b)[i];
		delete a; delete b;
//.........这里部分代码省略.........

/**
 * Reverse the 'src' list of RefRecords into the 'dst' list.  Don't
 * modify 'src'.
 */
void reverseRefRecords(
	const EList<RefRecord>& src,
	EList<RefRecord>& dst,
	bool recursive,
	bool verbose)
{
	dst.clear();
	{
		EList<RefRecord> cur;
		for(int i = (int)src.size()-1; i >= 0; i--) {
			bool first = (i == (int)src.size()-1 || src[i+1].first);
			// Clause after the || on next line is to deal with empty FASTA
			// records at the end of the 'src' list, which would be wrongly
			// omitted otherwise.
			if(src[i].len || (first && src[i].off == 0)) {
				cur.push_back(RefRecord(0, src[i].len, first));
				first = false;
			}
			if(src[i].off) cur.push_back(RefRecord(src[i].off, 0, first));
		}
		bool mergedLast;
		for(int i = 0; i < (int)cur.size(); i++) {
			mergedLast = false;
			assert(cur[i].off == 0 || cur[i].len == 0);
			if(i < (int)cur.size()-1 && cur[i].off != 0 && !cur[i+1].first) {
				dst.push_back(RefRecord(cur[i].off, cur[i+1].len, cur[i].first));
				i++;
				mergedLast = true;
			} else {
				dst.push_back(cur[i]);
			}
		}
	}
	//if(verbose) {
	//	cout << "Source: " << endl;
	//	printRecords(cout, src);
	//	cout << "Dest: " << endl;
	//	printRecords(cout, dst);
	//}
#ifndef NDEBUG
	size_t srcnfirst = 0, dstnfirst = 0;
	for(size_t i = 0; i < src.size(); i++) {
		if(src[i].first) {
			srcnfirst++;
		}
	}
	for(size_t i = 0; i < dst.size(); i++) {
		if(dst[i].first) {
			dstnfirst++;
		}
	}
	assert_eq(srcnfirst, dstnfirst);
	if(!recursive) {
		EList<RefRecord> tmp;
		reverseRefRecords(dst, tmp, true);
		assert_eq(tmp.size(), src.size());
		for(size_t i = 0; i < src.size(); i++) {
			assert_eq(src[i].len, tmp[i].len);
			assert_eq(src[i].off, tmp[i].off);
			assert_eq(src[i].first, tmp[i].first);
		}
	}
#endif
}

//.........这里部分代码省略.........
	cerr << "PASSED" << endl;

	cerr << "Test S2bDnaString ...";
	{
		const char *str =
			"ACGTACGTAC" "ACGTACGTAC" "ACGTACGTAC"
			"ACGTACGTAC" "ACGTACGTAC" "ACGTACGTAC";
		const char *gs =
			"GGGGGGGGGG" "GGGGGGGGGG" "GGGGGGGGGG"
			"GGGGGGGGGG" "GGGGGGGGGG" "GGGGGGGGGG";
		for(size_t i = 0; i < 60; i++) {
			S2bDnaString s(str, i, true);
			S2bDnaString sr;
			BTDnaString s2(str, i, true);
			assert(sstr_eq(s, s2));
			if(i >= 10) {
				BTDnaString s3;
				s.windowGetDna(s3, true, 3, 4);
				assert(sstr_eq(s3.toZBuf(), (const char*)"TACG"));
				s.windowGetDna(s3, false, 3, 4);
				assert(sstr_eq(s3.toZBuf(), (const char*)"CGTA"));
				assert_eq('A', s.toChar(0));
				assert_eq('G', s.toChar(2));
				assert_eq('A', s.toChar(4));
				assert_eq('G', s.toChar(6));
				assert_eq('A', s.toChar(8));
				
				s.reverseWindow(1, 8);
				s2.reverseWindow(1, 8);
				
				assert_eq('A', s.toChar(1));
				assert_eq('T', s.toChar(2));
				assert_eq('G', s.toChar(3));
				assert_eq('C', s.toChar(4));
				assert_eq('A', s.toChar(5));
				assert_eq('T', s.toChar(6));
				assert_eq('G', s.toChar(7));
				assert_eq('C', s.toChar(8));
				assert(sstr_eq(s, s2));

				s.reverseWindow(1, 8);
				s2.reverseWindow(1, 8);
				assert(sstr_eq(s, s2));
			}
			if(i > 1) {
				s.reverse();
				sr.installReverseChars(str, i);
				s2.reverse();
				assert(sstr_eq(s, s2));
				assert(sstr_eq(sr, s2));
				s.reverse();
				sr.reverse();
				assert(sstr_neq(s, s2));
				assert(sstr_neq(sr, s2));
				s.fill(2);
				s2.reverse();
				assert(sstr_leq(s, gs));
				assert(sstr_gt(s, s2));
				assert(sstr_gt(s, sr));
				s2.fill(2);
				sr.fill(2);
				assert(sstr_eq(s, s2));
				assert(sstr_eq(s, sr));
			}
		}
		S2bDnaString s(str, true);
		S2bDnaString sr;
		BTDnaString s2(str, true);
		assert(sstr_eq(s2.toZBuf(), str));
		assert(sstr_eq(s, s2));
		s.reverse();
		sr.installReverseChars(str);
		s2.reverse();
		assert(sstr_eq(s, s2));
		assert(sstr_eq(sr, s2));
		s.reverse();
		sr.reverse();
		assert(sstr_neq(s, s2));
		assert(sstr_neq(sr, s2));
	}
	cerr << "PASSED" << endl;

	cerr << "Test operator=() ...";
	{
		S2bDnaString s;
		s.installChars(string("gtcagtca"));
		assert(sstr_eq(s.toZBuf(), (const char *)"GTCAGTCA"));
	}
	cerr << "PASSED" << endl;
	
	cerr << "Conversions from string ...";
	{
		SStringExpandable<char> se(string("hello"));
		EList<SStringExpandable<char> > sel;
		sel.push_back(SStringExpandable<char>(string("hello")));
	}
	cerr << "PASSED" << endl;
	
	cerr << "PASSED" << endl;
}