Fix a bug with SV calling from clipped SFS.

Author: Parsoa

clipper.cpp

@@ -98,41 +98,34 @@ vector<Clip> Clipper::filter_tooclose_clips(const vector<Clip> &clips, interval_
     return fclips;
 }
 
-int binary_search_left(const vector<Clip>& clips, int begin, int end, const Clip& query) {
+// find smallest right that is larger than query
+int binary_search(const vector<Clip>& clips, int begin, int end, const Clip& query) {
+    //for (int i = 0; i < clips.size(); i++) {
+    //    if (query.p < clips[i].p) {
+    //        return i ;
+    //    }
+    //}
+    //return -1 ;
     if (begin > end || begin >= clips.size()) {
         return -1 ;
     }
     int m = (begin + end) / 2 ;
     if (clips[m].p == query.p) {
-        if (m != 0) {
-            return m - 1 ;
-        } else {
-            return -1 ;
-        }
-    } else if (clips[m].p < query.p) {
-        return binary_search_left(clips, m, end, query) ;
-    } else {
-        return binary_search_left(clips, begin, m - 1, query) ;
-    }
-}
-
-int binary_search_right(const vector<Clip>& clips, int begin, int end, const Clip& query) {
-    if (begin > end || begin >= clips.size()) {
-        return -1 ;
-    }
-    int m = (begin + end) / 2 ;
-    if (clips[m].p == query.p) {
-        if (m != 0) {
+        if (m + 1 < clips.size()) {
             return m + 1 ;
         } else {
-            return -1 ;
+            return m ;
         }
     } else if (clips[m].p > query.p) {
-        return binary_search_right(clips, begin, m, query) ;
+        if (m > 0 && clips[m - 1].p < query.p) {
+            return m ;
+        }
+        return binary_search(clips, begin, m - 1, query) ;
     } else {
-        return binary_search_right(clips, m + 1, end, query) ;
+        return binary_search(clips, m + 1, end, query) ;
     }
 }
+
 void Clipper::call(int threads, interval_tree_t<int>& vartree) {
     lprint({"Predicting SVS from", to_string(clips.size()), "clipped SFS on", to_string(threads), "threads.."}) ;
     vector<Clip> rclips;
@@ -144,9 +137,9 @@ void Clipper::call(int threads, interval_tree_t<int>& vartree) {
             rclips.push_back(clip);
         }
     }
-    lprint({"Preprocessing clipped SFS.."}) ;
     lprint({to_string(lclips.size()), "left clips."}) ;
     lprint({to_string(rclips.size()), "right clips."}) ;
+    lprint({"Preprocessing clipped SFS.."}) ;
     #pragma omp parallel for num_threads(2) schedule(static,1)
     for (int i = 0; i < 2; i++) {
         if (i == 0) {
@@ -178,7 +171,7 @@ void Clipper::call(int threads, interval_tree_t<int>& vartree) {
         int t = omp_get_thread_num() ;
         string chrom = lc.chrom ;
         // we get the closest right clip
-        int r = binary_search_right(rclips, 0, rclips.size() - 1, lc) ;
+        int r = binary_search(rclips, 0, rclips.size() - 1, lc) ;
         if (r == -1) {
             continue ;
         }
@@ -202,7 +195,7 @@ void Clipper::call(int threads, interval_tree_t<int>& vartree) {
         int t = omp_get_thread_num() ;
         string chrom = rc.chrom ;
         // we get the closest right clip
-        int l = binary_search_left(lclips, 0, lclips.size() - 1, rc) ;
+        int l = binary_search(lclips, 0, lclips.size() - 1, rc) ;
         if (l == -1) {
             continue ;
         }

extender.cpp

@@ -28,7 +28,7 @@ void Extender::run(int _threads) {
         }
         clips.insert(clips.begin(), _p_clips[i].begin(), _p_clips[i].end()) ;
     }
-    lprint({"Extension complete."}) ;
+    lprint({"Extension complete.", to_string(clips.size()), " clipped SFS."}) ;
     // build a separate interval tree for each chromosome
     cluster_no_interval_tree() ;
     // put all clusters in a single vector
@@ -60,6 +60,7 @@ void Extender::run(int _threads) {
     lprint({"Extracted", to_string(svs.size()), "SVs."});
     sam_close(bam_file) ;
     filter_sv_chains() ;
+    //lprint({"Error stats: ", to_string(skip_1), " can't be mapped, ", to_string(skip_2), " can't be extended ", to_string(skip_3), " anomalies.", to_string(skip_4), " clipped."}) ;
 }
 
 pair<int, int> Extender::get_unique_kmers(const vector<pair<int, int>> &alpairs, const uint k, const bool from_end, string chrom) {
@@ -207,9 +208,6 @@ bool Extender::load_batch_bam(int threads, int batch_size, int p) {
         if (aln->core.flag & BAM_FUNMAP || aln->core.flag & BAM_FSUPPLEMENTARY || aln->core.flag & BAM_FSECONDARY) {
             continue ;
         }
-        if (aln->core.qual != 0) {
-            continue ;
-        }
         char *qname = bam_get_qname(aln);
         if (SFSs->find(qname) == SFSs->end()) {
             continue ;
@@ -298,6 +296,7 @@ void Extender::extend_alignment(bam1_t* aln, int index) {
             uint l = bam_cigar_oplen(*(cigar + 0));
             if (op == BAM_CSOFT_CLIP) {
                 lclip = make_pair(aln->core.pos, l);
+                skip_4++ ;
             }
             // we skip this SFS
             continue ;
@@ -306,6 +305,7 @@ void Extender::extend_alignment(bam1_t* aln, int index) {
             uint l = bam_cigar_oplen(*(cigar + aln->core.n_cigar - 1));
             if (op == BAM_CSOFT_CLIP) {
                 rclip = make_pair(bam_endpos(aln), l);
+                skip_4++ ;
             }
             // we skip this SFS
             continue ;
@@ -382,12 +382,12 @@ void Extender::extend_alignment(bam1_t* aln, int index) {
         } else {
             _p_extended_sfs[index].push_back(ExtSFS(string(chrom), string(qname), prekmer.second, postkmer.second + kmer_size));
         }
-        if (lclip.second > 0) {
-            _p_clips[index].push_back(Clip(qname, chrom, lclip.first, lclip.second, true)) ;
-        }
-        if (rclip.second > 0) {
-            _p_clips[index].push_back(Clip(qname, chrom, rclip.first, rclip.second, false)) ;
-        }
+    }
+    if (lclip.second > 0) {
+        _p_clips[index].push_back(Clip(qname, chrom, lclip.first, lclip.second, true)) ;
+    }
+    if (rclip.second > 0) {
+        _p_clips[index].push_back(Clip(qname, chrom, rclip.first, rclip.second, false)) ;
     }
 }
 
@@ -607,7 +607,9 @@ void Extender::extract_sfs_sequences() {
         cerr << endl ;
     }
     for (int i = 0; i < threads; i++) {
-        free(seq[i]) ;
+        if (len[i] > 0) {
+            free(seq[i]) ;
+        }
         bam_destroy1(_p_aln[i]) ;
         sam_close(_p_bam_file[i]) ;
     }
@@ -784,6 +786,9 @@ void Extender::call() {
 
 void Extender::filter_sv_chains() {
     std::sort(svs.begin(), svs.end()) ;
+    if (svs.size() < 2) {
+        return ;
+    }
     lprint({to_string(svs.size()), "SVs before chain filtering."}) ;
     vector<SV> _svs ;
     auto& prev = svs[0] ;
@@ -796,7 +801,7 @@ void Extender::filter_sv_chains() {
         }
         auto& sv = svs[i] ;
         if (sv.chrom == prev.chrom && sv.s - prev.e < 2 * sv.l && prev.type == sv.type) {
-            cout << sv.chrom << " " << sv.s << " " << sv.type << " ---- " << prev.s << endl ; 
+            //cout << sv.chrom << " " << sv.s << " " << sv.type << " ---- " << prev.s << endl ; 
             // check for sequence similarity
             double w_r = max((double)sv.w, (double)prev.w) / min((double)sv.w, (double)prev.w) ;
             double l_r = min((double)sv.l, (double)prev.l) / max((double)sv.l, (double)prev.l) ;
@@ -807,7 +812,7 @@ void Extender::filter_sv_chains() {
                 } else {
                     d = rapidfuzz::fuzz::ratio(sv.altall, prev.altall) ;
                 }
-                cout << w_r << " " << l_r << " " << d << endl ;
+                //cout << w_r << " " << l_r << " " << d << endl ;
                 if (d > 70) {
                     if (sv.w > prev.w) {
                         _svs.push_back(sv) ;

extender.hpp

@@ -35,6 +35,7 @@ class Extender {
     uint skip_1 = 0 ;      // SFS skipped since no first/last base can be placed from read alignment (should be rare)
     uint skip_2 = 0 ;      // SFS skipped since it couldn't be extended
     uint skip_3 = 0 ;      // SFS skipped since reads starts/ends inside a cluster
+    uint skip_4 = 0 ;      // SFS skipped since reads starts/ends inside a cluster
     uint small_cl = 0 ;    // number of cluster (before clustering) with low support
     uint extcl = 0 ;       // number of extended clusters (after clustering)
     uint small_extcl = 0 ; // number of extended clusters (after clustering) with low support

ping_pong.cpp

@@ -57,9 +57,6 @@ bool PingPong::backward_search(rld_t *index, const uint8_t *P, int p2) {
 // However non-reconstructed reads are going to produce loads of crappy SFS, unless we filter them
 void PingPong::ping_pong_search(rld_t *index, uint8_t* P, int l, std::vector<sfs_type_t>& solutions, bool isreconstructed, bam1_t* aln) {
     //DEBUG(cerr << "Read Length: " << l << endl ;)
-    if (!isreconstructed) {
-        return ;
-    }
     // this is not needed anymore
     //int current_interval = -1 ;
     //vector<pair<int, int>> m_intervals ;
@@ -268,12 +265,16 @@ batch_type_t PingPong::process_batch(rld_t* index, int p, int i) {
     } else {
         for (int j = 0; j < read_seqs[p][i].size(); j++) {
             char *qname = bam_get_qname(bam_entries[p][i][j]) ;
+            bool isreconstructed = reconstructed_reads.find(qname) != reconstructed_reads.end() ;
             if (config->putative) {
                 if (ignored_reads.find(qname) != ignored_reads.end()) {
                     continue ;
                 }
+                // was not ignored, so either it's reconstructed or not:
+                if (!isreconstructed) {
+                    continue ;
+                }
             }
-            bool isreconstructed = reconstructed_reads.find(qname) != reconstructed_reads.end() ;
             //cout << qname << " " << isreconstructed << endl ;
             ping_pong_search(index, read_seqs[p][i][j], read_seq_lengths[p][i][j], solutions[qname], isreconstructed, bam_entries[p][i][j]) ;
         }

reconstructor.cpp

@@ -172,7 +172,7 @@ void Reconstructor::reconstruct_read(bam1_t* alignment, char* read_seq, string c
             new_cigar.push_back(cigar_offsets[m]) ;
         } else {
             cout << "Illegal Cigar OP" << endl ;
-            break ;            
+            break ;
             //if (cigar_offsets[m].second == BAM_CPAD || cigar_offsets[m].second == BAM_CHARD_CLIP || cigar_offsets[m].second == BAM_CBACK) {
         }
         m += 1 ;