libStatGen Software 1
PileupElementBaseQual.cpp
1/*
2 * Copyright (C) 2010 Regents of the University of Michigan
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18#include <stdexcept>
19
20#include "PileupElementBaseQual.h"
21
22PileupElementBaseQual::PileupElementBaseQual()
23 : PileupElement(),
24 myBases(NULL),
25 myQualities(NULL),
26 myAllocatedSize(0),
27 myIndex(-1),
28 myAddDelAsBase(false)
29{
30 myAllocatedSize = 1024;
31 myBases = (char*)malloc(myAllocatedSize + 1);
32 myQualities = (char*)malloc(myAllocatedSize + 1);
33 if((myBases == NULL ) || (myQualities == NULL))
34 {
35 // TODO, check for malloc failures.
36 std::cerr << "Failed Memory Allocation\n";
37 }
38}
39
40// NOTE that this method does not actually copy, it just resets.
41PileupElementBaseQual::PileupElementBaseQual(const PileupElementBaseQual& q)
42 : PileupElement(),
43 myBases(NULL),
44 myQualities(NULL),
45 myAllocatedSize(0),
46 myIndex(-1)
47{
48 myAllocatedSize = 1024;
49 myBases = (char*)malloc(myAllocatedSize + 1);
50 myQualities = (char*)malloc(myAllocatedSize + 1);
51 myAddDelAsBase = q.myAddDelAsBase;
52 if((myBases == NULL ) || (myQualities == NULL))
53 {
54 // TODO, check for malloc failures.
55 std::cerr << "Failed Memory Allocation\n";
56 }
57}
58
59
60PileupElementBaseQual::~PileupElementBaseQual()
61{
62 if(myBases != NULL)
63 {
64 free(myBases);
65 myBases = NULL;
66 }
67 if(myQualities != NULL)
68 {
69 free(myQualities);
70 myQualities = NULL;
71 }
72}
73
74
75// Add an entry to this pileup element.
77{
78 // Call the base class:
80
81 // Increment the index
82 ++myIndex;
83
84 // if the index has gone beyond the allocated space, double the size.
85 if(myIndex >= myAllocatedSize)
86 {
87 char* tempBuffer = (char*)realloc(myBases, myAllocatedSize * 2);
88 if(tempBuffer == NULL)
89 {
90 std::cerr << "Memory Allocation Failure\n";
91 // TODO
92 return;
93 }
94 myBases = tempBuffer;
95 tempBuffer = (char*)realloc(myQualities, myAllocatedSize * 2);
96 if(tempBuffer == NULL)
97 {
98 std::cerr << "Memory Allocation Failure\n";
99 // TODO
100 return;
101 }
102 myQualities = tempBuffer;
103 myAllocatedSize = myAllocatedSize * 2;
104 }
105
106 Cigar* cigar = record.getCigarInfo();
107
108 if(cigar == NULL)
109 {
110 throw std::runtime_error("Failed to retrieve cigar info from the record.");
111 }
112
113
114 int32_t readIndex =
116
117 // If the readPosition is N/A, this is a deletion.
118 if(readIndex != CigarRoller::INDEX_NA)
119 {
120 char base = record.getSequence(readIndex);
121 char qual = record.getQuality(readIndex);
122 if(qual == UNSET_QUAL)
123 {
124 qual = ' ';
125 }
126 myBases[myIndex] = base;
127 myQualities[myIndex] = qual;
128 }
129 else if(myAddDelAsBase)
130 {
131 // This version adds deletions as bases.
132 myBases[myIndex] = '-';
133 myQualities[myIndex] = '0';
134 }
135 else
136 {
137 // Do not add a deletion.
138 // Did not add any entries, so decrement the index counter since the
139 // index was not used.
140 --myIndex;
141 }
142}
143
145{
146 if(getRefPosition() != UNSET_POSITION)
147 {
148 myBases[myIndex+1] = '\0';
149 myQualities[myIndex+1] = '\0';
150 std::cout << getChromosome() << "\t" << getRefPosition() << "\tN\t" << myIndex+1 << "\t";
151 std::cout << myBases << "\t";
152 std::cout << myQualities;
153 std::cout << "\n";
154 }
155}
156
157void PileupElementBaseQual::reset(int refPosition)
158{
159 // Call the base class.
160 PileupElement::reset(refPosition);
161
162 myIndex = -1;
163}
164
This class represents the CIGAR without any methods to set the cigar (see CigarRoller for that).
Definition: Cigar.h:84
static const int32_t INDEX_NA
Value associated with an index that is not applicable/does not exist, used for converting between que...
Definition: Cigar.h:492
int32_t getQueryIndex(int32_t refOffset)
Return the query index associated with the specified reference offset or INDEX_NA based on this cigar...
Definition: Cigar.cpp:202
This class inherits from the base class and stores base and qualities.
virtual void reset(int32_t refPosition)
Resets the entry, setting the new position associated with this element.
virtual void analyze()
Perform the analysis associated with this class.
virtual void addEntry(SamRecord &record)
Add an entry to this pileup element.
This is a base class pileup component, representing the information for one reference position.
Definition: PileupElement.h:27
const char * getChromosome() const
Get the chromosome name stored in this element.
Definition: PileupElement.h:51
int32_t getRefPosition() const
Get the reference position stored in this element.
Definition: PileupElement.h:54
virtual void addEntry(SamRecord &record)
Add an entry to this pileup element.
virtual void reset(int32_t refPosition)
Resets the entry, setting the new position associated with this element.
Class providing an easy to use interface to get/set/operate on the fields in a SAM/BAM record.
Definition: SamRecord.h:52
Cigar * getCigarInfo()
Returns a pointer to the Cigar object associated with this record.
Definition: SamRecord.cpp:1836
int32_t get0BasedPosition()
Get the 0-based(BAM) leftmost position of the record.
Definition: SamRecord.cpp:1319
const char * getQuality()
Returns the SAM formatted quality string (QUAL).
Definition: SamRecord.cpp:1638
const char * getSequence()
Returns the SAM formatted sequence string (SEQ), translating the base as specified by setSequenceTran...
Definition: SamRecord.cpp:1568