1 // ==========================================================================
2 //                                  Gustaf
3 // ==========================================================================
4 // Copyright (c) 2011-2018, Kathrin Trappe, FU Berlin
5 // All rights reserved.
6 //
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are met:
9 //
10 //     * Redistributions of source code must retain the above copyright
11 //       notice, this list of conditions and the following disclaimer.
12 //     * Redistributions in binary form must reproduce the above copyright
13 //       notice, this list of conditions and the following disclaimer in the
14 //       documentation and/or other materials provided with the distribution.
15 //     * Neither the name of Knut Reinert or the FU Berlin nor the names of
16 //       its contributors may be used to endorse or promote products derived
17 //       from this software without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE
23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
29 // DAMAGE.
30 //
31 // ==========================================================================
32 // Author: Kathrin Trappe <kathrin.trappe@fu-berlin.de>
33 // ==========================================================================
34 
35 #ifndef SEQAN_APPS_GUSTAF_GUSTAF_MATEPAIRS_H_
36 #define SEQAN_APPS_GUSTAF_GUSTAF_MATEPAIRS_H_
37 
38 #include <seqan/basic.h>
39 #include <seqan/sequence.h>
40 #include <seqan/seq_io.h>
41 #include "stellar_routines.h"
42 
43 using namespace seqan;
44 
45 // Imports mate pairs from two files, joins them, stores the joining position in
46 // String readJoinPositions, and stores the sequences in the StringSet seqs and
47 // their identifiers in the StringSet ids
48 template <typename TSequence, typename TId>
49 inline bool
_importSequences(CharString const & fileNameL,CharString const & fileNameR,bool revCompl,StringSet<TSequence> & seqs,StringSet<TId> & ids,StringSet<CharString> & quals,String<unsigned> & readJoinPositions)50 _importSequences(CharString const & fileNameL,
51                  CharString const & fileNameR,
52                  bool revCompl,
53                  StringSet<TSequence> & seqs,
54                  StringSet<TId> & ids,
55                  StringSet<CharString> & quals,
56                  String<unsigned> & readJoinPositions
57                  )
58 {
59     typedef typename Iterator<StringSet<TId>, Standard>::Type TIdSetIterator;
60 
61     seqan::SeqFileIn l, r;
62     if (!open(l, toCString(fileNameL)) || !open(r, toCString(fileNameR)))
63     {
64         std::cerr << "Failed to open file." << std::endl;
65         return false;
66     }
67 
68     TSequence seq;
69     TSequence seqL;
70     TSequence seqR;
71     TId id;
72     TId sId;
73     CharString qual;
74     CharString qualL;
75     CharString qualR;
76     while (!atEnd(l) || !atEnd(r))
77     {
78         if (readRecord(id, seqL, qualL, l) != 0)
79         {
80             std::cerr << "Problem reading from first input file." << std::endl;
81             return false;
82         }
83         // Note: overwrites id with ID of right mate, which is done also by gustaf_join_mates
84         if (readRecord(id, seqR, qualR, r) != 0)
85         {
86             std::cerr << "Problem reading from first input file." << std::endl;
87             return false;
88         }
89 
90         appendValue(readJoinPositions, length(seqL));
91         if (revCompl)
92         {
93             reverseComplement(seqR);
94             reverse(qualR);
95         }
96         append(seq, seqL);
97         append(seq, seqR);
98         append(qual, qualL);
99         append(qual, qualR);
100         appendValue(seqs, seq, Generous());
101         appendValue(quals, qual, Generous());
102 
103         _getShortId(sId, id);
104         appendValue(ids, sId, Generous());
105         clear(seq);
106         clear(qual);
107     }
108 
109     // Check for dupliacte id entries.
110     StringSet<TId> uniqueIds = ids;
111     std::sort(begin(uniqueIds, Standard()), end(uniqueIds, Standard()), IdComparator<TId>());  // O(n*log(n))
112     TIdSetIterator itOldEnd = end(uniqueIds, Standard());
113     TIdSetIterator itNewEnd = std::unique(begin(uniqueIds, Standard()), end(uniqueIds, Standard()), IdComparator<TId>());  // O(n)
114 
115     --itNewEnd;
116     unsigned diff = itOldEnd - itNewEnd;
117     if (length(ids) - diff > 0)
118     {
119         std::cout << "Found nonunique sequence IDs" << std::endl;
120         return false;
121     }
122     std::cout << "Loaded " << length(ids) << " paired sequences" << std::endl;
123     return true;
124 }
125 
126 // /////////////////////////////////////////////////////////////////////////////
127 // Imports mate pairs from two files, joins them, stores the joining position in
128 // String readJoinPositions, and stores the sequences in the StringSet seqs and
129 // their identifiers in the StringSet ids
130 template <typename TSequence, typename TId>
131 inline bool
_importSequences(CharString const & fileNameL,CharString const & fileNameR,bool revCompl,StringSet<TSequence> & seqs,StringSet<TId> & ids,String<unsigned> & readJoinPositions)132 _importSequences(CharString const & fileNameL,
133                  CharString const & fileNameR,
134                  bool revCompl,
135                  StringSet<TSequence> & seqs,
136                  StringSet<TId> & ids,
137                  String<unsigned> & readJoinPositions
138                  )
139 {
140     typedef typename Iterator<StringSet<TId>, Standard>::Type TIdSetIterator;
141 
142     seqan::SeqFileIn leftMates, rightMates;
143     if (!open(leftMates, toCString(fileNameL)))
144     {
145         std::cerr << "Failed to open " << fileNameL << " file." << std::endl;
146         return false;
147     }
148     if (!open(rightMates, toCString(fileNameR)))
149     {
150         std::cerr << "Failed to open " << fileNameR << " file." << std::endl;
151         return false;
152     }
153 
154     TSequence seq;
155     TSequence seqL;
156     TSequence seqR;
157     TId idL;
158     TId idR;
159     TId sId;
160     unsigned seqCount = 0;
161     for (unsigned i = 0; !atEnd(leftMates) && !atEnd(rightMates); ++i, ++seqCount)
162     {
163         try
164         {
165             readRecord(idL, seqL, leftMates);
166             readRecord(idR, seqR, rightMates);
167         }
168         catch (seqan::IOError const & ioErr)
169         {
170             std::cerr << "Problem reading mates (" << ioErr.what() << ")\n";
171             return false;
172         }
173 
174         if (revCompl)
175             reverseComplement(seqR);
176         appendValue(readJoinPositions, length(seqL));
177         append(seq, seqL);
178         append(seq, seqR);
179         // Note: saving ID of right(!) mate per default
180         appendValue(seqs, seq, Generous());
181 
182         _getShortId(sId, idR);
183         appendValue(ids, sId, Generous());
184         clear(seq);
185     }
186     if (!atEnd(leftMates) || !atEnd(rightMates))
187     {
188         std::cerr << "Unequal number of left and right mates!\n";
189         return false;
190     }
191 
192     std::cout << "Loaded " << seqCount << " mate pair sequence" << ((seqCount > 1) ? "s." : ".") << std::endl;
193 
194     StringSet<TId> uniqueIds = ids;
195     // Check for dupliacte id entries.
196     std::sort(begin(uniqueIds, Standard()), end(uniqueIds, Standard()), IdComparator<TId>());  // O(n*log(n))
197     TIdSetIterator itOldEnd = end(uniqueIds, Standard());
198     TIdSetIterator itNewEnd = std::unique(begin(uniqueIds, Standard()), end(uniqueIds, Standard()), IdComparator<TId>());  // O(n)
199 
200     --itNewEnd;
201     unsigned diff = itOldEnd - itNewEnd;
202     if (length(ids) - diff > 0)
203     {
204         std::cout << "Found nonunique sequence IDs" << std::endl;
205         return false;
206     }
207     std::cout << "Loaded " << length(ids) << " paired sequences" << std::endl;
208     return true;
209 }
210 template <typename TSequence, typename TId, typename TPos>
211 inline bool
_isLeftMate(StellarMatch<TSequence,TId> const & sMatch,TPos const & joiningPos)212 _isLeftMate(StellarMatch<TSequence, TId> const & sMatch, TPos const & joiningPos)
213 {
214     if (sMatch.begin2 > joiningPos)
215         return false;
216     if (sMatch.end2 < joiningPos)
217         return true;
218     if ((joiningPos - sMatch.begin2) > (sMatch.end2 - joiningPos))
219         return true;
220     return false;
221 }
222 
223 template <typename TString, typename TId, typename TMatch>
224 inline bool
_countMateMatches(TString const & confMateMatches,String<TMatch> const & queryMatches,TId & dbID,unsigned support)225 _countMateMatches(TString const & confMateMatches, String<TMatch> const & queryMatches, TId & dbID, unsigned support)
226 {
227     unsigned count = 0;
228     // StellarMatches in confMateMatches have the correct mate distance, but the IntervalTree does not check for correct
229     // database IDs, this is done here (i.e. only matches on the same reference are counted)
230     // Note: matches on the opposite strand in the correct distance are considered valid
231     for (unsigned i = 0; i < length(confMateMatches); ++i)
232     {
233         unsigned index = confMateMatches[i];
234         if (queryMatches[index].id == dbID)
235             ++count;
236     }
237     // Support must/should dependend on the size of mateIntervalTrees!
238     return count >= support;
239 }
240 
241 // sMatch is part of the left mate, valid (right) mate matches can only be upstream but cannot exceed database length
242 template <typename TMatch, typename TMSplazerChain>
243 inline bool
_checkRightMateMatches(TMatch const & sMatch,String<TMatch> const & queryMatches,TMSplazerChain const & gustafChain,MSplazerOptions const & options)244 _checkRightMateMatches(TMatch const & sMatch, String<TMatch> const & queryMatches, TMSplazerChain const & gustafChain, MSplazerOptions const & options)
245 {
246     // typedef typename TMatch::TPos TPos;
247     typedef unsigned TPos;
248     String<unsigned> confMateMatches;
249     // Begin and end position of valid region for mate matches
250     unsigned dbLength = length(source(sMatch.row1));
251     TPos mateIntervalBegin = _min(dbLength, sMatch.begin1 + options.libSize - gustafChain.mateJoinPosition - options.libError);
252     TPos mateIntervalEnd = _min(dbLength, sMatch.end1 + options.libSize - gustafChain.mateJoinPosition + options.libError);
253     SEQAN_ASSERT_GEQ_MSG(dbLength, mateIntervalBegin, "Calculated valid region for right mate exceeds database length");
254     SEQAN_ASSERT_GEQ_MSG(dbLength, mateIntervalEnd, "Calculated valid region for right mate exceeds database length");
255     findIntervals(confMateMatches, gustafChain.rightMateTree, mateIntervalBegin, mateIntervalEnd);
256     /*
257     std::cout << "checkRightMatches " << mateIntervalBegin << "," << mateIntervalEnd << ";";
258     for (unsigned i = 0; i < length(confMateMatches); ++i)
259         std::cout << confMateMatches[i] << ",";
260     std::cout << std::endl;
261     */
262     return _countMateMatches(confMateMatches, queryMatches, sMatch.id, options.mateSupport);
263 }
264 
265 // sMatch is part of the right mate, valid (left) mate matches can only be downstream but cannot be below 0
266 template <typename TMatch, typename TMSplazerChain>
267 inline bool
_checkLeftMateMatches(TMatch const & sMatch,String<TMatch> const & queryMatches,TMSplazerChain const & gustafChain,MSplazerOptions const & options)268 _checkLeftMateMatches(TMatch const & sMatch, String<TMatch> const & queryMatches, TMSplazerChain const & gustafChain, MSplazerOptions const & options)
269 {
270     // typedef typename TMatch::TPos TPos;
271     typedef unsigned TPos;
272     String<unsigned> confMateMatches;
273     // Begin and end position of valid region for mate matches
274     TPos mateIntervalBegin = _max(0u, sMatch.begin1 - options.libSize + gustafChain.mateJoinPosition - options.libError);
275     TPos mateIntervalEnd = _max(0u, sMatch.end1 - options.libSize + gustafChain.mateJoinPosition + options.libError);
276     SEQAN_ASSERT_GEQ_MSG(mateIntervalBegin, 0u, "Calculated valid region for left mate is < 0");
277     SEQAN_ASSERT_GEQ_MSG(mateIntervalEnd, 0u, "Calculated valid region for left mate is < 0");
278     findIntervals(confMateMatches, gustafChain.leftMateTree, mateIntervalBegin, mateIntervalEnd);
279     /*
280     std::cout << "checkLeftMatches " << mateIntervalBegin << "," << mateIntervalEnd << ";";
281     for (unsigned i = 0; i < length(confMateMatches); ++i)
282         std::cout << confMateMatches[i] << ",";
283     std::cout << std::endl;
284     */
285     return _countMateMatches(confMateMatches, queryMatches, sMatch.id, options.mateSupport);
286 }
287 
288 // Check for one match if it is confirmed by ANY of the other (mate) matches
289 template <typename TMatch, typename TMSplazerChain>
290 inline bool
_checkMateMatches(TMatch const & sMatch,String<TMatch> const & queryMatches,TMSplazerChain const & gustafChain,MSplazerOptions const & options)291 _checkMateMatches(TMatch const & sMatch, String<TMatch> const & queryMatches, TMSplazerChain const & gustafChain, MSplazerOptions const & options)
292 {
293     if (_isLeftMate(sMatch, gustafChain.mateJoinPosition))
294             return _checkRightMateMatches(sMatch, queryMatches, gustafChain, options);
295     return _checkLeftMateMatches(sMatch, queryMatches, gustafChain, options);
296 }
297 
298 
299 // Check for two matches, each from a different mate, if they both apply to the libSize+sd, i.e. the BP is artificial
300 // Assumptions: sMatch1 < sMatch2, valid order and valid gap between matches regarding read sequence
301 template <typename TMatch, typename TMSplazerChain>
302 inline bool
_artificialBP(TMatch const & sMatch1,TMatch const & sMatch2,TMSplazerChain const & gustafChain)303 _artificialBP(TMatch const & sMatch1, TMatch const & sMatch2, TMSplazerChain const & gustafChain)//, MSplazerOptions const & options)
304 {
305     // Check if both matches are from different mates
306     if (_isLeftMate(sMatch1, gustafChain.mateJoinPosition) && _isLeftMate(sMatch2, gustafChain.mateJoinPosition))
307         return false;
308     if (!_isLeftMate(sMatch1, gustafChain.mateJoinPosition) && !_isLeftMate(sMatch2, gustafChain.mateJoinPosition))
309         return false;
310     // Check the distance between inner match position on the reference
311     /*
312     typedef typename TMatch::TPos TPos;
313     TPos dist = sMatch2.begin1 - sMatch1.end1;
314     if (dist < (options.libSize - options.libError))
315         return false;
316     if (dist > (options.libSize + options.libError))
317         return false;
318         */
319     return true;
320 }
321 
322 
323 // Intitialisation of graph structure for combinable StellarMatches of a read
324 template <typename TSequence, typename TId, typename TMSplazerChain>
_initialiseGraphMatePairsNoBreakends(QueryMatches<StellarMatch<TSequence,TId>> & queryMatches,TMSplazerChain & chain,MSplazerOptions const & options)325 void _initialiseGraphMatePairsNoBreakends(QueryMatches<StellarMatch<TSequence, TId> > & queryMatches,
326                       TMSplazerChain & chain,
327                       MSplazerOptions const & options)
328 {
329     // std::cerr << " Initialising graph structure " << std::endl;
330     typedef typename TMSplazerChain::TGraph TGraph;
331     typedef typename EdgeDescriptor<TGraph>::Type TEdgeDescriptor;
332     typedef typename Iterator<String<StellarMatch<TSequence, TId> > >::Type TIterator;
333     typedef typename TMSplazerChain::TInterval TInterval;
334 
335     String<TInterval> leftMateMatches;
336     String<TInterval> rightMateMatches;
337     TIterator itStellarMatches = begin(queryMatches.matches);
338     TIterator itEndStellarMatches = end(queryMatches.matches);
339     // The default vertex descriptor is an integer. So if inserted in the same order, the vertex descriptor value is
340     // the same as the position of the corresponding vertex within the QueryMatches --> since we can easily iterate
341     // through the QueryMatches and use the iterator, we wont keep track of the vertex descriptors.
342     unsigned index = 0;
343     for (; itStellarMatches < itEndStellarMatches; goNext(itStellarMatches))
344     {
345         addVertex(chain.graph);
346         // Collect intervals for mate IntervalTrees
347         if (_isLeftMate(*itStellarMatches, chain.mateJoinPosition))
348             appendValue(leftMateMatches, TInterval((*itStellarMatches).begin1, (*itStellarMatches).end1, index));
349         else
350             appendValue(rightMateMatches, TInterval((*itStellarMatches).begin1, (*itStellarMatches).end1, index));
351         ++index;
352     }
353 
354     /*
355     std::cout << "Left " << std::endl;
356     for (unsigned i = 0; i < length(leftMateMatches); ++i)
357         std::cout << leftMateMatches[i].i1 << "," << leftMateMatches[i].i2 << "," << leftMateMatches[i].cargo << " ";
358     std::cout << std::endl;
359     std::cout << "Right " << std::endl;
360     for (unsigned i = 0; i < length(rightMateMatches); ++i)
361         std::cout << rightMateMatches[i].i1 << "," << rightMateMatches[i].i2 << "," << rightMateMatches[i].cargo << " ";
362     std::cout << std::endl;
363     */
364 
365     // Create IntervalTree for mates
366     createIntervalTree(chain.rightMateTree, rightMateMatches);
367     createIntervalTree(chain.leftMateTree, leftMateMatches);
368 
369     // std::cerr << " Created graph " << std::endl;
370     // Add start and end to graph and property map
371     chain.startVertex = addVertex(chain.graph);
372     chain.endVertex = addVertex(chain.graph);
373 
374     int cargo = 0;
375     resize(chain.breakpoints.slotLookupTable, 2 * length(queryMatches.matches));
376     // Adding edges to start and end vertices
377     for (unsigned i = 0; i < length(queryMatches.matches); ++i)
378     {
379         cargo = static_cast<int>(queryMatches.matches[i].begin2);
380         if (cargo < (options.initGapThresh + 1))
381         {
382             cargo += chain.matchDistanceScores[i];
383             // Note that the penalty for no mate matches is also always on the ingoing edge,
384             // so for the graph initialsation, we only check this for the start vertex.
385             if (!_checkMateMatches(queryMatches.matches[i], queryMatches.matches, chain, options))
386                 cargo += options.noMateMatchesPen;
387             TEdgeDescriptor edge = addEdge(chain.graph, chain.startVertex, i, cargo);
388             resizeEdgeMap(chain.breakpoints.slotLookupTable, chain.graph);
389             assignProperty(chain.breakpoints, edge);
390         }
391         cargo = static_cast<int>(length(source(queryMatches.matches[i].row2))) -
392                 static_cast<int>(queryMatches.matches[i].end2);
393         if (cargo < (options.initGapThresh + 1))
394         {
395             TEdgeDescriptor edge = addEdge(chain.graph, i, chain.endVertex, cargo);
396             resizeEdgeMap(chain.breakpoints.slotLookupTable, chain.graph);
397             assignProperty(chain.breakpoints, edge);
398         }
399     }
400 }
401 
402 // Intitialisation of graph structure for combinable StellarMatches of a read
403 // including breakends
404 template <typename TSequence, typename TId, typename TMSplazerChain>
_initialiseGraphMatePairs(QueryMatches<StellarMatch<TSequence,TId>> & queryMatches,TId & queryId,TMSplazerChain & chain,MSplazerOptions const & options)405 void _initialiseGraphMatePairs(QueryMatches<StellarMatch<TSequence, TId> > & queryMatches,
406                       TId & queryId,
407                       TMSplazerChain & chain,
408                       MSplazerOptions const & options)
409 {
410 
411     // std::cerr << " Initialising graph structure " << std::endl;
412     typedef typename TMSplazerChain::TGraph TGraph;
413     typedef typename EdgeDescriptor<TGraph>::Type TEdgeDescriptor;
414     typedef typename Iterator<String<StellarMatch<TSequence, TId> > >::Type TIterator;
415     typedef typename TMSplazerChain::TInterval TInterval;
416     typedef Breakpoint<TSequence, TId> TBreakpoint;
417 
418     String<TInterval> leftMateMatches;
419     String<TInterval> rightMateMatches;
420     TIterator itStellarMatches = begin(queryMatches.matches);
421     TIterator itEndStellarMatches = end(queryMatches.matches);
422     // The default vertex descriptor is an integer. So if inserted in the same order the vertex descriptor value is
423     // the same as the position of the corresponding vertex within the QueryMatches --> since we can easily iterate
424     // through the QueryMatches and use the iterator we wont keep track of the vertex descriptors
425     unsigned index = 0;
426     for (; itStellarMatches < itEndStellarMatches; goNext(itStellarMatches))
427     {
428         addVertex(chain.graph);
429         // Collect intervals for mate IntervalTrees
430         if (_isLeftMate(*itStellarMatches, chain.mateJoinPosition))
431             appendValue(leftMateMatches, TInterval((*itStellarMatches).begin1, (*itStellarMatches).end1, index));
432         else
433             appendValue(rightMateMatches, TInterval((*itStellarMatches).begin1, (*itStellarMatches).end1, index));
434         ++index;
435     }
436 
437     // Create IntervalTree for mates
438     createIntervalTree(chain.rightMateTree, rightMateMatches);
439     createIntervalTree(chain.leftMateTree, leftMateMatches);
440 
441     // std::cerr << " Created graph " << std::endl;
442     // Add start and end to graph and property map
443     chain.startVertex = addVertex(chain.graph);
444     chain.endVertex = addVertex(chain.graph);
445 
446     int cargo = 0;
447     resize(chain.breakpoints.slotLookupTable, 2 * length(queryMatches.matches));
448     // Adding edges to start and end vertices
449     for (unsigned i = 0; i < length(queryMatches.matches); ++i)
450     {
451         // start vertex
452         cargo = static_cast<int>(queryMatches.matches[i].begin2);
453         if (cargo < (options.breakendThresh + 1))
454         {
455             cargo += chain.matchDistanceScores[i];
456             // Note that the penalty for no mate matches is also always on the ingoing edge,
457             // so for the graph initialsation, we only check this for the start vertex.
458             if (!_checkMateMatches(queryMatches.matches[i], queryMatches.matches, chain, options))
459                 cargo += options.noMateMatchesPen;
460             TEdgeDescriptor edge = addEdge(chain.graph, chain.startVertex, i, cargo);
461             resizeEdgeMap(chain.breakpoints.slotLookupTable, chain.graph);
462             if (cargo < (options.initGapThresh + 1))
463                 assignProperty(chain.breakpoints, edge);
464             else
465             {
466                 // TODO(ktrappe): needs trimming of x-drop/sloppy end
467                 TBreakpoint bp(queryMatches.matches[i].id,
468                                queryMatches.matches[i].id,
469                                queryMatches.matches[i].orientation,
470                                queryMatches.matches[i].orientation,
471                                queryMatches.matches[i].begin1,
472                                queryMatches.matches[i].begin1,
473                                queryMatches.matches[i].begin2,
474                                queryMatches.matches[i].begin2,
475                                queryId);
476                 bp.svtype = TBreakpoint::BREAKEND;
477                 // bp.imprecise = true;
478                 bp.breakend = 0; // left breakend
479                 assignProperty(chain.breakpoints, edge, bp);
480             }
481         }
482         // end vertex
483         cargo = static_cast<int>(length(source(queryMatches.matches[i].row2))) -
484                 static_cast<int>(queryMatches.matches[i].end2);
485         if (cargo < (options.breakendThresh + 1))
486         {
487             TEdgeDescriptor edge = addEdge(chain.graph, i, chain.endVertex, cargo);
488             resizeEdgeMap(chain.breakpoints.slotLookupTable, chain.graph);
489             if (cargo < (options.initGapThresh + 1))
490                 assignProperty(chain.breakpoints, edge);
491             else
492             {
493                 TBreakpoint bp(queryMatches.matches[i].id,
494                                queryMatches.matches[i].id,
495                                queryMatches.matches[i].orientation,
496                                queryMatches.matches[i].orientation,
497                                queryMatches.matches[i].end1,
498                                queryMatches.matches[i].end1,
499                                queryMatches.matches[i].end2,
500                                queryMatches.matches[i].end2,
501                                queryId);
502                 bp.svtype = TBreakpoint::BREAKEND;
503                 // bp.imprecise = true;
504                 bp.breakend = 1; // right breakend
505                 assignProperty(chain.breakpoints, edge, bp);
506             }
507         }
508     }
509 }
510 
511 #endif  // #ifndef _APPS_GUSTAF_GUSTAF_MATEPAIRS_H_
512