@@ -126,15 +126,7 @@ void Reconstructor::reconstruct_read(bam1_t* alignment, char* read_seq, string c
126126 memcpy (new_qual + n, qual + match_offset + ins_offset + soft_clip_offset, cigar_offsets[m].first ) ;
127127 n += cigar_offsets[m].first ;
128128 for (int j = 0 ; j < cigar_offsets[m].first ; j++) {
129- <<<<<<< HEAD
130- new_seq[n] = chromosome_seqs[chrom][ref_offset + j] ;
131- new_qual[n] = qual[soft_clip_offset + match_offset + ins_offset + j] ;
132- num_match += 1 ? chromosome_seqs[chrom][ref_offset + j] == read_seq[match_offset + ins_offset + soft_clip_offset + j] : 0 ;
133- num_mismatch += 1 ? chromosome_seqs[chrom][ref_offset + j] != read_seq[match_offset + ins_offset + soft_clip_offset + j] : 0 ;
134- n++ ;
135- =======
136129 num_mismatch += 1 ? ref_seq[ref_offset + j] != read_seq[match_offset + ins_offset + soft_clip_offset + j] : 0 ;
137- >>>>>>> poa-genotyper
138130 }
139131 num_match += cigar_offsets[m].first ;
140132 ref_offset += cigar_offsets[m].first ;
@@ -200,10 +192,8 @@ void Reconstructor::reconstruct_read(bam1_t* alignment, char* read_seq, string c
200192 // cout << bam_get_qname(alignment) << endl ;
201193 // only do this on first processing thread
202194 if (omp_get_thread_num () == 2 ) {
203- global_num_bases += num_match + num_mismatch + ins_offset + del_offset ;
204- global_num_match += num_match ;
195+ global_num_bases += num_match ;
205196 global_num_mismatch += num_mismatch ;
206- global_num_indel += ins_offset + del_offset ;
207197 }
208198 // how many errors and SNPs do we expect? 1/1000 each, so say if we see more than twice that then don't correct
209199 if (config->selective ) {
@@ -214,12 +204,9 @@ void Reconstructor::reconstruct_read(bam1_t* alignment, char* read_seq, string c
214204 return ;
215205 }
216206 // if we have so many deletions and insertions, then abort
217- // if ((ins_offset + del_offset) / (num_match + num_mismatch) > 3 * expected_indel_rate) {
218- // if (omp_get_thread_num() == 3) {
219- // num_ignored_reads += 1 ;
220- // }
221- // return ;
222- // }
207+ if (ins_offset + del_offset > 0.7 * strlen (read_seq)) {
208+ return ;
209+ }
223210 }
224211 if (should_ignore) {
225212 // check mismatch rate
@@ -292,16 +279,8 @@ void Reconstructor::run() {
292279 }
293280 }
294281 }
295- <<<<<<< HEAD
296- int p = 0 ;
297- int b = 0 ;
298- int batch_size = (10000 / config->threads ) * config->threads ;
299- lprint ({" Loading first batch.." });
300- for (int i = 0 ; i < 2 ; i++) {
301- =======
302282 for (int i = 0 ; i < modulo; i++) {
303283 bam_entries.push_back (vector<vector<bam1_t *>>(config->threads )) ;
304- >>>>>>> poa-genotyper
305284 for (int j = 0 ; j < config->threads ; j++) {
306285 for (int k = 0 ; k < batch_size / config->threads ; k++) {
307286 bam_entries[i][j].push_back (bam_init1 ()) ;
@@ -386,18 +365,11 @@ void Reconstructor::run() {
386365 if (s - t == 0 ) {
387366 s += 1 ;
388367 }
389- <<<<<<< HEAD
390- cerr << " [I] Processed batch " << std::left << std::setw (10 ) << b << " . Reads so far " << std::right << std::setw (12 ) << u << " . Reads per second: " << u / (s - t) << " . Time: " << std::setw (8 ) << std::fixed << s - t << " \n " ;
391- cerr << " [I] Process bases: " << std::left << std::setw (16 ) << uint64_t (global_num_bases) << " , num mismatch: " << std::setw (16 ) << uint64_t (global_num_mismatch) << " , mismatch rate: " << global_num_mismatch / global_num_bases << " , indel rate: " << global_num_indel / (global_num_match + global_num_mismatch) << " , ignored reads: " << num_ignored_reads << endl ;
392- expected_mismatch_rate = global_num_mismatch / global_num_bases ;
393- expected_indel_rate = global_num_indel / global_num_bases ;
394- =======
395368 cerr << " [I] Processed batch " << b << " . Reads so far " << reads_processed << " . Reads per second: " << reads_processed / (s - t) << " . Time: " << s - t << " \n " ;
396369 cerr << " [I] Processed bases: " << uint64_t (global_num_bases) << " , num mismatch: " << uint64_t (global_num_mismatch) << " , mismatch rate: " << global_num_mismatch / global_num_bases << " , ignored reads: " << num_ignored_reads << " \n " ;
397370 expected_mismatch_rate = global_num_mismatch / global_num_bases ;
398371 cerr << " \x1b [A" ;
399372 cerr << " \x1b [A" ;
400- >>>>>>> poa-genotyper
401373 }
402374 cerr << endl ;
403375 cerr << endl ;
@@ -490,4 +462,3 @@ bool Reconstructor::load_batch_bam(int threads, int batch_size, int p) {
490462 // lprint({"Loaded", to_string(n), "BAM reads.."});
491463 return n != 0 ? true : false ;
492464}
493-
0 commit comments