Skip to content

Commit 2baa8f5

Browse files
committed
ignore SV type
1 parent 198be12 commit 2baa8f5

File tree

1 file changed

+22
-5
lines changed

1 file changed

+22
-5
lines changed

src/compvcf.h

Lines changed: 22 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,7 @@ namespace sansa
7373
bool filterForPass;
7474
bool checkID;
7575
bool checkCT;
76+
bool checkSVT;
7677
int32_t qualthres;
7778
int32_t bpdiff;
7879
int32_t minsize;
@@ -98,10 +99,12 @@ namespace sansa
9899
typename TCompSVType::const_iterator itsv = std::lower_bound(compsv.begin(), compsv.end(), CompSVRecord(basesv[i].tid, earliestStart));
99100
for(uint32_t j = (itsv - compsv.begin()); j < compsv.size(); ++j) {
100101
if (basesv[i].tid < compsv[j].tid) break; // Sorted by tid
101-
if (c.checkCT) {
102-
if (basesv[i].svt != compsv[j].svt) continue;
103-
} else {
104-
if (_addID(basesv[i].svt) != _addID(compsv[j].svt)) continue; // Compare SV types (BND, INV,...) but not CT
102+
if (c.checkSVT) {
103+
if (c.checkCT) {
104+
if (basesv[i].svt != compsv[j].svt) continue;
105+
} else {
106+
if (_addID(basesv[i].svt) != _addID(compsv[j].svt)) continue; // Compare SV types (BND, INV,...) but not CT
107+
}
105108
}
106109
if (basesv[i].mtid != compsv[j].mtid) continue;
107110
if (std::abs(basesv[i].svStart - compsv[j].svStart) > c.bpdiff) {
@@ -277,6 +280,7 @@ namespace sansa
277280
}
278281

279282
// BNDs
283+
int32_t svStartVal = rec->pos;
280284
if (svtVal == "BND") {
281285
int32_t pos2val = -1;
282286
if ((_isKeyPresent(hdr, "CHR2")) && (_isKeyPresent(hdr, "POS2"))) {
@@ -306,6 +310,14 @@ namespace sansa
306310
if ((!chr2Name.empty()) && (pos2val != -1)) {
307311
svEndVal = pos2val;
308312
svLenVal = 0;
313+
// Some SV callers use BND for intra-chromosomal, make sure POS < END
314+
if (chr2Name == std::string(bcf_hdr_id2name(hdr, rec->rid))) {
315+
if (svStartVal > svEndVal) {
316+
int32_t tmpVal = svStartVal;
317+
svStartVal = svEndVal;
318+
svEndVal = tmpVal;
319+
}
320+
}
309321
} else {
310322
success=false;
311323
std::cerr << "Error: Could not parse BND ALT allele " << rec->d.allele[1] << std::endl;
@@ -365,7 +377,7 @@ namespace sansa
365377
if (c.chrmap.find(chr2Name) == c.chrmap.end()) c.chrmap.insert(std::make_pair(chr2Name, c.chrmap.size()));
366378
sv.mtid = c.chrmap[chr2Name];
367379
}
368-
sv.svStart = rec->pos;
380+
sv.svStart = svStartVal;
369381
sv.svEnd = svEndVal;
370382
// Use some canonical ordering, swap if necessary
371383
if ((svtVal == "BND") && (sv.tid > sv.mtid)) {
@@ -537,6 +549,7 @@ namespace sansa
537549
("sizeratio,s", boost::program_options::value<float>(&c.sizeratio)->default_value(0.5), "min. SV size ratio")
538550
("divergence,d", boost::program_options::value<float>(&c.divergence)->default_value(0.3), "max. SV allele divergence")
539551
("outprefix,o", boost::program_options::value<std::string>(&c.outprefix)->default_value("out"), "output prefix")
552+
("nosvt,t", "Ignore the SV type")
540553
("pass,p", "Filter sites for PASS")
541554
("ignore,i", "Ignore duplicate IDs")
542555
("ct,c", "Require matching CT value in addition to SV type")
@@ -580,6 +593,10 @@ namespace sansa
580593
if (vm.count("ct")) c.checkCT = true;
581594
else c.checkCT = false;
582595

596+
// Ignore SV type
597+
if (vm.count("nosvt")) c.checkSVT = false;
598+
else c.checkSVT = true;
599+
583600
// Check base VCF file
584601
std::set<std::string> baseSamples;
585602
if (vm.count("base")) {

0 commit comments

Comments
 (0)