1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
// Author: Mark Chaisson
#pragma once
#include <alignment/datastructures/alignment/AlignmentCandidate.hpp>
#include <pbdata/SMRTSequence.hpp>
#include <iostream>
#include <string>
#include <vector>
class ReadAlignments
{
public:
/*
This class stores the alignments from a read. A read may be
aligned in several different modes:
1. Fullread - Treat the read as a unit from start to end
2. Subread - Align each subread independently
3. CCSDeNovo - Only align the CCS sequence from a read
4. CCSAllPass - Align the de novo ccs sequences and then the
subreads to where the denovo ccs aligned.
5. CCSFullPass - Same as allpass, except using only complete
subreads.
6. ZmwSubreads - Align subreads of each zmw to where the longest
subread of the zmw aligned to.
The alignments are a raggad array of n sequences; n is 1 for cases
1 and 3, the number of subreads for cases 2 and 4, and the number
of full length passes for case 5.
A ReadAligments class must only have alignments for a single type
of read in it.
*/
std::vector<std::vector<T_AlignmentCandidate *> > subreadAlignments;
std::vector<SMRTSequence> subreads;
AlignMode alignMode;
SMRTSequence read;
inline int GetNAlignedSeq();
inline bool AllSubreadsHaveAlignments();
inline void Clear();
inline void Resize(int nSeq);
inline void CheckSeqIndex(int seqIndex);
inline void SetSequence(int seqIndex, SMRTSequence &seq);
inline void AddAlignmentForSeq(int seqIndex, T_AlignmentCandidate *alignmentPtr);
inline void AddAlignmentsForSeq(int seqIndex,
std::vector<T_AlignmentCandidate *> &seqAlignmentPtrs);
// Copy all T_AlignmentCandidate objects (to which subreadAlignment[seqIndex]
// is pointing) to newly created objects, and then return pointers to the new
// objects.
inline std::vector<T_AlignmentCandidate *> CopySubreadAlignments(int seqIndex);
inline void Print(std::ostream &out = std::cout);
inline ~ReadAlignments();
};
inline int ReadAlignments::GetNAlignedSeq() { return subreadAlignments.size(); }
inline bool ReadAlignments::AllSubreadsHaveAlignments()
{
int i, nAlignedSeq;
nAlignedSeq = subreadAlignments.size();
for (i = 0; i < nAlignedSeq; i++) {
if (subreadAlignments[i].size() == 0) {
return false;
}
}
return true;
}
inline void ReadAlignments::Clear()
{
int i;
int nAlignedSeq;
for (i = 0, nAlignedSeq = subreadAlignments.size(); i < nAlignedSeq; i++) {
int nAlignments;
int a;
for (a = 0, nAlignments = subreadAlignments[i].size(); a < nAlignments; a++) {
delete subreadAlignments[i][a];
}
subreadAlignments[i].clear();
}
for (i = 0, nAlignedSeq = subreads.size(); i < nAlignedSeq; i++) {
subreads[i].Free();
}
subreadAlignments.clear();
read.Free();
}
inline void ReadAlignments::Resize(int nSeq)
{
subreadAlignments.resize(nSeq);
subreads.resize(nSeq);
}
inline void ReadAlignments::CheckSeqIndex(int seqIndex)
{
if (seqIndex < 0 or seqIndex >= int(subreads.size())) {
std::cout << "ERROR, adding a sequence to an unallocated position." << std::endl;
assert(0);
}
}
inline void ReadAlignments::SetSequence(int seqIndex, SMRTSequence &seq)
{
CheckSeqIndex(seqIndex);
subreads[seqIndex] = seq;
}
inline void ReadAlignments::AddAlignmentForSeq(int seqIndex, T_AlignmentCandidate *alignmentPtr)
{
CheckSeqIndex(seqIndex);
subreadAlignments[seqIndex].push_back(alignmentPtr);
}
inline void ReadAlignments::AddAlignmentsForSeq(
int seqIndex, std::vector<T_AlignmentCandidate *> &seqAlignmentPtrs)
{
CheckSeqIndex(seqIndex);
subreadAlignments[seqIndex].insert(subreadAlignments[seqIndex].end(), seqAlignmentPtrs.begin(),
seqAlignmentPtrs.end());
}
inline std::vector<T_AlignmentCandidate *> ReadAlignments::CopySubreadAlignments(int seqIndex)
{
std::vector<T_AlignmentCandidate *> ret;
for (int i = 0; i < int(subreadAlignments[seqIndex].size()); i++) {
T_AlignmentCandidate *q = new T_AlignmentCandidate();
*q = *(subreadAlignments[seqIndex][i]);
ret.push_back(q);
}
return ret;
}
inline void ReadAlignments::Print(std::ostream &out)
{
out << "A ReadAlignments object with " << subreadAlignments.size()
<< " groups of subread alignments." << std::endl;
for (int i = 0; i < int(subreadAlignments.size()); i++) {
out << " subreadAlignment group [" << i << "/" << subreadAlignments.size() << "] has "
<< subreadAlignments[i].size() << " alignments." << std::endl;
for (int j = 0; j < int(subreadAlignments[i].size()); j++) {
out << " [" << i << "][" << j << "/" << subreadAlignments[i].size() << "]"
<< std::endl;
subreadAlignments[i][j]->Print(out);
}
}
out << " read: ";
read.Print(out);
out << std::endl << std::endl;
}
inline ReadAlignments::~ReadAlignments() { read.Free(); }
|
