Slice Tools Home SourceForge Page libSlice Home Modules Files Examples Additional Information

SliceData.cc Go to the documentation of this file. #include "SliceData.hh" #include <stdlib.h> // for atol namespace libSlice { using std::string; SliceData::SliceData(const ReadManager & readManager, const char * index, const char * gindex, const char * existing, long sliceRangeIndex, long sliceRangeGindex) : m_readManager(readManager) { m_slice = NULL; m_highQualityThreshold = 0; m_nucData = ""; m_readIdData = ""; m_qualValData = ""; m_errorData = ""; m_errorCode = SLICE_NOERROR; m_errorOperation = STROP_NONE; m_doAmbiguity = 0; m_scoringModel = 0; m_index = 0; m_gindex = 0; m_indexstr = "Unknown"; m_existing = (existing) ? existing : "-"; if (index) { m_index = atol(index); m_index += sliceRangeIndex; m_indexstr = index; } if (gindex) { m_gindex = atol(gindex); m_gindex += sliceRangeGindex; } double normal = 0.2; setDistributionTable(normal, normal, normal, normal, normal); } SliceData::~SliceData() { if (m_slice) { if (m_slice->rc) delete [] m_slice->rc; if (m_slice->qv) delete [] m_slice->qv; if (m_slice->bc) delete [] m_slice->bc; delete m_slice; m_slice = NULL; } } void SliceData::setDistributionTable(double freqA, double freqC, double freqG, double freqT, double freqGap) { m_distribution.freqA = freqA; m_distribution.freqC = freqC; m_distribution.freqG = freqG; m_distribution.freqT = freqT; m_distribution.freqGap = freqGap; } void SliceData::setHighQualityThreshold(int highQualityThreshold) { m_highQualityThreshold = highQualityThreshold; } void SliceData::addNuc(const char * nucData) { if (nucData) m_nucData.append(nucData); } void SliceData::addQualVal(const char * qualVal) { if (qualVal) m_qualValData.append(qualVal); } void SliceData::addReadId(const char * readId) { if (readId) m_readIdData.append(readId); } void SliceData::setConsensus(char consensus) { m_existing = consensus; if (m_slice) { m_slice->c = consensus; } } // TODO: Investigate strtok char * SliceData::foreach(const string & data, int numresults, StringOperation operation, const char * underflowError, const char * overflowError) { const char * databuf = data.c_str(); char * retval = new char [numresults+1]; retval[numresults] = '\0'; string::size_type startpos = 0; string::size_type endpos = 0; for (int i=0; i < numresults; i++) { endpos = data.find(' ', startpos); if (endpos == data.npos) { if (i != (numresults-1)) { m_errorOperation = operation; m_errorCode = SLICE_ERR_UNDERFLOW; throw Exception(underflowError, "SliceData::foreach"); } endpos = data.length() + 1; } char elementResult = '\0'; if (operation == STROP_GET_QUAL_VAL) { elementResult = (char) (atoi(databuf + startpos) & 0xff); } else if (operation == STROP_IS_REVERSE_COMPLIMENT) { const string & element = data.substr(startpos, (endpos - startpos)); const ReadData * rd = m_readManager.find(element.c_str()); if (!rd) { m_errorOperation = operation; m_errorCode = SLICE_ERR_NO_READID; m_errorData = element; throw Exception("ReadID " + element + " not found", "SliceData::foreach"); } elementResult = rd->isReverseCompliment(); } retval[i] = elementResult; startpos = endpos + 1; } if (numresults && (endpos != data.length()+1)) { m_errorOperation = operation; m_errorCode = SLICE_ERR_OVERFLOW; throw Exception(overflowError, "SliceData::foreach"); } return retval; } char * SliceData::getQualVals(int len) { // There better only be number data in m_qualValData if (m_qualValData.find_first_not_of("0123456789 ") != m_qualValData.npos) { m_errorCode = SLICE_ERR_QV_NONNUMERIC; m_errorOperation = STROP_GET_QUAL_VAL; throw Exception("Invalid Quality Value in slice " + m_indexstr, "SliceData::getQualVals"); } return foreach(m_qualValData, len, STROP_GET_QUAL_VAL, "Number of Bases > Number of QualVals", "Number of QualVals > Number of Bases"); } char * SliceData::getReverseCompliments(int len) { return foreach(m_readIdData, len, STROP_IS_REVERSE_COMPLIMENT, "Number of Bases > Number of ReadIds", "Number of ReadIds > Number of Bases"); } char SliceData::getExisting() const { char retval = '-'; if (m_existing.length() == 1) { retval = m_existing[0]; } else { //TODO: Decide what validation should be done (if any) } return retval; } unsigned short SliceData::getDepth() const { return m_nucData.length(); } char * SliceData::getBases(int len) const { char * retval = new char [len + 1]; strcpy(retval, m_nucData.c_str()); return retval; } void SliceData::flattenSlice() { if (!m_slice) { m_slice = new libSlice_Slice; memset (m_slice, 0, sizeof(libSlice_Slice)); int depth = getDepth(); m_slice->dcov = depth; m_slice->c = getExisting(); m_slice->bc = getBases(depth); m_slice->qv = getQualVals(depth); m_slice->rc = getReverseCompliments(depth); } } char SliceData::getBase(unsigned short index) { flattenSlice(); if (index > getDepth()) { throw(Exception("Index > Depth", "getBase")); } return m_slice->bc[index]; } char SliceData::getQualVal(unsigned short index) { flattenSlice(); if (index > getDepth()) { throw(Exception("Index > Depth", "getQualVal")); } return m_slice->qv[index]; } //TODO: Refactor to utilize foreach ReadData * SliceData::getRead(unsigned short index) const { if (index > getDepth()) { throw(Exception("Index > Depth", "getRead")); } string::size_type startpos = 0; string::size_type endpos = 0; while (index) { startpos = m_readIdData.find(' ', startpos); startpos++; if (startpos == m_readIdData.npos) { throw Exception("Index > Number of ReadIds", "getReadId"); } index--; } endpos = m_readIdData.find(' ', startpos); if (endpos == m_readIdData.npos) { endpos = m_readIdData.length(); } const string & element = m_readIdData.substr(startpos, (endpos - startpos)); ReadData * rd = m_readManager.find(element.c_str()); if (!rd) { throw Exception("ReadID " + element + " not found", "SliceData::getRead"); } return rd; } /* * This turned out to be slower than getRead * const ReadData * SliceData::getRead2(unsigned short index) { int depth = getDepth(); if (index > depth) { throw(Exception("Index > Depth", "getRead2")); } if (m_readVector.empty()) { if (depth > 1) { int len = m_readIdData.length(); char buffer[len+1]; strcpy(buffer, m_readIdData.c_str()); char * start = NULL; for (int i = 0; i < len; i++) { if (buffer[i] != ' ' && start == NULL) { start = &buffer[i]; } else if (buffer[i] == ' ') { buffer[i] = '\0'; if (start) { // Start points to the beginning of the null terminated string to // extract m_readVector.push_back(m_readManager.find(start)); start = NULL; } } } } else { m_readVector.push_back(m_readManager.find(m_readIdData.c_str())); } if (m_readVector.size() != depth) { throw Exception("Number of reads doesn't match number of bases", "SliceData::getRead2"); } } return m_readVector[index]; } */ long SliceData::getIndex() const { return m_index; } long SliceData::getGindex() const { return m_gindex; } char SliceData::getConsensus() { flattenSlice(); return m_slice->c; } ConsensusData SliceData::doCalculations(bool recallEmpty) { flattenSlice(); int errCode; libSlice_Consensus resultConsensus; // This should probably never be set if (recallEmpty) { libSlice_setRecallEmpty(recallEmpty); } // Don't do ambiguity here, get that separately if ((errCode = libSlice_getConsensusParam(m_slice, &resultConsensus, &m_distribution, m_highQualityThreshold, 0))) { m_errorCode = SLICE_ERR_GETCONSENSUS; throw Exception("Error Calculating consensus", "_getConsensus"); } if (m_doAmbiguity == 2) { char origFlags = resultConsensus.ambiguityFlags; libSlice_updateAmbiguity(m_slice, &resultConsensus, m_highQualityThreshold); if (0 && origFlags != resultConsensus.ambiguityFlags) { cerr << "Annotation: Updating Ambiguity " << libSlice_convertAmbiguityFlags(origFlags) << " to " << libSlice_convertAmbiguityFlags(resultConsensus.ambiguityFlags) << " at gindex " << m_gindex << endl; } } else if (m_doAmbiguity == 3) { char origFlags = resultConsensus.ambiguityFlags; libSlice_updateAmbiguityConic(m_slice, &resultConsensus, m_highQualityThreshold); if (0 && origFlags != resultConsensus.ambiguityFlags) { cerr << "Conic: Updating Ambiguity " << libSlice_convertAmbiguityFlags(origFlags) << " to " << libSlice_convertAmbiguityFlags(resultConsensus.ambiguityFlags) << " at gindex " << m_gindex << endl; } } char ambiguityCode; if (m_slice->dcov == 0 && !recallEmpty) { // Don't update the slice ambiguityCode = m_slice->c; } else { ambiguityCode = libSlice_convertAmbiguityFlags(resultConsensus.ambiguityFlags); if (m_doAmbiguity) { m_slice->c = ambiguityCode; } else { m_slice->c = resultConsensus.consensus; } } m_existing = (string) "" + m_slice->c; char qualityClass = libSlice_getConsQC(m_slice, m_highQualityThreshold); return ConsensusData(&resultConsensus, ambiguityCode, qualityClass, m_doAmbiguity); } void SliceData::print(ostream & out, const string & prefix) { flattenSlice(); out << prefix << "<Slice Index=\"" << m_index << "\" Gindex=\"" << m_gindex << "\" Existing=\"" << m_existing << "\" Depth=\"" << getDepth() << "\">" << endl; for (int i=0; i < m_slice->dcov; i++) { out << prefix << " <Base i=\"" << i << "\" bc=\"" << m_slice->bc[i] << "\" qv=\"" << (int) m_slice->qv[i] << "\" rc=\"" << (int) m_slice->rc[i] << "\" />" << endl; } out << prefix << "</Slice>" << endl; } void SliceData::printXML(ostream & out, const string & prefix) { out << prefix << "<Slice Index=\"" << m_index << "\" Gindex=\"" << m_gindex << "\" Existing=\"" << m_existing << "\">" << endl; out << prefix << " <Nuc>" << m_nucData << "</Nuc>" << endl; out << prefix << " <Qualval>" << m_qualValData << "</Qualval>" << endl; out << prefix << " <ReadID>" << m_readIdData << "</ReadID>" << endl; out << prefix << "</Slice>" << endl; } void SliceData::resetReadId() { m_readIdData = ""; if (m_slice) { delete m_slice; m_slice = NULL; } } void SliceData::resetQualVal() { m_qualValData = ""; if (m_slice) { delete m_slice; m_slice = NULL; } } void SliceData::resetNuc() { m_nucData = ""; if (m_slice) { delete m_slice; m_slice = NULL; } } void SliceData::resetIndex(long index, long gindex) { m_index = index; m_gindex = gindex; } const string & SliceData::getQualValData() const { return m_qualValData; } const string & SliceData::getReadIdData() const { return m_readIdData; } const string & SliceData::getNucData() const { return m_nucData; } void SliceData::printError(ostream & out) const { if (m_errorCode != SLICE_NOERROR) { out << " Error Slice index=" << m_index << " code=" << m_errorCode << " msg="; switch (m_errorCode) { case SLICE_NOERROR: out << "No Error"; break; case SLICE_ERR_QV_NONNUMERIC: out << "Non Numberic values in Quality Value"; break; case SLICE_ERR_UNDERFLOW: out << "Less data than required in operation"; break; case SLICE_ERR_OVERFLOW: out << "More data than required in operation"; break; case SLICE_ERR_NO_READID: out << "ReadId not found" << " readid=" << m_errorData; break; case SLICE_ERR_GETCONSENSUS: out << "Internal error in getConsensus"; break; case SLICE_ERR_CALCULATION_NOT_DONE: out << "Results requested before calculations performed"; break; }; if (m_errorCode == SLICE_ERR_OVERFLOW || m_errorCode == SLICE_ERR_UNDERFLOW) { switch (m_errorOperation) { case STROP_GET_QUAL_VAL: out << " operation=" << m_errorOperation << " op_str=" << "getQualval"; break; case STROP_IS_REVERSE_COMPLIMENT: out << " operation=" << m_errorOperation << " op_str=" << "getReverseCompliment"; break; case STROP_NONE: break; }; } // out << endl; } } int SliceData::getErrorCode() const { return m_errorCode; } int SliceData::getErrorOperation() const { return m_errorOperation; } void SliceData::setDoAmbiguity(int doAmbiguity) { m_doAmbiguity = doAmbiguity; } char SliceData::getBaseCompliment(char base) { return libSlice_getCompliment(base); } void SliceData::setScoringModel(int scoringModel) { m_scoringModel = scoringModel; } } ostream & operator<< (ostream & os, libSlice::SliceData & slice) { slice.print(os); return os; } // Convienent way to print *slices out ostream & operator<< (ostream & os, libSlice::SliceData * slice) { slice->print(os); return os; }