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 
22 PileupElementBaseQual::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.
41 PileupElementBaseQual::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 
60 PileupElementBaseQual::~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 =
115  cigar->getQueryIndex(getRefPosition(), record.get0BasedPosition());
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 
157 void PileupElementBaseQual::reset(int refPosition)
158 {
159  // Call the base class.
160  PileupElement::reset(refPosition);
161 
162  myIndex = -1;
163 }
164 
Cigar
This class represents the CIGAR without any methods to set the cigar (see CigarRoller for that).
Definition: Cigar.h:84
SamRecord::get0BasedPosition
int32_t get0BasedPosition()
Get the 0-based(BAM) leftmost position of the record.
Definition: SamRecord.cpp:1307
PileupElementBaseQual::reset
virtual void reset(int32_t refPosition)
Resets the entry, setting the new position associated with this element.
Definition: PileupElementBaseQual.cpp:157
PileupElement::getChromosome
const char * getChromosome() const
Get the chromosome name stored in this element.
Definition: PileupElement.h:51
PileupElementBaseQual
This class inherits from the base class and stores base and qualities.
Definition: PileupElementBaseQual.h:26
PileupElementBaseQual::analyze
virtual void analyze()
Perform the analysis associated with this class.
Definition: PileupElementBaseQual.cpp:144
Cigar::INDEX_NA
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
PileupElement::reset
virtual void reset(int32_t refPosition)
Resets the entry, setting the new position associated with this element.
Definition: PileupElement.cpp:67
PileupElement::getRefPosition
int32_t getRefPosition() const
Get the reference position stored in this element.
Definition: PileupElement.h:54
SamRecord::getQuality
const char * getQuality()
Returns the SAM formatted quality string (QUAL).
Definition: SamRecord.cpp:1626
SamRecord::getSequence
const char * getSequence()
Returns the SAM formatted sequence string (SEQ), translating the base as specified by setSequenceTran...
Definition: SamRecord.cpp:1556
PileupElement::addEntry
virtual void addEntry(SamRecord &record)
Add an entry to this pileup element.
Definition: PileupElement.cpp:44
Cigar::getQueryIndex
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
SamRecord
Class providing an easy to use interface to get/set/operate on the fields in a SAM/BAM record.
Definition: SamRecord.h:52
SamRecord::getCigarInfo
Cigar * getCigarInfo()
Returns a pointer to the Cigar object associated with this record.
Definition: SamRecord.cpp:1824
PileupElement
This is a base class pileup component, representing the information for one reference position.
Definition: PileupElement.h:27
PileupElementBaseQual::addEntry
virtual void addEntry(SamRecord &record)
Add an entry to this pileup element.
Definition: PileupElementBaseQual.cpp:76