Merge pull request #33 from gatoravi/cis-ase-leftout

Left out a commit in #24

src/cis-ase/cis_ase_identifier.cc

@@ -111,6 +111,8 @@ void CisAseIdentifier::parse_options(int argc, char* argv[]) {
     cerr << "\nReference fasta file: " << ref_;
     cerr << "\nAnnotation file: " << gtf_;
     cerr << "\nMinimum read-depth for variants: " << min_depth_;
+    cerr << "\nWindow around somatic-variants to look for transcripts: " <<
+            transcript_variant_window_;
     cerr << endl;
 }
 
@@ -275,7 +277,7 @@ bool CisAseIdentifier::process_rna_hom(bcf_hdr_t* bcf_hdr, int tid,
 }
 
 //Get the window pertinent to this variant
-//Get the transcripts within a certain window
+//Get the transcripts within a certain distance from the somatic variant
 //Return the window that encompasses all these transcripts.
 string CisAseIdentifier::get_relevant_window(const char* chr, int pos) {
     CHRPOS min_start = pos;
@@ -405,30 +407,17 @@ void CisAseIdentifier::process_snps_in_window(string region) {
 
 //ASE identification starts here
 void CisAseIdentifier::identify_ase() {
-    htsFile *test_bcf = NULL;
-    bcf_hdr_t *test_header = NULL;
-    bcf1_t *test_record = bcf_init();
-    test_bcf = bcf_open(somatic_vcf_.c_str(), "r");
-    if(test_bcf == NULL) {
-        throw std::runtime_error("Unable to open file.");
-    }
-    test_header = bcf_hdr_read(test_bcf);
-    if(test_header == NULL) {
-        throw std::runtime_error("Unable to read header.");
-    }
-    while(bcf_read(test_bcf, test_header, test_record) == 0) {
-        string somatic_region = common::create_region_string(bcf_hdr_id2name(test_header, test_record->rid),
-                                                             test_record->pos+1, test_record->pos+1);
+    while(bcf_read(somatic_vcf_fh_,
+                   somatic_vcf_header_, somatic_vcf_record_) == 0) {
+        string somatic_region = common::create_region_string(bcf_hdr_id2name(somatic_vcf_header_, somatic_vcf_record_->rid),
+                                                             somatic_vcf_record_->pos+1, somatic_vcf_record_->pos+1);
         cerr << endl << "somatic region is " << somatic_region << endl;
         set_mpileup_conf_region(somatic_conf_, somatic_region);
         mpileup_run(&somatic_conf_,
            &CisAseIdentifier::process_somatic_het,
            somatic_dna_mmc_);//The workhorse
         free_mpileup_conf(somatic_conf_);
     }
-    bcf_hdr_destroy(test_header);
-    bcf_destroy(test_record);
-    bcf_close(test_bcf);
 }
 
 //Free relevant pointers