c++/EEDB/Tools/Paraclu.cpp
/////////////////////////////////////////////////////////////////////////////////////
//
// original ported version. close to direct interpretation of the original ParaClu algorithm
// but flipped recursion logic so that parent decides based on children and self
// also working with the sliding input buffer and the ZENBU datamodel objects
//
/////////////////////////////////////////////////////////////////////////////////////
original ported version. close to direct interpretation of the original ParaClu algorithm
but flipped recursion logic so that parent decides based on children and self
also working with the sliding input buffer and the ZENBU datamodel objects
original ported version . close to direct interpretation of the original ParaClu algorithm but flipped recursion logic so that parent decides based on children and self also working with the sliding input buffer and the ZENBU datamodel objects
void  EEDB::Tools::Paraclu::_process_cluster_range_full(FeatureLink *beg, FeatureLink *end, double parentDensity) {
  if((beg==NULL) || (end==NULL)) { return; }
  long   cluster_length = end->feature->chrom_start() - beg->feature->chrom_start() + 1;
  double totalValue     = _calc_total_value(beg, end);
  double density        = totalValue / cluster_length;
  if((totalValue < _min_cutoff_value) || (density < _min_density)) { 
    if(cluster_length > _max_cluster_length) { 
      _discard_region(beg, end);
    }
    return; 
  } 
  if(beg==end) {  
    double stability = density/parentDensity; 
    if((_selection_mode==FULL_HIERARCHY) || ((_selection_mode==STABILITY_CUT) && (stability >= _min_stability))) {
      EEDB::Feature *feature = _generate_cluster(beg, end);
      if(feature) {
        feature->metadataset()->add_tag_data("density", d_to_string(density));
        feature->metadataset()->add_tag_data("maxDensity", d_to_string(density));  
        feature->metadataset()->add_tag_data("parentDensity", d_to_string(parentDensity));
        feature->metadataset()->add_tag_data("stability", d_to_string(stability));
        _output_buffer->insert_feature(feature);
      } 
    }
    return;
  }
   find best mid-cut-point
  _weakestPrefix(beg, end);
  _weakestSuffix(beg, end);
  double maxDensity;
  FeatureLink *mid = NULL;
  if(_minPrefixDensity < _minSuffixDensity) {
    maxDensity = _minPrefixDensity;
    mid = _minPrefix;
    fprintf(stderr, "choose prefix %f %d\n", _minPrefixDensity, _minPrefix->feature->chrom_start());
  } else {
    maxDensity = _minSuffixDensity;
    mid = _minSuffix;
    fprintf(stderr, "choose suffix %f %d\n", _minSuffixDensity, _minSuffix->feature->chrom_start());
  }  
  if no mid point pop out
  if(mid==NULL) { 
    did not find a cut so bounce out
    fprintf(stderr, "Paraclu PROBLEM mid is null\n");
    return;
  }
  printf("  maxDensity %f\n", maxDensity);
  printf("  parentDensity %f\n", parentDensity);
  printf("  mid link %s\n", mid->feature->chrom_location().c_str());
  double child_stability = maxDensity / parentDensity;
  if((cluster_length <= _max_cluster_length) && (maxDensity>parentDensity)) {
    if((_selection_mode==FULL_HIERARCHY) || ((_selection_mode==STABILITY_CUT) && (child_stability >= _min_stability))) {
      EEDB::Feature *feature = _generate_cluster(beg, end);
      if(feature) {
        feature->metadataset()->add_tag_data("density", d_to_string(density));
        feature->metadataset()->add_tag_data("maxDensity", d_to_string(maxDensity));
        feature->metadataset()->add_tag_data("parentDensity", d_to_string(parentDensity));
        feature->metadataset()->add_tag_data("child_stability", d_to_string(child_stability));
        _output_buffer->insert_feature(feature);
        fprintf(stderr, "make cluster %ld .. %ld  [ %ld ] [total=%f density=%f]\n", beg->feature->chrom_start(), end->feature->chrom_start(), cluster_length, totalValue, density);
      }
    }
    if((_selection_mode==STABILITY_CUT) && (child_stability >= _min_stability)) {
       no need to dig deeper, martin's paraclu-cut picks biggest cluster in hierarchy
       Remove clusters that are contained in larger clusters.  This assumes
       that the input is sorted by start position (ascending) then end
       position (descending).
      fprintf(stderr, "stop stability_cut\n");
      return;
    }
  }
  double newMinDensity = max(parentDensity, maxDensity);
  double newMinDensity = parentDensity;
  if(cluster_length <= _max_cluster_length) { newMinDensity = max(parentDensity, maxDensity); }
  FeatureLink *mid1 = mid->prev;
  long   left_length  = mid1->feature->chrom_start() - beg->feature->chrom_start() + 1;
  fprintf(stderr, "left\n");
  _process_cluster_range_full(beg, mid1, newMinDensity);
  only dig down right side when inside left most cluster <= max_cluster_length
  if(cluster_length <= _max_cluster_length) { 
    fprintf(stderr, "right\n");
    _process_cluster_range_full(mid, end, newMinDensity);
  }
  special transition, from above max_cluster into left-most cluster
  once finished here we poping completely out so discard this left-most region
  if((left_length <= _max_cluster_length) && (cluster_length > _max_cluster_length)) {
    fprintf(stderr, "discard left-side %ld .. %ld  [ %ld ] [total=%f density=%f]\n", beg->feature->chrom_start(), end->feature->chrom_start(), cluster_length, totalValue, density);
    _discard_region(beg, mid1);
  }
}
void EEDB :: Tools :: Paraclu :: _process_cluster_range_full ( FeatureLink * beg , FeatureLink * end , double parentDensity ) { if ( ( beg == NULL ) || ( end == NULL ) ) { return ; } long cluster_length = end -> feature -> chrom_start ( ) - beg -> feature -> chrom_start ( ) + 1 ; double totalValue = _calc_total_value ( beg , end ) ; double density = totalValue / cluster_length ; if ( ( totalValue < _min_cutoff_value ) || ( density < _min_density ) ) { if ( cluster_length > _max_cluster_length ) { _discard_region ( beg , end ) ; } return ; } if ( beg == end ) { double stability = density / parentDensity ; if ( ( _selection_mode == FULL_HIERARCHY ) || ( ( _selection_mode == STABILITY_CUT ) && ( stability >= _min_stability ) ) ) { EEDB :: Feature * feature = _generate_cluster ( beg , end ) ; if ( feature ) { feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( density ) ) ; feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( density ) ) ; feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( parentDensity ) ) ; feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( stability ) ) ; _output_buffer -> insert_feature ( feature ) ; } } return ; } _weakestPrefix ( beg , end ) ; _weakestSuffix ( beg , end ) ; double maxDensity ; FeatureLink * mid = NULL ; if ( _minPrefixDensity < _minSuffixDensity ) { maxDensity = _minPrefixDensity ; mid = _minPrefix ; } else { maxDensity = _minSuffixDensity ; mid = _minSuffix ; } if ( mid == NULL ) { return ; } double child_stability = maxDensity / parentDensity ; if ( ( cluster_length <= _max_cluster_length ) && ( maxDensity > parentDensity ) ) { if ( ( _selection_mode == FULL_HIERARCHY ) || ( ( _selection_mode == STABILITY_CUT ) && ( child_stability >= _min_stability ) ) ) { EEDB :: Feature * feature = _generate_cluster ( beg , end ) ; if ( feature ) { feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( density ) ) ; feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( maxDensity ) ) ; feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( parentDensity ) ) ; feature -> metadataset ( ) -> add_tag_data ( " " , d_to_string ( child_stability ) ) ; _output_buffer -> insert_feature ( feature ) ; } } if ( ( _selection_mode == STABILITY_CUT ) && ( child_stability >= _min_stability ) ) { return ; } } double newMinDensity = max ( parentDensity , maxDensity ) ; FeatureLink * mid1 = mid -> prev ; long left_length = mid1 -> feature -> chrom_start ( ) - beg -> feature -> chrom_start ( ) + 1 ; _process_cluster_range_full ( beg , mid1 , newMinDensity ) ; if ( cluster_length <= _max_cluster_length ) { _process_cluster_range_full ( mid , end , newMinDensity ) ; } if ( ( left_length <= _max_cluster_length ) && ( cluster_length > _max_cluster_length ) ) { _discard_region ( beg , mid1 ) ; } }
original ported version.
original ported version .
//OK this is a closer attempt at martin's original algorithm, full hierarchy //recursive algorithm which subdivides sites //original paraclu called the parameter minDensity, but parentDensity is better name //if(beg==end) {  return; } //fprintf(stderr, "_process_cluster_range_full  %ld .. %ld  [ %ld ] [total=%f density=%f]\n", beg->feature->chrom_start(), end->feature->chrom_start(), cluster_length, totalValue, density); //weak region so can fall out //if >max_length then we are walking down the left-side, so can discard and fall out //fprintf(stderr, "discard large-weak region %ld .. %ld  [ %ld ] [total=%f density=%f]\n", beg->feature->chrom_start(), end->feature->chrom_start(), cluster_length, totalValue, density); //my stability //fprintf(stderr, "_process_cluster_range  beg == end density[%f] parentDensity[%f] stability[%f] minstab[%f]\n", density, parentDensity, stability, _min_stability); //this is part of the hierarchy //fprintf(stderr, "single 1bp stable cluster [%ld] & STABLE [%f]\n", cluster_length, stability); //1bp so no children // // find best mid-cut-point // //fprintf(stderr, "choose prefix %f %d\n", _minPrefixDensity, _minPrefix->feature->chrom_start()); //fprintf(stderr, "choose suffix %f %d\n", _minSuffixDensity, _minSuffix->feature->chrom_start()); //if no mid point pop out //did not find a cut so bounce out //fprintf(stderr, "Paraclu PROBLEM mid is null\n"); //printf("  maxDensity %f\n", maxDensity); //printf("  parentDensity %f\n", parentDensity); //printf("  mid link %s\n", mid->feature->chrom_location().c_str()); //this is part of the hierarchy /*
     std::cerr << beg->feature->chrom_name() << "\t" << beg->feature->strand() << "\t"
     << std::setprecision(20)
     << beg->feature->chrom_start() << "\t" << end->feature->chrom_start() << "\t"
     << cluster_length << "\t" << totalValue << "\t"
     << std::setprecision(3)
     << parentDensity << "\t" << maxDensity ;
     std::cerr << "\t" << stability;
     if(stability >= _min_stability) { std::cerr << "\tSTABLE"; }
     std::cerr << "\n";
     */ //fprintf(stderr, "make cluster %ld .. %ld  [ %ld ] [total=%f density=%f]\n", beg->feature->chrom_start(), end->feature->chrom_start(), cluster_length, totalValue, density); // no need to dig deeper, martin's paraclu-cut picks biggest cluster in hierarchy // Remove clusters that are contained in larger clusters.  This assumes // that the input is sorted by start position (ascending) then end // position (descending). //fprintf(stderr, "stop stability_cut\n"); //double newMinDensity = parentDensity; //if(cluster_length <= _max_cluster_length) { newMinDensity = max(parentDensity, maxDensity); } //fprintf(stderr, "left\n"); //only dig down right side when inside left most cluster <= max_cluster_length //fprintf(stderr, "right\n"); //special transition, from above max_cluster into left-most cluster //once finished here we poping completely out so discard this left-most region //fprintf(stderr, "discard left-side %ld .. %ld  [ %ld ] [total=%f density=%f]\n", beg->feature->chrom_start(), end->feature->chrom_start(), cluster_length, totalValue, density);