InputStream.cpp

/* sdm: simple demultiplexer
Copyright (C) 2013 Falk Hildebrand
email: Falk.Hildebrand@gmail.com

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/


#ifndef _CRT_SECURE_NO_DEPRECATE
#define _CRT_SECURE_NO_DEPRECATE
#endif

#include <mutex>
#include <thread>
#include <atomic>
#include "InputStream.h"
#include "Common.h"
#include "FastxReader.h"

// Utilities are implemented in Common.cpp (defined once in Common.cpp)

namespace {
constexpr int MAX_IO_WARNINGS_TO_PRINT = 20;
std::atomic<int> g_ioWarningsSeen{ 0 };
std::atomic<bool> g_useClassicSeedSelection{ false };

std::string normalizedRecordId(const std::string& header) {
	if (header.empty()) { return std::string(); }
	std::string id = header;
	if (id[0] == '@' || id[0] == '>') {
		id.erase(0, 1);
	}
	remove_paired_info(id);
	return id;
}

std::string exactRecordId(const std::string& header) {
	if (header.empty()) { return std::string(); }
	const size_t start = (header[0] == '@' || header[0] == '>') ? 1u : 0u;
	const size_t end = header.find_first_of(" \t", start);
	return header.substr(start, end == std::string::npos ? std::string::npos : end - start);
}

bool recordBufferHasData(const std::vector<std::string>& record) {
	return std::any_of(record.begin(), record.end(), [](const std::string& value) {
		return !value.empty();
	});
}

void validateSynchronizedHeaders(const std::vector<std::string>& first,
	const std::vector<std::string>& second, const char* secondLabel) {
	const std::string firstId = first.empty() ? std::string() : normalizedRecordId(first[0]);
	const std::string secondId = second.empty() ? std::string() : normalizedRecordId(second[0]);
	if (firstId.empty() || secondId.empty() || firstId != secondId) {
		throw std::runtime_error("Input records are not synchronized: R1 id '" + firstId +
			"' does not match " + secondLabel + " id '" + secondId + "'");
	}
}

bool readFastaRecord(ifbufstream* fasta, ifbufstream* quality, std::vector<std::string>& record) {
	record.assign(3, std::string());
	if (!fasta || fasta->eof()) {
		if (quality && !quality->eof()) {
			throw std::runtime_error("QUAL input contains records after its FASTA input ended");
		}
		return false;
	}

	const bool gotHeader = fasta->getlines(record[0], false);
	const bool gotSequence = gotHeader && fasta->getlines(record[1], true);
	if (!gotHeader && record[0].empty()) {
		return false;
	}
	if (!gotHeader || !gotSequence || record[0].empty() || record[0][0] != '>') {
		throw std::runtime_error("Truncated or malformed FASTA record near header '" + record[0] + "'");
	}

	if (quality) {
		std::string qualityHeader;
		const bool gotQualityHeader = quality->getlines(qualityHeader, false);
		const bool gotQualities = gotQualityHeader && quality->getlines(record[2], true, true);
		if (!gotQualityHeader || !gotQualities) {
			throw std::runtime_error("QUAL input ended before FASTA record '" + normalizedRecordId(record[0]) + "'");
		}
		const std::string fastaId = exactRecordId(record[0]);
		const std::string qualityId = exactRecordId(qualityHeader);
		if (fastaId.empty() || qualityId.empty() || fastaId != qualityId) {
			throw std::runtime_error("FASTA/QUAL headers are not synchronized: '" + fastaId +
				"' versus '" + qualityId + "'");
		}
	}
	return true;
}
}

void setDNAseedSelectionClassic(bool enabled) {
	g_useClassicSeedSelection.store(enabled, std::memory_order_relaxed);
}

bool getDNAseedSelectionClassic() {
	return g_useClassicSeedSelection.load(std::memory_order_relaxed);
}


bool whoIsBetter_old(shared_ptr<DNA> d1, shared_ptr<DNA> d2, shared_ptr<DNA> dM,
	shared_ptr<DNA> r1, shared_ptr<DNA> r2, shared_ptr<DNA> rM,
	float& ever_best, bool forSeed) {

	//if (forSeed) {		return false;	}
	//check if two primers present
	if (d2 == nullptr) {//hard reason .. only for PacBio etc reads
		if (d1->has2PrimersDetected() && !r1->has2PrimersDetected()) { return true; }
		if (!d1->has2PrimersDetected() && r1->has2PrimersDetected()) { return false; }
	}
	else {
		if (d2->getRevPrimDetect() && !r1->getFwdPrimDetect() && !r2->getRevPrimDetect()) { return true; }
	}
	//check if at least 1 primers present
	if (d1->getFwdPrimDetect() && !r1->getFwdPrimDetect()) { return true; }//hard reason
	if (!d1->getFwdPrimDetect() && r1->getFwdPrimDetect()) { return false; }

	double d1pid(1.), refpid(1.);
	if (ever_best >= 0) {
		d1pid = (double)d1->getTempFloat(); refpid = (double)r1->getTempFloat();
		if (d1pid > ever_best) {
			ever_best = (float)d1pid;
		}
		//everbest is likely 100.f (ref OTUs)
		if (d1pid < (refpid - 0.3f) || d1pid < (ever_best - 0.4f)) { return false; }
	}

	bool dMerge(true), rMerg(true);
	double curL = (double)d1->getMergeLength();
	if (curL < 0) {
		curL = (double)d1->length();
		//if (d2 != NULL) { curL += (double) d2->length(); }
		dMerge = false;
	}
	double refL = (double)r1->getMergeLength();
	if (refL < 0) {
		refL = (double)r1->length();
		//if (r2 != NULL) { refL += (double) r2->length(); }
		rMerg = false;
	}

	//first check if d1 has merged, but ref did not.. clearly go for d, hard filter
	//if (d1->getMergeLength() != -1 && r1->getMergeLength() == -1) { return true; }

	//hard check on length ratios.. too drastically small , don't use d1
	if ((curL) / (refL) < BestLengthRatio) { return false; }
	//at least 90% length of "good" hit
	//if (r1->getMergeErrors() < 0) {//no merge, can look at read1 only
	//	if (d1->length() / r1->length() < RefLengthRatio) { return false; }
	//}


	float dmergErrSco = 0.f;	float refMergErrSco = 0.f;
	//only check further if both comparisons did merge
	if (r1->getMergeLength() >= 0 && d1->getMergeLength() >= 0) {
		if (d1->getMergeErrors() > 0) {
			dmergErrSco = (float)d1->getMergeErrors();// *log10((maxQErr - (float)d1->getMergeErrorsQual()));
			dmergErrSco += d1->getMergeErrorsQual() / 30;
		}
		if (r1->getMergeErrors() > 0) {
			refMergErrSco = (float)r1->getMergeErrors() + r1->getMergeErrorsQual() / 30;// *log10((maxQErr - (float)r1->getMergeErrorsQual()));
		}
		//scale to 1
		float maxMerr = max(refMergErrSco, dmergErrSco);
		dmergErrSco /= maxMerr;	refMergErrSco /= maxMerr;
		dmergErrSco = 1.f - dmergErrSco; refMergErrSco = 1.f - refMergErrSco;
	}

	//choose merged DNA if possible; d2 no longer needed then
	shared_ptr<DNA> dx, rx;
	bool allowMergeGuide = false;
	if (allowMergeGuide && dM != nullptr) { dx = dM;	d2 = nullptr; }
	else { dx = d1; }
	if (allowMergeGuide && rM != nullptr) { rx = rM; r2 = nullptr; }
	else { rx = r1; }
	float thScore = dx->getAvgQual(); //*(d1pid / 100)* log((float)curL);
	float rScore = rx->getAvgQual();// *(refpid / 100)* log((float)refL);//r1->length()
	//if (thScore > rScore) {
		//also check for stable lowest score
		// if (d1->minQual() > r1->minQual() - MinQualDiff) { return true; }
	//}
	float maxScore = max(thScore, rScore);
	thScore /= maxScore;	rScore /= maxScore;


	//checks if the new DNA has a better overall quality
	double dAcSc = dx->getAccumError();	double dLen = dx->mem_length();
	double tAcSc = rx->getAccumError();	double tLen = rx->mem_length();
	if (d2 != nullptr) { dAcSc += d2->getAccumError(); dLen += d2->mem_length(); }
	if (r2 != nullptr) { tAcSc += r2->getAccumError(); tLen += r2->mem_length(); }
	dAcSc /= dLen;	tAcSc /= tLen;

	//calculate ratios to compare more easily among metrics
	double ratAccErr = log(dAcSc) / log(tAcSc); //smaller better, log changes terms around
	double ratLength = double(curL) / (double)refL; //higher better
	double ratId = d1pid / refpid * 10 - 9.; //higher better //10-fold weighting
	if ((ratAccErr * ratLength * ratId) > 1.) {
		return true;
	}
	//normalize and invert
	//normalize to gene length, and invert to convert to positive score system
	/*maxScore = max(dAcSc, tAcSc);
	dAcSc /= maxScore;	tAcSc /= maxScore;
	dAcSc = 1.f - dAcSc; tAcSc = 1.f - tAcSc;
	//norm to 1, to compare to other terms
	double maxEr = max(dAcSc, tAcSc);
	dAcSc /= maxEr; tAcSc /= maxEr;
	//dmergErrSco refMergErrSco dAcSc + tAcSc + thScore rScore
	if ( (dAcSc) *(d1pid / 100) * log((float)curL)
		> (tAcSc) *(refpid / 100) * log((float)refL )) {

		if (dx->minQual() > rx->minQual() - MinQualDiff) { return true; }
//		return true;

	}
	*/

	return false;
}


//compares two DNA entries, decides which one has overall better stats
bool whoIsBetter(shared_ptr<DNA> d1, shared_ptr<DNA> d2, shared_ptr<DNA> dM, 
			shared_ptr<DNA> r1, shared_ptr<DNA> r2, shared_ptr<DNA> rM, 
			float& ever_best, bool forSeed, DNAunique* merge_stats_owner,
			uint64_t dSiz , uint64_t rSiz) {

	if (forSeed && getDNAseedSelectionClassic()) {
		return whoIsBetter_old(d1, d2, dM, r1, r2, rM, ever_best, forSeed);
	}
	
	if (d1 == nullptr || r1 == nullptr) {
		return false;
	}

	//if (forSeed) { 		return false; 	}
	//check if two primers present
	if (d2 == nullptr) {//hard reason .. only for PacBio etc reads
		if (d1->has2PrimersDetected() && !r1->has2PrimersDetected()) { return true; }
		if (!d1->has2PrimersDetected() && r1->has2PrimersDetected()) { return false; }
 } else {
		//if (r2 != nullptr && d2->getRevPrimDetect() && !r1->getFwdPrimDetect() && !r2->getRevPrimDetect()) {	return true;}
	}
	//check if at least 1 primers present
	if (d1->getFwdPrimDetect() && !r1->getFwdPrimDetect()) { return true; }//hard reason
	if (!d1->getFwdPrimDetect() && r1->getFwdPrimDetect()) { return false; }

	double d1pid(1.), refpid(1.);
	if (ever_best >=0) {
		d1pid = (double) d1->getTempFloat(); refpid = (double)r1->getTempFloat();
		if (d1pid > ever_best) {
			ever_best = (float) d1pid;
		}
		//everbest is likely 100.f (ref OTUs)
		//weakened a bit from 0.3/0.4 to 0.6/0.8, to allow for more variability in seed selection, which is not as critical as final selection
		//if (d1pid < (refpid - 0.6f) || d1pid < (ever_best - 0.8f)) { return false; }
		//reversed
		if (d1pid < (refpid - 0.3f) || d1pid < (ever_best - 0.5f)) { return false; }
	}
	double ratId = 1.;
	if (refpid > 0.) {
		if (forSeed) {
			// Stronger identity discrimination for seed selection (legacy-like behavior).
			ratId = (double(d1pid) / double(refpid)) * 50 - 49.0;
		}
		else {
			ratId = pointFiveRatio(double(d1pid), (double)refpid, 2.0);
		}
	}


	bool dMerge(true), rMerg(true);
	double curL = (double)d1->getMergeLength();
	if (curL < 0) {
		curL = (double) d1->length();
		//if (d2 != NULL) { curL += (double) d2->length(); }
		dMerge = false;
	}
	double refL = (double)r1->getMergeLength();
	if (refL < 0) {
		refL = (double) r1->length();
		//if (r2 != NULL) { refL += (double) r2->length(); }
		rMerg = false;
	}
 if (merge_stats_owner != nullptr) {
		merge_stats_owner->recordWhoIsBetterCompare(d1->getMergeLength(), r1->getMergeLength());
	}
 double mergedFraction = 0.0;
	double avgMergeLen = 0.0;
	if (merge_stats_owner != nullptr) {
		const auto& ms = merge_stats_owner->getWhoIsBetterMergeStats();
		const uint64_t totalCompared = ms.merged_compares + ms.non_merged_compares;
		if (totalCompared > 0) {
			mergedFraction = static_cast<double>(ms.merged_compares) / static_cast<double>(totalCompared);
		}
		avgMergeLen = ms.merged_compares > 0 ? static_cast<double>(ms.accumulated_merge_length) / static_cast<double>(ms.merged_compares) : 0.0;
	}
	bool mergAdv = dMerge && !rMerg;//if d merged but ref did not, go for d, hard filter

	//first check if d1 has merged, but ref did not.. clearly go for d, hard filter
	//if (d1->getMergeLength() != -1 && r1->getMergeLength() == -1) { return true; }

	//hard check on length ratios.. too drastically small , don't use d1
 if (refL <= 0.) { return false; }
	if ((curL) / (refL) < BestLengthRatio && !mergAdv) { return false; }
	//at least 90% length of "good" hit
	//if (r1->getMergeErrors() < 0) {//no merge, can look at read1 only
	//	if (d1->length() / r1->length() < RefLengthRatio) { return false; }
	//}

	//also check for number of N in sequence
	int d1Ncnt = d1->numNonCanonicalDNA(true); int r1Ncnt = r1->numNonCanonicalDNA(true);

	float CanoScore = 1.f;
	if (d1Ncnt > 0 || r1Ncnt > 0) {
		CanoScore = 1.f - (float)(d1Ncnt - r1Ncnt) / (float)(d1Ncnt + r1Ncnt + 10);
	}

	float dmergErrSco = 0.f; 	float refMergErrSco = 0.f;
	//only check further if both comparisons did merge
	if (r1->getMergeLength() >=0 && d1->getMergeLength() >=0) {
		if (d1->getMergeErrors() > 0) {
			dmergErrSco = (float)d1->getMergeErrors();// *log10((maxQErr - (float)d1->getMergeErrorsQual()));
			dmergErrSco += d1->getMergeErrorsQual()/30;
		}
		if (r1->getMergeErrors() > 0) {
			refMergErrSco = (float)r1->getMergeErrors() + r1->getMergeErrorsQual() / 30;// *log10((maxQErr - (float)r1->getMergeErrorsQual()));
		}
		//scale to 1
 float maxMerr = max(refMergErrSco, dmergErrSco);
		if (maxMerr > 0.f) {
			dmergErrSco /= maxMerr; 	refMergErrSco /= maxMerr;
			dmergErrSco = 1.f - dmergErrSco; refMergErrSco = 1.f - refMergErrSco;
		}
	}

	//choose merged DNA if possible; d2 no longer needed then
	shared_ptr<DNA> dx, rx;
	bool allowMergeGuide = dMerge && rMerg;
	if (allowMergeGuide && dM != nullptr) { dx = dM; 	d2 = nullptr; } 	else { dx = d1; }
	if (allowMergeGuide && rM != nullptr) { rx = rM; r2 = nullptr;}else { rx = r1; }
	float thScore = dx->getAvgQual(); //*(d1pid / 100)* log((float)curL);
	float rScore = rx->getAvgQual();// *(refpid / 100)* log((float)refL);//r1->length()
	//if (thScore > rScore) {
		//also check for stable lowest score
		// if (d1->minQual() > r1->minQual() - MinQualDiff) { return true; }
	//}
	float maxScore = max(thScore, rScore);
	if (maxScore > 0.f) {
		thScore /= maxScore; 	rScore /= maxScore;
	}


	/*
	* accumulated error was overall too harsh.. switch to qual ratio instead..
	//checks if the new DNA has a better overall quality
	double dAcSc = dx->getAccumError();	double dLen = dx->mem_length();
	double tAcSc = rx->getAccumError();	double tLen = rx->mem_length();
	if (d2 != nullptr) { dAcSc += d2->getAccumError(); dLen += d2->mem_length(); }
	if (r2 != nullptr) { tAcSc += r2->getAccumError(); tLen += r2->mem_length(); }
	if (dLen <= 0. || tLen <= 0.) {
		return false;
	}
	//dAcSc /= dLen; 	tAcSc /= tLen;
	
	//calculate ratios to compare more easily among metrics
	double ratAccErr = dAcSc / tAcSc;///logRatio(dAcSc, tAcSc);
	*/

	double dAccuE = dx->getAccumError(); double rAccuE = rx->getAccumError();
	double AccErrRatio = pointFiveRatio(dAccuE, rAccuE,0.25);

	double dQ=dx->getAvgQual(); double rQ=rx->getAvgQual();
	if (d2 != nullptr) { dQ += d2->getAvgQual(); }
	if (r2 != nullptr) { rQ += r2->getAvgQual(); }

	double qualRatio = pointFiveRatio(dQ, rQ, 0.5); //higher better


	//double ratLength = double(curL) / (double)refL; //higher better
	double ratLength = pointFiveRatio(double(curL), (double)refL); //higher better, but not too much better, as this can be due to chimeras etc



	//preference based on merge status, but modulated by how often merged DNAs have been better in the past (mergedFraction)
	double mergePreference = 1.0;
	if (dMerge != rMerg) {
		if (forSeed) {
			// In seed mode, disable historical merge-bias but still favor merged reads.
			mergePreference = dMerge ? 1.35 : 0.75;
		}
		else if (dMerge) {
			mergePreference = 1.0 + mergedFraction;
		}
		else {
			mergePreference = 1.0 - (0.5 * mergedFraction);
		}
	}

	//also take into account the merge length and if this is the averaged norm:
	double mergeLenDevRatio = 1.;
	if (!forSeed && dMerge && rMerg) {
		double dMLDIFF = (double) abs(avgMergeLen - d1->getMergeLength());
		double rMLDIFF = (double) abs(avgMergeLen - r1->getMergeLength());
		mergeLenDevRatio = (dMLDIFF+ avgMergeLen) / (rMLDIFF+ avgMergeLen);
		//mergeLenDeviation = (double)(rMLDIFF - dMLDIFF) / (double)(avgMergeLen + 1e-12); //positive if d is closer to average merge length, negative if r is closer
	}

	//also capture the size difference.. as higher abudant reads might be more trustworthy, but only if other metrics are not too different
	double logsizeRatio = pointFiveRatio(double(dSiz), double(rSiz),0.5); //higher better for r


	double thresh(1.01f);
	auto clampScore = [](double v) {
		if (!std::isfinite(v)) {
			return 1.0;
		}
		if (v < 0.2) {
			return 0.2;
		}
		if (v > 2.5) {
			return 2.5;
		}
		return v;
	};

	double wCano = 0.10;
	double wAccErr = 0.20;
	double wQual = 0.25;
	double wLen = 0.10;
	double wId = 0.15;
	double wMergePref = 0.08;
	double wSize = 0.06;
	double wMergeLenDev = 0.06;

	if (forSeed) {
		// Seed selection: prioritize identity and sequence quality.
		wCano = 0.08; //some preference for fewer Ns, but not too much, as this can be due to chimeras etc
		wAccErr = 0.44; //some preference for lower accumulated error, but not too much, as this can be affected by length and other factors
		wQual = 0.14; //strong preference for higher average quality, as this is a strong indicator of good seeds
		wLen = 0.08;//some preference for longer length, but not too much, as this can be due to chimeras etc
		wId = 0.42;//strong preference for higher identity to ref OTU, as this is a strong indicator of good seeds
		wMergePref = 0.10;//some preference for merged reads, but not too much, as this can be affected by the specific dataset and merge performance
		wSize = 0.02;//small preference for higher abundance, but not too much, as this can be affected by the specific dataset and amplification biases
		wMergeLenDev = 0.00;//no preference for merge length deviation in seed selection, as this is more relevant for representative selection where we want to avoid outliers, but for seeds we want to allow some diversity in merge lengths
	}
	else {
		// Derep representative selection: prioritize merged-read behavior first,
		// then merge-length consistency and sequence quality.
		wCano = 0.08;
		wAccErr = 0.32;
		wQual = 0.16;
		wLen = 0.08;
		wId = 0.0;
		wMergePref = 0.34;
		wSize = 0.04;
		wMergeLenDev = 0.38;
	}

	double cumulativeQualityScore =
		std::pow(clampScore(CanoScore), wCano) *
		std::pow(clampScore(AccErrRatio), wAccErr) *
		std::pow(clampScore(qualRatio), wQual) *
		std::pow(clampScore(ratLength), wLen) *
		std::pow(clampScore(ratId), wId) *
		std::pow(clampScore(mergePreference), wMergePref) *
		std::pow(clampScore(logsizeRatio), wSize) *
		std::pow(clampScore(mergeLenDevRatio), wMergeLenDev);

	// Previous formula (kept for reference):
	// if ((CanoScore * AccErrRatio * qualRatio * ratLength * ratId * mergePreference * logsizeRatio * mergeLenDevRatio) > thresh) {
	if (cumulativeQualityScore > thresh) {
		return true;
	}
	//normalize and invert
	//normalize to gene length, and invert to convert to positive score system
	/*maxScore = max(dAcSc, tAcSc);
	dAcSc /= maxScore;	tAcSc /= maxScore;
	dAcSc = 1.f - dAcSc; tAcSc = 1.f - tAcSc;
	//norm to 1, to compare to other terms
	double maxEr = max(dAcSc, tAcSc); 
	dAcSc /= maxEr; tAcSc /= maxEr;
	//dmergErrSco refMergErrSco dAcSc + tAcSc + thScore rScore
	if ( (dAcSc) *(d1pid / 100) * log((float)curL)
		> (tAcSc) *(refpid / 100) * log((float)refL )) {
		
		if (dx->minQual() > rx->minQual() - MinQualDiff) { return true; }
//		return true;

	}
	*/

	return false;
}


bool whoIsBetter(shared_ptr<DNAunique> d1,
	shared_ptr<DNAunique> r1,
	float& ever_best, bool forSeed) {
	

	bool ret = whoIsBetter(d1, d1->getPair(), d1->getMerge(),
		r1, r1->getPair(), r1->getMerge(), ever_best, forSeed, d1.get(),
		d1->totalSum(), r1->totalSum());

	return ret;

}

double pointFiveRatio(double v1, double v2, double scale ) {
	
	double ret = 1.0;
	
	if (scale <= 0.) { cerr << "pointFiveRatio::scale not ok: "<<scale<<endl; scale = 0.5; }
	if ((v1 + v2) <= 0.) { return 1.0; }
	
	ret = (v1 / (v1 + v2)) + 0.5;
	if (scale != 0.5) {
		ret = (-1. * (1. - ret) * scale) + 1;
	}
	return ret;
}

double logRatio(double v1, double v2) {
	double ret = 1.f;
	if (v1 > 0. && v2 > 0.) {
		double logRef = std::log(v2);
		if (std::isfinite(logRef) && std::abs(logRef) > 1e-12) {
			ret = std::log(v1) / logRef; //smaller better, log changes terms around
		}
		else {
			ret = (v1 / (v1+v2)) + 0.5;
		}
	}
	if (ret > 0.) {
		return ret;
	}
	return -1. * ret;
}

double loglogRatio(double v1, double v2) {
	double ret = 1.f;
	if (v1 > 0. && v2 > 0.) {
		double logRef = std::log(std::log(v2));
		if (std::isfinite(logRef) && std::abs(logRef) > 1e-12) {
			ret = std::log(std::log(v1)) / logRef; //smaller better, log changes terms around
		}
		else {
			ret = v1 / v2;
		}
	}
	if (ret > 0.) {
		return ret;
	}
	return -1. * ret;
}





///////////////////////////////////////////////////////////////
//INPUT STREAMER


///////////////////////////////////////////////////////////////
//INPUT STREAMER

InputStreamer::~InputStreamer(){
	allStreamClose();
	for (size_t i = 0; i < pairs_read.size(); i++) {
		cdbg("~InputStreamer: pairsRead" + to_string(i) + ": " + to_string(pairs_read[i]) + " ");
	}
	cdbg("Returned DNA objects: "+to_string(N_DNAreturned) + " pairs: " + to_string(N_DNApairsReturned) + " ");
}

bool InputStreamer::read_fasta_entry(ifbufstream* fasta_is, ifbufstream* quality_is, shared_ptr<DNA> tdn1, shared_ptr<DNA> tdn2, int &cnt){
	if (fasta_is->eof()) return false;
	string line;
	string tseq;
	string tqual;
	string lineQ;
	fasta_is->getlines(line,false);
	if (!qualAbsent) {
		quality_is->getlines(lineQ,false);
	}
	cnt++;

	if (cnt == 1) {
		if (line[0] != '>' && fasta_is) {
			cerr << "ERROR: Line 1 in fasta file does not start with \">\" \n";
			exit(23);
		}
		tdn1->setHeader(line.substr(1));
		if (!fasta_is->getlines(line,true))
			return false;
		if (!qualAbsent) {
			quality_is->getlines(lineQ,true);
		}
		cnt++;
	}

	tdn1->setSequence(tseq);
	size_t lsize = tseq.size();
	vector<qual_score> Iqual(lsize, 11);
	if (!qualAbsent) {
		rtrim(tqual);
		const char* lQ = tqual.c_str();
		uint ii(0);
		qual_score nn(0);

		for (; ii < lsize; ii++) {
			nn = (qual_score)parseInt(&lQ);
			Iqual[ii] = nn;
			if (*lQ == '\0') break;
		}

		if (Iqual.size() != lsize) { cerr << "Unequal fasta/qual length in read_fasta_entry\n"; exit(923); }
	}
	tdn1->setQual(move(Iqual));
	if (!line.empty() && line[0] == '>') {
		tdn2->setHeader(line.substr(1));
	} else {
		tdn2->setHeader("");
	}
	return true;
}

void InputStreamer::IO_Error(string x) {
	fqverMTX.lock();
	int seen = g_ioWarningsSeen.fetch_add(1, std::memory_order_relaxed);
	if (seen < MAX_IO_WARNINGS_TO_PRINT) {
		cerr << x << endl;
	}
	else if (seen == MAX_IO_WARNINGS_TO_PRINT) {
		cerr << "Further I/O parse warnings suppressed after " << MAX_IO_WARNINGS_TO_PRINT << " messages." << endl;
	}
	if (DieOnError) exit(632);
	if ((int)ErrorLog.size() < MAX_IO_WARNINGS_TO_PRINT) {
		ErrorLog.push_back(x);
	}
	fqverMTX.unlock();
}

shared_ptr<DNA> InputStreamer::read_fastq_entry_fast(istream & fna, int& lnCnt, bool& corrupt) {
	string line;
	if (!fna) { return NULL; }
	if (!safeGetline(fna, line)) { return NULL; }
	shared_ptr<DNA> tdn = make_shared<DNA>("", "");
	if (line.length() == 0) { return tdn; }
	while (line[0] != '@') {
		IO_Error("ERROR on line " + itos(lnCnt) + ": Could not find '@' when expected (file likely corrupt, trying to recover):\n" + line);
		corrupt = true;
		if (!safeGetline(fna, line)) { corrupt = true; return NULL; }
	}
	tdn->setHeader(line.substr(1));
	if (!safeGetline(fna, tdn->getSequence())) { corrupt = true; return NULL; }
	if (!safeGetline(fna, line)) { corrupt = true; return NULL; }
	if (line[0] != '+') { IO_Error("Error input line " + itos(lnCnt + 2) + ": Could not find '+' when expected (file likely corrupt, aborting):\n" + line); corrupt = true; return tdn; }

	vector<qual_score> Iqual(tdn->mem_length(), 0);
	if (!safeGetline(fna, line)) { corrupt = true; return NULL; }
	uint qcnt(0); uint lline = (uint)line.length();
	for (; qcnt < lline; qcnt++) {
		qual_score q = (qual_score)line[qcnt] - fastQver;
		Iqual[qcnt] = q;
	}

	if (qcnt == tdn->mem_length()) {
		tdn->setQual(move(Iqual));
	} else if (line.length() + qcnt != tdn->length()) {
		corrupt = true;
		IO_Error("Error input line " + itos(lnCnt + 3) + ": More quality positions than nucleotides detected for sequence\n " + tdn->getId());
	}
	lnCnt+=4;
	corrupt = false;
	return tdn;
}

void InputStreamer::minmaxQscore(shared_ptr<DNA> t, bool& print) {
	fqverMTX.lock();
	vector<qual_score> Qs = t->getQual();
	for (size_t i = 0; i < Qs.size(); i++) { minmaxQscore(Qs[i], print); }
	fqverMTX.unlock();
}

qual_score InputStreamer::minmaxQscore(qual_score t,bool& print) {
	if (t < 0) {
		if (fqSolexaFmt ){
 if (t < -5 && print) { cerr << "Unusually low Solexa quality score (" << t << "); setting to 0.\n"; print = false; }
		} else {
 if (fastQver == 64 && t >= -5) {
				fqSolexaFmt = true;
				if (print) { cerr << "Detected negative qualities in +64 range; enabling Solexa interpretation.\n"; print = false; }
			}
			else if (fastQver == 64 && t < -5) {
				fastQver = 33;
				fqSolexaFmt = false;
				if (print) { cerr << "Detected qualities incompatible with +64; switching to Phred+33.\n"; print = false; }
			}
			else { if (print) { cerr << "Unusually low quality score (" << t << "); setting to 0.\n"; print = false; } }
		}
		t = 0;
	}
	if (minQScore > t) { minQScore = t; }
	else if (maxQScore < t) { maxQScore = t; }
	return t;
}

bool InputStreamer::checkInFileStatus() {
	for (uint i = 0; i < 3; i++) {
		if (isFasta) {
			if (fasta_istreams[i] != NULL && !fasta_istreams[i]->eof()) { return true; }
		} else {
			if (fastq_istreams[i] != NULL && !fastq_istreams[i]->eof()) { return true; }
		}
	}
	return false;
}

void InputStreamer::allStreamReset() {
	std::unique_lock<std::shared_mutex> lock(protect);
	resetStats();
	cdbg("Resetting input streams");
	for (uint i = 0; i < 3; i++) {
		if (fasta_istreams[i] != NULL) { fasta_istreams[i]->reset(); }
		if (quality_istreams[i] != NULL ) { quality_istreams[i]->reset(); }
		if (fastq_istreams[i] != NULL ) { fastq_istreams[i]->reset(); }
	}
	cdbg("Done");
}

void InputStreamer::allStreamClose(){
	for (uint i = 0; i < 3; i++){
		if (fasta_istreams[i] != NULL) { delete fasta_istreams[i]; fasta_istreams[i] = NULL;}
		if (quality_istreams[i] != NULL) { delete quality_istreams[i]; quality_istreams[i] = NULL; }
		if (fastq_istreams[i] != NULL) { delete fastq_istreams[i]; fastq_istreams[i] = NULL; }
	}
	hasMIDs = false;

 if (!isFasta && minQScore < SCHAR_MAX) { /* maxminQualWarns_fq() moved/removed with utilities */ }
}

void InputStreamer::jumpToNextDNA(bool& stillMore, int pos) {
	if (isFasta) {
		stillMore = read_fasta_entry(fasta_istreams[pos], quality_istreams[pos], dnaTemp1[pos], dnaTemp2[pos], lnCnt[pos]);
		dnaTemp1[pos] = dnaTemp2[pos];
		dnaTemp2[pos] = make_shared <DNA>("", "");
	} else {
		fastq_istreams[pos]->jumpLines(4);
		if (fastq_istreams[pos]->eof()) { stillMore = false; }
	}
}

vector<shared_ptr<DNA>> InputStreamer::getDNApairs() {
	// One critical section owns the full R1/R2/MID set so concurrent consumers
	// cannot interleave mates from adjacent records.
	std::unique_lock<std::shared_mutex> lock(protect);
	vector<shared_ptr<DNA>> ret(3, nullptr);
	ret[0] = getDNAUnlocked(0);
	if (numPairs > 1) {
		ret[1] = getDNAUnlocked(1);
	}
	if (hasMIDs) {
		ret[2] = getDNAUnlocked(2);
	}

	const bool allEnded = ret[0] == nullptr && (numPairs <= 1 || ret[1] == nullptr) &&
		(!hasMIDs || ret[2] == nullptr);
	if (allEnded) { return ret; }
	if (ret[0] == nullptr || (numPairs > 1 && ret[1] == nullptr) || (hasMIDs && ret[2] == nullptr)) {
		throw std::runtime_error("Input streams ended at different records while reading a synchronized read set");
	}
	if (!ret[0]->seal() || ret[0]->isEmpty()) { ret[0] = nullptr; }
	if (numPairs > 1) {
		if (!ret[1]->seal() || ret[1]->isEmpty()) { ret[1] = nullptr; }
		if (ret[0] && ret[1] && !ret[0]->sameHeadPosFree(ret[1])) {
			throw std::runtime_error("Paired input headers do not match: '" + ret[0]->getShortId() +
				"' versus '" + ret[1]->getShortId() + "'");
		}
	}
	if (hasMIDs) {
		if (!ret[2]->seal() || ret[2]->isEmpty()) { ret[2] = nullptr; }
		if (ret[0] && ret[2] && !ret[0]->sameHeadPosFree(ret[2])) {
			throw std::runtime_error("R1 and MID input headers do not match: '" + ret[0]->getShortId() +
				"' versus '" + ret[2]->getShortId() + "'");
		}
		if (ret[2]) { ret[2]->setMIDseq(true); }
	}
	if (!ret[0] || (numPairs > 1 && !ret[1]) || (hasMIDs && !ret[2])) {
		throw std::runtime_error("Malformed record encountered in a synchronized input set");
	}

	N_DNApairsReturned++;
	
	return ret;
}


string InputStreamer::current_infiles() {
	string ret;
	for (int pos = 0; pos < 3; ++pos) {
		ifbufstream* stream = isFasta ? fasta_istreams[pos] : fastq_istreams[pos];
		if (!stream) { continue; }
		if (!ret.empty()) { ret += ", "; }
		ret += stream->getInFile();
	}
	return ret;
}

bool InputStreamer::getDNAlines(multi_tmp_lines* tmpO, int blocks, bool MIDuse,bool safe) {

	


	std::unique_lock<std::shared_mutex> inputGuard(protect, std::defer_lock);
	if (safe) { inputGuard.lock(); }

	assert(tmpO->size() == blocks);
	size_t k(0);
	for (k = 0; k < blocks; k++) {
		// MID/index reads are synchronized metadata; read limits count biological mates only.
		const uint64_t recordsNeeded = 1u + (numPairs > 1 ? 1u : 0u);
		if (_globalMaxRdsRead > 0 &&
			(_globalRdsRead >= _globalMaxRdsRead ||
			 _localRdsRead > _globalMaxRdsRead - _globalRdsRead ||
			 recordsNeeded > _globalMaxRdsRead - _globalRdsRead - _localRdsRead)) {
			tmpO->setSize(k);
			return false;
		}
		const bool b1 = this->getDNAlines(tmpO->tmp[k][0], 0);
		bool b2 = true;
		if (numPairs > 1) {
			b2 = this->getDNAlines(tmpO->tmp[k][1], 1);
		}
		bool b3 = true;
		if (MIDuse) {
			b3 = this->getDNAlines(tmpO->tmp[k][2], 2);
		}

		const bool partial1 = !b1 && recordBufferHasData(tmpO->tmp[k][0]);
		const bool partial2 = numPairs > 1 && !b2 && recordBufferHasData(tmpO->tmp[k][1]);
		const bool partial3 = MIDuse && !b3 && recordBufferHasData(tmpO->tmp[k][2]);
		const bool allEnded = !b1 && (numPairs <= 1 || !b2) && (!MIDuse || !b3);
		if (allEnded && !partial1 && !partial2 && !partial3) {
			tmpO->setSize(k);
			return false;
		}
		if (!b1 || (numPairs > 1 && !b2) || (MIDuse && !b3) || partial1 || partial2 || partial3) {
			tmpO->setSize(k);
			const std::string message =
				"Input streams are truncated or ended at different records: " + current_infiles();
			if (DieOnError) {
				throw std::runtime_error(message);
			}
			IO_Error(message + "; stopping after recovered records");
			return false;
		}

		if (numPairs > 1) {
			validateSynchronizedHeaders(tmpO->tmp[k][0], tmpO->tmp[k][1], "R2");
		}
		if (MIDuse) {
			validateSynchronizedHeaders(tmpO->tmp[k][0], tmpO->tmp[k][2], "MID");
		}
		_localRdsRead += recordsNeeded;
	}
	return true;
}
bool InputStreamer::getDNAlines(vector<string>& ret, int pos) {
	
	//vector<string> ret(4,"");
	bool stillMore = true;
	if (pos < 0 || pos > 2) {
		return false;
	}
	if (pos == 1 && numPairs <= 1) {
		return false;
	}
	bool repairInStream(false);
	if (isFasta) {
		if (ret.size() != 3) { ret.resize(3); }
	}
	else {
		if (ret.size() < 4) { ret.resize(4); }
	}
	std::fill(ret.begin(), ret.end(), std::string());
	//vector<string>tmpLines(4, "");
	if (isFasta) {//get DNA from fasta + qual_ files
		stillMore = readFastaRecord(fasta_istreams[pos], quality_istreams[pos], ret);
		if (stillMore) { ++lnCnt[pos]; }
	} else {// fastq format
		if (fastq_istreams[pos] == nullptr || fastq_istreams[pos]->eof()) {
			return false;
		}
		//fqRead
		stillMore = fastq_istreams[pos]->get4lines(ret);
		if (!stillMore && recordBufferHasData(ret)) {
			const std::string message = "Truncated FASTQ record in " + fastq_istreams[pos]->getInFile();
			if (DieOnError) {
				throw std::runtime_error(message);
			}
			IO_Error(message + "; discarding incomplete record");
			ret.assign(4, std::string());
		}
		if (stillMore) {
			std::string malformedMessage;
			if (ret[0].empty() || ret[0][0] != '@') {
				malformedMessage = "Malformed FASTQ header";
			}
			else if (ret[2].empty() || ret[2][0] != '+') {
				malformedMessage = "Malformed FASTQ plus line for " + ret[0];
			}
			else if (ret[1].size() != ret[3].size()) {
				malformedMessage = "FASTQ sequence/quality length mismatch for " + ret[0];
			}
			if (!malformedMessage.empty()) {
				if (DieOnError) {
					throw std::runtime_error(malformedMessage);
				}
				IO_Error(malformedMessage + "; discarding malformed record");
				ret.assign(4, std::string());
				stillMore = false;
			}
		}
		/*for (int xx = 0; xx < 4; xx++) {
			stillMore = fastq_istreams[pos]->getline(ret[xx]);
		}*/
 if (stillMore) {
			lnCnt[pos] += 4;
		}
 if (stillMore && fastQver == 0) {
		 fastQver = auto_fq_version(ret);
		 QverSet = true;			
		}
	}
	if (stillMore) {
		pairs_read[pos]++;
	}
	return stillMore;
}

shared_ptr<DNA> InputStreamer::getDNA(int pos) {
	std::unique_lock<std::shared_mutex> lock(protect);
	return getDNAUnlocked(pos);
	}

shared_ptr<DNA> InputStreamer::getDNAUnlocked(int pos){
	//if (sync) {
	//	while (desync(pos)) {
	//		jumpToNextDNA(stillMore, pos);
	//	}
	//}
	bool stillMore = true;
	if (pos == 1 && numPairs <= 1) {
		return nullptr;
	}
	shared_ptr<DNA> ret(nullptr);
	bool corrupt(false); //corrupt state isn't implemented for fnaread
	bool repairInStream(false);
	if (isFasta) {//get DNA from fasta + qual_ files
		vector<string>tmpStr(3, "");
		stillMore = readFastaRecord(fasta_istreams[pos], quality_istreams[pos], tmpStr);
		if (!stillMore) { return nullptr; }
		++lnCnt[pos];
		++pairs_read[pos];

		ret = str2DNA(tmpStr, keepPairHD, fastQver,pos);

		if (!stillMore || fasta_istreams[pos]->eof() || (!*(fasta_istreams[pos]))) {
			if (ret != nullptr) { if (!ret->seal() || ret->isEmpty()) { ret = nullptr; } } //delete ret;
			stillMore = false; 
		}
		else if (ret == nullptr || !ret->seal() || ret->isEmpty()) {
			corrupt = true;
		}
	}
	else { //fqRead
		if (fastq_istreams[pos] == nullptr || fastq_istreams[pos]->eof()) {
			stillMore = false;
			return nullptr;
		}
		vector<string>tmpLines(4, "");
		stillMore = getDNAlines(tmpLines,pos);
		if (stillMore) {
			ret = str2DNA(tmpLines, keepPairHD, fastQver,pos);
		}

		if (!stillMore || fastq_istreams[pos]->eof()) {
			if (ret != nullptr) { if (!ret->seal() || ret->isEmpty()) { ret = nullptr; } }
			stillMore = false;
		} else if (ret == nullptr || !ret->seal() || ret->isEmpty()) {
			corrupt = true;
		}

		if (ret != nullptr && !keepPairHD) {//better cut in early
			ret->setHeader(ret->getPositionFreeId());
		}
		if (corrupt) {
			ret = nullptr;
			repairInStream = true;
		}
	}
	N_DNAreturned++;

	return ret;
}

void InputStreamer::openMIDseqs(string p,string in){
	if (in==""){return;}
	
	if (fastq_istreams[2] != NULL){
		throw std::runtime_error("MID file was already initialized");
	}
#ifdef DEBUG
	cerr << "Open Mid sequence_ file" << endl;
#endif

	string file_type = "MID specific fastq";
	string tmp = (p + in);

	fastq_istreams[2] = DBG_NEW ifbufstream(tmp, (size_t)round(INPUT_BUFFER_SIZE*0.6),doTIO);
	if (fastq_istreams[2]->eof()) {
		cerr << "\nCouldn't find " << file_type << " " << tmp << "!\n Aborting..\n";		exit(4);
	}
	hasMIDs=true;
}

bool InputStreamer::detectAndSwitchFileFormat(ifbufstream* stream, const std::string& filepath) {
	if (stream == nullptr || stream->eof()) {
		return false; // Can't detect format if stream is not available
	}

	// Peek at the first character to determine actual format
	int firstChar = stream->peekFirstChar();

	if (firstChar == EOF) {
		return false; // Empty file
	}

	bool detectedFasta = (firstChar == '>');
	bool detectedFastq = (firstChar == '@');

	// Check if detected format differs from assumed format
	if (detectedFasta && !isFasta) {
		throw std::runtime_error("Input format mismatch: FASTA data supplied through FASTQ input: " + filepath);
	}
	else if (detectedFastq && isFasta) {
		throw std::runtime_error("Input format mismatch: FASTQ data supplied through FASTA input: " + filepath);
	}

	return false; // No format switch needed
}

bool InputStreamer::setupFastq(string path, string fileS, int& pairs, string subsPairs,
	bool simu, bool verbose) {
	allStreamClose();
	minQScore = SCHAR_MAX;
	maxQScore = -1;
	fqSolexaFmt = false;
	resetStats();
	vector<string> tfas = splitByCommas(fileS);
	if (pairs == -1) {
		pairs = (int)tfas.size();
		if (pairs > 1) { cerr << "Paired input (\",\" separated) detected\n"; }
	}
	numPairs = pairs;
	string p1(""), p2(""), midp("");
	string xtraMsg = "";
	if (getCurFileN() > 0) {
		xtraMsg = " " + itos(getCurFileN()) + " of " + itos(totalFileNumber);
		if (BCnumber > 1) {
			xtraMsg = ", looking for " + itos(BCnumber) + "BCs.\n";
		}
		else {
			xtraMsg = "";
		}
	}

	if (tfas.size() != (uint)pairs && subsPairs == "") {
		cerr << "Unequal number of files (" << tfas.size() << ") and option-set paired files (" << pairs << ").\n Aborting...\n"; exit(52);
	}
	string file1("");
	if (tfas.size() == 3) {
		if (tfas.size() != 3) { cerr << "Could not detect 3 input files in string\n" << fileS << "\n Aborting.." << endl; exit(76); }
		midp = path + tfas[1];
		p1 = path + tfas[0];
		p2 = path + tfas[2];
		file1 = tfas[0];
		//		cerr << p1 << " + " << p2 << " and " << midp << endl;
	}
	else if (tfas.size() == 2) {
		p1 = path + tfas[0];
		p2 = path + tfas[1];
		file1 = tfas[0];
	}
	else if (tfas.size() == 1) {
		p1 = path + fileS;
		file1 = fileS;
		//cerr << "Reading fastq " << p1 << endl;
	}
	if (subsPairs == "1") {
		p2 = "";
	}
	else 	if (subsPairs == "2") {
		p1 = p2; p2 = "";
	}
	if (!subsPairs.empty()) {
		pairs = 1;
		numPairs = 1;
		midp.clear();
	}

	if (simu) {

		return fileExists(p1) && fileExists(p2) && fileExists(midp);
	}


	if (verbose && !p1.empty() && !p2.empty()) {
		if (!midp.empty()) {
			//cerr << "Reading paired fastq + MID file" << xtraMsg<<"."<<endl;
		}
		else {
			//cerr << "Reading paired fastq" << xtraMsg << "." << endl;
			//cerr << p1 << " + " << p2 << endl;
		}
	}
	else {
		//cerr << "Reading fastq" << xtraMsg << "." << endl;
		//cerr << p1 << endl;
	}
	if (verbose) {
		cerr << "At " << file1 << ":\n";
	}


	bool suc = setupFastq_2(p1, p2, midp);
	//read progress bar setup
	//_measure(*fastq_istreams[0]);
	//allStreamClose();
	//setupFastq_2(p1, p2, midp);
	return suc;
}


bool InputStreamer::setupFastaQual2(string p1, string p2, string file_type) {
	if (p1.empty()) return false;
	int pos = 0;
	if (fasta_istreams[0] != nullptr) pos = 1;

	fasta_istreams[pos] = DBG_NEW ifbufstream(p1, INPUT_BUFFER_SIZE, doTIO);
	if (fasta_istreams[pos]->eof()) {
		cerr << "\nWarning: Could not open or empty " << file_type << " " << p1 << " — skipping this input.\n";
		fasta_istreams[pos] = nullptr;
		return false;
	}

	// Detect and switch format if needed
	detectAndSwitchFileFormat(fasta_istreams[pos], p1);

	if (!p2.empty()) {
		quality_istreams[pos] = DBG_NEW ifbufstream(p2, INPUT_BUFFER_SIZE, doTIO);
		if (quality_istreams[pos]->eof()) {
			cerr << "\nWarning: Could not open or empty quality file " << p2 << " — continuing without quality stream.\n";
			quality_istreams[pos] = nullptr;
			qualAbsent = true;
		}
		else {
			qualAbsent = false;
		}
	}
	else {
		quality_istreams[pos] = nullptr;
		qualAbsent = true;
	}

	return true;
}
// Prints out the progress bar
inline void InputStreamer::_print(int cur, float prog) {

	std::cerr << std::fixed << std::setprecision(2)
		<< "\r [" << std::string(cur, '#')
		<< std::string(_max + 1 - cur, ' ') << "] " << 100 * prog << "%";

	if (prog == 1) std::cerr << std::endl;
	else std::cerr.flush();

}
inline bool InputStreamer::_drawbar(istream & tar) {
 if (_fileLength <= 0) { return false; }
 int pos((int) tar.tellg());
 float prog(pos / float(_fileLength)); // percentage of infile already read
 if (pos == -1 || prog > 1.f) {
 _print(_max + 1, 1);
 _fileLength = 0;
 return false;
 }
 
 // Number of #'s as function of current progress "prog"
 int cur((int) ceil(prog * (float) _max));
 if (_last != cur) _last = cur, _print(cur, prog);
 return prog == 1.f;
}

inline void InputStreamer::_measure(istream& tar) {
	tar.seekg(0, ios_base::end);
	_fileLength = (int) tar.tellg();
	tar.seekg(0, ios_base::beg);
	tar.clear();
}

bool InputStreamer::setupFastq_2(string p1, string p2, string midp) {
	string file_type = "";
	size_t bufS = INPUT_BUFFER_SIZE;

	if (!p1.empty()) {
		file_type = "fastq file 1";
		fastq_istreams[0] = DBG_NEW ifbufstream(p1, (size_t)round(bufS * 0.8), doTIO);
		if (fastq_istreams[0]->eof()) {
			cerr << "\nWarning: Could not open or empty " << file_type << " " << p1 << " — skipping this input.\n";
			fastq_istreams[0] = nullptr;
			return false;
		}
		// Detect and switch format if needed
		detectAndSwitchFileFormat(fastq_istreams[0], p1);
	}
	if (!p2.empty()) {
		file_type = "fastq file 2";
		fastq_istreams[1] = DBG_NEW ifbufstream(p2, (size_t)round(bufS * 1.2), doTIO);
		if (fastq_istreams[1]->eof()) {
			cerr << "\nWarning: Could not open or empty " << file_type << " " << p2 << " — skipping this input.\n";
			fastq_istreams[1] = nullptr;
			return false;
		}
		// Detect and switch format if needed
		detectAndSwitchFileFormat(fastq_istreams[1], p2);
	}
	if (!midp.empty()) {
		this->openMIDseqs("", midp);
	}
	return true;
}

string InputStreamer::setupInput(string path, int t, const string& uniqueFastxFile,
	filesStr& files, int& paired, string onlyPair,
	string& mainFilename, bool simulate) {
	string mainFilepath("");
	vector<string> fastqFiles = files.FastqF;
	vector<string> fastaFiles = files.FastaF;
	vector<string> qualityFiles = files.QualF;
	vector<string> midFiles = files.MIDfq;

	if (isFasta) {
		if (fastaFiles[t] != uniqueFastxFile) {
			cerr << "Error in matching FASTA target filenames.\n";
			exit(11);
		}
		if (!this->setupFastaQual(path, fastaFiles[t], qualityFiles[t], paired, onlyPair)) {
			mainFilename = fastaFiles[t];
			return "";
		}
		mainFilepath = path + fastaFiles[t];
		mainFilename = fastaFiles[t];
	}
	else {
		if (fastqFiles[t] != uniqueFastxFile) {
			cerr << "Error in matching target filenames.\n";
			exit(11);
		}
		if (!this->setupFastq(path, fastqFiles[t], paired, onlyPair, simulate, !simulate)) {
			mainFilename = fastqFiles[t];
			return "";
		}
		mainFilepath = path + fastqFiles[t];
		mainFilename = fastqFiles[t];
	}
	if ((size_t)t < midFiles.size()) {
		this->openMIDseqs(path, midFiles[t]);
	}
	return mainFilepath;
}

void InputStreamer::setupFna(string in) {
	int paired = 1;
	string fmt = detectSeqFmt(in);
	if (fmt == "-i_fastq") {
		setupFastq("", in, paired, "", false, false);
	}
	else {
		setupFastaQual("", in, "", paired, "", false);
	}
}

qual_score InputStreamer::auto_fq_version() {
	if (fastQver != 0) {
		return fastQver;
	}

	return (qual_score)auto_fq_version(minQScore, maxQScore);
}

qual_score InputStreamer::auto_fq_version(qual_score minQScoreIn, qual_score maxQScoreIn) {
 // minQScoreIn/maxQScoreIn are raw ASCII quality characters from FASTQ lines.
	// Prefer Phred+33 unless there is strong evidence for legacy +64.
	if ((int)minQScoreIn < 59) { return 33; }
	if ((int)maxQScoreIn > 74) { return 64; }
	return 33;
}
qual_score InputStreamer::auto_fq_version(const vector<string>& ret) {
	qual_score minAscii = SCHAR_MAX;
	qual_score maxAscii = 0;
	for (unsigned char ch : ret[3]) {
		qual_score q = static_cast<qual_score>(ch);
		if (q < minAscii) { minAscii = q; }
		if (q > maxAscii) { maxAscii = q; }
	}

	if (minAscii < minQScore) { minQScore = minAscii; }
	if (maxAscii > maxQScore) { maxQScore = maxAscii; }

	const bool modernIlluminaHeader =
		(ret[0].find(" 1:N:") != string::npos) ||
		(ret[0].find(" 2:N:") != string::npos) ||
		(ret[0].find(" 1:Y:") != string::npos) ||
		(ret[0].find(" 2:Y:") != string::npos);

	if (modernIlluminaHeader) {
		fastQver = 33;
		fqSolexaFmt = false;
	}
	else {
		fastQver = static_cast<qual_score>(auto_fq_version(minQScore, maxQScore));
	}
	return fastQver;
}


bool InputStreamer::setupFastaQual(string path, string fasta, string qual, int& pairs, string subsPairs, bool simu) {
	allStreamClose();
	resetStats();
	vector<string> fna = splitByCommas(fasta);
	vector<string> qua = qual.empty() ? vector<string>(0) : splitByCommas(qual);

	if (pairs == -1) {
		pairs = (int)fna.size();
	}
	if (pairs > 1 && subsPairs.empty() && fna.size() < static_cast<size_t>(pairs)) {
		throw std::runtime_error("Paired FASTA input is missing one or more mate files");
	}
	if (!qua.empty() && pairs > 1 && subsPairs.empty() && qua.size() < static_cast<size_t>(pairs)) {
		throw std::runtime_error("Paired QUAL input is missing one or more mate files");
	}
	numPairs = pairs;

	fastaFilepathTemp[0] = (fna.empty() ? "" : path + fna[0]);
	fastaFilepathTemp[1] = (fna.size() > 1 ? path + fna[1] : "");
	qualityFilepathTemp[0] = (qua.size() > 0 ? path + qua[0] : "");
	qualityFilepathTemp[1] = (qua.size() > 1 ? path + qua[1] : "");

	if (subsPairs == "1") {
		fastaFilepathTemp[1].clear();
		qualityFilepathTemp[1].clear();
	}
	else if (subsPairs == "2") {
		fastaFilepathTemp[0] = fastaFilepathTemp[1];
		qualityFilepathTemp[0] = qualityFilepathTemp[1];
	}
	if (!subsPairs.empty()) {
		pairs = 1;
		numPairs = 1;
		fastaFilepathTemp[1].clear();
		qualityFilepathTemp[1].clear();
	}

	if (simu) {
		bool ok = !fastaFilepathTemp[0].empty() && fileExists(fastaFilepathTemp[0]);
		if (!qualityFilepathTemp[0].empty()) ok = ok && fileExists(qualityFilepathTemp[0]);
		if (pairs == 2 && !fastaFilepathTemp[1].empty()) ok = ok && fileExists(fastaFilepathTemp[1]);
		if (pairs == 2 && !qualityFilepathTemp[1].empty()) ok = ok && fileExists(qualityFilepathTemp[1]);
		return ok;
	}

	bool suc = setupFastaQual2(fastaFilepathTemp[0], qualityFilepathTemp[0]);
	if (pairs == 2 && !fastaFilepathTemp[1].empty()) {
		suc = suc && setupFastaQual2(fastaFilepathTemp[1], qualityFilepathTemp[1], "fasta file 2");
	}
	return suc;
}