cmd_cram_dsc_pileup.cpp

#include "cramore.h"
#include "bcf_filtered_reader.h"
#include "sam_filtered_reader.h"
#include "sc_drop_seq.h"
#include "genomeLoci.h"

#include <cassert>

// TODO : Record strand info
// TODO : Reduce memory footprint
int32_t cmdCramDscPileup(int32_t argc, char** argv) {
  // input/output files
  SAMFilteredReader sr;
  BCFFilteredReader vr;
  std::string outPrefix;
  // tags for droplet barcodes and UMI barcodes
  std::string tagGroup("CB");
  std::string tagUMI("UB");
  int32_t capBQ = 40;
  int32_t minBQ = 13;
  int32_t minTD = 0;
  sr.filt.exclude_flag = 0x0f04;
  sr.filt.minMQ = 20;
  // input options for VCFs (sites)
  std::vector<std::string> smIDs;
  vr.vfilt.maxAlleles = 2;    
  // output options
  vr.verbose = 10000;
  sr.verbose = 1000000;
  // options to filter droplets
  std::string groupList;
  int32_t minTotalReads = 0;
  int32_t minUniqReads = 0;
  int32_t minCoveredSNPs = 0;
  bool skipUmiFlag = false;
  int32_t excludeFlag = 0x0704;

  paramList pl;

  BEGIN_LONG_PARAMS(longParameters)
    LONG_PARAM_GROUP("Options for input SAM/BAM/CRAM", NULL)
    LONG_STRING_PARAM("sam",&sr.sam_file_name, "Input SAM/BAM/CRAM file. Must be sorted by coordinates and indexed")
    LONG_STRING_PARAM("tag-group",&tagGroup, "Tag representing readgroup or cell barcodes, in the case to partition the BAM file into multiple groups. For 10x genomics, use CB")
    LONG_STRING_PARAM("tag-UMI",&tagUMI, "Tag representing UMIs. For 10x genomiucs, use UB")
    LONG_INT_PARAM("exclude-flag",&excludeFlag, "SAM/BAM flag to exclude")

    LONG_PARAM_GROUP("Options for input VCF/BCF", NULL)
    LONG_STRING_PARAM("vcf",&vr.bcf_file_name, "Input VCF/BCF file, containing the AC and AN field")
    LONG_MULTI_STRING_PARAM("sm",&smIDs, "List of sample IDs to compare to (default: use all)")
    LONG_STRING_PARAM("sm-list",&vr.sample_id_list, "File containing the list of sample IDs to compare")        

    LONG_PARAM_GROUP("Output Options", NULL)
    LONG_STRING_PARAM("out",&outPrefix,"Output file prefix")
    LONG_INT_PARAM("sam-verbose",&sr.verbose, "Verbose message frequency for SAM/BAM/CRAM")
    LONG_INT_PARAM("vcf-verbose",&vr.verbose, "Verbose message frequency for VCF/BCF")
    LONG_PARAM("skip-umi", &skipUmiFlag, "Do not generate [prefix].umi.gz file, which stores the regions covered by each barcode/UMI pair")

    LONG_PARAM_GROUP("SNP-overlapping Read filtering Options", NULL)
    LONG_INT_PARAM("cap-BQ", &capBQ, "Maximum base quality (higher BQ will be capped)")
    LONG_INT_PARAM("min-BQ", &minBQ, "Minimum base quality to consider (lower BQ will be skipped)")
    LONG_INT_PARAM("min-MQ", &sr.filt.minMQ, "Minimum mapping quality to consider (lower MQ will be ignored)")
    LONG_INT_PARAM("min-TD", &minTD, "Minimum distance to the tail (lower will be ignored)")
    LONG_INT_PARAM("excl-flag", &sr.filt.exclude_flag, "SAM/BAM FLAGs to be excluded")

    LONG_PARAM_GROUP("Cell/droplet filtering options", NULL)
    LONG_STRING_PARAM("group-list",&groupList, "List of tag readgroup/cell barcode to consider in this run. All other barcodes will be ignored. This is useful for parallelized run")    
    LONG_INT_PARAM("min-total", &minTotalReads, "Minimum number of total reads for a droplet/cell to be considered")
    LONG_INT_PARAM("min-uniq", &minUniqReads, "Minimum number of unique reads (determined by UMI/SNP pair) for a droplet/cell to be considered")
    LONG_INT_PARAM("min-snp", &minCoveredSNPs, "Minimum number of SNPs with coverage for a droplet/cell to be considered")
  END_LONG_PARAMS();

  pl.Add(new longParams("Available Options", longParameters));
  pl.Read(argc, argv);
  pl.Status();

  std::set<std::string> bcdSet;
  if ( !groupList.empty() ) {
    tsv_reader tsv_group_list(groupList.c_str());
    while( tsv_group_list.read_line() > 0 ) {
      bcdSet.insert(tsv_group_list.str_field_at(0));
    }
    notice("Finished loading %u droplet/cell barcodes to consider", bcdSet.size());
  }

  for(int32_t i=0; i < (int32_t)smIDs.size(); ++i) {
    vr.add_specified_sample(smIDs[i].c_str());
  }

  vr.unlimited_buffer = true;
  vr.vfilt.maxAlleles = 2;
  sr.set_buffer_size(1);
  //sr.unlimited_buffer = true;

  notice("Initializing BCF reader..");
  vr.init_params();
  notice("Initializing SAM reader..");  
  sr.init_params();

  int32_t n_warning_no_gtag = 0;
  int32_t n_warning_no_utag = 0;  
  
  if ( outPrefix.empty() )
    error("[E:%s:%d %s] --out parameter is missing",__FILE__,__LINE__,__PRETTY_FUNCTION__);

  char gtag[2] = {0,0};
  char utag[2] = {0,0};    

  if ( tagGroup.empty() ) { // do nothing
  }
  else if ( tagGroup.size() == 2 ) {
    gtag[0] = tagGroup.at(0);
    gtag[1] = tagGroup.at(1);    
  }
  else {
    error("[E:%s:%d %s] Cannot recognize group tag %s. It is suppose to be a length 2 string",__FILE__,__LINE__,__FUNCTION__,tagGroup.c_str());
  }

  if ( tagUMI.empty() ) { // do nothing
  }
  else if ( tagUMI.size() == 2 ) {
    utag[0] = tagUMI.at(0);
    utag[1] = tagUMI.at(1);    
  }
  else {
    error("[E:%s:%d %s] Cannot recognize UMI tag %s. It is suppose to be a length 2 string",__FILE__,__LINE__,__FUNCTION__,tagUMI.c_str());
  }    

  // scan VCF and CRAM simultaneously
  // read a variant first

  sc_dropseq_lib_t scl;

  std::vector<int32_t> snpids;
  //std::vector<int32_t> cellids;  
  
  /*
  if ( !vr.parse_posteriors(vr.cdr.hdr, vr.cursor(), field.c_str(), 0.01) )
    error("[E:%s] Cannot parse posterior probability at %s:%d", __PRETTY_FUNCTION__, bcf_hdr_id2name(vr.cdr.hdr,vr.cursor()->rid), vr.cursor()->pos+1);
  */

  // check if the chromosome names are in the same order between BCF and SAM
  std::map<int32_t,int32_t> rid2tids; // chromosome mapping VCF -> BAM
  std::map<int32_t,int32_t> tid2rids; // chromosome mapping BAM -> VCF
  std::vector<std::string> tchroms; // chromosomes appeared in SAM/BAM
  std::vector<std::string> rchroms; // chromosomes appeared in VCF/BCF
  int32_t ntids = bam_hdr_get_n_targets(sr.hdr);
  int32_t prevrid = -1;
  for(int32_t i=0; i < ntids; ++i) {
    const char* chrom = bam_get_chromi(sr.hdr, i);
    tchroms.push_back(chrom);
    int32_t rid = bcf_hdr_name2id(vr.cdr.hdr, chrom);
    if ( rid >= 0 ) {
      if ( prevrid >= rid ) {
	const char* prevchrom = bcf_hdr_id2name(vr.cdr.hdr, prevrid);
	error("[E:%s] Your VCF/BCF files and SAM/BAM/CRAM files have different ordering of chromosomes. SAM/BAM/CRAM file has %s before %s, but VCF/BCF file has %s after %s", __PRETTY_FUNCTION__, prevchrom, chrom, prevchrom, chrom);
      }
      rid2tids[rid] = i;
      tid2rids[i] = rid;
      rchroms.resize(rid+1);
      rchroms[rid] = chrom;
      prevrid = rid;
    }
  }

  if ( rid2tids.empty() || tid2rids.empty() || ( rid2tids.size() != tid2rids.size() ) ) {
    error("Your VCF/BCF files and SAM/BAM/CRAM files does not have any matching chromosomes, or some chromosome names are duplicated. This usually can be resolved by running 'bcftools reheader -f [ref.fasta.fai] [in.vcf.gz] > [out.vcf.gz]'");
  }
  notice("Identified %zu overlapping chromosomes between VCF and BAM", rid2tids.size());

  // read a single variant from VCF/BCF
  if ( !vr.read() )
    error("[E:%s Cannot read any single variant from %s]", __PRETTY_FUNCTION__, vr.bcf_file_name.c_str());  
  
  //int32_t nv = vr.get_nsamples();
  //double* gps = new double[nv*3];
  //for(int32_t i=0; i < nv * 3; ++i) 
  //  gps[i] = vr.get_posterior_at(i);
  int32_t snpid = scl.add_snp( vr.cursor()->rid, vr.cursor()->pos+1, vr.cursor()->d.allele[0][0], vr.cursor()->d.allele[1][0], vr.calculate_af(true), NULL);
  snpids.push_back(snpid);

  int32_t ibeg = 0;
  char base, qual;
  int32_t rpos;
  kstring_t readseq = {0,0,0};
  kstring_t readqual = {0,0,0};

  int32_t nReadsMultiSNPs = 0, nReadsSkipBCD = 0, nReadsPass = 0, nReadsRedundant = 0, nReadsN = 0, nReadsLQ = 0, nReadsTMP = 0, nReadsFilt = 0;

  // bcd -> umi -> bed : this will require a lot of memory
  std::vector<std::map<std::string, std::pair<genomeLoci,genomeLoci>> > umiLoci;  // fwd/rev pair
    
  while( sr.read() ) { // read SAM file, processing each individual read
    bam1_t* b = sr.cursor();    
    int32_t endpos = bam_endpos(b);
    const char* chrom = bam_get_chrom(sr.hdr, b);
    int32_t tid2rid = bcf_hdr_name2id(vr.cdr.hdr, chrom);
    bool noBCF = false;

    int32_t bamFlag = bam_get_flag(b);
    if ( bamFlag & excludeFlag )  {
      ++nReadsFilt;
      continue;
    }
    bool revStrand = (bamFlag & 0x0010) ? true : false;
    
    if ( tid2rid < 0 ) { // no matching BCF entry in the chromosome, skip;
      noBCF = true;
      //continue;
    }

    if ( vr.cursor()->rid >= rchroms.size() ) { // new contig absent in the header was introduced
      error("The contig %s was absent in the VCF header. This happens when contigs in VCF header do not match to actual records. Run 'bcftools reheader -f [ref.fasta.fai] [in.vcf.gz] > [out.vcf.gz] to resolve this problem", chrom);
    }

    int32_t n_cleared = vr.clear_buffer_before( bcf_hdr_id2name(vr.cdr.hdr, vr.cursor()->rid), b->core.pos );
    //for(int32_t i=ibeg; i < ibeg+n_cleared; ++i) {
    //  v_umis.clear();
    //}
    ibeg += n_cleared;
      
    // add new snps
    if ( !noBCF ) {
      while( ( !vr.eof ) && ( ( vr.cursor()->rid < tid2rid ) || ( ( vr.cursor()->rid == tid2rid ) && ( vr.cursor()->pos < endpos ) ) ) ) {
	if ( vr.read() ) {
	  double af = vr.calculate_af(true);
	  snpid = scl.add_snp( vr.cursor()->rid, vr.cursor()->pos+1, vr.cursor()->d.allele[0][0], vr.cursor()->d.allele[1][0], af, NULL);
	  snpids.push_back(snpid);
	}
	else {
	  //error("Cannot read new SNP");
	}
      }
    }
    
    // get barcode
    int32_t ibcd = 0;
    if ( tagGroup.empty() ) {
      ibcd = scl.add_cell(".");
    }
    else {
      uint8_t *bcd = (*gtag) ? (uint8_t*) bam_aux_get(b, gtag) : NULL;
      const char* sbcd = ".";
      if ( ( bcd != NULL ) && ( *bcd == 'Z' ) ) {
	sbcd = bam_aux2Z(bcd);
      }
      else {
	if ( n_warning_no_gtag < 10 ) {
	  notice("WARNING: Cannot find Droplet/Cell tag %s from %d-th read %s at %s:%d-%d. Treating all of them as a single group", tagUMI.c_str(), sr.n_read, bam_get_qname(b), chrom, b->core.pos, endpos);
	}
	else if ( n_warning_no_gtag == 10 ) {
	  notice("WARNING: Suppressing 10+ missing Droplet/Cell tag warnings...");
	}
	++n_warning_no_gtag;	
      }
      
      if ( bcdSet.empty() || ( bcdSet.find(sbcd) != bcdSet.end() ) ) {
	ibcd = scl.add_cell(sbcd);
      }
      else {
	++nReadsSkipBCD;
	continue;
      }
    }
    
    ++nReadsTMP;

    // get UMI
    std::string sumi(".");
    if ( tagUMI.empty() ) { // if UMI tag was not defined
      sumi.clear();
      catprintf(sumi,"%x",rand()); // give a random UMI
    }
    else {
      uint8_t *umi = (*utag) ? (uint8_t*) bam_aux_get(b, utag) : NULL;      
      if ( ( umi != NULL ) && ( *umi == 'Z' ) ) {
	sumi = bam_aux2Z(umi);
      }
      else {
	if ( n_warning_no_utag < 10 ) {
	  notice("WARNING: Cannot find UMI tag %s from %d-th read %s at %s:%d-%d. Treating all of them as a single UMI", tagUMI.c_str(), sr.n_read, bam_get_qname(b), bam_get_chrom(sr.hdr, b), b->core.pos, endpos);
	}
	else if ( n_warning_no_utag == 10 ) {
	  notice("WARNING: Suppressing 10+ UMI warnings...");
	}
	++n_warning_no_utag;
	//error("[E:%s] Cannot find UMI tag %d %d %x %s %s %x", __PRETTY_FUNCTION__, sr.nbuf, sr.ridx, sr.cursor(), bcd, utag, umi);
      }
    }

    
    ++scl.cell_totl_reads[ibcd];

    if ( !skipUmiFlag ) {  // count UMI
      if ( ibcd >= (int32_t)umiLoci.size() )
	umiLoci.resize((ibcd + 1) * 2);
      
      //if ( rand() % 100000 == 0 )
      //notice("ibcd = %d, sumi = %s", ibcd, sumi.c_str());
      
      genomeLoci& loci = revStrand ? umiLoci[ibcd][sumi].second : umiLoci[ibcd][sumi].first;
      if ( b->core.n_cigar ) { // go over CIGAR string to find 'M's
	//int32_t rlen = b->core.l_qseq;
	int32_t cpos = b->core.pos;
	int32_t rpos = 0;
	kstring_t str = {0, 0, 0};
	
	uint32_t *cigar = bam_get_cigar(b);
	for (uint32_t i = 0; i < b->core.n_cigar; ++i) {
	  char op = bam_cigar_opchr(cigar[i]);
	  str.l = 0;
	  kputw(bam_cigar_oplen(cigar[i]), &str);
	  char* stop;
	  int32_t len = strtol(str.s, &stop, 10);
	  assert(stop);
	  
	  // if M is observed, add the region to loci
	  if (op=='M') {
	    loci.add(chrom, cpos+1, cpos+len);
	    cpos += len;
	    rpos += len;
	  }
	  else if ( ( op=='D' ) || ( op=='N' ) ) {
	    cpos += len;
	  }
	  else if (op=='S' || op=='I') {
	    rpos += len;
	  }
	}
	loci.resolveOverlaps();
      }
      //umiLoci[ibcd][sumi].add( chrom, b->core.pos+1, endpos );
    }

    if ( noBCF ) continue;  // skip the part that reads variants overlaps info
    
    // genotype all reads together
    int32_t nv_pass = 0;
    int32_t nv_redundant = 0;
    int32_t nv_valid = 0;
    int32_t allele, bq;

    //if ( rand() % 10000 == 0 )
    //notice("Reading between %s:%d-%d at %s:%d to %d i=beg=%d, nbuf=%d, vidx=%d, size=%u prevpos=%d", bam_get_chrom(sr.hdr, sr.cursor()), sr.cursor()->core.pos+1, bam_endpos(sr.cursor()), bcf_hdr_id2name(vr.cdr.hdr, scl.snps[ibeg].rid), scl.snps[ibeg].pos, scl.snps[ibeg+vr.nbuf-1].pos, ibeg, vr.nbuf, vr.vidx, vr.vbufs.size(), scl.snps[ibeg-1].pos);
    
    for(int32_t i=ibeg; i < ibeg+vr.nbuf; ++i) {
      bam_get_base_and_qual_and_read_and_qual(b, (uint32_t)scl.snps[i].pos-1, base, qual, rpos, &readseq, &readqual);
      if ( rpos == BAM_READ_INDEX_NA ) {
	//if ( rand() % 1000 == 0 ) 
	//notice("Cannot find any informative read between %s:%d-%d at %s:%d", bam_get_chrom(sr.hdr, b), b->core.pos+1, bam_endpos(b), bcf_hdr_id2name(vr.cdr.hdr, scl.snps[i].rid), scl.snps[i].pos);
	continue;
      }
      if ( base == 'N' ) continue;

      ++nv_valid;

      if ( qual-33 < minBQ ) { continue; }
      if ( rpos < minTD-1 ) { continue; }
      if ( rpos + minTD > b->core.l_qseq ) { continue; }

      allele = ( base == scl.snps[i].ref ) ? 0 : ( ( base == scl.snps[i].alt ) ? 1 : 2 );
      bq = qual-33 > capBQ ? capBQ : qual-33;

      if ( scl.add_read(snpids[i], ibcd, sumi.c_str(), allele, bq) )
	++nv_pass;
      else
	++nv_redundant;
    }

    if ( nv_pass > 1 ) ++nReadsMultiSNPs;
    if ( nv_pass > 0 ) ++nReadsPass;
    else if ( nv_redundant > 0 ) ++nReadsRedundant;
    else if ( nv_valid > 0 ) ++nReadsLQ;
    else ++nReadsN;
  }

  if ( n_warning_no_utag > 10 ) 
    notice("WARNING: Suppressed a total of %d UMI warnings...", n_warning_no_utag);
    
  if ( n_warning_no_gtag > 10 ) 
    notice("WARNING: Suppressed a total of %d droplet/cell barcode warnings...", n_warning_no_gtag);
  
  notice("Finished reading %d markers from the VCF file", (int32_t)snpids.size());

  //notice("Finished processing %d reads across %d variants across %d barcodes", nReadsPass, (int32_t)v_poss.size(), (int32_t)bcMap.size(), (int32_t)bcMap.size());
  notice("Total number input reads : %d", sr.n_read);
  notice("Total number of unmapped or QC-failed reads : %d", nReadsFilt);
  notice("Total number of read-QC-passed reads : %d ", sr.n_read - sr.n_skip); //, nReadsN + nReadsUnique + nReadsLQ + nReadsPass);
  notice("Total number of skipped reads with ignored barcodes : %d", nReadsSkipBCD);
  notice("Total number of non-skipped reads with considered barcodes : %d", nReadsTMP);  
  notice("Total number of gapped/noninformative reads : %d", nReadsN);
  notice("Total number of base-QC-failed reads : %d", nReadsLQ);  
  notice("Total number of redundant reads : %d", nReadsRedundant);
  notice("Total number of pass-filtered reads : %d", nReadsPass);
  notice("Total number of pass-filtered reads overlapping with multiple SNPs : %d", nReadsMultiSNPs);

  sr.close();
  //vr.close();

  notice("Starting to prune out cells with too few reads...");
  int32_t nRemoved = 0;
  if ( minTotalReads + minUniqReads + minCoveredSNPs < 0 ) {
    for(int32_t i=0; i < scl.nbcs; ++i) {
      if ( ( scl.cell_totl_reads[i] < minTotalReads ) || ( scl.cell_uniq_reads[i] < minUniqReads) || ( (int32_t)scl.cell_umis[i].size() < minCoveredSNPs ) ) {
	for(std::map<int32_t,sc_snp_droplet_t*>::iterator it = scl.cell_umis[i].begin();
	    it != scl.cell_umis[i].end(); ++it) {
	  delete it->second;
	  scl.snp_umis[it->first].erase(i);
	}
	scl.cell_umis[i].clear();
	++nRemoved;
      }
    }
  }
  notice("Finishing pruning out %d cells with too few reads...", nRemoved);


  if ( (int32_t)snpids.size() != scl.nsnps )
    error("[E:%s snpids.size() = %u != scl.nsnps = %d",__PRETTY_FUNCTION__, snpids.size(), scl.nsnps);

  // Calculate average genotype probability
  /*
  double* gp0s = (double*) calloc(scl.nsnps * 3, sizeof(double)); 
  for(int32_t i=0; i < scl.nsnps; ++i) {
    for(int32_t j=0; j < nv; ++j) {
      gp0s[i*3] += scl.snps[i].gps[3*j];
      gp0s[i*3+1] += scl.snps[i].gps[3*j+1];
      gp0s[i*3+2] += scl.snps[i].gps[3*j+2];      
    }
    gp0s[i*3] /= nv;
    gp0s[i*3+1] /= nv;
    gp0s[i*3+2] /= nv;    
    } */

  //std::map<int64_t, size_t> snpcell2idx;
  //std::map<int64_t, size_t> cellsnp2idx;  
  //std::vector<double> gls;

  htsFile* wC = hts_open((outPrefix+".cel.gz").c_str(),"wg");
  htsFile* wV = hts_open((outPrefix+".var.gz").c_str(),"wg");
  htsFile* wP = hts_open((outPrefix+".plp.gz").c_str(),"wg");
  
  if ( ( wC == NULL ) || ( wV == NULL ) || ( wP == NULL ) )    
    error("[E:%s:%d %s] Cannot create %s.* file",__FILE__,__LINE__,__FUNCTION__,outPrefix.c_str());

  notice("Writing cell information");
  std::vector<std::string> v_bcs(scl.bc_map.size());
  hprintf(wC, "#DROPLET_ID\tBARCODE\tNUM.READ\tNUM.UMI\tNUM.UMIwSNP\tNUM.SNP\n");
  for(std::map<std::string,int32_t>::iterator it = scl.bc_map.begin(); it != scl.bc_map.end(); ++it) {
    if ( !v_bcs[it->second].empty() )
      error("Duplicate position for barcode %s at %d", it->first.c_str(), it->second);
    v_bcs[it->second] = it->first;
    //hprintf(wC, "%s\t%d\n", it->first.c_str(), it->second);
  }
  for(int32_t i=0; i < (int32_t)v_bcs.size(); ++i) {
    if ( !skipUmiFlag ) {
      hprintf(wC, "%d\t%s\t%d\t%u\t%d\t%u\n", i, v_bcs[i].c_str(), scl.cell_totl_reads[i], umiLoci[i].size(), scl.cell_uniq_reads[i], scl.cell_umis[i].size());
    }
    else {
      hprintf(wC, "%d\t%s\t%d\t.\t%d\t%u\n", i, v_bcs[i].c_str(), scl.cell_totl_reads[i], scl.cell_uniq_reads[i], scl.cell_umis[i].size());      
    }
  }
  v_bcs.clear();
  notice("Finished writing cell information");
  hts_close(wC);

  notice("Build indices for sparse matrix");
  //size_t k;
  int32_t i;

  if ( !skipUmiFlag ) {
    htsFile* wU = hts_open((outPrefix+".umi.gz").c_str(),"wg");
    if ( wU == NULL )
      error("[E:%s:%d %s] Cannot create %s.umi.gz file",__FILE__,__LINE__,__FUNCTION__,outPrefix.c_str());
    
    for(i=0; i < (int32_t)umiLoci.size(); ++i) {  // for each cell barcode, print the number of distinct UMIs
      for(std::map<std::string, std::pair<genomeLoci,genomeLoci> >::iterator itu = umiLoci[i].begin();
	  itu != umiLoci[i].end(); ++itu) {
	if ( !( itu->second.first.empty() && itu->second.second.empty() ) ) { // at least fwd or rev is not empty
	  hprintf(wU, "%d\t%s\t%u\t%u", i, itu->first.c_str(), itu->second.first.totalLength(), itu->second.second.totalLength());
	  for(itu->second.first.rewind(); !itu->second.first.isend(); itu->second.first.next()) {
	    hprintf(wU, "\t%s:%d:%d:+", itu->second.first.it->chrom.c_str(), itu->second.first.it->beg1, itu->second.first.it->end0 - itu->second.first.it->beg1 + 1);
	  }
	  for(itu->second.second.rewind(); !itu->second.second.isend(); itu->second.second.next()) {
	    hprintf(wU, "\t%s:%d:%d:-", itu->second.second.it->chrom.c_str(), itu->second.second.it->beg1, itu->second.second.it->end0 - itu->second.second.it->beg1 + 1);
	  }	  
	  hprintf(wU, "\n");
	}
      }
    }
    hts_close(wU);
  }
  
  //double tmp;
  //for(i=0, k=0; i < scl.nsnps; ++i) {
  hprintf(wV, "#SNP_ID\tCHROM\tPOS\tREF\tALT\tAF\n");
  hprintf(wP, "#DROPLET_ID\tSNP_ID\tALLELES\tBASEQS\n");  
  for(i=0; i < scl.nsnps; ++i) {    
    hprintf(wV, "%d\t%s\t%d\t%c\t%c\t%.5lf\n", i, rchroms[scl.snps[i].rid].c_str(), scl.snps[i].pos, scl.snps[i].ref, scl.snps[i].alt, scl.snps[i].af); //0.5*gp0s[i*3+1] + gp0s[i*3+2]);
    
    std::map<int32_t,sc_snp_droplet_t*>& cells = scl.snp_umis[i];
    if ( cells.empty() ) continue;
    std::map<int32_t,sc_snp_droplet_t*>::iterator it;

    //gls.resize(3*(k+cells.size()));
    //double GLs[3];
    for(it = cells.begin(); it != cells.end(); ++it) {
      //GLs[0] = GLs[1] = GLs[2] = 1.0;
      std::string s_als;
      std::string s_bqs;

      for(sc_snp_droplet_it_t it2=it->second->begin(); it2 != it->second->end(); ++it2) {
	uint8_t al = ( it2->second >> 24 ) & 0x00ff;
	uint8_t bq = ( it2->second >> 16 ) & 0x00ff;
	s_als += ('0' + al);
	s_bqs += (33 + bq);
      }
      hprintf(wP, "%d\t%d\t%s\t%s\n", it->first, i, s_als.c_str(), s_bqs.c_str());
    }
  }

  notice("Finished builing indices for sparse matrix");
  hts_close(wP);
  hts_close(wV);  
  
  return 0;
}