<commit_before>#
# script to calculate the gene density features on a database
# should work on any species database
#

#
# It will only run on databases with genes ...
#

# I think the right thing here is too generate densities on the biggest 50-100 toplevel slices ...
# The website will be happy with about 150 bins I think.


use strict;

use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::DensityType;
use Bio::EnsEMBL::DensityFeature;
use Getopt::Long;

my $bin_count = 150;
my $long_slice_count = 100;


my ( $host, $user, $pass, $port, $dbname  );

$port = 3306 ; 

my ( $block_count, $genome_size, $block_size );

GetOptions( "host=s", \$host,
	    "user=s", \$user,
	    "pass=s", \$pass,
	    "port=i", \$port,
	    "dbname=s", \$dbname
	  );



unless ($host || $user || $pass || $dbname) { 
	print "\n\nusage : perl gene_density.pl -host <HOST> -user <USER> -pass <PASS> -port <3306> -dbname <DATABASENAME> \n\n" ; 
	exit(0) ; 
	}

my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(-host => $host,
					    -user => $user,
					    -port => $port,
					    -pass => $pass,
					    -dbname => $dbname);

my $sth = $db->dbc()->prepare( "select count(*) from gene" );
$sth->execute();

my ( $gene_count )  = $sth->fetchrow_array();

if( ! $gene_count ) {
  print STDERR "No gene density for $dbname.\n";
  exit();
}

#
# Could be database without seq_regions
#  Then have to try and attach core db
#
$sth = $db->dbc()->prepare( "select count(*)  from seq_region" );
$sth->execute();
my ( $seq_region_count ) = $sth->fetchrow_array();
if( ! $seq_region_count ) {
  #
  # for the time being only do core dbs
  # no dbs with no seq regions
  #
  print STDERR "No gene density for $dbname, no seq_regions.\n";
  exit();

  if( ($dbname =~ /_est_/ ) || ( $dbname =~ /_vega_/ ) || ( $dbname =~ /_cdna_/ ) ) {
    my $dna_db_name = $dbname;
    $dna_db_name =~ s/(_estgene_|_vega_|_cdna_)/_core_/;
    my $dna_db =  new Bio::EnsEMBL::DBSQL::DBAdaptor
      (-host => $host,
       -user => $user,
       -port => $port,
       -pass => <PASSWORD>,
       -dbname => $dna_db_name );
    print STDERR "Attaching $dna_db_name to $dbname.\n";
    $db->dnadb( $dna_db );
  } else {
    print STDERR "No gene density for $dbname, no seq_regions.\n";
    exit();
  }
}

#
# Get the adaptors needed;
#

my $dfa = $db->get_DensityFeatureAdaptor();
my $dta = $db->get_DensityTypeAdaptor();
my $aa  = $db->get_AnalysisAdaptor();
my $slice_adaptor = $db->get_SliceAdaptor();

#
# block size estimation
#

my $top_slices = $slice_adaptor->fetch_all('toplevel');
my @sorted_slices = sort { $b->seq_region_length() <=> $a->seq_region_length()} @$top_slices;

	
my $analysis = new Bio::EnsEMBL::Analysis (-program     => "gene_density_calc.pl",
					   -database    => "ensembl",
					   -gff_source  => "gene_density_calc.pl",
					   -gff_feature => "density",
					   -logic_name  => "knownGeneDensity");

$aa->store( $analysis );

$analysis = new Bio::EnsEMBL::Analysis (-program     => "gene_density_calc.pl",
					-database    => "ensembl",
					-gff_source  => "gene_density_calc.pl",
					-gff_feature => "density",
					-logic_name  => "geneDensity");

$aa->store( $analysis );

#
# Now the actual feature calculation loop
#


foreach my $known (1, 0) {
  #
  # Create new analysis object for density calculation.
  #

  if($known) {
    $analysis = $aa->fetch_by_logic_name('knownGeneDensity');
  } else {
    $analysis = $aa->fetch_by_logic_name('geneDensity');
  }

  #
  # Create new density type.
  #

  my $dt = Bio::EnsEMBL::DensityType->new(-analysis   => $analysis,
					  -region_features => $bin_count,
					  -value_type => 'sum');

  $dta->store($dt);

  my $slice_count = 0;
  my ( $current, $current_start, $current_end  );
  
  foreach my $slice (@sorted_slices){

    $block_size = $slice->length() / $bin_count;

    my @density_features=();

    print "Gene densities for ".$slice->seq_region_name().
      " with block size $block_size\n";
    $current_end = 0;
    $current = 0;

    while($current_end < $slice->end()) {
      $current += $block_size;
      $current_start = $current_end+1;
      $current_end = int( $current + 1 );

      if( $current_end < $current_start ) { 
	$current_end = $current_start;
      }

      if( $current_end > $slice->end() ) {
	$current_end = $slice->end();
      }


      my $sub_slice = $slice->sub_Slice( $current_start , $current_end );

      my $count =0;

      #
      # Store info for genes (ignore pseudo genes)
      #

      foreach my $gene (@{$sub_slice->get_all_Genes()}){
	if($gene->biotype() !~ /pseudogene/i and $gene->start >=1 ) {
	  $count++ if(!$known || $gene->is_known());
	}
      }

      push @density_features, Bio::EnsEMBL::DensityFeature->new
	(-seq_region    => $slice,
	 -start         => $current_start,
	 -end           => $current_end,
	 -density_type  => $dt,
	 -density_value => $count);
    }

    $dfa->store(@density_features);
    print "Created ", scalar @density_features, " gene density features.\n";

    last if( $slice_count++ > $long_slice_count );
  }
}







#
# helper to draw an ascii representation of the density features
#
sub print_features {
  my $features = shift;

  return if(!@$features);

  my $sum = 0;
  my $length = 0;
#  my $type = $features->[0]->{'density_type'}->value_type();

  print("\n");
  my $max=0;
  foreach my $f (@$features) {
    if($f->density_value() > $max){
      $max=$f->density_value();
    }
  }
  if( !$max ) { $max = 1 };

  foreach my $f (@$features) {
    my $i=1;
    for(; $i< ($f->density_value()/$max)*40; $i++){
      print "*";
    }
    for(my $j=$i;$j<40;$j++){
      print " ";
    }
    print "  ".$f->density_value()."\t".$f->start()."\n";
  }
}




  


<commit_msg>Now delete old features before starting to avoid creating duplicates.<commit_after>#
# script to calculate the gene density features on a database
# should work on any species database
#

#
# It will only run on databases with genes ...
#

# I think the right thing here is too generate densities on the biggest 50-100 toplevel slices ...
# The website will be happy with about 150 bins I think.


use strict;

use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::DensityType;
use Bio::EnsEMBL::DensityFeature;
use Getopt::Long;

my $bin_count = 150;
my $long_slice_count = 100;


my ( $host, $user, $pass, $port, $dbname  );

$port = 3306 ; 

my ( $block_count, $genome_size, $block_size );

GetOptions( "host=s", \$host,
	    "user=s", \$user,
	    "pass=s", \$pass,
	    "port=i", \$port,
	    "dbname=s", \$dbname
	  );



unless ($host || $user || $pass || $dbname) {
	print "\n\nusage : perl gene_density.pl -host <HOST> -user <USER> -pass <<PASSWORD>> -port <3306> -dbname <DATABASENAME> \n\n" ;
	exit(0) ;
	}

my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(-host => $host,
					    -user => $user,
					    -port => $port,
					    -pass => $pass,
					    -dbname => $dbname);

my $sth = $db->dbc()->prepare( "select count(*) from gene" );
$sth->execute();

my ( $gene_count )  = $sth->fetchrow_array();

if( ! $gene_count ) {
  print STDERR "No gene density for $dbname.\n";
  exit();
}

#
# Could be database without seq_regions
#  Then have to try and attach core db
#
$sth = $db->dbc()->prepare( "select count(*)  from seq_region" );
$sth->execute();
my ( $seq_region_count ) = $sth->fetchrow_array();
if( ! $seq_region_count ) {
  #
  # for the time being only do core dbs
  # no dbs with no seq regions
  #
  print STDERR "No gene density for $dbname, no seq_regions.\n";
  exit();

  if( ($dbname =~ /_est_/ ) || ( $dbname =~ /_vega_/ ) || ( $dbname =~ /_cdna_/ ) ) {
    my $dna_db_name = $dbname;
    $dna_db_name =~ s/(_estgene_|_vega_|_cdna_)/_core_/;
    my $dna_db =  new Bio::EnsEMBL::DBSQL::DBAdaptor
      (-host => $host,
       -user => $user,
       -port => $port,
       -pass => $pass,
       -dbname => $dna_db_name );
    print STDERR "Attaching $dna_db_name to $dbname.\n";
    $db->dnadb( $dna_db );
  } else {
    print STDERR "No gene density for $dbname, no seq_regions.\n";
    exit();
  }
}

#
# Get the adaptors needed;
#

my $dfa = $db->get_DensityFeatureAdaptor();
my $dta = $db->get_DensityTypeAdaptor();
my $aa  = $db->get_AnalysisAdaptor();
my $slice_adaptor = $db->get_SliceAdaptor();

print "Deleting old knownGeneDensity and geneDensity features\n";
$sth = $db->dbc->prepare("DELETE df, dt, a FROM density_feature df, density_type dt, analysis a WHERE a.analysis_id=dt.analysis_id AND dt.density_type_id=df.density_type_id AND a.logic_name IN ('knownGeneDensity', 'geneDensity')");
$sth->execute();

#
# block size estimation
#

my $top_slices = $slice_adaptor->fetch_all('toplevel');
my @sorted_slices = sort { $b->seq_region_length() <=> $a->seq_region_length()} @$top_slices;

	
my $analysis = new Bio::EnsEMBL::Analysis (-program     => "gene_density_calc.pl",
					   -database    => "ensembl",
					   -gff_source  => "gene_density_calc.pl",
					   -gff_feature => "density",
					   -logic_name  => "knownGeneDensity");

$aa->store( $analysis );

$analysis = new Bio::EnsEMBL::Analysis (-program     => "gene_density_calc.pl",
					-database    => "ensembl",
					-gff_source  => "gene_density_calc.pl",
					-gff_feature => "density",
					-logic_name  => "geneDensity");

$aa->store( $analysis );

#
# Now the actual feature calculation loop
#


foreach my $known (1, 0) {
  #
  # Create new analysis object for density calculation.
  #

  if($known) {
    $analysis = $aa->fetch_by_logic_name('knownGeneDensity');
  } else {
    $analysis = $aa->fetch_by_logic_name('geneDensity');
  }

  #
  # Create new density type.
  #

  my $dt = Bio::EnsEMBL::DensityType->new(-analysis   => $analysis,
					  -region_features => $bin_count,
					  -value_type => 'sum');

  $dta->store($dt);

  my $slice_count = 0;
  my ( $current, $current_start, $current_end  );
  
  foreach my $slice (@sorted_slices){

    $block_size = $slice->length() / $bin_count;

    my @density_features=();

    print "Gene densities for ".$slice->seq_region_name().
      " with block size $block_size\n";
    $current_end = 0;
    $current = 0;

    while($current_end < $slice->end()) {
      $current += $block_size;
      $current_start = $current_end+1;
      $current_end = int( $current + 1 );

      if( $current_end < $current_start ) { 
	$current_end = $current_start;
      }

      if( $current_end > $slice->end() ) {
	$current_end = $slice->end();
      }


      my $sub_slice = $slice->sub_Slice( $current_start , $current_end );

      my $count =0;

      #
      # Store info for genes (ignore pseudo genes)
      #

      foreach my $gene (@{$sub_slice->get_all_Genes()}){
	if($gene->biotype() !~ /pseudogene/i and $gene->start >=1 ) {
	  $count++ if(!$known || $gene->is_known());
	}
      }

      push @density_features, Bio::EnsEMBL::DensityFeature->new
	(-seq_region    => $slice,
	 -start         => $current_start,
	 -end           => $current_end,
	 -density_type  => $dt,
	 -density_value => $count);
    }

    $dfa->store(@density_features);
    print "Created ", scalar @density_features, " gene density features.\n";

    last if( $slice_count++ > $long_slice_count );
  }
}







#
# helper to draw an ascii representation of the density features
#
sub print_features {
  my $features = shift;

  return if(!@$features);

  my $sum = 0;
  my $length = 0;
#  my $type = $features->[0]->{'density_type'}->value_type();

  print("\n");
  my $max=0;
  foreach my $f (@$features) {
    if($f->density_value() > $max){
      $max=$f->density_value();
    }
  }
  if( !$max ) { $max = 1 };

  foreach my $f (@$features) {
    my $i=1;
    for(; $i< ($f->density_value()/$max)*40; $i++){
      print "*";
    }
    for(my $j=$i;$j<40;$j++){
      print " ";
    }
    print "  ".$f->density_value()."\t".$f->start()."\n";
  }
}




  



{'fraction_non_alphanumeric': 0.1773565747165341, 'fraction_numerical': 0.0071078016584870535, 'mean_word_length': 3.044832306639288, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 168, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2488371', 'n_tokens_mistral': 4756, 'n_tokens_neox': 4361, 'n_words': 1024}