multiply file support

Author: panguangze

phaser.cpp

@@ -21,9 +21,6 @@ Phaser::Phaser(const std::string &fnvcf, const std::string &fnout, std::vector<s
             frbed = new BEDReader(fnbed.data());
         }
     }
-    for (const auto& item : fnfrags) {
-        frfrags.push_back(new FragmentReader(item.data()));
-    }
     frvcf = new VCFReader(fnvcf.data());
     fwvcf = new VCFWriter(frvcf->header, fnout.data());
 //    frfrag = new FragmentReader(fnfrag.data());

spectral.cpp

@@ -128,9 +128,6 @@ void Spectral::reset()
     for (auto item : frs) {
         item->set_window_info(start_variant_idx, end_variant_idx_overlap, end_variant_idx_intended);
     }
-    for (auto item : frs) {
-        item->set_window_info(start_variant_idx, end_variant_idx_overlap, end_variant_idx_intended);
-    }
 //    fr->set_window_info(start_variant_idx, end_variant_idx_overlap, end_variant_idx_intended);
     this->variant_graph.reset(variant_count);
     this->raw_graph = new double[n * n];
@@ -142,19 +139,27 @@ void Spectral::reset()
         this->raw_count[i] = 0;
     }
     this->frag_buffer.clear();
+//    read graph
+    ViewMap weighted_graph(raw_graph, n, n);
+    CViewMap count_graph(raw_count, n, n);
     for (int i = 0; i < OPERATIONS.size(); i++) {
         auto item = OPERATIONS[i];
         if (item == MODE_10X)
-            read_fragment_10x(i);
+            read_fragment_10x(i,weighted_graph,count_graph);
         else if (item == MODE_HIC)
-            read_fragment_hic(i);
+            read_fragment_hic(i,weighted_graph,count_graph);
         else if (item == MODE_PE)
-            read_fragment(i);
+            read_fragment(i,weighted_graph,count_graph);
         else if (item == MODE_NANOPORE)
-            read_fragment_nanopore(i);
+            read_fragment_nanopore(i,weighted_graph,count_graph);
         else if (item == MODE_PACBIO)
-            read_fragment_pacbio(i);
+            read_fragment_pacbio(i,weighted_graph,count_graph);
+    }
+//    cal prob graph
+    if (HAS_TENX){
+        add_snp_edge_barcode(weighted_graph, count_graph);
     }
+    cal_prob_matrix(weighted_graph, count_graph, nullptr, nullptr, nullptr);
 }
 
 GMatrix Spectral::slice_submat(std::set<uint> &variants_mat, GMatrix &adj_mat)
@@ -212,23 +217,23 @@ CMatrix Spectral::slice_submat(std::set<uint> &variants_mat, bool t)
 }
 
 // read fragment matrix
-void Spectral::read_fragment(int frIdx)
+void Spectral::read_fragment(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph)
 {
     auto fr = frs[frIdx];
 //    this->frag_buffer.clear();
     Fragment fragment;
-    ViewMap weighed_graph(raw_graph, n, n);
-    CViewMap count_graph(raw_count, n, n);
+//    ViewMap weighed_graph(raw_graph, n, n);
+//    CViewMap count_graph(raw_count, n, n);
     while (fr->get_next_pe(fragment))
     {
-        add_snp_edge(fragment, weighed_graph, count_graph);
+        add_snp_edge(fragment, weighted_graph, count_graph);
         this->frag_buffer.push_back(fragment);
         fragment.reset();
     }
-    cal_prob_matrix(weighed_graph, count_graph, nullptr, nullptr, nullptr);
+//    cal_prob_matrix(weighted_graph, count_graph, nullptr, nullptr, nullptr);
 }
 
-void Spectral::read_fragment_10x(int frIdx)
+void Spectral::read_fragment_10x(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph)
 {
     auto fr = frs[frIdx];
     this->region_frag_stats->clear();
@@ -238,8 +243,8 @@ void Spectral::read_fragment_10x(int frIdx)
     frbed->read_region_frag(this->chromo_phaser->chr_name.c_str(), start_pos, end_pos, this->region_frag_stats);
     //this->barcode_linker->regionFragStats = region_frag_stats;
     Fragment fragment;
-    ViewMap weighed_graph(raw_graph, n, n);
-    CViewMap count_graph(raw_count, n, n);
+//    ViewMap weighed_graph(raw_graph, n, n);
+//    CViewMap count_graph(raw_count, n, n);
 
     while (fr->get_next_tenx(fragment))
     {
@@ -250,61 +255,61 @@ void Spectral::read_fragment_10x(int frIdx)
             //add_barcode_info(fragment, barcode_linker);
         fragment.reset();
     }
-    add_snp_edge_barcode(weighed_graph, count_graph);
-    cal_prob_matrix(weighed_graph, count_graph, nullptr, nullptr, nullptr);
+//    add_snp_edge_barcode(weighted_graph, count_graph);
+//    cal_prob_matrix(weighted_graph, count_graph, nullptr, nullptr, nullptr);
 }
 
 //TODO: only unused HiC linker should be stored
-void Spectral::read_fragment_hic(int frIdx)
+void Spectral::read_fragment_hic(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph)
 {
     auto fr = frs[frIdx];
     Fragment fragment;
 //    this->frag_buffer.clear();
-    CViewMap count_graph(raw_count, n, n);
-    ViewMap weighed_graph(raw_graph, n, n);
+//    CViewMap count_graph(raw_count, n, n);
+//    ViewMap weighed_graph(raw_graph, n, n);
     while (fr->get_next_hic(fragment))
     {
-        add_snp_edge(fragment, weighed_graph, count_graph);
+        add_snp_edge(fragment, weighted_graph, count_graph);
         this->frag_buffer.push_back(fragment);
         if ( fragment.snps[0].first >= phasing_window->prev_window_start)
             if (fragment.insertion_size >= 5000 && fragment.insertion_size <= 40000000)
                 this->hic_linker_container.add_HiC_info(fragment);
         fragment.reset();
     }
-    cal_prob_matrix(weighed_graph, count_graph, nullptr, nullptr, nullptr);
+//    cal_prob_matrix(weighted_graph, count_graph, nullptr, nullptr, nullptr);
 }
 
-void Spectral::read_fragment_nanopore(int frIdx)
+void Spectral::read_fragment_nanopore(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph)
 {
     auto fr = frs[frIdx];
 //    this->frag_buffer.clear();
     Fragment fragment;
-    ViewMap weighed_graph(raw_graph, n, n);
-    CViewMap count_graph(raw_count, n, n);
+//    ViewMap weighed_graph(raw_graph, n, n);
+//    CViewMap count_graph(raw_count, n, n);
     while (fr->get_next_nanopore(fragment))
     {
-        add_snp_edge(fragment, weighed_graph, count_graph);
+        add_snp_edge(fragment, weighted_graph, count_graph);
         this->frag_buffer.push_back(fragment);
         fragment.reset();
     }
     //this->q_aver = q_sum / this->frag_buffer.size() / 6;
-    cal_prob_matrix(weighed_graph, count_graph, nullptr, nullptr, nullptr);
+//    cal_prob_matrix(weighted_graph, count_graph, nullptr, nullptr, nullptr);
 }
 
-void Spectral::read_fragment_pacbio(int frIdx)
+void Spectral::read_fragment_pacbio(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph)
 {
     auto fr = frs[frIdx];
 //    this->frag_buffer.clear();
     Fragment fragment;
-    ViewMap weighed_graph(raw_graph, n, n);
-    CViewMap count_graph(raw_count, n, n);
+//    ViewMap weighed_graph(raw_graph, n, n);
+//    CViewMap count_graph(raw_count, n, n);
     while (fr->get_next_pacbio(fragment))
     {
-        add_snp_edge(fragment, weighed_graph, count_graph);
+        add_snp_edge(fragment, weighted_graph, count_graph);
         this->frag_buffer.push_back(fragment);
         fragment.reset();
     }
-    cal_prob_matrix(weighed_graph, count_graph, nullptr, nullptr, nullptr);
+//    cal_prob_matrix(weighted_graph, count_graph, nullptr, nullptr, nullptr);
 }
 
 // aid function
@@ -625,6 +630,8 @@ void Spectral::solver()
     //for each block (variant) in phasing window
     for (auto i : phasing_window->current_window_idxes)
     {
+        if (i == 15)
+            int tmp = 9;
         mat_idx = phasing_window->var_idx2mat_idx(i);
         if (met_idx.find(mat_idx) == met_idx.end())
             met_idx.insert(mat_idx);
@@ -667,11 +674,11 @@ void Spectral::solver()
         }
     }
 //    TODO, if we need this
-    if (false)
-    {
-        for (auto start_idx: this->phased_block_starts)
-            barcode_aware_filter(start_idx.first);
-    }
+//    if (HAS_TENX)
+//    {
+//        for (auto start_idx: this->phased_block_starts)
+//            barcode_aware_filter(start_idx.first);
+//    }
 }
 
 void Spectral::add_snp_edge_barcode_subroutine(ViewMap &sub_weighted_graph, CViewMap &sub_count_graph, VariantGraph &sub_variant_graph, std::map<uint, int> &subroutine_map, std::map<uint, uint> & subroutine_blk_start)
@@ -880,7 +887,11 @@ void Spectral::find_connected_component_dfs(const Eigen::MatrixBase<Derived> &ad
         std::set<uint> &nxt_vars = sub_variant_graph.graph[var_idx];
 
         uint idx = phasing_window->mat_idx2var_idx(subroutine_blk_start[var_idx]);
+        //        Fixme, block_to_merge equals to starting block, this is a bug
+//        if(idx == starting_block->block_id) return;
         ptr_PhasedBlock block_to_merge = phasing_window->blocks[idx];
+        if (block_to_merge->block_id == starting_block->block_id)
+            return;
         if (adj_mat(2 * prev_var_idx, 2 * var_idx) > 0)
         {
             if (block_to_merge->size() == 1)
@@ -953,7 +964,11 @@ void Spectral::find_connected_component_dfs(const Eigen::MatrixBase<Derived> &ad
         std::set<uint> &nxt_vars = this->variant_graph.graph[var_idx];
 
         uint idx = phasing_window->mat_idx2var_idx(var_idx);
+        //        Fixme, block_to_merge equals to starting block, this is a bug
+//        if(idx == starting_block->block_id) return;
         ptr_PhasedBlock block_to_merge = phasing_window->blocks[idx];
+        if (block_to_merge->block_id == starting_block->block_id)
+            return;
         if (adj_mat(2 * prev_var_idx, 2 * var_idx) > 0)
         {
             if (block_to_merge->size() == 1)

spectral.h

@@ -82,7 +82,7 @@ class Spectral
 
     std::map<uint, uint> phased_block_starts;
     std::set<uint> block_tobe_split;
-    FragmentReader *fr;             //Fragment matrix reader
+//    FragmentReader *fr;             //Fragment matrix reader
     std::vector<FragmentReader* > frs;
     BEDReader *frbed;
     HiCLinkerContainer hic_linker_container;
@@ -130,11 +130,11 @@ class Spectral
     CMatrix slice_submat(std::set<uint> &variants_mat, bool t);
     CMatrix slice_submat(std::set<uint> &variants_mat, bool t, CMatrix &adj_mat);
     void filter_inconsistency();
-    void read_fragment_10x(int frIdx);
-    void read_fragment(int frIdx);
-    void read_fragment_hic(int frIdx);
-    void read_fragment_pacbio(int frIdx);
-    void read_fragment_nanopore(int frIdx);
+    void read_fragment_10x(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph);
+    void read_fragment(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph);
+    void read_fragment_hic(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph);
+    void read_fragment_pacbio(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph);
+    void read_fragment_nanopore(int frIdx, ViewMap &weighted_graph, CViewMap &count_graph);
     void poss_phase_error_correction(uint block_start_idx);
     void fragment_supported_flipping_score(ptr_PhasedBlock &phased_block, Fragment & fragment, int *supporting_reads_count, double *supporting_weight_count, std::map<uint, std::set<uint>> &connection_map);
     int locate_block_valid_start(const Eigen::VectorXd &vec);