1
- #include < queue>
1
+ #include < queue>
2
2
#include < vector>
3
3
#include < iostream>
4
4
#include < map>
@@ -16,7 +16,7 @@ using namespace std;
16
16
typedef uint32_t DepthType; // type for depth of coverage, kept it small
17
17
const char ref_char = ' >' ; // reference prefix for "counts" output
18
18
19
- const string VERSION = " %prog 1.1.3 "
19
+ const string VERSION = " %prog 1.1.4 "
20
20
" \n Copyright (C) 2014-2019 Giovanni Birolo and Andrea Telatin\n "
21
21
" https://github.com/telatin/covtobed - License MIT"
22
22
" .\n "
@@ -61,11 +61,17 @@ class Input {
61
61
bool get_next_alignment (BamAlignment & alignment) {
62
62
bool more_alignments, good_alignment;
63
63
do {
64
- debug cerr << " [M] " << alignment.Name << " :" << alignment.Position << " | Is mapped? " << alignment.IsMapped () << endl;
64
+ debug cerr << " [M] Read on ref#" << alignment.RefID << " pos:" << alignment.Position <<
65
+ " \n\t | Is mapped? " << alignment.IsMapped () << " | AlignmentFlag:" << alignment.AlignmentFlag << endl;
65
66
more_alignments = input_bams.GetNextAlignmentCore (alignment);
66
67
if (discard_invalid_alignments) {
67
68
good_alignment = alignment.IsMapped () && alignment.MapQuality >= min_mapq
68
- && !alignment.IsDuplicate () && !alignment.IsFailedQC () && alignment.IsPrimaryAlignment () ;
69
+ && !alignment.IsDuplicate () && !alignment.IsFailedQC () && alignment.IsPrimaryAlignment ();
70
+
71
+ // 1.2.0 ProperPair
72
+ if (alignment.IsPaired () && ! alignment.IsProperPair () ) {
73
+ good_alignment = false ;
74
+ }
69
75
} else {
70
76
good_alignment = alignment.IsMapped () && alignment.MapQuality >= min_mapq;
71
77
}
@@ -210,11 +216,6 @@ int main(int argc, char *argv[]) {
210
216
min_mapq = options.get (" min_mapq" );
211
217
}
212
218
213
- if (physical_coverage and only_valid) {
214
- // https://github.com/telatin/covtobed/issues/11
215
- cerr << " Parameters --physical-coverage and --discard-invalid-alignments are currently mutually exclusive." << endl;
216
- exit (0 );
217
- }
218
219
219
220
try {
220
221
// open input and output
@@ -253,6 +254,7 @@ int main(int argc, char *argv[]) {
253
254
while (more_alignments_for_ref && next_change == alignment.Position ) {
254
255
if (physical_coverage) {
255
256
if (alignment.InsertSize > 0 ) {
257
+ debug cerr << " [phy] pos:" << alignment.Position << " size:" << alignment.InsertSize << endl;
256
258
coverage_ends.push ({alignment.Position + alignment.InsertSize , alignment.IsReverseStrand ()});
257
259
coverage.inc (alignment.IsReverseStrand ());
258
260
}
@@ -272,9 +274,13 @@ int main(int argc, char *argv[]) {
272
274
debug cerr << " [<] Coverage is " << coverage << " from " << next_change << endl;
273
275
last_pos = next_change;
274
276
} while (last_pos != ref.RefLength );
275
- debug cerr << " [_] Completed at " << ref.RefName << " :" << last_pos << endl;
277
+ debug cerr << " [_] Completed at " << ref.RefName << " :" << last_pos << " [coverage: " << coverage_ends. size () << " ] " << endl;
276
278
// reference ended
277
- assert (coverage_ends.empty ());
279
+ if (!coverage_ends.empty ()) {
280
+ cerr << " Coverage is not zero at the end of " << ref.RefName << endl;
281
+ cerr << " Try samtools fixmate on the input file" << endl;
282
+ }
283
+ // 1.2.0 -- removed: assert(coverage_ends.empty());
278
284
}
279
285
// assert(!more_alignments);
280
286
if (more_alignments) {
0 commit comments