Skip to content

Commit

Permalink
Translocation bug fix.
Browse files Browse the repository at this point in the history
  • Loading branch information
@tobiasrausch
tobiasrausch committed Feb 5, 2014
1 parent 5f9c8c1 commit 7f5ab1d
Show file tree
Hide file tree
Showing 4 changed files with 40 additions and 21 deletions.
6 changes: 3 additions & 3 deletions src/delly.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1821,7 +1821,7 @@ inline int run(Config const& c, TSVType svType) {
for(; ((vecNext != vecEnd) && (abs(_minCoord(vecNext->Position, vecNext->MatePosition, svType) + vecNext->Length - minCoord) <= overallMaxISize)) ; ++vecNext) {
// Check that mate chr agree (only for translocations)
if (vecBeg->MateRefID!=vecNext->MateRefID) continue;

// Check combinability of pairs
if (_pairsDisagree(minCoord, maxCoord, vecBeg->Length, vecBeg->maxNormalISize, _minCoord(vecNext->Position, vecNext->MatePosition, svType), _maxCoord(vecNext->Position, vecNext->MatePosition, svType), vecNext->Length, vecNext->maxNormalISize, _getSpanOrientation(*vecBeg, vecBeg->libOrient, svType), _getSpanOrientation(*vecNext, vecNext->libOrient, svType), svType)) continue;

Expand Down Expand Up @@ -1905,7 +1905,7 @@ inline int run(Config const& c, TSVType svType) {
int32_t clusterMateRefID=g[itWEdge->source]->MateRefID;
_initClique(g[itWEdge->source], svStart, svEnd, wiggle, svType);
int connectionType = _getSpanOrientation(*g[itWEdge->source], g[itWEdge->source]->libOrient, svType);
if (svStart >= svEnd) continue;
if ((clusterRefID==clusterMateRefID) && (svStart >= svEnd)) continue;
clique.insert(itWEdge->source);

// Grow the clique from the seeding edge
Expand Down Expand Up @@ -1949,7 +1949,7 @@ inline int run(Config const& c, TSVType svType) {
for(;itC!=clique.end();++itC) mapQV.push_back(g[*itC]->MapQuality);
std::sort(mapQV.begin(), mapQV.end());
svRec.peMapQuality = mapQV[mapQV.size()/2];
if (svRec.svStartEnd > svRec.svEndBeg) {
if ((refIndex==clusterMateRefID) && (svRec.svStartEnd > svRec.svEndBeg)) {
unsigned int midPointDel = ((svRec.svEnd - svRec.svStart) / 2) + svRec.svStart;
svRec.svStartEnd = midPointDel -1;
svRec.svEndBeg = midPointDel;
Expand Down
51 changes: 35 additions & 16 deletions src/spanning.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,27 @@ namespace torali {
};


template<typename TPos, typename TCT>
struct SVSizeType {
TPos start;
TPos end;
TCT ct;

SVSizeType() {}

SVSizeType(TPos const s, TPos const e, TCT const c) : start(s), end(e), ct(c) {}
};

template<typename TSVSize>
struct SortSVSizes : public std::binary_function<TSVSize, TSVSize, bool>
{
inline bool operator()(TSVSize const& s1, TSVSize const& s2) {
return (s1.ct < s2.ct) || ((s1.ct==s2.ct) && (s1.start < s2.start)) || ((s1.ct==s2.ct) && (s1.start == s2.start) && (s1.end < s2.end));
}
};



template<typename TPos, typename TQual, typename TArrayType>
inline void
_addReadAndBpCounts(std::vector<HitInterval<TPos, TQual> > const& hit_vector, TArrayType* bp_count)
Expand Down Expand Up @@ -214,15 +235,16 @@ namespace torali {
++show_progress;

// Collect all SV sizes on this chromosome
typedef std::vector<std::pair<int, int> > TSVSizes;
typedef SVSizeType<int, int> TSVSizeType;
typedef std::vector<TSVSizeType> TSVSizes;
TSVSizes svSizes;
typename TSVs::const_iterator itSV = svs.begin();
typename TSVs::const_iterator itSVEnd = svs.end();
for(;itSV!=itSVEnd;++itSV) {
if (itSV->chr == refIndex) svSizes.push_back(std::make_pair(itSV->svStart, itSV->svEnd));
else if (itSV->chr2 == refIndex) svSizes.push_back(std::make_pair(itSV->svEnd, itSV->svStart));
if (itSV->chr == refIndex) svSizes.push_back(TSVSizeType(itSV->svStart, itSV->svEnd, itSV->ct));
else if (itSV->chr2 == refIndex) svSizes.push_back(TSVSizeType(itSV->svEnd, itSV->svStart, itSV->ct));
}
std::sort(svSizes.begin(), svSizes.end());
std::sort(svSizes.begin(), svSizes.end(), SortSVSizes<TSVSizeType>());
if (svSizes.empty()) continue;

// Iterate all samples
Expand Down Expand Up @@ -264,6 +286,7 @@ namespace torali {
if (libIt->second.median == 0) continue; // Single-end library
int outerISize = std::abs(al.Position - al.MatePosition) + al.Length;

// Abnormal paired-end
if ((getStrandIndependentOrientation(al) != libIt->second.defaultOrient) || (outerISize > libIt->second.maxNormalISize) || (al.RefID!=al.MateRefID)) {
if (_acceptedInsertSize(libIt->second.maxNormalISize, libIt->second.median, abs(al.InsertSize), svType)) continue; // Normal paired-end (for deletions only)
if (_acceptedOrientation(libIt->second.defaultOrient, getStrandIndependentOrientation(al), svType)) continue; // Orientation disagrees with SV type
Expand All @@ -272,19 +295,19 @@ namespace torali {
int32_t const minPos = _minCoord(al.Position, al.MatePosition, svType);
int32_t const maxPos = _maxCoord(al.Position, al.MatePosition, svType);

typename TSVSizes::const_iterator itSize = std::lower_bound(svSizes.begin(), svSizes.end(), std::make_pair(minPos, maxPos));
typename TSVSizes::const_iterator itSize = std::lower_bound(svSizes.begin(), svSizes.end(), TSVSizeType(minPos, maxPos, _getSpanOrientation(al, libIt->second.defaultOrient, svType)), SortSVSizes<TSVSizeType>());
bool validSize = false;
if (itSize != svSizes.end()) {
++itSize;
if (itSize != svSizes.end()) validSize = (!_pairsDisagree(minPos, maxPos, al.Length, libIt->second.median, itSize->first, itSize->second, al.Length, libIt->second.median, _getSpanOrientation(al, libIt->second.defaultOrient, svType), _getSpanOrientation(al, libIt->second.defaultOrient, svType), svType));
if (itSize != svSizes.end()) validSize = (!_pairsDisagree(minPos, maxPos, al.Length, libIt->second.maxNormalISize, itSize->start, itSize->end, al.Length, libIt->second.maxNormalISize, _getSpanOrientation(al, libIt->second.defaultOrient, svType), itSize->ct, svType));
--itSize;
}
if ((!validSize) && (itSize != svSizes.end())) {
validSize = (!_pairsDisagree(minPos, maxPos, al.Length, libIt->second.median, itSize->first, itSize->second, al.Length, libIt->second.median, _getSpanOrientation(al, libIt->second.defaultOrient, svType), _getSpanOrientation(al, libIt->second.defaultOrient, svType), svType));
validSize = (!_pairsDisagree(minPos, maxPos, al.Length, libIt->second.maxNormalISize, itSize->start, itSize->end, al.Length, libIt->second.maxNormalISize, _getSpanOrientation(al, libIt->second.defaultOrient, svType), itSize->ct, svType));
}
if ((!validSize) && (!svSizes.empty()) && (itSize != svSizes.begin())) {
--itSize;
validSize = (!_pairsDisagree(itSize->first, itSize->second, al.Length, libIt->second.median, minPos, maxPos, al.Length, libIt->second.median, _getSpanOrientation(al, libIt->second.defaultOrient, svType), _getSpanOrientation(al, libIt->second.defaultOrient, svType), svType));
validSize = (!_pairsDisagree(itSize->start, itSize->end, al.Length, libIt->second.maxNormalISize, minPos, maxPos, al.Length, libIt->second.maxNormalISize, itSize->ct, _getSpanOrientation(al, libIt->second.defaultOrient, svType), svType));
}
if (!validSize) continue;
}
Expand Down Expand Up @@ -316,18 +339,14 @@ namespace torali {
} else if ((getStrandIndependentOrientation(al) != libIt->second.defaultOrient) || (outerISize > libIt->second.maxNormalISize) || (al.RefID!=al.MateRefID)) {
// Missing spanning coverage
if (_mateIsUpstream(libIt->second.defaultOrient, (al.AlignmentFlag & 0x0040), (al.AlignmentFlag & 0x0010)))
missingSpan.push_back(THitInterval(al.Position, al.Position + libIt->second.median, pairQuality));
//missingSpan.push_back(THitInterval(al.Position, al.Position + libIt->second.maxNormalISize, pairQuality));
missingSpan.push_back(THitInterval(al.Position, al.Position + libIt->second.maxNormalISize, pairQuality));
else
missingSpan.push_back(THitInterval(std::max(0, al.Position + al.Length - libIt->second.median), al.Position + al.Length, pairQuality));
//missingSpan.push_back(THitInterval(std::max(0, al.Position + al.Length - libIt->second.maxNormalISize), al.Position + al.Length, pairQuality));
missingSpan.push_back(THitInterval(std::max(0, al.Position + al.Length - libIt->second.maxNormalISize), al.Position + al.Length, pairQuality));
if (al.RefID==al.MateRefID) {
if (_mateIsUpstream(libIt->second.defaultOrient, !(al.AlignmentFlag & 0x0040), (al.AlignmentFlag & 0x0020)))
missingSpan.push_back(THitInterval(al.MatePosition, al.MatePosition + libIt->second.median, pairQuality));
//missingSpan.push_back(THitInterval(al.MatePosition, al.MatePosition + libIt->second.maxNormalISize, pairQuality));
missingSpan.push_back(THitInterval(al.MatePosition, al.MatePosition + libIt->second.maxNormalISize, pairQuality));
else
missingSpan.push_back(THitInterval(std::max(0, al.MatePosition + al.Length - libIt->second.median), al.MatePosition + al.Length, pairQuality));
//missingSpan.push_back(THitInterval(std::max(0, al.MatePosition + al.Length - libIt->second.maxNormalISize), al.MatePosition + al.Length, pairQuality));
missingSpan.push_back(THitInterval(std::max(0, al.MatePosition + al.Length - libIt->second.maxNormalISize), al.MatePosition + al.Length, pairQuality));
}
}
++libIt->second.unique_pairs;
Expand Down
2 changes: 1 addition & 1 deletion src/tags.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -409,7 +409,7 @@ namespace torali {
_pairsDisagree(TSize const pair1Min, TSize const pair1Max, TSize const pair1ReadLength, TISize const pair1maxNormalISize, TSize const pair2Min, TSize const pair2Max, TSize const pair2ReadLength, TISize const pair2maxNormalISize, int const ct1, int const ct2, SVType<TranslocationTag>) {
// Do both pairs support the same translocation type
if (ct1 != ct2) return true;

// Check read offsets
if (ct1%2==0) {
if ((pair2Min + pair2ReadLength - pair1Min) > pair1maxNormalISize) return true;
Expand Down
2 changes: 1 addition & 1 deletion src/version.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ namespace torali
std::cout << "certain conditions (GPL); for license details use '-l'." << std::endl;
std::cout << "This program comes with ABSOLUTELY NO WARRANTY; for details use '-w'." << std::endl;
std::cout << std::endl;
std::cout << title << " (Version: 0.3.1)" << std::endl;
std::cout << title << " (Version: 0.3.2)" << std::endl;
std::cout << "Contact: Tobias Rausch ([email protected])" << std::endl;
std::cout << "**********************************************************************" << std::endl;
std::cout << std::endl;
Expand Down

0 comments on commit 7f5ab1d

Please sign in to comment.