GenomeSequence.cpp
/*
* Copyright (C) 2010-2012 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/>.
*/
#include "assert.h"
#include "ctype.h"
#include "stdio.h"
#include "Error.h"
#include "Generic.h"
#include "GenomeSequence.h"
#include <algorithm>
#include <istream>
#include <fstream>
#include <sstream>
#include <stdexcept>
#if defined(_WIN32)
#include <io.h>
#ifndef R_OK
#define R_OK 4
#endif
#endif
// not general use:
#include "CSG_MD5.h"
//
// given a read in a string, pack it into the vector of
// bytes coded as two bases per byte.
//
// The goal is to allow us to more rapidly compare against
// the genome, which is itself packed 2 bases per byte.
//
// Unfortunately, the match position may be odd or even,
// so as a result, we also need to be able to prepad
// with 'N' bases to get the byte alignment the same, so
// padWithNCount may be a positive number indicating how
// many N bases to prepend.
//
void PackedRead::set(const char *rhs, int padWithNCount)
{
clear();
// pad this packed read with 'N' bases two at a time
while (padWithNCount>1)
{
packedBases.push_back(
BaseAsciiMap::base2int[(int) 'N'] << 4 |
BaseAsciiMap::base2int[(int) 'N']
);
padWithNCount -= 2;
length+=2;
}
// when we have only one base, pack one 'N' base with
// the first base in rhs if there is one.
if (padWithNCount)
{
// NB: *rhs could be NUL, which is ok here - just keep
// the length straight.
packedBases.push_back(
BaseAsciiMap::base2int[(int) *rhs] << 4 |
BaseAsciiMap::base2int[(int) 'N']
);
// two cases - have characters in rhs or we don't:
if (*rhs)
{
length+=2; // pad byte plus one byte from rhs
rhs++;
}
else
{
length++;
}
padWithNCount--; // should now be zero, so superfluous.
assert(padWithNCount==0);
}
// pad pairs of bases from rhs, two at a time:
while (*rhs && *(rhs+1))
{
packedBases.push_back(
BaseAsciiMap::base2int[(int) *(rhs+1)] << 4 |
BaseAsciiMap::base2int[(int) *(rhs+0)]
);
rhs+=2;
length+=2;
}
// if there is an odd base left at the end, put it
// in a byte all its own (low 4 bits == 0):
if (*rhs)
{
packedBases.push_back(
BaseAsciiMap::base2int[(int) *(rhs+0)]
);
length++;
}
return;
}
std::string GenomeSequence::IntegerToSeq(unsigned int n, unsigned int wordsize) const
{
std::string sequence("");
for (unsigned int i = 0; i < wordsize; i ++)
sequence += "N";
unsigned int clearHigherBits = ~(3U << (wordsize<<1)); // XXX refactor - this appears several places
if (n > clearHigherBits)
error("%d needs to be a non-negative integer < clearHigherBits\n", n);
for (unsigned int i = 0; i < wordsize; i++)
{
sequence[wordsize-1-i] = BaseAsciiMap::int2base[n & 3];
n >>= 2;
}
return sequence;
}
GenomeSequence::GenomeSequence()
{
constructorClear();
}
void GenomeSequence::constructorClear()
{
_debugFlag = 0;
_progressStream = NULL;
_colorSpace = false;
_createOverwrite = false;
}
void GenomeSequence::setup(const char *referenceFilename)
{
setReferenceName(referenceFilename);
if (_progressStream) *_progressStream << "open and prefetch reference genome " << referenceFilename << ": " << std::flush;
if (open(false))
{
std::cerr << "Failed to open reference genome " << referenceFilename << std::endl;
std::cerr << errorStr << std::endl;
exit(1);
}
prefetch();
if (_progressStream) *_progressStream << "done." << std::endl << std::flush;
}
GenomeSequence::~GenomeSequence()
{
// free up resources:
_umfaFile.close();
}
//
// mapped open.
//
// if the file exists, map in into memory, and fill in a few useful
// fields.
//
bool GenomeSequence::open(bool isColorSpace, int flags)
{
bool rc;
if (isColorSpace)
{
_umfaFilename = _baseFilename + "-cs.umfa";
}
else
{
_umfaFilename = _baseFilename + "-bs.umfa";
}
if(access(_umfaFilename.c_str(), R_OK) != 0)
{
// umfa file doesn't exist, so try to create it.
if(create(isColorSpace))
{
// Couldon't access or create the umfa.
std::cerr << "GenomeSequence::open: failed to open file "
<< _umfaFilename
<< " also failed creating it."
<< std::endl;
return true;
}
}
rc = genomeSequenceArray::open(_umfaFilename.c_str(), flags);
if (rc)
{
std::cerr << "GenomeSequence::open: failed to open file "
<< _umfaFilename
<< std::endl;
return true;
}
_colorSpace = header->_colorSpace;
return false;
}
void GenomeSequence::sanityCheck(MemoryMap &fasta) const
{
unsigned int i;
unsigned int genomeIndex = 0;
for (i=0; i<fasta.length(); i++)
{
switch (fasta[i])
{
case '>':
while (fasta[i]!='\n' && fasta[i]!='\r') i++;
break;
case '\n':
case '\r':
break;
default:
assert(BaseAsciiMap::base2int[(int)(*this)[genomeIndex]] ==
BaseAsciiMap::base2int[(int) fasta[i]]);
genomeIndex++;
break;
}
}
}
#define HAS_SUFFIX(str, suffix) ((strlen(suffix) < str.size()) && (str.substr(str.size() - strlen(suffix)) == suffix))
//
// referenceFilename is either a fasta or a UM fasta (.fa or .umfa)
// filename. In both cases, the suffix gets removed and the
// base name is kept for later use depending on context.
// @return always return false
//
bool GenomeSequence::setReferenceName(std::string referenceFilename)
{
if (HAS_SUFFIX(referenceFilename, ".fa"))
{
_referenceFilename = referenceFilename;
_baseFilename = _referenceFilename.substr(0, referenceFilename.size() - 3);
}
else if (HAS_SUFFIX(referenceFilename, ".umfa"))
{
_baseFilename = referenceFilename.substr(0, referenceFilename.size() - 5);
}
else if (HAS_SUFFIX(referenceFilename, "-cs.umfa"))
{
_baseFilename = referenceFilename.substr(0, referenceFilename.size() - 8);
}
else if (HAS_SUFFIX(referenceFilename, "-bs.umfa"))
{
_baseFilename = referenceFilename.substr(0, referenceFilename.size() - 8);
}
else
{
_baseFilename = referenceFilename;
}
_fastaFilename = _baseFilename + ".fa";
if (HAS_SUFFIX(referenceFilename, ".fasta"))
{
_referenceFilename = referenceFilename;
_baseFilename = _referenceFilename.substr(0, referenceFilename.size() - 6);
_fastaFilename = _baseFilename + ".fasta";
}
return false;
}
//
// this works in lockstep with ::create to populate
// the per chromosome header fields size and md5
// checksum.
//
// It relies on header->elementCount being set to
// the data length loaded so far ... not the ultimate
// reference length.
//
bool GenomeSequence::setChromosomeMD5andLength(uint32_t whichChromosome)
{
if (whichChromosome>=header->_chromosomeCount) return true;
ChromosomeInfo *c = &header->_chromosomes[whichChromosome];
c->size = header->elementCount - c->start;
MD5_CTX md5Context;
uint8_t md5Signature[MD5_DIGEST_LENGTH];
//
// it's easier to ensure we do this right if we just do it
// in one big chunk:
//
char *md5Buffer = (char *) malloc(c->size);
MD5Init(&md5Context);
for (genomeIndex_t i = 0; i < c->size; i ++)
{
md5Buffer[i] = (*this)[c->start + i];
}
MD5Update(&md5Context, (unsigned char *) md5Buffer, c->size);
MD5Final((unsigned char *) &md5Signature, &md5Context);
free(md5Buffer);
for (int i=0; i<MD5_DIGEST_LENGTH; i++)
{
// cheesy but works:
sprintf(c->md5+2*i, "%02x", md5Signature[i]);
}
// redundant, strictly speaking due to sprintf NUL terminating
// it's output strings, but put it here anyway.
c->md5[2*MD5_DIGEST_LENGTH] = '\0';
return false;
}
//
// Given a buffer with a fasta format contents, count
// the number of chromosomes in it and return that value.
//
static bool getFastaStats(const char *fastaData, size_t fastaDataSize, uint32_t &chromosomeCount, uint64_t &baseCount)
{
chromosomeCount = 0;
baseCount = 0;
bool atLineStart = true;
//
// loop over the fasta file, essentially matching for the
// pattern '^>.*$' and counting them.
//
for (size_t fastaIndex = 0; fastaIndex < fastaDataSize; fastaIndex++)
{
switch (fastaData[fastaIndex])
{
case '\n':
case '\r':
atLineStart = true;
break;
case '>':
{
if (!atLineStart) break;
chromosomeCount++;
//
// eat the rest of the line
//
while (fastaIndex < fastaDataSize && fastaData[fastaIndex]!='\n' && fastaData[fastaIndex]!='\r')
{
fastaIndex++;
}
break;
}
default:
baseCount++;
atLineStart = false;
break;
}
}
return false;
}
class PackedSequenceData : public std::vector<uint8_t>
{
std::vector<uint8_t> m_packedBases;
size_t m_baseCount;
void set(size_t index, uint8_t value) {
m_packedBases[index>>1] =
(m_packedBases[index>>1] // original value
& ~(7<<((index&0x01)<<2))) // logical AND off the original value
| ((value&0x0f)<<((index&0x1)<<2)); // logical OR in the new value
}
public:
void reserve(size_t baseCount) {m_packedBases.reserve(baseCount/2);}
size_t size() {return m_baseCount;}
void clear() {m_packedBases.clear(); m_baseCount = 0;}
uint8_t operator [](size_t index)
{
return (m_packedBases[index>>1] >> ((index&0x1)<<2)) & 0xf;
}
void push_back(uint8_t base);
};
//
// 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.
//
bool loadFastaFile(const char *filename,
std::vector<PackedSequenceData> &sequenceData,
std::vector<std::string> &chromosomeNames)
{
InputFile inputStream(filename, "r", InputFile::DEFAULT);
if(!inputStream.isOpen()) {
std::cerr << "Failed to open file " << filename << "\n";
return true;
}
int whichChromosome = -1;
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++;
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.
break;
}
}
return false;
}
//
// recreate the umfa file from a reference fasta format file
//
// The general format of a FASTA file is best described
// on wikipedia at http://en.wikipedia.org/wiki/FASTA_format
//
// The format parsed here is a simpler subset, and is
// described here http://www.ncbi.nlm.nih.gov/blast/fasta.shtml
//
bool GenomeSequence::create(bool isColor)
{
setColorSpace(isColor);
if (_baseFilename=="")
{
std::cerr << "Base reference filename is empty." << std::endl;
return true;
}
if (isColorSpace())
{
_umfaFilename = _baseFilename + "-cs.umfa";
}
else
{
_umfaFilename = _baseFilename + "-bs.umfa";
}
if (!_createOverwrite && access(_umfaFilename.c_str(), R_OK) == 0)
{
std::cerr << "Output file '" << _umfaFilename << "' exists or is not writable - please remove." << std::endl;
return true;
}
MemoryMap fastaFile;
if (fastaFile.open(_fastaFilename.c_str()))
{
std::cerr << "failed to open input fasta file '" << _fastaFilename << "'." << std::endl;
return true;
}
std::cerr << "Creating FASTA "
<< (isColorSpace() ? "color space " : "")
<< "binary cache file '"
<< _umfaFilename
<< "'."
<< std::endl;
std::cerr << std::flush;
//
// simple ptr to fasta data -- just treat the memory map
// as an array of fastaDataSize characters...
//
const char *fasta = (const char *) fastaFile.data;
size_t fastaDataSize = fastaFile.length();
uint32_t chromosomeCount = 0;
uint64_t baseCount = 0;
getFastaStats(fasta, fastaDataSize, chromosomeCount, baseCount);
if (genomeSequenceArray::create(_umfaFilename.c_str(), baseCount, chromosomeCount))
{
std::cerr << "failed to create '"
<< _umfaFilename
<< "'."
<< std::endl;
perror("");
return true;
}
header->elementCount = 0;
header->_colorSpace = isColorSpace();
header->setApplication(_application.c_str());
header->_chromosomeCount = chromosomeCount;
//
// clear out the variable length chromosome info array
//
for (uint32_t i=0; i<header->_chromosomeCount; i++) header->_chromosomes[i].constructorClear();
std::string chromosomeName;
//
// for converting the reference to colorspace, the first base is always 5 (in base space it is 'N')
signed char lastBase = BaseAsciiMap::base2int[(int) 'N'];
bool terminateLoad = false;
int percent = -1, newPercent;
uint32_t whichChromosome = 0;
for (uint64_t fastaIndex = 0; fastaIndex < fastaDataSize; fastaIndex++)
{
if (_progressStream)
{
newPercent = (int) (1.0 * fastaIndex / fastaDataSize) * 100;
if (newPercent>percent)
{
*_progressStream << "\r" << newPercent << "% ";
*_progressStream << std::flush;
percent = newPercent;
}
}
switch (fasta[fastaIndex])
{
case '\n':
case '\r':
break;
case '>':
{
chromosomeName = "";
fastaIndex++; // skip the > char
//
// pull out the chromosome new name
//
while (!isspace(fasta[fastaIndex]))
{
chromosomeName += fasta[fastaIndex++]; // slow, but who cares
}
//
// eat the rest of the line
//
while (fasta[fastaIndex]!='\n' && fasta[fastaIndex]!='\r')
{
fastaIndex++;
}
//
// save the Chromosome name and index into our
// header so we can use them later.
//
ChromosomeInfo *c = &header->_chromosomes[whichChromosome];
c->setChromosomeName(chromosomeName.c_str());
c->start = header->elementCount;
// c->size gets computed at the next '>' line or at the EOF
if (whichChromosome>0)
{
//
// compute md5 checksum for the chromosome that we just
// loaded (if there was one) - note that on the last
// chromosome, we have to duplicate this code after
// the end of this loop
//
setChromosomeMD5andLength(whichChromosome - 1);
}
whichChromosome++;
if (whichChromosome > header->_chromosomeCount)
{
std::cerr << "BUG: Exceeded computed chromosome count ("
<< header->_chromosomeCount
<< ") - genome is now truncated at chromosome "
<< header->_chromosomes[header->_chromosomeCount-1].name
<< " (index "
<< header->_chromosomeCount
<< ")."
<< std::endl;
terminateLoad = true;
}
break;
}
default:
// save the base pair value
// Note: invalid characters come here as well, but we
// let ::set deal with mapping them.
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 =
BaseAsciiMap::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++,
BaseAsciiMap::int2base[(int) color]);
}
lastBase = thisBase;
}
else
{
set(header->elementCount++, toupper(fasta[fastaIndex]));
}
break;
}
//
// slightly awkward exit handling when we exceed the fixed
// number of chromosomes
//
if (terminateLoad) break;
}
//
// also slightly awkward code to handle the last dangling chromosome...
// all we should need to do is compute the md5 checksum
//
if (whichChromosome==0)
{
fastaFile.close();
throw std::runtime_error("No chromosomes found - aborting!");
}
else
{
setChromosomeMD5andLength(whichChromosome-1);
}
fastaFile.close();
if (_progressStream) *_progressStream << "\r";
std::cerr << "FASTA binary cache file '"
<< _umfaFilename
<< "' created."
<< std::endl;
//
// leave the umfastaFile open in case caller wants to use it
//
return false;
}
int GenomeSequence::getChromosomeCount() const
{
return header->_chromosomeCount;
}
//return chromosome index: 0, 1, ... 24;
int GenomeSequence::getChromosome(genomeIndex_t position) const
{
if (position == INVALID_GENOME_INDEX) return INVALID_CHROMOSOME_INDEX;
if (header->_chromosomeCount == 0)
return INVALID_CHROMOSOME_INDEX;
int start = 0;
int stop = header->_chromosomeCount - 1;
// eliminate case where position is in the last chromosome, since the loop
// below falls off the end of the list if it in the last one.
if (position > header->_chromosomes[stop].start)
return (stop);
while (start <= stop)
{
int middle = (start + stop) / 2;
if (position >= header->_chromosomes[middle].start && position < header->_chromosomes[middle + 1].start)
return middle;
if (position == header->_chromosomes[middle + 1].start)
return (middle + 1);
if (position > header->_chromosomes[middle + 1].start)
start = middle + 1;
if (position < header->_chromosomes[middle].start)
stop = middle - 1;
}
return -1;
}
//
// Given a chromosome name and 1-based chromosome index, return the
// genome index (0 based) into sequence for it.
//
// NB: the header->chromosomes array contains zero based genome positions
//
genomeIndex_t GenomeSequence::getGenomePosition(
const char *chromosomeName,
unsigned int chromosomeIndex) const
{
genomeIndex_t i = getGenomePosition(chromosomeName);
if (i == INVALID_GENOME_INDEX) return INVALID_GENOME_INDEX;
return i + chromosomeIndex - 1;
}
genomeIndex_t GenomeSequence::getGenomePosition(
int chromosome,
unsigned int chromosomeIndex) const
{
if (chromosome<0 || chromosome >= (int) header->_chromosomeCount) return INVALID_GENOME_INDEX;
genomeIndex_t i = header->_chromosomes[chromosome].start;
if (i == INVALID_GENOME_INDEX) return INVALID_GENOME_INDEX;
return i + chromosomeIndex - 1;
}
//
// return the genome index (0 based) of the start of the named
// chromosome. If none is found, INVALID_GENOME_INDEX is returned.
//
// XXX may need to speed this up - and smarten it up with some
// modest chromosome name parsing.... e.g. '%d/X/Y' or 'chr%d/chrX/chrY' or
// other schemes.
//
genomeIndex_t GenomeSequence::getGenomePosition(const char *chromosomeName) const
{
int chromosome = getChromosome(chromosomeName);
if (chromosome==INVALID_CHROMOSOME_INDEX) return INVALID_GENOME_INDEX;
return header->_chromosomes[chromosome].start;
}
int GenomeSequence::getChromosome(const char *chromosomeName) const
{
unsigned int i;
for (i=0; i<header->_chromosomeCount; i++)
{
if (strcmp(header->_chromosomes[i].name, chromosomeName)==0)
{
return i;
}
}
return INVALID_CHROMOSOME_INDEX;
}
//
// Given a read, reverse the string and swap the base
// pairs for the reverse strand equivalents.
//
void GenomeSequence::getReverseRead(std::string &read)
{
std::string newRead;
if (read.size()) for (int32_t i=(int) read.size() - 1; i>=0; i--)
{
newRead.push_back(BasePair(read[i]));
}
read = newRead;
}
void GenomeSequence::getReverseRead(String& read)
{
int i = 0;
int j = read.Length()-1;
char temp;
while (i < j)
{
temp = read[j];
read[j] = read[i];
read[i] = temp;
}
}
#define ABS(x) ( (x) > 0 ? (x) : -(x) )
int GenomeSequence::debugPrintReadValidation(
std::string &read,
std::string &quality,
char direction,
genomeIndex_t readLocation,
int sumQuality,
int mismatchCount,
bool recurse
)
{
int validateSumQ = 0;
int validateMismatchCount = 0;
int rc = 0;
std::string genomeData;
for (uint32_t i=0; i<read.size(); i++)
{
if (tolower(read[i]) != tolower((*this)[readLocation + i]))
{
validateSumQ += quality[i] - '!';
// XXX no longer valid:
if (direction=='F' ? i<24 : (i >= (read.size() - 24))) validateMismatchCount++;
genomeData.push_back(tolower((*this)[readLocation + i]));
}
else
{
genomeData.push_back(toupper((*this)[readLocation + i]));
}
}
assert(validateSumQ>=0);
if (validateSumQ != sumQuality && validateMismatchCount == mismatchCount)
{
printf("SUMQ: Original Genome: %s test read: %s : actual sumQ = %d, test sumQ = %d\n",
genomeData.c_str(),
read.c_str(),
validateSumQ,
sumQuality
);
rc++;
}
else if (validateSumQ == sumQuality && validateMismatchCount != mismatchCount)
{
printf("MISM: Original Genome: %s test read: %s : actual mismatch %d test mismatches %d\n",
genomeData.c_str(),
read.c_str(),
validateMismatchCount,
mismatchCount
);
rc++;
}
else if (validateSumQ != sumQuality && validateMismatchCount != mismatchCount)
{
printf("BOTH: Original Genome: %s test read: %s : actual sumQ = %d, test sumQ = %d, actual mismatch %d test mismatches %d\n",
genomeData.c_str(),
read.c_str(),
validateSumQ,
sumQuality,
validateMismatchCount,
mismatchCount
);
rc++;
}
if (recurse && ABS(validateMismatchCount - mismatchCount) > (int) read.size()/2)
{
printf("large mismatch difference, trying reverse strand: ");
std::string reverseRead = read;
std::string reverseQuality = quality;
getReverseRead(reverseRead);
reverse(reverseQuality.begin(), reverseQuality.end());
rc = debugPrintReadValidation(reverseRead, reverseQuality, readLocation, sumQuality, mismatchCount, false);
}
return rc;
}
#undef ABS
bool GenomeSequence::wordMatch(unsigned int index, std::string &word) const
{
for (uint32_t i = 0; i<word.size(); i++)
{
if ((*this)[index + i] != word[i]) return false;
}
return true;
}
bool GenomeSequence::printNearbyWords(unsigned int index, unsigned int deviation, std::string &word) const
{
for (unsigned int i = index - deviation; i < index + deviation; i++)
{
if (wordMatch(i, word))
{
std::cerr << "word '"
<< word
<< "' found "
<< i - index
<< " away from position "
<< index
<< "."
<< std::endl;
}
}
return false;
}
void GenomeSequence::dumpSequenceSAMDictionary(std::ostream &file) const
{
for (unsigned int i=0; i<header->_chromosomeCount; i++)
{
file
<< "@SQ"
<< "\tSN:" << header->_chromosomes[i].name // name
<< "\tLN:" << header->_chromosomes[i].size // number of bases
<< "\tAS:" << header->_chromosomes[i].assemblyID // e.g. NCBI36.3
<< "\tM5:" << header->_chromosomes[i].md5
<< "\tUR:" << header->_chromosomes[i].uri
<< "\tSP:" << header->_chromosomes[i].species // e.g. Homo_sapiens
<< std::endl;
}
}
void GenomeSequence::dumpHeaderTSV(std::ostream &file) const
{
file << "# Reference: " << _baseFilename << std::endl;
file << "# SN: sample name - must be unique" << std::endl;
file << "# AS: assembly name" << std::endl;
file << "# SP: species" << std::endl;
file << "# LN: chromosome/contig length" << std::endl;
file << "# M5: chromosome/contig MD5 checksum" << std::endl;
file << "# LN and M5 are only printed for informational purposes." << std::endl;
file << "# Karma will only set those values when creating the index." << std::endl;
file << "SN" << "\t" << "AS" << "\t" << "SP" << "\t" << "UR" << "\t" << "LN" << "\t" << "M5" << std::endl;
for (unsigned int i=0; i<header->_chromosomeCount; i++)
{
file
<< header->_chromosomes[i].name // name
<< "\t" << header->_chromosomes[i].assemblyID // e.g. NCBI36.3
<< "\t" << header->_chromosomes[i].uri
<< "\t" << header->_chromosomes[i].species // e.g. Homo_sapiens
<< "\t" << header->_chromosomes[i].size // number of bases
<< "\t" << header->_chromosomes[i].md5
<< std::endl;
}
}
void GenomeSequence::getString(std::string &str, int chromosome, uint32_t index, int baseCount) const
{
//
// calculate the genome index for the lazy caller...
//
genomeIndex_t genomeIndex = header->_chromosomes[chromosome].start + index - 1;
getString(str, genomeIndex, baseCount);
}
void GenomeSequence::getString(String &str, int chromosome, uint32_t index, int baseCount) const
{
std::string string;
this-> getString(string, chromosome, index, baseCount);
str = string.c_str();
}
void GenomeSequence::getString(std::string &str, genomeIndex_t index, int baseCount) const
{
str.clear();
if (baseCount > 0)
{
for (int i=0; i<baseCount; i++)
{
str.push_back((*this)[index + i]);
}
}
else
{
// if caller passed negative basecount, give them
// the read for the 3' end
for (int i=0; i< -baseCount; i++)
{
str.push_back(BaseAsciiMap::base2complement[(int)(*this)[index + i]]);
}
}
}
void GenomeSequence::getString(String &str, genomeIndex_t index, int baseCount) const
{
std::string string;
getString(string, index, baseCount);
str = string.c_str();
}
void GenomeSequence::getHighLightedString(std::string &str, genomeIndex_t index, int baseCount, genomeIndex_t highLightStart, genomeIndex_t highLightEnd) const
{
str.clear();
if (baseCount > 0)
{
for (int i=0; i<baseCount; i++)
{
char base = (*this)[index + i];
if (in(index+i, highLightStart, highLightEnd))
base = tolower(base);
str.push_back(base);
}
}
else
{
// if caller passed negative basecount, give them
// the read for the 3' end
for (int i=0; i< -baseCount; i++)
{
char base = BaseAsciiMap::base2complement[(int)(*this)[index + i]];
if (in(index+i, highLightStart, highLightEnd))
base = tolower(base);
str.push_back(base);
}
}
}
void GenomeSequence::print30(genomeIndex_t index) const
{
std::cout << "index: " << index << "\n";
for (genomeIndex_t i=index-30; i<index+30; i++)
std::cout << (*this)[i];
std::cout << "\n";
for (genomeIndex_t i=index-30; i<index; i++)
std::cout << " ";
std::cout << "^";
std::cout << std::endl;
}
void GenomeSequence::getMismatchHatString(std::string &result, const std::string &read, genomeIndex_t location) const
{
result = "";
for (uint32_t i=0; i < read.size(); i++)
{
if (read[i] == (*this)[location+i])
result.push_back(' ');
else
result.push_back('^');
}
}
void GenomeSequence::getMismatchString(std::string &result, const std::string &read, genomeIndex_t location) const
{
result = "";
for (uint32_t i=0; i < read.size(); i++)
{
if (read[i] == (*this)[location+i])
result.push_back(toupper(read[i]));
else
result.push_back(tolower(read[i]));
}
}
genomeIndex_t GenomeSequence::simpleLocalAligner(std::string &read, std::string &quality, genomeIndex_t index, int windowSize) const
{
int bestScore = 1000000; // either mismatch count or sumQ
genomeIndex_t bestMatchLocation = INVALID_GENOME_INDEX;
for (int i=-windowSize; i<windowSize; i++)
{
int newScore;
if (i<0 && ((uint32_t) -i) > index) continue;
if (index + i + read.size() >= getNumberBases()) continue;
if (quality=="")
{
newScore = this->getMismatchCount(read, index + i);
}
else
{
newScore = this->getSumQ(read, quality, index + i);
}
if (newScore < bestScore)
{
bestScore = newScore;
bestMatchLocation = index + i;
}
}
return bestMatchLocation;
}
std::ostream &operator << (std::ostream &stream, genomeSequenceMmapHeader &h)
{
stream << (MemoryMapArrayHeader &) h;
stream << "chromosomeCount: " << h._chromosomeCount << std::endl;
stream << "isColorSpace: " << h._colorSpace << std::endl;
stream << "chromosomeCount: " << h._chromosomeCount << std::endl;
uint64_t totalSize = 0;
for (uint32_t i=0; i < h._chromosomeCount; i++)
{
totalSize += h._chromosomes[i].size;
stream << "Chromosome Index " << i << " name: " << h._chromosomes[i].name << std::endl;
stream << "Chromosome Index " << i << " whole genome start: " << h._chromosomes[i].start << std::endl;
stream << "Chromosome Index " << i << " whole genome size: " << h._chromosomes[i].size << std::endl;
stream << "Chromosome Index " << i << " md5 checksum: " << h._chromosomes[i].md5 << std::endl;
stream << "Chromosome Index " << i << " assemblyID: " << h._chromosomes[i].assemblyID << std::endl;
stream << "Chromosome Index " << i << " species: " << h._chromosomes[i].species << std::endl;
stream << "Chromosome Index " << i << " URI: " << h._chromosomes[i].uri << std::endl;
}
stream << "Total Genome Size: " << totalSize << " bases."<< std::endl;
if (totalSize != h.elementCount)
{
stream << "Total Genome Size: does not match elementCount!\n";
}
stream << std::endl;
return stream;
}
void GenomeSequence::getChromosomeAndIndex(std::string &s, genomeIndex_t i) const
{
int whichChromosome = 0;
whichChromosome = getChromosome(i);
if (whichChromosome == INVALID_CHROMOSOME_INDEX)
{
s = "invalid genome index"; // TODO include the index in error
}
else
{
std::ostringstream buf;
genomeIndex_t chromosomeIndex = i - getChromosomeStart(whichChromosome) + 1;
buf << header->_chromosomes[whichChromosome].name << ":" << chromosomeIndex;
#if 0
buf << " (GenomeIndex " << i << ")";
#endif
s = buf.str();
}
return;
}
void GenomeSequence::getChromosomeAndIndex(String& s, genomeIndex_t i) const
{
std::string ss;
getChromosomeAndIndex(ss, i);
s = ss.c_str();
return;
}
//
// This is intended to be a helper routine to get dbSNP files
// loaded. In some cases, we will load into an mmap() file (ie
// when we are creating it), in others, we will simply be loading
// an existing dbSNP file into RAM (when the binary file does not
// exist or when we are running with useMemoryMapFlag == false.
//
// Assume that dbSNP exists, is writable, and is the right size.
//
// Using the dbSNPFilename given, mark each dbSNP position
// with a bool true.
//
// Return value:
// True: if populateDBSNP() succeed
// False: if not succeed
bool GenomeSequence::populateDBSNP(
mmapArrayBool_t &dbSNP,
IFILE inputFile) const
{
assert(dbSNP.getElementCount() == getNumberBases());
if(inputFile == NULL)
{
// FAIL, file not opened.
return(false);
}
std::string chromosomeName;
std::string position;
genomeIndex_t chromosomePosition1; // 1-based
uint64_t ignoredLineCount = 0;
// Read til the end of the file.
char* postPosPtr = NULL;
while(!inputFile->ifeof())
{
chromosomeName.clear();
position.clear();
// Read the chromosome
if(inputFile->readTilTab(chromosomeName) <= 0)
{
// hit either eof or end of line, check if
// it is a header.
if(chromosomeName.size()>0 && chromosomeName[0]=='#')
{
// header, so just continue.
continue;
}
// Not the header, so this line is poorly formatted.
++ignoredLineCount;
// Continue to the next line.
continue;
}
// Check if it is a header line.
if(chromosomeName.size()>0 && chromosomeName[0]=='#')
{
// did not hit eof or end of line,
// so discard the rest of the line.
inputFile->discardLine();
continue;
}
// Not a header, so read the position.
if(inputFile->readTilTab(position) > 0)
{
// Additional data on the line, so discard it.
inputFile->discardLine();
}
// Convert the position to a string.
chromosomePosition1 = strtoul(position.c_str(), &postPosPtr, 0);
if(postPosPtr == position.c_str())
{
++ignoredLineCount;
continue;
}
// 1-based genome index.
genomeIndex_t genomeIndex =
getGenomePosition(chromosomeName.c_str(), chromosomePosition1);
// if the genome index is invalid, ignore it
if((genomeIndex == INVALID_GENOME_INDEX) ||
(genomeIndex > getNumberBases()))
{
ignoredLineCount++;
continue;
}
dbSNP.set(genomeIndex, true);
}
if (ignoredLineCount > 0)
{
std::cerr << "GenomeSequence::populateDBSNP: ignored " << ignoredLineCount << " SNP positions due to invalid format of line." << std::endl;
return false;
}
return true;
}
bool GenomeSequence::loadDBSNP(
mmapArrayBool_t &dbSNP,
const char *inputFileName) const
{
//
// the goal in this section of code is to allow the user
// to either specify a valid binary version of the SNP file,
// or the original text file that it gets created from.
//
// To do this, we basically open, sniff the error message,
// and if it claims it is not a binary version of the file,
// we go ahead and treat it as the text file and use the
// GenomeSequence::populateDBSNP method to load it.
//
// Further checking is really needed to ensure users don't
// mix a dbSNP file for a different reference, since it is really
// easy to do.
//
if (strlen(inputFileName)!=0)
{
std::cerr << "Load dbSNP file '" << inputFileName << "': " << std::flush;
if (dbSNP.open(inputFileName, O_RDONLY))
{
//
// failed to open, possibly due to bad magic.
//
// this is really awful ... need to have a return
// code that is smart enough to avoid this ugliness:
//
if (dbSNP.getErrorString().find("wrong type of file")==std::string::npos)
{
std::cerr << "Error: " << dbSNP.getErrorString() << std::endl;
exit(1);
}
//
// we have a file, assume we can load it as a text file
//
IFILE inputFile = ifopen(inputFileName, "r");
if(inputFile == NULL)
{
std::cerr << "Error: failed to open " << inputFileName << std::endl;
exit(1);
}
std::cerr << "(as text file) ";
// anonymously (RAM resident only) create:
dbSNP.create(getNumberBases());
// now load it into RAM
populateDBSNP(dbSNP, inputFile);
ifclose(inputFile);
}
else
{
std::cerr << "(as binary mapped file) ";
}
std::cerr << "DONE!" << std::endl;
return false;
}
else
{
return true;
}
}
#if defined(TEST)
void simplestExample(void)
{
GenomeSequence reference;
genomeIndex_t index;
// a particular reference is set by:
// reference.setFastaName("/usr/cluster/share/karma/human_g1k_v37_12CS.fa")
//
// In the above example, the suffix .fa is stripped and replaced with .umfa,
// which contains the actual file being opened.
//
if (reference.open())
{
perror("GenomeSequence::open");
exit(1);
}
index = 1000000000; // 10^9
//
// Write the base at the given index. Here, index is 0 based,
// and is across the whole genome, as all chromosomes are sequentially
// concatenated, so the allowed range is
//
// 0.. (reference.getChromosomeStart(last) + reference.getChromosomeSize(last))
//
// (where int last = reference.getChromosomeCount() - 1;)
//
std::cout << "base[" << index << "] = " << reference[index] << std::endl;
//
// Example for finding chromosome and one based chromosome position given
// and absolute position on the genome in 'index':
//
int chr = reference.getChromosome(index);
genomeIndex_t chrIndex = index - reference.getChromosomeStart(chr) + 1; // 1-based
std::cout << "genome index " << index << " corresponds to chromosome " << chr << " position " << chrIndex << std::endl;
//
// Example for finding an absolute genome index position when the
// chromosome name and one based position are known:
//
const char *chromosomeName = "5";
chr = reference.getChromosome(chromosomeName); // 0-based
chrIndex = 100000; // 1-based
index = reference.getChromosomeStart(chr) + chrIndex - 1;
std::cout << "Chromosome '" << chromosomeName << "' position " << chrIndex << " corresponds to genome index position " << index << std::endl;
reference.close();
}
void testGenomeSequence(void)
{
GenomeSequence reference;
#if 0
std::string referenceName = "someotherreference";
if (reference.setFastaName(referenceName))
{
std::cerr << "failed to open reference file "
<< referenceName
<< std::endl;
exit(1);
}
#endif
std::cerr << "open and prefetch the reference genome: ";
// open it
if (reference.open())
{
exit(1);
}
std::cerr << "done!" << std::endl;
//
// For the human genome, genomeIndex ranges from 0 to 3.2x10^9
//
genomeIndex_t genomeIndex; // 0 based
unsigned int chromosomeIndex; // 1 based
unsigned int chromosome; // 0..23 or so
std::string chromosomeName;
//
// Here we'll start with a chromosome name, then obtain the genome
// index, and use it to find the base we want:
//
chromosomeName = "2";
chromosomeIndex = 1234567;
// this call is slow (string search for chromsomeName):
genomeIndex = reference.getGenomePosition(chromosomeName.c_str(), chromosomeIndex);
assert(genomeIndex!=INVALID_GENOME_INDEX);
std::cout << "Chromosome " << chromosomeName << ", index ";
std::cout << chromosomeIndex << " contains base " << reference[genomeIndex];
std::cout << " at genome index position " << genomeIndex << std::endl;
//
// now reverse it - given a genomeIndex from above, find the chromosome
// name and index:
//
// slow (binary search on genomeIndex):
chromosome = reference.getChromosome(genomeIndex);
unsigned int newChromosomeIndex;
// not slow:
newChromosomeIndex = genomeIndex - reference.getChromosomeStart(chromosome) + 1;
assert(chromosomeIndex == newChromosomeIndex);
// more testing... at least test and use PackedRead:
//
PackedRead pr;
pr.set("ATCGATCG", 0);
assert(pr.size()==8);
assert(pr[0]==BaseAsciiMap::base2int[(int) 'A']);
assert(pr[1]==BaseAsciiMap::base2int[(int) 'T']);
assert(pr[2]==BaseAsciiMap::base2int[(int) 'C']);
assert(pr[3]==BaseAsciiMap::base2int[(int) 'G']);
pr.set("ATCGATCG", 1);
assert(pr.size()==9);
pr.set("", 0);
assert(pr.size()==0);
pr.set("", 1);
assert(pr.size()==1);
pr.set("", 2);
assert(pr.size()==2);
pr.set("", 3);
assert(pr.size()==3);
assert(pr[0]==BaseAsciiMap::base2int[(int) 'N']);
assert(pr[1]==BaseAsciiMap::base2int[(int) 'N']);
assert(pr[2]==BaseAsciiMap::base2int[(int) 'N']);
pr.set("C", 1);
assert(pr.size()==2);
assert(pr[0]==BaseAsciiMap::base2int[(int) 'N']);
assert(pr[1]==BaseAsciiMap::base2int[(int) 'C']);
}
//
// After I build libcsg, I compile and run this test code using:
//
// g++ -DTEST -o try GenomeSequence.cpp -L. -lcsg -lm -lz -lssl
// you also may need -fno-rtti
// ./try
//
int main(int argc, const char **argv)
{
#if 1
simplestExample();
#else
testGenomeSequence();
#endif
exit(0);
}
#endif
ifopen
IFILE ifopen(const char *filename, const char *mode, InputFile::ifileCompression compressionMode=InputFile::DEFAULT)
Open a file with the specified name and mode, using a filename of "-" to indicate stdin/stdout.
Definition: InputFile.h:562
operator<<
InputFile & operator<<(InputFile &stream, const std::string &str)
Write to a file using streaming.
Definition: InputFile.h:736
ifclose
int ifclose(IFILE &file)
Close the file.
Definition: InputFile.h:580
BaseAsciiMap::base2complement
static unsigned char base2complement[]
This table maps 5' base space to the 3' complement base space values, as well as 5' color space value...
Definition: BaseAsciiMap.h:41
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
GenomeSequence
Create/Access/Modify/Load Genome Sequences stored as binary mapped files.
Definition: GenomeSequence.h:100
GenomeSequence::loadDBSNP
bool loadDBSNP(mmapArrayBool_t &dbSNP, const char *inputFileName) const
user friendly dbSNP loader.
Definition: GenomeSequence.cpp:1301
GenomeSequence::getChromosomeCount
int getChromosomeCount() const
Return the number of chromosomes in the genome.
Definition: GenomeSequence.cpp:731
GenomeSequence::getSumQ
int getSumQ(std::string &read, std::string &qualities, genomeIndex_t location) const
brute force sumQ - no sanity checking
Definition: GenomeSequence.h:501
GenomeSequence::getMismatchCount
int getMismatchCount(std::string &read, genomeIndex_t location, char exclude='\0') const
Return the mismatch count, disregarding CIGAR strings.
Definition: GenomeSequence.h:488
GenomeSequence::getChromosomeStart
genomeIndex_t getChromosomeStart(int chromosomeIndex) const
given a chromosome, return the genome base it starts in
Definition: GenomeSequence.h:246
GenomeSequence::getGenomePosition
genomeIndex_t getGenomePosition(const char *chromosomeName, unsigned int chromosomeIndex) const
given a chromosome name and position, return the genome position
Definition: GenomeSequence.cpp:779
GenomeSequence::~GenomeSequence
~GenomeSequence()
Close the file if open and destroy the object.
Definition: GenomeSequence.cpp:169
GenomeSequence::isColorSpace
bool isColorSpace() const
tell us if we are a color space reference or not
Definition: GenomeSequence.h:209
GenomeSequence::getChromosome
int getChromosome(genomeIndex_t position) const
given a whole genome index, get the chromosome it is located in
Definition: GenomeSequence.cpp:737
GenomeSequence::setReferenceName
bool setReferenceName(std::string referenceFilename)
set the reference name that will be used in open()
Definition: GenomeSequence.cpp:254
GenomeSequence::getNumberBases
genomeIndex_t getNumberBases() const
return the number of bases represented in this reference
Definition: GenomeSequence.h:216
GenomeSequence::GenomeSequence
GenomeSequence()
Simple constructor - no implicit file open.
Definition: GenomeSequence.cpp:139
GenomeSequence::open
bool open(bool isColorSpace=false, int flags=O_RDONLY)
open the reference specified using GenomeSequence::setReferenceName
Definition: GenomeSequence.cpp:182
InputFile
Class for easily reading/writing files without having to worry about file type (uncompressed,...
Definition: InputFile.h:37
InputFile::ifeof
int ifeof() const
Check to see if we have reached the EOF.
Definition: InputFile.h:386
InputFile::discardLine
int discardLine()
Read until the end of the line, discarding the characters, returning -1 returned for EOF and returnin...
Definition: InputFile.cpp:95
InputFile::DEFAULT
@ DEFAULT
Check the extension, if it is ".gz", treat as gzip, otherwise treat it as UNCOMPRESSED.
Definition: InputFile.h:45
InputFile::readTilTab
int readTilTab(std::string &field)
Read, appending the characters into the specified string until tab, new line, or EOF is found,...
Definition: InputFile.cpp:133
MemoryMapArrayHeader
Definition: MemoryMapArray.h:65
MemoryMapArray
Definition: MemoryMapArray.h:142
MemoryMapArray::open
bool open(const char *file, int flags=O_RDONLY)
open a previously created mapped vector
Definition: MemoryMapArray.h:269
MemoryMapArray::create
int create(const char *file, indexT elementCount, int optionalHeaderCount=0)
Create a vector with elementCount memebers.
Definition: MemoryMapArray.h:208
MemoryMap
There are a pair of related data structures in the operating system, and also a few simple algorithms...
Definition: MemoryMap.h:156
MemoryMap::open
virtual bool open(const char *file, int flags=O_RDONLY)
open a previously created mapped vector
Definition: MemoryMap.cpp:156
PackedRead
Definition: GenomeSequenceHelpers.h:142
PackedSequenceData
Definition: GenomeSequence.cpp:382
String
Definition: StringBasics.h:39
genomeSequenceMmapHeader
Definition: GenomeSequenceHelpers.h:80
ChromosomeInfo
Definition: GenomeSequenceHelpers.h:41
MD5_CTX
Definition: CSG_MD5.h:60