Loading...
Searching...
No Matches
ReferenceSequence.h
/*
* Copyright (C) 2010 Regents of the University of Michigan
*
* 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
* (at your option) any later version.
*
* 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 __REFERENCESEQUENCE_H
#define __REFERENCESEQUENCE_H
#include "BaseAsciiMap.h"
#include "Generic.h"
#include "PackedVector.h"
#include <algorithm>
#include <iostream>
//
// The namespace Sequence is for templated algorithms that are by and large
// independent of the underlying storage mechanism for the bases.
//
// They are written in such a way as to assume that the array operator []
// will return an ASCII representation of a base space nucleotide.
//
// In theory, this set of templates will work with a variety of combinations
// of means for representing bases - String, std::string, and others.
//
// The containers are expected to allow for an overidden [] operator,
// and provide a size() method to return the number of bases in the
// container.
//
// In containers where sequence data is placed, they must in addition
// have a clear() method as well as a push_back() method as done
// in std::string containers.
//
namespace Sequence {
//
// wordMatch(Sequnece, Index, Word) - compare a word to a sequence of bases
// Sequence is a generic container with a large set of bases
// Index is the starting point to start the comparison
// Word is the small sequence of bases to match
//
template<typename sequenceType, typename sequenceIndexType, typename wordType>
bool wordMatch(sequenceType &sequence, sequenceIndexType index, wordType &word)
{
if( (index + word.size()) >= sequence.size() ) return false;
for(size_t i = 0; i < word.size(); i++) {
if( sequence[index + i] != word[i]) return false;
}
return true;
}
//
// printNearbyWords(output, sequence, index, word, deviation) searches
// for 'deviation' bases on either side of the index into sequence
// and prints all occurrences where word appears.
//
template<typename sequenceType, typename sequenceIndexType, typename wordType>
void printNearbyWords(std::ostream &output, sequenceType &sequence, sequenceIndexType index, wordType &word, int deviation)
{
for (sequenceIndexType i = index - deviation; i < index + deviation; i++)
{
if (wordMatch(sequence, i, word))
{
output << "word '"
<< word
<< "' found "
<< i - index
<< " away from position "
<< index
<< "."
<< std::endl;
}
}
}
//
// getString(sequence, index, baseCount, word) - populate word with the 'baseCount'
// bases that occur at the 'index' starting position in sequence.
//
template<typename sequenceType, typename sequenceIndexType, typename wordType>
void getString(sequenceType &sequence, sequenceIndexType index, int baseCount, wordType &word)
{
word.clear();
for (sequenceIndexType i=0; i < (sequenceIndexType) baseCount; i++)
{
word.push_back(sequence[index + i]);
}
}
//
// getHighLightedString() is a debugging aid for printing "highlighted"
// subsets of bases, where the highlighting is done via turning the
// base into a lower case ASCII equivalent.
//
template<typename sequenceType, typename sequenceIndexType, typename wordType>
void getHighLightedString(
sequenceType &sequence,
sequenceIndexType index,
int baseCount,
wordType &word,
sequenceIndexType highLightStart,
sequenceIndexType highLightEnd)
{
word.clear();
for (sequenceIndexType i=0; i < (sequenceIndexType) baseCount; i++)
{
char base = sequence[index+i];
if(in(index+i, highLightStart, highLightEnd))
base = tolower(base);
word.push_back(base);
}
}
//
// printBaseContext() outputs a base at location 'index' along with 'baseCount'
// bases on either side of that base (default 30).
//
template<typename sequenceType, typename sequenceIndexType>
void printBaseContext(std::ostream &output, sequenceType &sequence, sequenceIndexType index, int baseCount = 30)
{
output << "index: " << index << std::endl;
for (sequenceIndexType i=index-baseCount; i<=index+baseCount; i++)
output << sequence[i];
output << std::endl;
for (sequenceIndexType i=index-baseCount; i<index; i++)
std::cout << " ";
output << "^";
output << std::endl;
}
template<typename sequenceType, typename sequenceIndexType>
void getMismatchHatString(sequenceType &sequence, sequenceIndexType location, std::string &result, std::string &read)
{
result = "";
for (uint32_t i=0; i < read.size(); i++)
{
if (read[i] == sequence[location+i])
result.push_back(' ');
else
result.push_back('^');
}
}
template<typename sequenceType, typename sequenceIndexType>
void getMismatchString(sequenceType &sequence, sequenceIndexType location, std::string &result, std::string &read)
{
result = "";
for (uint32_t i=0; i < read.size(); i++)
{
if (read[i] == sequence[location+i])
result.push_back(toupper(read[i]));
else
result.push_back(tolower(read[i]));
}
}
/// Return the mismatch count, disregarding CIGAR strings
///
/// \param read is the sequence we're counting mismatches in
/// \param location is where in the genmoe we start comparing
/// \param exclude is a wildcard character (e.g. '.' or 'N')
///
/// \return number of bases that don't match the reference, except those that match exclude
template<typename sequenceType, typename sequenceIndexType, typename readType>
int getMismatchCount(sequenceType &sequence, sequenceIndexType location, readType &read, char exclude='\0')
{
int mismatchCount = 0;
for (uint32_t i=0; i<read.size(); i++)
if (read[i]!=exclude) mismatchCount += read[i]!=sequence[location + i];
return mismatchCount;
}
/// brute force sumQ - no sanity checking
///
/// \param read shotgun sequencer read string
/// \param qualities phred quality string of same length
/// \param location the alignment location to check sumQ
template<typename sequenceType, typename sequenceIndexType, typename readType, typename qualityType>
int getSumQ(sequenceType &sequence, sequenceIndexType location, readType &read, qualityType &qualities)
{
int sumQ = 0;
for (uint32_t i=0; i<read.size(); i++)
sumQ += (read[i]!=sequence[location + i] ? (qualities[i]-33) : 0);
return sumQ;
};
template<typename sequenceType, typename sequenceIndexType, typename readType, typename qualityType>
sequenceIndexType simpleLocalAligner(
sequenceType &sequence,
sequenceIndexType index,
readType &read,
qualityType &quality,
int windowSize)
{
int bestScore = 1000000; // either mismatch count or sumQ
sequenceIndexType bestMatchLocation = -1;
for (int i=-windowSize; i<windowSize; i++)
{
int newScore;
// 'in' is a template, and types of all three args have to match
if( not in((size_t) index + i, (size_t) 0, sequence.size())) continue;
if (quality.size() == 0)
{
newScore = getMismatchCount(sequence, index + i, read);
}
else
{
newScore = getSumQ(sequence, index + i, read, quality);
}
if (newScore < bestScore)
{
bestScore = newScore;
bestMatchLocation = index + i;
}
}
return bestMatchLocation;
}
}
typedef uint32_t PackedVectorIndex_t;
//
// this is actually a base space implementation, since
// the load/set methods do indirection from the symbol
// value through a symbol value to symbol lookup table.
//
class PackedSequenceData : public PackedVector4Bit_t
{
public:
inline char operator[](PackedVectorIndex_t baseIndex)
{
return BaseAsciiMap::int2base[ (static_cast<PackedVector4Bit_t>(*this))[baseIndex]];
}
inline void set(PackedVectorIndex_t baseIndex, char value)
{
this->PackedVector4Bit_t::set(baseIndex, BaseAsciiMap::base2int[(uint32_t) value]);
}
inline void push_back(char value)
{
this->PackedVector4Bit_t::push_back(BaseAsciiMap::base2int[(uint32_t) value]);
}
};
std::ostream &operator << (std::ostream &stream, PackedSequenceData &v)
{
for(size_t i=0; i<v.size(); i++) {
stream << i << ": " << v[i] << std::endl;
}
return stream;
}
//
// Load a fasta format file from filename into the buffer
// provided by the caller.
// While parsing the fasta file, record each chromosome name,
// its start location, and its size.
//
// NB: the caller must implement the logic to determine how
// large the sequence data is. There is no correct way to do
// this, because we can't reliably estimate here how much sequence
// data is contained in a compressed file.
//
// To safely pre-allocate space in sequenceData, use the reserve() method
// before calling this function.
//
template<typename SequenceDataType, typename ChromosomeNameType>
bool loadFastaFile(const char *filename,
std::vector<SequenceDataType> &sequenceData,
std::vector<ChromosomeNameType> &chromosomeNames)
{
InputFile inputStream(filename, "r", InputFile::DEFAULT);
if(!inputStream.isOpen()) {
std::cerr << "Failed to open file " << filename << "\n";
return true;
}
int whichChromosome = -1;
//
// chromosomeNames is cheap to clear, so do it here.
//
// NB: I explicitly choose not to clear the sequence data
// container, this allows the caller to pre-allocate based
// on their knowledge of the size of the expected genome.
//
chromosomeNames.clear();
char ch;
while((ch = inputStream.ifgetc()) != EOF) {
switch (ch)
{
case '\n':
case '\r':
break;
case '>':
{
std::string chromosomeName = "";
//
// pull out the chromosome new name
//
while (!isspace((ch = inputStream.ifgetc())) && ch != EOF)
{
chromosomeName += ch; // slow, but who cares
}
//
// eat the rest of the line
//
do {
ch = inputStream.ifgetc();
} while(ch != EOF && ch != '\n' && ch != '\r');
//
// save the Chromosome name and index into our
// header so we can use them later.
//
chromosomeNames.push_back(chromosomeName);
whichChromosome++;
sequenceData.resize(whichChromosome+1);
break;
}
default:
// we get here for sequence data.
//
// save the base value
// Note: invalid characters come here as well, but we
// let ::set deal with mapping them.
sequenceData[whichChromosome].push_back(toupper(ch));
#if 0
if (isColorSpace())
{
//
// anything outside these values represents an invalid base
// base codes: 0-> A, 1-> C, 2-> G, 3-> T
// colorspace: 0-> blue, 1-> green, 2-> oragne, 3->red
//
const char fromBase2CS[] =
{
/* 0000 */ 0, // A->A
/* 0001 */ 1, // A->C
/* 0010 */ 2, // A->G
/* 0011 */ 3, // A->T
/* 0100 */ 1, // C->A
/* 0101 */ 0, // C->C
/* 0110 */ 3, // C->G
/* 0111 */ 2, // C->T
/* 1000 */ 2, // G->A
/* 1001 */ 3, // G->C
/* 1010 */ 0, // G->G
/* 1011 */ 1, // G->T
/* 1100 */ 3, // T->A
/* 1101 */ 2, // T->C
/* 1110 */ 1, // T->G
/* 1111 */ 0, // T->T
};
//
// we are writing color space values on transitions,
// so we don't write a colorspace value when we
// get the first base value.
//
// On second and subsequent bases, write based on
// the index table above
//
char thisBase = base2int[(int)(fasta[fastaIndex])];
if (lastBase>=0)
{
char color;
if (lastBase>3 || thisBase>3) color=4;
else color = fromBase2CS[(int)(lastBase<<2 | thisBase)];
// re-use the int to base, because ::set expects a base char (ATCG), not
// a color code (0123). It should only matter on final output.
set(header->elementCount++, int2base[(int) color]);
}
lastBase = thisBase;
}
else
{
set(header->elementCount++, toupper(fasta[fastaIndex]));
}
#endif
break;
}
}
return false;
}
#endif
BaseAsciiMap::base2int
static unsigned char base2int[256+1]
Map ASCII values to a 2 (or 3) bit encoding for the base pair value for just base space (ACTGNactgn).
Definition BaseAsciiMap.h:61
BaseAsciiMap::int2base
static const char int2base[]
Convert from int representation to the base.
Definition BaseAsciiMap.h:38
InputFile
Class for easily reading/writing files without having to worry about file type (uncompressed,...
Definition InputFile.h:37
InputFile::DEFAULT
@ DEFAULT
Check the extension, if it is ".gz", treat as gzip, otherwise treat it as UNCOMPRESSED.
Definition InputFile.h:45
PackedSequenceData
Definition GenomeSequence.cpp:382