Add output VCF file option

src/cis-ase/cis_ase_identifier.cc

@@ -125,6 +125,9 @@ void CisAseIdentifier::parse_options(int argc, char* argv[]) {
     if(use_binomial_model_) {
         cerr << "\nUsing the binomial model for modeling RNAseq ASE";
     }
+    if(output_file_ != "NA") {
+        cerr << "\nWriting VCF output to " << output_file_;
+    }
     cerr << endl;
 }
 
@@ -431,7 +434,7 @@ void CisAseIdentifier::process_snps_in_window(string somatic_region, BED region)
                                 &CisAseIdentifier::process_germline_het,
                                 germline_dna_mmc_)) {
                         cerr << "DNA is het. potential ASE " << snp_region << endl;
-                        vcf_op_.print_line();
+                        vcf_op_.print_line(ofs_);
                     } else {
                         cerr << "DNA not het" << endl;
                     }
@@ -507,11 +510,12 @@ void CisAseIdentifier::run() {
     mmc_init_all(); //load all the mmcs
     load_reference(); //load reference genome
     gtf_parser_.load(); //load gene annotations
+    set_ostream(); //Set the output stream
     annotate_exonic_polymorphisms();
     open_somatic_vcf();
     open_poly_vcf();
     mpileup_init_all();
-    vcf_op_.print_header();
+    vcf_op_.print_header(ofs_);
     identify_ase();//Start running the pileups and looking at GTs
     cleanup();//Cleanup file handles
 }

src/cis-ase/cis_ase_identifier.h

@@ -118,13 +118,13 @@ struct VcfRecord {
         p_hom_rna = -1;
     }
     //Print the variant line
-    void print_line() {
+    void print_line(ostream& out = std::cout) {
         string id = ".";
         string qual = ".";
         //This variant satisfies ASE criterion
         string filter = "PASS";
         string info = construct_info();
-        cout << chr << "\t" <<
+        out << chr << "\t" <<
                 pos << "\t" <<
                 id << "\t" <<
                 ref << "\t" <<
@@ -140,18 +140,18 @@ struct VcfRecord {
                "P_HOM_RNA=" + common::num_to_str(p_hom_rna);
     }
     //VCF header
-    void print_header() {
-        cout << "##fileformat=VCFv4.2" << endl;
-        cout << "##INFO=<ID=SOMATIC_VARIANT,Number=1,Type=String,"
+    void print_header(ostream& out = std::cout) {
+        out << "##fileformat=VCFv4.2" << endl;
+        out << "##INFO=<ID=SOMATIC_VARIANT,Number=1,Type=String,"
                 "Description=\"Somatic variant proximal to ASE variant.\"";
-        cout << endl;
-        cout << "##INFO=<ID=P_HET_DNA,Number=1,Type=Float,"
+        out << endl;
+        out << "##INFO=<ID=P_HET_DNA,Number=1,Type=Float,"
                 "Description=\"Posterior probability of het in the DNA at ASE site.\"";
-        cout << endl;
-        cout << "##INFO=<ID=P_HOM_RNA,Number=1,Type=Float,"
+        out << endl;
+        out << "##INFO=<ID=P_HOM_RNA,Number=1,Type=Float,"
                 "Description=\"Posterior probability of hom in the RNA at ASE site.\"";
-        cout << endl;
-        cout << "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO" << endl;
+        out << endl;
+        out << "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO" << endl;
     }
     //Reset all fields
     void reset() {
@@ -355,7 +355,7 @@ class CisAseIdentifier {
         string output_file_;
         //Which polymorphisms to look at
         string relevant_poly_annot_;
-        //output stream to output annotated junctions file
+        //output stream to output ASE variants in VCF format
         ofstream ofs_;
         //Somatic VCF file handle
         htsFile *somatic_vcf_fh_;
@@ -460,6 +460,19 @@ class CisAseIdentifier {
         //create the map, where list of exonic variants are
         //indexed by "chr:bin"
         void annotate_exonic_polymorphisms();
+        //If output file is not empty, attempt to open
+        //If output file is empty, set to cout
+        void set_ostream() {
+            if(output_file_ == "NA") {
+                common::copy_stream(cout, ofs_);
+                return;
+            }
+            ofs_.open(output_file_.c_str());
+            if(!ofs_.is_open()) {
+                throw runtime_error("Unable to open " +
+                        output_file_);
+            }
+        }
 };
 
 #endif //CIS_ASE_IDENTIFIER_