GenomeSequenceHelpers.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 _GENOME_SEQUENCE_HELPERS_H
#define _GENOME_SEQUENCE_HELPERS_H
#if !defined(MD5_DIGEST_LENGTH)
#define MD5_DIGEST_LENGTH 16
#endif
#include "MemoryMapArray.h"
#include <stdint.h>
//
// ChromosomeInfo represents the per chromosome information
// necessary to write out SAM/BAM records. In addition, it
// contains a single internal index used to point to the vector
// offset where the chromosome bases start.
//
// This is mildly non-optimal for larger collections of chromosomes
// or contigs - one use case described having millions of contigs,
// in which case, this structure alone would take a gigabyte or more.
//
struct ChromosomeInfo
{
static const int MAX_GENOME_INFO_STRING=128;
void constructorClear()
{
memset(this,0, sizeof(*this));
}
void setChromosomeName(const char *n)
{
strncpy(name, n, sizeof(name)-1);
name[sizeof(name)-1] = '\0';
}
genomeIndex_t start; // internal offset to combined genome vector
genomeIndex_t size; // SAM SQ:LN value
char md5[2*MD5_DIGEST_LENGTH + 1]; // 32 chars plus NUL, SAM SQ:M5 value
char name[MAX_GENOME_INFO_STRING]; // SAM SQ:SN value
char assemblyID[MAX_GENOME_INFO_STRING]; // SAM SQ:AS value
char uri[MAX_GENOME_INFO_STRING]; // SAM SQ:UR value
char species[MAX_GENOME_INFO_STRING]; // SAM SQ:SP value
// handy setting methods:
void setAssemblyID(const char *newID)
{
strncpy(assemblyID, newID, sizeof(assemblyID)-1);
name[sizeof(name)-1] = '\0';
}
void setSpecies(const char *newSpecies)
{
strncpy(species, newSpecies, sizeof(species)-1);
species[sizeof(species)-1] = '\0';
}
void setURI(const char *newURI)
{
strncpy(uri, newURI, sizeof(uri)-1);
uri[sizeof(uri)-1] = '\0';
}
};
class genomeSequenceMmapHeader : public MemoryMapArrayHeader
{
friend class GenomeSequence;
friend std::ostream &operator << (std::ostream &, genomeSequenceMmapHeader &);
private:
uint32_t _chromosomeCount;
bool _colorSpace;
ChromosomeInfo _chromosomes[0];
public:
//
// getHeaderSize is special in that it must not access any
// member variables, since it is called when the header has
// not been created yet.
//
static size_t getHeaderSize(int chromosomeCount)
{
return sizeof(genomeSequenceMmapHeader) + sizeof(ChromosomeInfo[1]) * chromosomeCount;
}
//
// below methods return TRUE if it failed, false otherwise (primarily
// a length check).
//
};
std::ostream &operator << (std::ostream &stream, genomeSequenceMmapHeader &h);
//
// define the genomeSequence array type:
//
// NB the access/set routines use the encoded base values in the range
// 0-15, not the corresponding base pair character.
//
inline uint8_t genomeSequenceAccess(void *base, genomeIndex_t index)
{
if ((index&1)==0)
{
return ((uint8_t *) base)[index>>1] & 0xf;
}
else
{
return (((uint8_t *) base)[index>>1] & 0xf0) >> 4;
}
};
inline void genomeSequenceSet(void *base, genomeIndex_t index, uint8_t v)
{
if ((index&1)==0)
{
((uint8_t *) base)[index>>1] = (((uint8_t *) base)[index>>1] & 0xf0) | v;
}
else
{
((uint8_t *) base)[index>>1] = (((uint8_t *) base)[index>>1] & 0x0f) | v<<4;
}
}
inline size_t mmapGenomeSequenceElementCount2Bytes(genomeIndex_t i)
{
return sizeof(uint8_t) * i / 2;
}
class PackedRead
{
void set(int index, int val)
{
packedBases[index>>1] =
(packedBases[index>>1] // original value
& ~(7<<((index&0x01)<<2))) // logical AND off the original value
| ((val&0x0f)<<((index&0x1)<<2)); // logical OR in the new value
}
public:
std::vector<uint8_t> packedBases;
uint32_t length;
int size()
{
return length;
}
void clear()
{
packedBases.clear();
length=0;
}
void set(const char *rhs, int padWithNCount = 0);
uint8_t operator [](int index)
{
return (packedBases[index>>1] >> ((index&0x1)<<2)) & 0xf;
}
};
#endif
operator<<
InputFile & operator<<(InputFile &stream, const std::string &str)
Write to a file using streaming.
Definition: InputFile.h:736
GenomeSequence
Create/Access/Modify/Load Genome Sequences stored as binary mapped files.
Definition: GenomeSequence.h:100
MemoryMapArrayHeader
Definition: MemoryMapArray.h:65
PackedRead
Definition: GenomeSequenceHelpers.h:142
genomeSequenceMmapHeader
Definition: GenomeSequenceHelpers.h:80
ChromosomeInfo
Definition: GenomeSequenceHelpers.h:41