Class SequenceData

blacklist.cpp:24–141

class SequenceData {
	public:
		std::vector<int> binsInput;
		std::vector<int> binsSpikes;
		std::vector<int> binsMultimapping;
		std::vector<int> binsTemp;
		int totalReads;
		int zeroMulti = 0;

		SequenceData(std::string bamFile, std::string bamIndexFile) {
			this->bamFile = bamFile;
			this->bamIndexFile = bamIndexFile;
		}

		SequenceData() {
		}

		void setBamFile(std::string bamFile) {
			this->bamFile = bamFile;
		}

		std::string getBamFile() {
			return bamFile;
		}

		void setBamIndexFile(std::string bamIndexFile) {
			this->bamIndexFile = bamIndexFile;
		}

		void getInputBins(std::vector<uint8_t>& mappability, int binSize, int binOverlap, std::string refName)
		{
			this->binsInput.clear();
			this->binsSpikes.clear();
			this->binsMultimapping.clear();
			this->binsTemp.clear();

			BamTools::BamReader reader;
			BamTools::BamAlignment al;
			std::vector<int> tempCounts (mappability.size(), 0);
			std::vector<int> multiCounts (mappability.size(), 0);
			totalReads = 0;
			int testCntr = 0;
			int zmultiCntr = 0;

			reader.Open(this->bamFile);
			reader.OpenIndex(this->bamIndexFile); // Significant speedup here

 		  	// This restricts output to a specific chromosome
			if(reader.GetReferenceID(refName) == -1) {
				//cout << "No Reads to load!\n"; //DEBUG
			} else {
			  	reader.SetRegion(reader.GetReferenceID(refName), 0, reader.GetReferenceID(refName), mappability.size());

 				 // Keep vector of locations and counts of reads
 			  	// filter on mappability based on each read length
				while(reader.GetNextAlignmentCore(al)) {
			        	if(mappability[al.Position] > 0 && mappability[al.Position] <= al.Length) { // Filter mappability
			        		tempCounts[al.Position] = tempCounts[al.Position] + 1; // Count reads at each position that are uniquely mappable