enable multiply input

Author: panguangze

Diff for: CMakeLists.txt

@@ -24,8 +24,8 @@ aux_source_directory(./barcode_extractor BARCODE_SRC)
 
 add_executable(SpecHap ${SPEC_SRC})
 add_executable(BarcodeExtract ${BARCODE_SRC})
-target_link_libraries(SpecHap ${HTSlib_LIBRARIES} ${ARPACK_LIBRARIES})
-target_link_libraries(BarcodeExtract ${HTSlib_LIBRARIES})
+target_link_libraries(SpecHap -lhts -larpack)
+target_link_libraries(BarcodeExtract -lhts)
 
 #add_custom_target(extractHAIRS
 #        ALL COMMAND make

Diff for: main.cpp

@@ -2,8 +2,16 @@
 #include "phaser.h"
 #include "results.h"
 #include "type.h"
+#include <string>
+#include <sstream>
 
 
+#define PROTOCOL_HIC "hic"
+#define PROTOCOL_NGS "ngs"
+#define PROTOCOL_TENX "tenx"
+#define PROTOCOL_PACBIO "pacbio"
+#define PROTOCOL_NANOPORE "nanopore"
+#define PROTOCOL_HYBRID "hybrid"
 bool HYBRID = false; 
 bool NEW_FORMAT = false;
 bool KEEP_PS = false;
@@ -14,14 +22,18 @@ int MAX_BARCODE_SPANNING = 60000;
 int MAX_HIC_INSERTION= 40000000;
 int RECURSIVE_LIMIT = 15;
 int OPERATION = MODE_PE;
+bool HAS_HIC;
+bool HAS_TENX;
+std::vector<int> OPERATIONS;
 
 
 enum optionIndex
 {
     UNKNOWN, HELP, VCF, FRAGMENT, OUT, TENX, HIC, _WINDOW_SIZE, COVERAGE, _RECURSIVE_LIMIT, NANOPORE, PACBIO, NOSORT,
-    _MAX_BARCODE_SPANNING_LENGTH, _WINDOW_OVERLAP, STATS, _NEWFORMAT, USESECONDARY, _KEEP_PHASING_INFO, _BASEOFFSET, _HYBRID
+    _MAX_BARCODE_SPANNING_LENGTH, _WINDOW_OVERLAP, STATS, _NEWFORMAT, USESECONDARY, _KEEP_PHASING_INFO, _BASEOFFSET, _HYBRID,PROTOCOLS
 };
 
+
 struct Arg : public option::Arg
 {
     static void printError(const char *msg1, const option::Option &opt, const char *msg2)
@@ -62,7 +74,7 @@ const option::Descriptor usage[] =
                 {UNKNOWN,           0, "",  "",                       Arg::Unknown,  "Usage: "},
                 {HELP,              0, "h", "help",                   Arg::None,     "\t--help\tPrint this message and exit."},
                 {VCF,               0, "v", "vcf",                    Arg::Required, "\t-v <arg>,\t--vcf=<arg>\tSorted Heterozygous VCF File, gzipped, tbi or csi index required."},
-                {FRAGMENT,          0, "f", "frag",                   Arg::Required, "\t-f <arg>,\t--frag=<arg>\tFragment File Generated by ExtractHairs, sort by SNP position is required."},
+                {FRAGMENT,          0, "f", "frag",                   Arg::Required, "\t-f <arg>,\t--frag=<arg>\tFragment File Generated by ExtractHairs, sort by SNP position is required, multiply files with comma split is permitted"},
                 {STATS,             0, "s", "frag_stat",              Arg::Required,  "\t-s <arg>,\t--frag_stat=<arg>\tFragment file statistic, in bed format, required for 10x."},
                 {OUT,               0, "o", "out",                    Arg::Required, "\t-o <arg>,\t--out=<arg>\tOutput VCF File."},
                 {_WINDOW_SIZE,      0, "", "window_size",             Arg::Numeric,  "\t-w [<arg>],\t--window_size [<arg>]\tPhasing Window Size, default=250."},
@@ -71,6 +83,7 @@ const option::Descriptor usage[] =
                 {HIC,               0, "H", "hic",                    Arg::None,     "\t-H,\t--hic\tSpecified for HiC data."},
                 {PACBIO,            0, "P", "pacbio",                 Arg::None,     "\t-P,\t--pacbio\tSpecified for Pacbio data."},
                 {NANOPORE,          0, "N", "nanopore",               Arg::None,     "\t-N,\t--nanopore\tSpecified for Nanopore data."},
+                {PROTOCOLS,         0,"p","protocols",                Arg::Required, "\t-p,\t--protocols\t Sequence protocols for corresponding to frags, split with comma, hic,ngs, tenx, pacbio, nanopore is supported"},
                 {_HYBRID,           0, "", "hybrid",                  Arg::None,     "\t--hybrid\tSpecified for hybrid data type."},
                 {_NEWFORMAT,        0, "", "new_format",              Arg::None,     "\t--new_format\tSpecified when using new_format with extractHair"},
                 {_BASEOFFSET,       0, "", "base_offset",             Arg::Numeric,  "\t--base_offset\tQuality of set for read base, default is 33."},
@@ -138,30 +151,54 @@ int main(int argc, char *argv[])
         WINDOW_SIZE = int(atoi(options[_WINDOW_SIZE].arg));
 
     std::string fnbed; 
-    if (options[TENX])
-    {    
-        OPERATION = MODE_10X;
-        fnbed = options[STATS].arg;
-    }
-    else if (options[HIC])
-        OPERATION = MODE_HIC;
-    else if (options[PACBIO])
-        OPERATION = MODE_PACBIO;
-    else if (options[NANOPORE])
-        OPERATION = MODE_NANOPORE;
-
-    if (options[_HYBRID])   //hybrid mode, using pacbio/oxford nanopore, ngs and hic. 
-    {    
-        HYBRID = true;
-        OPERATION = MODE_HYBRID;        
-    }
-    
+//    if (options[TENX])
+//    {
+//        OPERATION = MODE_10X;
+//        fnbed = options[STATS].arg;
+//    }
+//    else if (options[HIC])
+//        OPERATION = MODE_HIC;
+//    else if (options[PACBIO])
+//        OPERATION = MODE_PACBIO;
+//    else if (options[NANOPORE])
+//        OPERATION = MODE_NANOPORE;
+//
+//    if (options[_HYBRID])   //hybrid mode, using pacbio/oxford nanopore, ngs and hic.
+//    {
+//        HYBRID = true;
+//        OPERATION = MODE_HYBRID;
+//    }
+//
     std::string invcf = options[VCF].arg;
     std::string out = options[OUT].arg;
     std::string frag = options[FRAGMENT].arg;
-    
-
-    Phaser *phaser = new Phaser(invcf, out, frag, fnbed);
+//    multiply frags file support
+    std::vector<std::string> frags;
+    std::string token;
+    std::istringstream iss(options[FRAGMENT].arg);
+    while(std::getline(iss, token, ','))
+        frags.push_back(token);
+//    protocols
+    std::istringstream protocols(options[PROTOCOLS].arg);
+    while(std::getline(protocols, token, ',')) {
+        if (token == PROTOCOL_HIC){
+            OPERATIONS.push_back(MODE_HIC);
+            HAS_HIC = true;
+        }
+        else if (token == PROTOCOL_NGS) OPERATIONS.push_back(MODE_PE);
+        else if (token == PROTOCOL_TENX) {
+            HAS_TENX = true;
+            OPERATIONS.push_back(MODE_10X);
+            fnbed = options[STATS].arg;
+        } else if (token == PROTOCOL_NANOPORE) OPERATIONS.push_back(MODE_NANOPORE);
+        else if (token == PROTOCOL_PACBIO) OPERATIONS.push_back(MODE_PACBIO);
+        else if (token == PROTOCOL_HYBRID) OPERATIONS.push_back(MODE_HYBRID);
+        else {
+            std::cerr << "SpecHap: Error, --protocols muste be split with comma, and only hic, ngs, tenx, pacbio and nanopore are supported\n";
+            exit(1);
+        }
+    }
+    auto *phaser = new Phaser(invcf, out, frags, fnbed);
     phaser->phasing();
     delete phaser;
     return 0;

Diff for: phaser.cpp

@@ -12,32 +12,46 @@
 
 // TODO clarify between variant count and block count
 
-Phaser::Phaser(const std::string &fnvcf, const std::string &fnout, const std::string &fnfrag, const std::string &fnbed)
+Phaser::Phaser(const std::string &fnvcf, const std::string &fnout, std::vector<std::string> & fnfrags, const std::string &fnbed)
 {
-    
+    for (int i = 0; i < fnfrags.size(); i++) {
+        auto item = fnfrags[i];
+        frfrags.push_back(new FragmentReader(item.data()));
+        if(OPERATIONS[i] == MODE_10X) {
+            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());
+//    frfrag = new FragmentReader(fnfrag.data());
     frbed = nullptr;
     coverage = 30;  //deprecated
 
-    if (OPERATION == MODE_10X)
+    if (HAS_TENX)
         frbed = new BEDReader(fnbed.data());
 
     bool use_secondary = false;
     threshold = 1e-5;
 
-    spectral = new Spectral(frfrag, frbed, threshold, coverage, use_secondary);
+    spectral = new Spectral(frfrags, frbed, threshold, coverage, use_secondary);
 }
 
 Phaser::~Phaser()
 {
     delete frvcf;
     delete fwvcf;
-    delete frfrag;
+//    delete frfrag;
+    for (auto & item : frfrags) {
+        delete item;
+        item = nullptr;
+    }
+    frfrags.clear();
+    frfrags.shrink_to_fit();
     delete spectral;
-    if (frbed != nullptr)
-        delete frbed;
+    delete frbed;
 }
 
 
@@ -127,7 +141,10 @@ void Phaser::phasing()
         else
             load_contig_records(chromo_phaser);
         chromo_phaser->construct_phasing_window_initialize();
-        frfrag->set_prev_chr_var_count(prev_variant_count);
+        for (auto item : frfrags) {
+            item->set_prev_chr_var_count(prev_variant_count);
+        }
+//        frfrag->set_prev_chr_var_count(prev_variant_count);
         spectral->set_chromo_phaser(chromo_phaser);
         phasing_by_chrom(chromo_phaser->variant_count, chromo_phaser);
 
@@ -142,13 +159,16 @@ void Phaser::phasing()
 
 void Phaser::phasing_by_chrom(uint var_count, ChromoPhaser *chromo_phaser)
 {
-    frfrag->set_curr_chr_var_count(var_count);
+    for (auto item : frfrags) {
+        item->set_curr_chr_var_count(var_count);
+    }
+//    frfrag->set_curr_chr_var_count(var_count);
 
     while (chromo_phaser->phased->rest_blk_count > 0)
     {
         if (chromo_phaser->phased->rest_blk_count > chromo_phaser->init_block_count)
             break;
-        if (OPERATION == MODE_10X)
+        if (HAS_TENX)
             chromo_phaser->phased->update_phasing_info(MAX_BARCODE_SPANNING);
         else
             {
@@ -160,7 +180,7 @@ void Phaser::phasing_by_chrom(uint var_count, ChromoPhaser *chromo_phaser)
         spectral->solver();
     }
         //std::cout << chromo_phaser->phased->phased_blk_count << std::endl;
-    if (OPERATION == MODE_HIC)
+    if (HAS_HIC)
     {
         phase_HiC_poss(chromo_phaser);
     }

Diff for: phaser.h

@@ -24,7 +24,7 @@ class Phaser
 {
 public:
     Phaser() = default;
-    explicit Phaser(const std::string & fnvcf, const std::string & fnout, const std::string & fnfrag, const std::string &fnbed);
+    explicit Phaser(const std::string & fnvcf, const std::string & fnout, std::vector<std::string> & fnfrags, const std::string &fnbed);
     ~Phaser();
     void phasing();
 
@@ -33,7 +33,8 @@ class Phaser
     double threshold;
     VCFReader *frvcf;
     VCFWriter *fwvcf;
-    FragmentReader *frfrag;
+    std::vector<FragmentReader* > frfrags;
+//    FragmentReader *frfrag;
     BEDReader *frbed;
     Spectral *spectral;
     int coverage;