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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
|
/*++
Module Name:
AlignmentAdjuster.cpp
Abstract:
Code to adjust alignments to avoid indels at ends, or overhanging contig boundaries
Authors:
Bill Bolosky, May, 2016
Environment:
User mode service.
Revision History:
Pulled out of other places in SNAP
--*/
#include "stdafx.h"
#include "AlignmentAdjuster.h"
AlignmentAdjuster::AlignmentAdjuster(const Genome *_genome)
{
genome = _genome;
}
void
AlignmentAdjuster::AdjustAlignment(Read *read, SingleAlignmentResult *result)
{
int additionalFrontClipping = 0;
if (result->status == NotFound) {
//
// No adjustment required on unaligned reads.
//
return;
}
read->setAdditionalFrontClipping(0);
int cumulativeAddFrontClipping = 0;
unsigned nAdjustments = 0;
const int cigarBufLen = 1000;
char cigarBuf[cigarBufLen + 1];
char dataBuffer[MAX_READ_LENGTH];
const char *data = dataBuffer;
const char *clippedData;
int netIndel = 0;
GenomeDistance dataLength;
const Genome::Contig *contig;
GenomeDistance extraBasesClippedAfter;
const char *reference;
int cigarBufUsed = 0;
for (;;) {
// Write the data and quality strings. If the read is reverse complemented, these need to
// be backwards from the original read. Also, both need to be unclipped.
GenomeDistance fullLength = read->getUnclippedLength();
GenomeDistance basesClippedBefore;
GenomeDistance basesClippedAfter;
dataLength = read->getDataLength();
if (fullLength > MAX_READ_LENGTH) {
WriteErrorMessage("AlignmentAdjuster: read too long (how did it get past the reader?), %d > %d\n", fullLength, MAX_READ_LENGTH);
soft_exit(1);
}
if (result->direction == RC) {
for (unsigned i = 0; i < fullLength; i++) {
dataBuffer[fullLength - 1 - i] = COMPLEMENT[read->getUnclippedData()[i]];
}
clippedData = &data[fullLength - dataLength - read->getFrontClippedLength()];
basesClippedBefore = fullLength - dataLength - read->getFrontClippedLength();
basesClippedAfter = read->getFrontClippedLength();
} else {
memcpy(dataBuffer, read->getUnclippedData(), read->getUnclippedLength());
clippedData = read->getData();
basesClippedBefore = read->getFrontClippedLength();
basesClippedAfter = fullLength - dataLength - basesClippedBefore;
}
GenomeDistance extraBasesClippedBefore = 0;
contig = genome->getContigForRead(result->location, read->getDataLength(), &extraBasesClippedBefore);
_ASSERT(NULL != contig && contig->length > genome->getChromosomePadding());
GenomeLocation genomeLocation = result->location + extraBasesClippedBefore;
clippedData += extraBasesClippedBefore;
dataLength -= extraBasesClippedBefore;
if (genomeLocation + dataLength > contig->beginningLocation + contig->length - genome->getChromosomePadding()) {
//
// The read hangs off the end of the contig. Soft clip it at the end. This is a tentative amount that assumes no net indels in the
// mapping, we'll refine it later if needed.
//
extraBasesClippedAfter = genomeLocation + dataLength - (contig->beginningLocation + contig->length - genome->getChromosomePadding());
} else {
extraBasesClippedAfter = 0;
}
reference = genome->getSubstring(genomeLocation, dataLength - extraBasesClippedAfter);
_ASSERT(NULL != reference); // Since we just clipped it to fit.
if (NULL == reference) { // This shouldn't happen, but the old code checked it, so I'll leave it in.
result->status = NotFound;
result->location = InvalidGenomeLocation;
return;
}
result->score = lv.computeEditDistanceNormalized(
reference,
(int)(dataLength - extraBasesClippedAfter + MAX_K), // Add space incase of indels. We know there's enough, because the reference is padded.
clippedData,
(int)(dataLength - extraBasesClippedAfter),
MAX_K - 1,
cigarBuf,
cigarBufLen,
true, // useM (doesn't matter, since we're throwing the cigar string away anyway)
COMPACT_CIGAR_STRING,
&cigarBufUsed,
&additionalFrontClipping,
&netIndel);
if (0 == additionalFrontClipping) {
break;
}
nAdjustments++;
// redo if read modified (e.g. to add soft clipping, or move alignment for a leading I.
const Genome::Contig *originalContig = genome->getContigAtLocation(result->location);
const Genome::Contig *newContig = genome->getContigAtLocation(result->location + additionalFrontClipping);
if (newContig == NULL || newContig != originalContig || result->location + additionalFrontClipping > originalContig->beginningLocation + originalContig->length - genome->getChromosomePadding() ||
nAdjustments > read->getDataLength()) {
//
// Altering this would push us over a contig boundary, or we're stuck in a loop. Just give up on the read.
//
result->status = NotFound;
result->location = InvalidGenomeLocation;
return;
} else {
cumulativeAddFrontClipping += additionalFrontClipping;
read->setAdditionalFrontClipping(__max(0,cumulativeAddFrontClipping));
result->clippingForReadAdjustment = __max(0, cumulativeAddFrontClipping);
result->location = result->location + additionalFrontClipping;
}
} // for ever
//
// Code copied from SAM.cpp:computeCigar()
//
// Normally, we'd be done. However, if the amount that we would clip at the end of the read because of hanging off of the end
// of the contig changed, then we need to recompute. In some cases this is an iterative processess as we add or remove bits
// of read.
//
GenomeDistance newExtraBasesClippedAfter = __max(0, result->location + dataLength + netIndel - (contig->beginningLocation + contig->length - genome->getChromosomePadding()));
for (GenomeDistance pass = 0; pass < dataLength; pass++) {
if (newExtraBasesClippedAfter == extraBasesClippedAfter) {
return;
}
extraBasesClippedAfter = newExtraBasesClippedAfter;
result->score = lv.computeEditDistanceNormalized(
reference,
(int)(dataLength - extraBasesClippedAfter + MAX_K), // Add space incase of indels. We know there's enough, because the reference is padded.
data,
(int)(dataLength - extraBasesClippedAfter),
MAX_K - 1,
cigarBuf,
cigarBufLen,
true, // useM (doesn't matter, since we're throwing the cigar string away anyway)
COMPACT_CIGAR_STRING,
&cigarBufUsed,
&additionalFrontClipping,
&netIndel);
newExtraBasesClippedAfter = __max(0, result->location + dataLength + netIndel - (contig->beginningLocation + contig->length - genome->getChromosomePadding()));
}
_ASSERT(!"cigar computation didn't converge");
}
void
AlignmentAdjuster::AdjustAlignments(Read **reads, PairedAlignmentResult *result)
{
for (int i = 0; i < NUM_READS_PER_PAIR; i++) {
SingleAlignmentResult singleResult;
singleResult.direction = result->direction[i];
singleResult.location = result->location[i];
singleResult.mapq = result->mapq[i];
singleResult.score = result->score[i];
singleResult.status = result->status[i];
singleResult.clippingForReadAdjustment = 0;
AdjustAlignment(reads[i], &singleResult);
result->direction[i] = singleResult.direction;
result->location[i] = singleResult.location;
result->mapq[i] = singleResult.mapq;
result->score[i] = singleResult.score;
result->status[i] = singleResult.status;
result->clippingForReadAdjustment[i] = singleResult.clippingForReadAdjustment;
}
}
|
