Major modification on demuxlet and freemuxlet, including faster versi…

…on of freemuxlet
statgen · Sep 20, 2019 · 7b141e3 · 7b141e3
1 parent 97e351b
commit 7b141e3
Show file tree

Hide file tree

Showing 25 changed files with 2,234 additions and 3,556 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -25,7 +25,6 @@ include_directories(${HTS_INCLUDE_DIRS})
 add_executable(popscle
         allele.cpp
         allele.h
-        bam_endian.h
         bam_ordered_reader.cpp
         bam_ordered_reader.h
         bcf_chunked_reader.cpp
@@ -35,23 +34,13 @@ add_executable(popscle
         bcf_filtered_reader.h
         bcf_variant_key.h
         bgzf.h
-        bh_tsne.cpp
-        bh_tsne.h
         cmd_cram_demuxlet.cpp
         cmd_cram_digital_pileup.cpp
         cmd_cram_freemuxlet.cpp
-        cmd_cram_gtf_util.cpp
-        cmd_cram_mux_pileup.cpp
-        cmd_cram_simuxlet.cpp
-        cmd_plp_find_mux_clust.cpp
+	cmd_cram_freemux2.cpp
         cmd_plp_make_dge_matrix.cpp
-        cmd_sc_kallisto_count.cpp
-        cmd_sc_map_stamps.cpp
-        cmd_sc_multinom_em.cpp
-        cmd_sc_multinom_gibbs.cpp
         commands.cpp
         commands.h
-        config.h
         Constant.h
         cramore.cpp
         cramore.h
@@ -78,9 +67,6 @@ add_executable(popscle
         interval_tree.h
         IO.cpp
         IO.h
-        knetfile.h
-        ksort.h
-        kstring.h
         log_tool.cpp
         log_tool.h
         louvain.h
@@ -97,9 +83,6 @@ add_executable(popscle
         sam_ordered_writer.h
         sc_drop_seq.cpp
         sc_drop_seq.h
-        sptree.cpp
-        sptree.h
-        tcf.h
         tsv_reader.cpp
         tsv_reader.h
         utils.cpp
@@ -110,7 +93,6 @@ add_executable(popscle
         variant_manip.h
         vntr.cpp
         vntr.h
-        vptree.h
         )
 
 find_library(ZLIB z HINTS /usr/lib/x86_64-linux-gnu/ /usr/lib/ /usr/lib64/)

diff --git a/bam_endian.h b/bam_endian.h
diff --git a/bam_ordered_reader.cpp b/bam_ordered_reader.cpp
@@ -100,7 +100,13 @@ BAMOrderedReader::BAMOrderedReader(std::string file_name, std::vector<GenomeInte
     interval_index = 0;
 
     random_access_enabled = intervals_present && index_loaded;
-};
+}
+
+BAMOrderedReader::~BAMOrderedReader()
+{
+    close();
+    //fprintf(stderr,"-- BAMOrderedReader object was deleted\n");
+}
 
 /**
  * Jump to interval. Returns false if not successful.
@@ -188,5 +194,17 @@ bool BAMOrderedReader::read(bam1_t *s)
  */
 void BAMOrderedReader::close()
 {
-    sam_close(file);
+    if ( hdr ) bam_hdr_destroy(hdr);
+    if ( itr ) bam_itr_destroy(itr);
+    if ( idx ) hts_idx_destroy(idx);
+    if ( s ) bam_destroy1(s);
+    if ( file ) sam_close(file);
+    if ( str.s ) free(str.s);
+
+    hdr = nullptr;
+    itr = nullptr;
+    idx = nullptr;
+    s = nullptr;
+    file = nullptr;
+    str.s = nullptr;
 }
diff --git a/bam_ordered_reader.h b/bam_ordered_reader.h
@@ -86,14 +86,7 @@ class BAMOrderedReader
      * @ref_fasta_file   reference FASTA file for CRAM
      */
     BAMOrderedReader(std::string file_name, std::vector<GenomeInterval>& intervals, std::string ref_fasta_file="");
-    ~BAMOrderedReader() {
-      if ( hdr ) bam_hdr_destroy(hdr);
-      if ( itr ) bam_itr_destroy(itr);
-      if ( idx ) hts_idx_destroy(idx);
-      if ( s ) bam_destroy1(s);
-      if ( str.s ) free(str.s);
-      //fprintf(stderr,"-- BAMOrderedReader object was deleted\n");
-    }
+    ~BAMOrderedReader();
 
     /**
      * Jump to interval. Returns false if not successful.

diff --git a/bcf_chunked_reader.cpp b/bcf_chunked_reader.cpp
@@ -264,23 +264,32 @@ bool BCFChunkedReader::initialize_current_interval() {
  */
 bool BCFChunkedReader::read(bcf1_t *v) {
   //bcf_clear(v);
-  if ( itr ) {
-    while( true ) {
+  //notice("BCFChunkedReader::read() called");
+  if ( itr ) {      // when reading by indices
+    while( true ) { // loop forever..
+      // first, 
       if ( itr && ( ( ( ftype.format == bcf ) && ( bcf_itr_next(file,itr,v) >= 0 ) ) ||
 		    ( ( ftype.format == vcf ) && ( tbx_itr_next(file,tbx,itr,&s) >= 0 ) ) ) ) {
-	if ( ( v->pos + 1 < chunk.chunk_intervals.it->beg1 ) || ( v->pos + 1 > chunk.chunk_intervals.it->end0 ) )
-	  continue;
-	if ( ftype.format == vcf ) {
+	if ( ftype.format == vcf ) 
 	  vcf_parse1(&s, hdr, v);
-	}
+
+	notice("Iterated a variant itr = %x, v->pos = %d, chr = %s, beg1 = %d, end0 = %d", itr, v->pos, chunk.chunk_intervals.it->chrom.c_str(), chunk.chunk_intervals.it->beg1, chunk.chunk_intervals.it->end0);
+
+	if ( ( v->pos + 1 < chunk.chunk_intervals.it->beg1 ) || ( v->pos + 1 > chunk.chunk_intervals.it->end0 ) )
+	    continue;
+
 	return true;
       }
       else {
+	notice("looking at next target intervals. itr = %x", itr);	
 	if ( target_intervals.next() ) {
-	  initialize_current_interval();
+	  notice("Initializing current interval");		  
+	  return initialize_current_interval();
 	}
 	else {
+	  notice("checking next file");	  
 	  if ( chunk.setNextFileName() ) {
+	    notice("setNextFileName returned true");	  	    
 	    if ( open_current_file() ) {
 	      // do nothing; current interval will be initialized
 	    }
@@ -292,6 +301,7 @@ bool BCFChunkedReader::read(bcf1_t *v) {
     }
   }
   else {
+    //notice("iterator absent");        
     if (bcf_read(file, hdr, v)==0)
       return true;
     else if ( chunk.setNextFileName() ) {