From 1339de5fc088d78e44a982f4a90ce5d3353bb86b Mon Sep 17 00:00:00 2001
From: Magali Ruffier <mr6@ebi.ac.uk>
Date: Thu, 20 Jun 2013 10:49:40 +0000
Subject: [PATCH] density scripts moved to
ensembl-production/modules/Bio/EnsEMBL/Production/Pipeline/Production
removing ensembl-core versions
---
misc-scripts/density_feature/DensityCalc.sh | 146 -----
.../density_feature/densities_multi_db.pl | 119 ----
.../density_feature/density_feature_fetch.pl | 97 ----
.../density_feature/gene_density_calc.pl | 379 -------------
.../density_feature/percent_gc_calc.pl | 292 ----------
.../density_feature/repeat_coverage_calc.pl | 394 --------------
.../density_feature/seq_region_stats.pl | 346 ------------
.../submit_density_features.pl | 510 ------------------
.../density_feature/variation_density.pl | 233 --------
9 files changed, 2516 deletions(-)
delete mode 100755 misc-scripts/density_feature/DensityCalc.sh
delete mode 100644 misc-scripts/density_feature/densities_multi_db.pl
delete mode 100644 misc-scripts/density_feature/density_feature_fetch.pl
delete mode 100644 misc-scripts/density_feature/gene_density_calc.pl
delete mode 100644 misc-scripts/density_feature/percent_gc_calc.pl
delete mode 100644 misc-scripts/density_feature/repeat_coverage_calc.pl
delete mode 100644 misc-scripts/density_feature/seq_region_stats.pl
delete mode 100644 misc-scripts/density_feature/submit_density_features.pl
delete mode 100644 misc-scripts/density_feature/variation_density.pl
diff --git a/misc-scripts/density_feature/DensityCalc.sh b/misc-scripts/density_feature/DensityCalc.sh
deleted file mode 100755
index 612cfe8ea6..0000000000
--- a/misc-scripts/density_feature/DensityCalc.sh
+++ /dev/null
@@ -1,146 +0,0 @@
-#!/bin/ksh
-
-# This script submits jobs to the farm to calculate the various density
-# features for a particular core database.
-
-# If an output dir is not specified, output from bjobs will be written to PWD.
-
-# Default values for command line switches:
-
-host='ens-staging' # -h
-port='3306' # -P
-user='ensadmin' # -u
-pass= # -p
-dbname= # -d
-species= # -s
-outdir= # -o
-
-
-while getopts 'h:P:d:s:u:p:o:' opt; do
- case ${opt} in
- h) host=${OPTARG} ;;
- P) port=${OPTARG} ;;
- d) dbname=${OPTARG} ;;
- s) species=${OPTARG};;
- u) user=${OPTARG} ;;
- p) pass=${OPTARG} ;;
- o) outdir=${OPTARG} ;;
- esac
-done
-
-if [[ -z ${host} || -z ${port} || -z ${dbname} || -z ${species} || -z ${user} || -z ${pass} ]]
-then
- print -u2 "Usage:\n\t$0 -h host -P port -d database -s species -u user -p password"
- exit 1
-fi
-
-# Enter an output dir so results don't get written in the same place as your script
-if [[ -z ${outdir} || ! -d ${outdir} ]]
-then
- outdir=$PWD
-fi
-
-print "Output dir is '${outdir}'"
-
-# Make sure this is a core database.
-if [[ -n ${dbname##*_core_*} ]]; then
- print -u2 "The database '${dbname}' is not a core database"
- exit 1
-fi
-
-print "Submitting percent GC calculation to queue 'normal'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_gc.out'"
-bsub -q normal -J gc_calc \
- -oo ${outdir}/${dbname}_gc.out \
- -eo ${outdir}/${dbname}_gc.err \
- perl ./percent_gc_calc.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -d ${dbname}
-
-print "Submitting gene density calculation to queue 'normal'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_gene.out'"
-bsub -q normal -J gene_density \
- -oo ${outdir}/${dbname}_gene.out \
- -eo ${outdir}/${dbname}_gene.err \
- perl ./gene_density_calc.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -d ${dbname}
-
-print "Submitting repeat coverage calculation to queue 'long'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_repeat.out'"
-bsub -q long -J repeat_cov \
- -oo ${outdir}/${dbname}_repeat.out \
- -eo ${outdir}/${dbname}_repeat.err \
- perl ./repeat_coverage_calc.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -d ${dbname}
-
-print "Submitting variation density calculation to queue 'normal'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_var.out'"
-bsub -q normal -J var_density \
- -oo ${outdir}/${dbname}_var.out \
- -eo ${outdir}/${dbname}_var.err \
- perl ./variation_density.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -s ${species}
-
-print "Submitting seq region gene stats calculation to queue 'normal'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_seqreg.out'"
-bsub -q normal -J seqreg_stats_gene \
- -oo ${outdir}/${dbname}_seqreg_gene.out \
- -eo ${outdir}/${dbname}_seqreg_gene.err \
- perl ./seq_region_stats.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -d ${dbname} \
- -s gene
-
-print "Submitting seq region snp stats calculation to queue 'normal'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_seqreg.out'"
-bsub -q normal -J seqreg_stats_snp \
- -oo ${outdir}/${dbname}_seqreg_snp.out \
- -eo ${outdir}/${dbname}_seqreg_snp.err \
- perl ./seq_region_stats.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -d ${dbname} \
- -s snp
-
-print "Submitting gene_gc content calculation to queue 'normal'"
-print "\tThe output from this job goes to the file"
-print "\t'${dbname}_genegc.out'"
-bsub -q normal -J genegc_stats \
- -oo ${outdir}/${dbname}_genegc.out \
- -eo ${outdir}/${dbname}_genegc.err \
- perl ../gene_gc.pl \
- -h ${host} \
- -port ${port} \
- -u ${user} \
- -p ${pass} \
- -pattern ${dbname}
-
-print "All jobs submitted."
-
-# $Id$
diff --git a/misc-scripts/density_feature/densities_multi_db.pl b/misc-scripts/density_feature/densities_multi_db.pl
deleted file mode 100644
index 12cae2b7bb..0000000000
--- a/misc-scripts/density_feature/densities_multi_db.pl
+++ /dev/null
@@ -1,119 +0,0 @@
-#
-# These script runs all density creation scripts on all
-# databases that match input pattern. It calls all configured
-# scripts of this directory with host, user,pass,port and dbname
-# parameter.
-#
-use strict;
-use Getopt::Long;
-use DBI;
-
-#
-# configure this hash and command line parameter section for
-# adding new scripts
-#
-
-my %density_scripts =
- ( "gene" => "gene_density_calc.pl",
- "gc" => "percent_gc_calc.pl",
- "repeat" => "repeat_coverage_calc.pl",
- "stats" => "seq_region_stats.pl",
- "snp" => "snp_lite_density.pl" );
-
-my @density_scripts = ();
-
-my ( $host, $user, $pass, $port, $dbpattern, $snpstats, $genestats, $stats );
-
-GetOptions( "host=s", \$host,
- "user=s", \$user,
- "pass=s", \$pass,
- "port=i", \$port,
- "genestats", \$genestats,
- "snpstats", \$snpstats,
- "gene", sub { push( @density_scripts, $density_scripts{'gene'} ); },
- "gc", sub { push( @density_scripts, $density_scripts{'gc'} ); },
- "repeat", sub { push( @density_scripts, $density_scripts{'repeat'} ); },
- "stats", \$stats,
- "snp", sub { push( @density_scripts, $density_scripts{'snp'} ); },
- "dbpattern=s", \$dbpattern
- );
-
-if( !$host || !$user ) {
- usage();
-}
-
-
-if($stats) {
- if($genestats || $snpstats) {
- print STDERR "Only specify one of -stats, -snpstats or -genestats\m";
- usage();
- }
- push @density_scripts, $density_scripts{'stats'};
-} elsif($genestats || $snpstats) {
- my $arg = '';
- $arg .= ' -genestats ' if($genestats);
- $arg .= ' -snpstats ' if($snpstats);
- my $line = $density_scripts{'stats'} . $arg;
- push @density_scripts, $line;
-}
-
-
-my $dsn = "DBI:mysql:host=$host";
-my $args = "-host $host ";
-if( $user ) {
- $args .= "-user $user ";
-}
-
-if( $port ) {
- $args .= "-port $port ";
- $dsn .= ";port=$port";
-}
-
-if( $pass ) {
- $args .= "-pass $pass ";
-}
-
-if( $host ) {
- $args .= "-host $host ";
-}
-
-
-my $db = DBI->connect( $dsn, $user, $pass );
-
-my @dbnames = map {$_->[0] } @{ $db->selectall_arrayref( "show databases" ) };
-
-for my $dbname ( @dbnames ) {
- if( $dbpattern ) {
- if( $dbname !~ /$dbpattern/ ) {
- next;
- }
- } else {
- print STDERR "$dbname\n";
- next;
- }
-
- if( ! @density_scripts ) {
- print "Matches $dbname.\n";
- } else {
- for my $scriptname ( @density_scripts ) {
- print( "/usr/local/ensembl/bin/perl $scriptname $args -dbname $dbname\n" );
- system( "/usr/local/ensembl/bin/perl $scriptname $args -dbname $dbname" );
- }
- }
-}
-
-
-sub usage {
- print STDERR "densities_multi_db.pl runs \n ";
- print STDERR <<EOF;
-
- with parameters -host [hostname] -user [username with write permit]
- -pass [password] -port [portnumber]
- on all databases that match -dbpattern [some_db_regexp]
- on -gene, -repeat, -snp, -stats/-genestats/-snpstats, -gc.
-
- Without -dbpattern it lists all available databases and does nothing.
-
-EOF
- exit();
-}
diff --git a/misc-scripts/density_feature/density_feature_fetch.pl b/misc-scripts/density_feature/density_feature_fetch.pl
deleted file mode 100644
index c0871384ad..0000000000
--- a/misc-scripts/density_feature/density_feature_fetch.pl
+++ /dev/null
@@ -1,97 +0,0 @@
-#!/usr/local/ensembl/bin/perl -w
-
-#
-# Ignore file anme this should calculate the gene count for a chromosome (X)
-#
-
-
-use strict;
-use lib '/tmp_mnt/nfs6/vol_vol1_homes/ianl/CVS/ensembl/modules';
-use lib '/tmp_mnt/nfs6/vol_vol1_homes/ianl/read_only/bioperl-live';
-
-
-use Data::Dumper;
-
-#use dbi;
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
-use Bio::EnsEMBL::DBSQL::SliceAdaptor;
-use Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor;
-use Bio::EnsEMBL::DBSQL::DensityTypeAdaptor;
-use Bio::EnsEMBL::Slice;
-use Bio::EnsEMBL::Analysis;
-use Bio::EnsEMBL::DensityType;
-use Bio::EnsEMBL::DensityFeature;
-print "hello\n";
-
-my $host = '127.0.0.1';
-my $user = 'ensro';
-#my $pass = 'ensembl';
-my $dbname = 'homo_sapiens_core_20_34';
-my $port = '5050';
-
-my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(-host => $host,
- -user => $user,
- -port => $port,
-# -pass => $pass,
- -dbname => $dbname);
-
-
-my $slice_adaptor = $db->get_SliceAdaptor();
-my $dfa = $db->get_DensityFeatureAdaptor();
-
-
-my $num_of_blocks = 20;
-my $interpolate = 1;
-
-
-#comment out as appropriate.
-#my $type="GeneDensity";
-#my $type="PercentGC";
-my $type ="PercentageRepeat";
-
-foreach my $chrom (qw(X Y)){
- if($interpolate){
- print "getting density feature for chromosome $chrom with $num_of_blocks blocks\n";
- }
- else{
- print "getting ALL density feature for chromosome $chrom\n";
- }
- my $slice = $slice_adaptor->fetch_by_region('chromosome',$chrom);
-
-
- my @features = @{$dfa->fetch_all_by_Slice($slice,$type, $num_of_blocks, $interpolate)};
-
- print_features(\@features);
-}
-
-
-#
-# 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();
- }
- }
- 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 " ".int($f->density_value())."\t".$f->start()."\t".$f->end()."\n";
- }
-}
diff --git a/misc-scripts/density_feature/gene_density_calc.pl b/misc-scripts/density_feature/gene_density_calc.pl
deleted file mode 100644
index 3994db3bb7..0000000000
--- a/misc-scripts/density_feature/gene_density_calc.pl
+++ /dev/null
@@ -1,379 +0,0 @@
-#!/usr/bin/env perl
-
-#
-# 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 to generate densities on the longest
-# 125 toplevel slices... The website will be happy with about 150 bins I
-# think.
-
-use strict;
-use warnings;
-
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
-use Bio::EnsEMBL::DensityType;
-use Bio::EnsEMBL::DensityFeature;
-use Bio::EnsEMBL::Utils::ConversionSupport;
-use Bio::EnsEMBL::Utils::CliHelper;
-
-my $bin_count = 150;
-my $max_slices = 100;
-
-my $cli_helper = Bio::EnsEMBL::Utils::CliHelper->new();
-my $optsd =
- [ @{ $cli_helper->get_dba_opts() }, @{ $cli_helper->get_dba_opts('m') } ];
-# add the print option
-push(@{$optsd},"print");
-my $opts = $cli_helper->process_args( $optsd, \&usage );
-
-if(!defined $opts->{mdbname} && !defined $opts->{mdbpattern}) {
- $opts->{mdbname} = 'ensembl_production';
-}
-
-if (! defined $opts->{mhost} ) {
- ( $opts->{mhost}, $opts->{mport}, $opts->{mdbname}, $opts->{muser} ) =
- ( 'ens-staging1', '3306', 'ensembl_production', 'ensro' );
-}
-
-my ($prod_dba) = @{ $cli_helper->get_dbas_for_opts( $opts, 1, 'm' ) };
-
-if ( !defined $prod_dba ) {
- usage();
-}
-
-my ( $block_count, $genome_size, $block_size );
-for my $db_args ( @{ $cli_helper->get_dba_args_for_opts($opts) } ) {
- print STDOUT "Connecting to " . $db_args->{-DBNAME} . "\n";
-
- my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(%$db_args);
- my $dbname = $db->dbc()->dbname();
-
- my $helper = $db->dbc()->sql_helper();
-
- my $gene_count =
- $helper->execute_single_result(
--SQL=>"select count(*) from gene join seq_region using (seq_region_id) join coord_system using (coord_system_id) where species_id=?",
- -PARAMS=>[$db->species_id()] );
-
- 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
- #
- my $seq_region_count =
- $helper->execute_single_result(
--SQL=>"select count(*) from seq_region join coord_system using (coord_system_id) where species_id=?",
- -PARAMS=>[$db->species_id()] );
- 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_/;
- $db_args->{dbname} = $dna_db_name;
- my $dna_db = new Bio::EnsEMBL::DBSQL::DBAdaptor($db_args);
- print STDOUT "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;
- #
-
- print STDOUT "Deleting old knownGeneDensity and geneDensity features\n";
-
- $helper->execute_update(
- -SQL=>qq/DELETE df, dt FROM density_feature df, density_type dt, analysis a,
- seq_region s, coord_system cs WHERE
- df.seq_region_id=s.seq_region_id AND s.coord_system_id=cs.coord_system_id AND cs.species_id=?
- AND a.analysis_id=dt.analysis_id AND dt.density_type_id=df.density_type_id
- AND a.logic_name IN ('knowngenedensity', 'genedensity')/, -PARAMS=>[$db->species_id()] );
-
- my $dfa = $db->get_DensityFeatureAdaptor();
- my $dta = $db->get_DensityTypeAdaptor();
- my $aa = $db->get_AnalysisAdaptor();
- my $slice_adaptor = $db->get_SliceAdaptor();
-
- my $support = 'Bio::EnsEMBL::Utils::ConversionSupport';
- my $analysis1 = $aa->fetch_by_logic_name('knowngenedensity');
- my $analysis2 = $aa->fetch_by_logic_name('genedensity');
-
- if ( !defined($analysis1) ) {
-
- my ( $display_label, $description ) =
- @{$prod_dba->dbc()->sql_helper()->execute(
--SQL=>"select distinct display_label, description from analysis_description where is_current = 1 and logic_name = 'knowngenedensity'"
- ) };
-
- $analysis1 =
- new Bio::EnsEMBL::Analysis( -program => "gene_density_calc.pl",
- -database => "ensembl",
- -gff_source => "gene_density_calc.pl",
- -gff_feature => "density",
- -logic_name => "knowngenedensity",
- -description => $description,
- -display_label => $display_label,
- -displayable => 1 );
-
- $aa->store($analysis1);
-
- } else {
- $analysis1->created( $support->date() );
- $aa->update($analysis1);
- }
-
- if ( !defined($analysis2) ) {
-
- my ( $display_label, $description ) =
- @{$prod_dba->dbc()->sql_helper()->execute(
--SQL=>"select distinct display_label, description from analysis_description where is_current = 1 and logic_name = 'genedensity'"
- ) };
-
- $analysis2 =
- new Bio::EnsEMBL::Analysis( -program => "gene_density_calc.pl",
- -database => "ensembl",
- -gff_source => "gene_density_calc.pl",
- -gff_feature => "density",
- -logic_name => "genedensity",
- -description => $description,
- -display_label => $display_label,
- -displayable => 1 );
-
- $aa->store($analysis2);
-
- } else {
- $analysis2->created( $support->date() );
- $aa->update($analysis2);
- }
-
- #
- # Now the actual feature calculation loop
- #
- my $total_count = 0;
-
- foreach my $known ( 1, 0 ) {
- #
- # Create new analysis object for density calculation.
- #
- my $analysis;
-
- if ($known) {
- $analysis = $aa->fetch_by_logic_name('knowngenedensity');
- } else {
- $analysis = $aa->fetch_by_logic_name('genedensity');
- }
-
- # Sort slices by coordinate system rank, then by length.
- my @sorted_slices =
- sort( { $a->coord_system()->rank() <=> $b->coord_system()->rank()
- || $b->seq_region_length() <=> $a->seq_region_length()
- } @{ $slice_adaptor->fetch_all('toplevel') } );
-
- #
- # 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 );
-
- while ( my $slice = shift @sorted_slices ) {
-
- $block_size = $slice->length()/$bin_count;
-
- my @density_features; #sf7
-
- print STDOUT "Gene densities for "
- . $slice->seq_region_name()
- . " with block size $block_size\n";
- $current_end = 0;
- $current = 0;
-
- while ( $current_end < $slice->length ) {
- $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 );
-
- if ($count > 0) {
- #density features with value = 0 are not stored
- $total_count++;
-}
- } ## end while ( $current_end < $slice...)
-
- $dfa->store(@density_features);
-
- last if ( $slice_count++ > $max_slices );
- } ## end while ( my $slice = shift...)
-
- } ## end foreach my $known ( 1, 0 )
- print STDOUT "Created $total_count gene density features.\n";
- print STDOUT "Finished with $dbname";
-} ## end for my $db_args ( @{ $cli_helper...})
-
-#
-# 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";
- }
-} ## end sub print_features
-
-sub usage {
- my $indent = ' ' x length($0);
- print <<EOF; exit(0);
-
-
-What does it do?
-
-Populates known gene density features in a database as well as gene density
-features from genes of any status.
-
-First it needs gene and seq_region tables to be populated. It then
-deletes all knownGeneDensity and geneDensity entries from the
-analysis, density_type and density_feature tables. All toplevel
-slices are fetched and sorted from longest to shortest. Each slice
-is divided into 150 bins (blocks). For each of the blocks or
-sub-slices, the number of genes is counted to give geneDensity and
-the number of genes with status KNOWN status is counted to give
-knownGeneDensity. All biotypes except pseudogene are counted.
-
-Input data: genes, top level seq regions, xrefs
-Output tables: analysis (logic_name: knownGeneDensity and geneDensity),
- analysis description, density_type, density_feature
-
-
-When to run it in the release cycle?
-
-It can be run after compara have handed over homologies and core have
-done xref projections.
-
-
-Which databases to run it on?
-
-It needs to be run across all core databases for every release.
-
-
-How long does it take?
-
-It takes about 10 mins to run for a database in normal queue,
-
-
-Usage:
-
- $0 -h host [-port port] -u user -p password \\
- $indent -d database | -pattern pattern \\
- $indent [-mhost ensembl_production host] [-mport ensembl_production host] [-muser ensembl_production host] \\
- $indent [-help] \\
-
-
- -h|host Database host to connect to
-
- -port Database port to connect to (default 3306)
-
- -u|user Database username
-
- -p|pass Password for user
-
- -d|dbname Database name
-
-
- -mhost ensembl_production database host to connect to
-
- -mport ensembl_production database port to connect to
-
- -muser ensembl_production database username
-
-
- -pattern Database name regexp
-
- -help This message
-
-
-EOF
-} ## end sub usage
-
diff --git a/misc-scripts/density_feature/percent_gc_calc.pl b/misc-scripts/density_feature/percent_gc_calc.pl
deleted file mode 100644
index b306c4f78b..0000000000
--- a/misc-scripts/density_feature/percent_gc_calc.pl
+++ /dev/null
@@ -1,292 +0,0 @@
-#!/usr/bin/env perl
-
-#
-# Calculate the GC content for top level seq_regions
-# small regions 500bp to be able to display on contigview
-# big regions genomesize / 4000 for 4000 features on the genome
-
-use strict;
-use warnings;
-#use dbi;
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
-use Bio::EnsEMBL::DBSQL::SliceAdaptor;
-use Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor;
-use Bio::EnsEMBL::DBSQL::DensityTypeAdaptor;
-use Bio::EnsEMBL::Slice;
-use Bio::EnsEMBL::Analysis;
-use Bio::EnsEMBL::DensityType;
-use Bio::EnsEMBL::DensityFeature;
-use Bio::EnsEMBL::Utils::ConversionSupport;
-use Bio::EnsEMBL::Utils::CliHelper;
-use Data::Dumper;
-
-my $cli_helper = Bio::EnsEMBL::Utils::CliHelper->new();
-my $optsd =
- [ @{ $cli_helper->get_dba_opts() }, @{ $cli_helper->get_dba_opts('m') } ];
-push( @{$optsd}, "print" );
-my $opts = $cli_helper->process_args( $optsd, \&usage );
-if ( !defined $opts->{mdbname} && !defined $opts->{mdbpattern} ) {
- $opts->{mdbname} = 'ensembl_production';
-}
-
-if ( !defined $opts->{mhost} ) {
- ( $opts->{mhost}, $opts->{mport}, $opts->{mdbname}, $opts->{muser} ) =
- ( 'ens-staging1', '3306', 'ensembl_production', 'ensro' );
-}
-
-my ($prod_dba) = @{ $cli_helper->get_dbas_for_opts( $opts, 1, 'm' ) };
-
-if ( !defined $prod_dba ) {
- usage();
-}
-
-for my $db_args ( @{ $cli_helper->get_dba_args_for_opts($opts) } ) {
- print STDOUT "Connecting to " . $db_args->{-DBNAME} . "\n";
-
- my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(%$db_args);
- my $dbname = $db->dbc()->dbname();
-
- my $helper = $db->dbc()->sql_helper();
-
- my $bin_count = 150;
- my $max_slices = 100;
-
- #
- # Check wether the script should run on given database
- #
- my $dna_count = $helper->execute_single_result(
- -SQL => q/select count(*) from dna join seq_region
-using (seq_region_id) join coord_system using (coord_system_id) where species_id=?/,
- -PARAMS => [ $db->species_id() ] );
-
- if ( !$dna_count ) {
- print STDERR "No dna, no gc content for $dbname.\n";
- exit();
- }
- print STDOUT "Deleting old PercentGC features\n";
- $helper->execute_update(
- -SQL => qq(
-DELETE df, dt
-FROM density_feature df, density_type dt, analysis a,
-seq_region sr, coord_system c
-WHERE
-sr.seq_region_id=df.seq_region_id
-AND sr.coord_system_id=c.coord_system_id
-AND c.species_id=?
-AND a.analysis_id=dt.analysis_id
-AND dt.density_type_id=df.density_type_id
-AND a.logic_name='percentgc'),
- -PARAMS => [ $db->species_id() ] );
-
- #
- # Get the adaptors needed;
- #
-
- my $slice_adaptor = $db->get_SliceAdaptor();
- my $dfa = $db->get_DensityFeatureAdaptor();
- my $dta = $db->get_DensityTypeAdaptor();
- my $aa = $db->get_AnalysisAdaptor();
-
- # Sort slices by coordinate system rank, then by length.
- my @sorted_slices =
- sort( { $a->coord_system()->rank() <=> $b->coord_system()->rank()
- || $b->seq_region_length() <=> $a->seq_region_length()
- } @{ $slice_adaptor->fetch_all('toplevel') } );
-
- my $analysis = $aa->fetch_by_logic_name('percentgc');
-
- if ( !defined($analysis) ) {
-
- my ( $display_label, $description ) =
- @{$prod_dba->dbc()->sql_helper()->execute(
- -SQL =>
-"select distinct display_label, description from analysis_description where is_current = 1 and logic_name = 'percentgc'"
- ) };
-
- $analysis =
- new Bio::EnsEMBL::Analysis( -program => "percent_gc_calc.pl",
- -database => "ensembl",
- -gff_source => "percent_gc_calc.pl",
- -gff_feature => "density",
- -logic_name => "percentgc",
- -description => $description,
- -display_label => $display_label,
- -displayable => 1 );
-
- $aa->store($analysis);
-
- } else {
-
- my $support = 'Bio::EnsEMBL::Utils::ConversionSupport';
- $analysis->created( $support->date() );
- $aa->update($analysis);
-
- }
-
- #
- # Create new density type.
- #
-
- my $density_type =
- Bio::EnsEMBL::DensityType->new( -analysis => $analysis,
- -region_features => $bin_count,
- -value_type => 'ratio' );
-
- $dta->store($density_type);
-
- my ( $current_start, $current_end, $current );
- my $slice_count = 0;
- my $block_size;
- my $total_count;
- while ( my $slice = shift @sorted_slices ) {
-
- $block_size = $slice->length()/$bin_count;
-
- my @density_features = ();
-
- print STDOUT "GC percentage for "
- . $slice->seq_region_name()
- . " with block size $block_size\n";
-
- $current_end = 0;
- $current = 0;
-
- while ( $current_end < $slice->length ) {
-
- $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 $gc = $sub_slice->get_base_count()->{'%gc'};
- my $df =
- Bio::EnsEMBL::DensityFeature->new( -seq_region => $slice,
- -start => $current_start,
- -end => $current_end,
- -density_type => $density_type,
- -density_value => $gc );
-#print join ("\t", $slice, $current_start, $current_end, $density_type, $gc, "\n");
- push( @density_features, $df );
- if ( $gc > 0 ) {
- $dfa->store($df);
- #density features with value = 0 are not stored
- $total_count++;
- }
- } ## end while ( $current_end < $slice...)
-
- last if ( $slice_count++ > $max_slices );
- } ## end while ( my $slice = shift...)
-
- print STDOUT "Created $total_count percent gc density features\n";
-} ## end for my $db_args ( @{ $cli_helper...})
-
-#
-# 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();
- }
- }
- 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";
- }
- # my $avg = undef;
- # $avg = $sum/$length if($length < 0);
- # print("Sum=$sum, Length=$length, Avg/Base=$sum");
-} ## end sub print_features
-
-sub usage {
- my $indent = ' ' x length($0);
- print <<EOF; exit(0);
-
-
-What does it do?
-
-First, it needs the dna table to be populated. It then deletes
-all PercentGC entries from the analysis, density_type and
-density_feature tables. All toplevel slices are fetched and sorted
-from longest to shortest. Each slice is divided into 150 bins
-(blocks). For each of the blocks or sub-slices, the %gc is
-calculated.
-
-Input data: dna sequence, top level seq regions
-Output tables: updates analysis creation date,
- density_type, density_feature
-
-
-When to run it in the release cycle?
-
-It can be run after genebuilders have finished their Xrefs
-(script not affected by projected Xrefs).
-
-
-Which databases to run it on?
-
-Run on core databases for new species or if one of the following changed:
- - dna sequence
- - assembly
-
-
-How long does it take?
-
-It takes about 15 mins to run for a database in normal queue,
-
-
-Usage:
-
- $0 -h host [-port port] -u user -p password \\
- $indent -d database \\
- $indent [-mhost ensembl_production host] [-mport ensembl_production host] [-muser ensembl_production host] \\
- $indent [-help] \\
-
- -h|host Database host to connect to
-
- -port Database port to connect to (default 3306)
-
- -u|user Database username
-
- -p|pass Password for user
-
- -pattern Database name regexp
-
- -mhost ensembl_production database host to connect to
-
- -mport ensembl_production database port to connect to
-
- -muser ensembl_production database username
-
- -help This message
-
-
-EOF
-
-} ## end sub usage
-
diff --git a/misc-scripts/density_feature/repeat_coverage_calc.pl b/misc-scripts/density_feature/repeat_coverage_calc.pl
deleted file mode 100644
index 30979c3fe6..0000000000
--- a/misc-scripts/density_feature/repeat_coverage_calc.pl
+++ /dev/null
@@ -1,394 +0,0 @@
-#!/usr/bin/env perl
-#
-# Calculate the repeat coverage for given database.
-# condition: 1k blocks to show contigview displays
-# 4000 blocks for a whole genome
-#
-# checks wether database contains repeats before doing anything
-
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
-use Bio::EnsEMBL::DBSQL::SliceAdaptor;
-use Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor;
-use Bio::EnsEMBL::DBSQL::DensityTypeAdaptor;
-use Bio::EnsEMBL::Slice;
-use Bio::EnsEMBL::Analysis;
-use Bio::EnsEMBL::DensityType;
-use Bio::EnsEMBL::DensityFeature;
-use Bio::EnsEMBL::Mapper::RangeRegistry;
-use Bio::EnsEMBL::Utils::Exception qw(warning throw);
-use Bio::EnsEMBL::Utils::ConversionSupport;
-
-use POSIX;
-
-use strict;
-use warnings;
-use Getopt::Long;
-use Bio::EnsEMBL::Utils::CliHelper;
-
-my $cli_helper = Bio::EnsEMBL::Utils::CliHelper->new();
-my $optsd =
- [ @{ $cli_helper->get_dba_opts() }, @{ $cli_helper->get_dba_opts('m') } ];
-# add the print option
-push(@{$optsd},"print");
-my $opts = $cli_helper->process_args( $optsd, \&usage );
-if(!defined $opts->{mdbname} && !defined $opts->{mdbpattern}) {
- $opts->{mdbname} = 'ensembl_production';
-}
-
-if (! defined $opts->{mhost} ) {
- ( $opts->{mhost}, $opts->{mport}, $opts->{mdbname}, $opts->{muser} ) =
- ( 'ens-staging1', '3306', 'ensembl_production', 'ensro' );
-}
-
-my ($prod_dba) = @{ $cli_helper->get_dbas_for_opts( $opts, 1, 'm' ) };
-
-if ( !defined $prod_dba ) {
- usage();
-}
-
-my $chunksize = 1_000_000;
-my $small_blocksize = 1_000;
-my $bin_count = 150;
-my $max_top_slice = 100;
-
-for my $db_args ( @{ $cli_helper->get_dba_args_for_opts($opts) } ) {
- print STDOUT "Connecting to " . $db_args->{-DBNAME} . "\n";
-
- my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(%$db_args);
- my $dbname = $db->dbc()->dbname();
-
- my $helper = $db->dbc()->sql_helper();
-
- my $repeat_count = $helper->execute_simple(
- -SQL=>q/select count(*) from repeat_feature join seq_region
-using (seq_region_id) join coord_system using (coord_system_id) where species_id=?/,
- -PARAMS=>[$db->species_id()] );
-
- if ( !$repeat_count ) {
- print STDERR "No repeat density for $dbname.\n";
- exit();
- }
-
- #
- # Get the adaptors needed;
- #
-
-#
-# Clean up old features first. Also remove analysis and density type entry as these are recreated.
-#
-
- print STDOUT "Deleting old PercentageRepeat features\n";
- $helper->execute_update(
--SQL=>q/DELETE df, dt FROM density_feature df, density_type dt, analysis a, seq_region sr, coord_system c
-WHERE
-sr.seq_region_id=df.seq_region_id AND c.coord_system_id=sr.coord_system_id AND c.species_id=?
-AND a.analysis_id=dt.analysis_id AND dt.density_type_id=df.density_type_id AND a.logic_name='percentagerepeat'/,
- -PARAMS=>[$db->species_id()] );
-
- my $slice_adaptor = $db->get_SliceAdaptor();
- my $dfa = $db->get_DensityFeatureAdaptor();
- my $dta = $db->get_DensityTypeAdaptor();
- my $aa = $db->get_AnalysisAdaptor();
-
- my $analysis = $aa->fetch_by_logic_name('percentagerepeat');
-
- if ( !defined($analysis) ) {
-
- my ( $display_label, $description ) =
- @{$prod_dba->dbc()->sql_helper()->execute(
--SQL=>"select distinct display_label, description from analysis_description where is_current = 1 and logic_name = 'percentagerepeat'"
- ) };
-
- $analysis =
- new Bio::EnsEMBL::Analysis( -program => "repeat_coverage_calc.pl",
- -database => "ensembl",
- -gff_source => "repeat_coverage_calc.pl",
- -gff_feature => "density",
- -logic_name => "percentagerepeat",
- -description => $description,
- -display_label => $display_label,
- -displayable => 1 );
-
- $aa->store($analysis);
- } else {
-
- my $support = 'Bio::EnsEMBL::Utils::ConversionSupport';
- $analysis->created( $support->date() );
- $aa->update($analysis);
-
- }
-
- my $slices = $slice_adaptor->fetch_all("toplevel");
- my @sorted_slices =
- sort { $b->seq_region_length() <=> $a->seq_region_length() } @$slices;
-
- my $small_density_type = Bio::EnsEMBL::DensityType->new(
- -analysis => $analysis,
- -block_size => $small_blocksize,
- -value_type => 'ratio' );
-
- my $variable_density_type =
- Bio::EnsEMBL::DensityType->new( -analysis => $analysis,
- -region_features => $bin_count,
- -value_type => 'ratio' );
-
- $dta->store($small_density_type);
- $dta->store($variable_density_type);
-
- my $slice_count = 0;
-
- while ( my $slice = shift @sorted_slices ) {
-
- #
- # do it for small and large blocks
- #
- print STDOUT ( "Working on seq_region "
- . $slice->seq_region_name()
- . " length "
- . $slice->seq_region_length() );
-
- my $rr = Bio::EnsEMBL::Mapper::RangeRegistry->new();
- my $chunk_end = 0;
- my $variable_end = 0;
- my $small_end = 0;
- my $small_start = 0;
- my $repeat_size;
- my $variable_start = 0;
- my $variable_blocksize =
- POSIX::ceil( $slice->seq_region_length()/
- $variable_density_type->region_features() );
- $slice_count++;
-
- while ( $chunk_end < $slice->length() ) {
- my $chunk_start = $chunk_end + 1;
- $chunk_end += $chunksize;
- $chunk_end = $slice->length() if $chunk_end > $slice->length();
-
- register( $rr, $slice, $chunk_start, $chunk_end );
-
- my @dfs = ();
-
- if ( $slice_count < $max_top_slice ) {
- while ( $variable_end + $variable_blocksize <= $chunk_end ) {
- # here we can do the variable sized repeat densities
- $variable_start = $variable_end + 1;
- $variable_end += $variable_blocksize;
-
- $repeat_size =
- $rr->overlap_size( "1", $variable_start, $variable_end );
- my $percentage_repeat =
- $repeat_size/$variable_blocksize*100;
-
- push( @dfs,
- Bio::EnsEMBL::DensityFeature->new(
- -seq_region => $slice,
- -start => $variable_start,
- -end => $variable_end,
- -density_type => $variable_density_type,
- -density_value => $percentage_repeat )
- );
-
- }
- }
-
- while ( $small_end + $small_blocksize <= $chunk_end ) {
- # here we can do the small sized density features
- $small_start = $small_end + 1;
- $small_end += $small_blocksize;
-
- $repeat_size =
- $rr->overlap_size( "1", $small_start, $small_end );
- my $percentage_repeat = $repeat_size/$small_blocksize*100;
-
- push( @dfs,
- Bio::EnsEMBL::DensityFeature->new(
- -seq_region => $slice,
- -start => $small_start,
- -end => $small_end,
- -density_type => $small_density_type,
- -density_value => $percentage_repeat
- ) );
- }
-
- if (@dfs) {
- $dfa->store(@dfs);
- } else {
- warning( "No repeat density calculated for "
- . $slice->name
- . " (chunk start $chunk_start, chunk end $chunk_end)." );
- }
-
- my $used_lower_limit =
- $small_start < $variable_start ? $small_start : $variable_start;
-
- # here some rr cleanup
- $rr->check_and_register( "1", 0, $used_lower_limit );
- } ## end while ( $chunk_end < $slice...)
-
- # missing the last bits
- if ( $small_end < $slice->length() ) {
- $small_start = $small_end + 1;
- $small_end = $slice->length();
-
- $repeat_size = $rr->overlap_size( "1", $small_start, $small_end );
- my $percentage_repeat =
- $repeat_size/( $small_end - $small_start + 1 )*100;
-
- $dfa->store( Bio::EnsEMBL::DensityFeature->new(
- -seq_region => $slice,
- -start => $small_start,
- -end => $small_end,
- -density_type => $small_density_type,
- -density_value => $percentage_repeat
- ) );
- }
-
- if ( $variable_end < $slice->length() && $slice_count < $max_top_slice )
- {
- $variable_start = $variable_end + 1;
- $variable_end = $slice->length();
-
- $repeat_size =
- $rr->overlap_size( "1", $variable_start, $variable_end );
- my $percentage_repeat =
- $repeat_size/( $variable_end - $variable_start + 1 )*100;
-
- $dfa->store( Bio::EnsEMBL::DensityFeature->new(
- -seq_region => $slice,
- -start => $variable_start,
- -end => $variable_end,
- -density_type => $variable_density_type,
- -density_value => $percentage_repeat )
- );
- }
-
- print STDOUT " DONE.\n";
- } ## end while ( my $slice = shift...)
-} ## end for my $db_args ( @{ $cli_helper...})
-
-sub register {
- my ( $rr, $slice, $start, $end ) = @_;
-
- my $subSlice = $slice->sub_Slice( $start, $end );
- my $repeats = $subSlice->get_all_RepeatFeatures();
- for my $repeat (@$repeats) {
- $rr->check_and_register( "1",
- $repeat->seq_region_start(),
- $repeat->seq_region_end() );
- }
-}
-
-#
-# 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();
- }
- }
- 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";
- }
- # my $avg = undef;
- # $avg = $sum/$length if($length < 0);
- # print("Sum=$sum, Length=$length, Avg/Base=$sum");
-} ## end sub print_features
-
-sub usage {
- my $indent = ' ' x length($0);
- print <<EOF; exit(0);
-
-What does it do?
-
-Calculates the percentage of repetetive elements for top level seq_regions.
-
-First it needs repeat_feature table to be populated. It then
-deletes all PercentageRepeat entries from the analysis, density_type
-and density_feature tables. All toplevel slices are fetched and sorted
-from longest to shortest. For each toplevel slice, both the small and
-the variable density repeats are calculated.
-
-Small repeats are done like this: Move along the toplevel slice 1 KB
-at a time. Find the %repeat for each of these 1 KB blocks.
-
-Variable repeats are done like this: Divide each slice into 150
-sub_slices. Move along the toplevel slice 1 MB at a time. Foreach
-sub_slice within the 1 MB, calculate the %repeat for that sub_slice.
-Variable repeats are only found for the 100 longest toplevel slices.
-
-Input data: repeat features, top level seq regions
-Output tables: updates analysis creation date,
- density_type (two entries, one for small_density type of
- length 1 KB and one for variable_density_type of length 1MB),
- density_feature
-
-
-When to run it in the release cycle?
-
-It can be run after genebuilders have finished their Xrefs
-(script not affected by projected Xrefs).
-
-
-Which databases to run it on?
-
-Run on core databases for new species or if one of the following changed:
- - dna sequence
- - assembly
- - repeats
-
-
-How long does it take?
-
-It takes about 1 hour to run for a database in the long queue. The script is
-slowed down considerably for very fragmented genomes with thousands of toplevel
-seqs.
-
-
-Usage:
- $0 -h host [-port port] -u user -p password \\
- $indent -d database \\
- $indent [-mhost ensembl_production host] [-mport ensembl_production host] [-muser ensembl_production host] \\
- $indent [-help] \\
-
-
- -h|host Database host to connect to
-
- -port Database port to connect to (default 3306)
-
- -u|user Database username
-
- -p|pass Password for user
-
- -mhost ensembl_production database host to connect to
-
- -mport ensembl_production database port to connect to
-
- -muser ensembl_production database username
-
- -pattern Database name regexp
-
- -help This message
-
-
-EOF
-
-} ## end sub usage
-
diff --git a/misc-scripts/density_feature/seq_region_stats.pl b/misc-scripts/density_feature/seq_region_stats.pl
deleted file mode 100644
index 5ec75579bf..0000000000
--- a/misc-scripts/density_feature/seq_region_stats.pl
+++ /dev/null
@@ -1,346 +0,0 @@
-#!/usr/bin/env perl
-
-use strict;
-use warnings;
-
-use Bio::EnsEMBL::Registry;
-use Bio::EnsEMBL::DBSQL::DBAdaptor;
-use Bio::EnsEMBL::Variation::DBSQL::DBAdaptor;
-use Bio::EnsEMBL::Utils::CliHelper;
-use Data::Dumper;
-
-my $cli_helper = Bio::EnsEMBL::Utils::CliHelper->new();
-my $optsd =
- [ @{ $cli_helper->get_dba_opts() }, @{ $cli_helper->get_dba_opts('m') } ];
-# add the print option
-push( @{$optsd}, "stats|s:s" );
-push(@{$optsd},"print");
-my $opts = $cli_helper->process_args( $optsd, \&usage );
-
-usage() if ( !$opts->{stats} || $opts->{stats} !~ /^(gene|snp)$/ );
-if(!defined $opts->{mdbname} && !defined $opts->{mdbpattern}) {
- $opts->{mdbname} = 'ensembl_production';
-}
-
-if ( !defined $opts->{mhost} ) {
- ( $opts->{mhost}, $opts->{mport}, $opts->{mdbname}, $opts->{muser} ) =
- ( 'ens-staging1', '3306', 'ensembl_production', 'ensro' );
-}
-
-my ($prod_dba) = @{ $cli_helper->get_dbas_for_opts( $opts, 1, 'm' ) };
-
-if ( !defined $prod_dba ) {
- usage();
-}
-
-my %attrib_codes = %{
- $prod_dba->dbc()->sql_helper()->execute_into_hash(
- -SQL =>
-"select distinct b.name, code from biotype b join attrib_type using(attrib_type_id) where is_current = 1 and db_type like '%core%' and object_type = 'gene' order by b.name"
- ) };
-
-#add known and novel protein coding attrib types
-$attrib_codes{'known protein_coding'} = 'GeneNo_knwCod';
-$attrib_codes{'novel protein_coding'} = 'GeneNo_novCod';
-
-my $genestats = 1 if ( $opts->{stats} eq 'gene' );
-my $snpstats = 1 if ( $opts->{stats} eq 'snp' );
-for my $db_args ( @{ $cli_helper->get_dba_args_for_opts($opts) } ) {
-
- my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(%$db_args);
- my $dbname = $db->dbc()->dbname();
-
- my $total_count = 0;
- # delete old attributes before starting
- if ($genestats) {
- foreach my $code ( values %attrib_codes ) {
- $db->dbc()->sql_helper()->execute_update(
- -SQL =>
-"DELETE sa FROM seq_region_attrib sa, attrib_type at, seq_region s, coord_system c WHERE s.seq_region_id=sa.seq_region_id AND c.coord_system_id=s.coord_system_id AND at.attrib_type_id=sa.attrib_type_id AND at.code=? AND c.species_id=?",
- -PARAMS => [ $code, $db->species_id() ] );
- }
- }
-
- if ($snpstats) {
- $db->dbc()->sql_helper()->execute_update(
- -SQL =>
-"DELETE sa FROM seq_region_attrib sa, attrib_type at, seq_region s, coord_system c WHERE s.seq_region_id=sa.seq_region_id AND c.coord_system_id=s.coord_system_id AND at.attrib_type_id=sa.attrib_type_id AND at.code=? AND c.species_id=?",
- -PARAMS => [ "SNPCount", $db->species_id() ] );
- }
-
- #
- # Only run on database with genes
- #
-
- my $genes_present;
-
- if ($genestats) {
- my $gene_count =
- $db->dbc()->sql_helper()->execute_single_result(
- -SQL =>
-"select count(*) from gene join seq_region using (seq_region_id) join coord_system using (coord_system_id) where species_id=?",
- -PARAMs => [ $db->species_id() ] );
- $genes_present = ($gene_count) ? 1 : 0;
- } else {
- $genes_present = 0;
- }
-
- #
- # and seq_regions
- #
- my $seq_region_count =
- $db->dbc()->sql_helper()->execute_single_result(
- -SQL =>
-"select count(*) from seq_region join coord_system using (coord_system_id) where species_id=?",
- -PARAMS => [ $db->species_id() ] );
- if ( !$seq_region_count ) {
- print STDERR "No seq_regions for $dbname\n";
- exit();
- }
-
- my $snps_present;
- my $snp_db;
-
- if ($snpstats) {
- $snp_db = variation_attach($db);
- if ( defined $snp_db ) { $snps_present = 1; }
- }
-
- my $slice_adaptor = $db->get_SliceAdaptor();
- my $attrib_adaptor = $db->get_AttributeAdaptor();
-
- # Do not include non-reference sequences ie. haplotypes for human
- #my $top_slices = $slice_adaptor->fetch_all( "toplevel" , undef, 1);
- my $top_slices = $slice_adaptor->fetch_all("toplevel");
-
- while ( my $slice = shift( @{$top_slices} ) ) {
- # print STDERR "Processing seq_region ", $slice->seq_region_name(), "\n";
-
- my @attribs;
-
- if ($genes_present) {
-
- my %attrib_counts;
- my %counts;
-
- my $genes = $slice->get_all_Genes();
-
- while ( my $gene = shift( @{$genes} ) ) {
-
- my $biotype = $gene->biotype();
- if ( $biotype =~ /coding/i && $biotype !~ /non_/i ) {
- if ( $gene->is_known() ) {
- $biotype = "known " . $biotype;
- } else {
- $biotype = "novel " . $biotype;
- }
- }
-
- $counts{$biotype}++;
-
- }
-
- for my $biotype ( keys %counts ) {
- my $attrib_code = $attrib_codes{$biotype};
- if ( !$attrib_code ) {
- print STDERR
- "Unspecified biotype \"$biotype\" in database $dbname.\n";
- next;
- }
-
- $attrib_counts{$attrib_code} += $counts{$biotype};
-
- }
-
- foreach my $attrib_code ( keys %attrib_counts ) {
- my $attrib_code_desc = $attrib_code;
- $attrib_code_desc =~ s/GeneNo_//;
- push @attribs,
- Bio::EnsEMBL::Attribute->new(
- -NAME => $attrib_code_desc . ' Gene Count',
- -CODE => $attrib_code,
- -VALUE => $attrib_counts{$attrib_code},
- -DESCRIPTION => 'Number of ' . $attrib_code_desc . ' Genes'
- );
-
- }
-
- } ## end if ($genes_present)
-
- if ($snps_present) {
- my $count =
- $snp_db->dbc()->sql_helper()->execute_single_result(
- -SQL =>
-"SELECT COUNT(*) FROM variation_feature WHERE seq_region_id = ?",
- -PARAMS => [ $slice->get_seq_region_id ] );
-
- push @attribs,
- Bio::EnsEMBL::Attribute->new(
- -NAME => 'SNP Count',
- -CODE => 'SNPCount',
- -VALUE => $count,
- -DESCRIPTION => 'Total Number of SNPs'
- );
- }
-
- $attrib_adaptor->store_on_Slice( $slice, \@attribs );
- my $slice_attrib_count = @attribs;
- $total_count += $slice_attrib_count;
- # print_chromo_stats([$slice]);
- } ## end while ( my $slice = shift...)
-
- print STDOUT
-"Written $total_count seq reqion attributes to database $dbname on server "
- . $db->dbc()->host() . "\n";
-
-} ## end for my $db_args ( @{ $cli_helper...})
-
-sub print_chromo_stats {
- my $chromosomes = shift;
-
- foreach my $chr (@$chromosomes) {
- print "\nchromosome: ", $chr->seq_region_name(), "\n";
- foreach my $attrib ( @{ $chr->get_all_Attributes() } ) {
- print " ", $attrib->name(), ": ", $attrib->value(), "\n";
- }
- }
-}
-
-#
-# tries to attach variation database.
-#
-
-sub variation_attach {
- my $db = shift;
-
- my $core_db_name;
- $core_db_name = $db->dbc->dbname();
- if ( $core_db_name !~ /_core_/ ) {
- return 0;
- }
- #
- # get a lost of all databases on that server
- #
- my $sth = $db->dbc->prepare("show databases");
- $sth->execute();
- my $all_db_names = $sth->fetchall_arrayref();
- my %all_db_names = map { ( $_->[0], 1 ) } @$all_db_names;
- my $snp_db_name = $core_db_name;
- $snp_db_name =~ s/_core_/_variation_/;
-
- if ( !exists $all_db_names{$snp_db_name} ) {
- return 0;
- }
-
- # this should register the dbadaptor with the Registry
- my $snp_db =
- Bio::EnsEMBL::Variation::DBSQL::DBAdaptor->new(
- -host => $db->dbc()->host(),
- -user => $db->dbc()->username(),
- -pass => $db->dbc()->password(),
- -port => $db->dbc()->port(),
- -dbname => $snp_db_name,
- -group => "variation",
- -species => "DEFAULT" );
-
- return $snp_db;
-} ## end sub variation_attach
-
-sub usage {
- my $indent = ' ' x length($0);
- print <<EOF; exit(0);
-
-
-For each toplevel slice, count the number of genes for each biotype
-(gene stats) or count the number of SNPs (snp stats).
-
-gene stats
-What does it do?
-
-Deletes all seq_region_attrib that have attrib_type code
-with a prefix 'GeneNo_'. All toplevel slices are fetched.
-
-Input data: dna seqence, genes, xrefs, xref projections
-Output tables: seq_region_attrib (attrib_type code with prefix 'GeneNo')
-
-
-When to run it in the release cycle?
-
-After core have finished xref projections
-
-
-Which databases to run it on?
-
-Run on all core databases (including otherfeatures, cdna etc) for each release.
-
-
-How long does it take?
-
-It takes about 10 mins to run for a database in normal queue,
-
-
-snp stats
-
-What does it do?
-
-Deletes out all seq_region_attrib that have attrib_type code of 'SNPCount'.
-Attach variation db if exists. All toplevel slices are fetched.
-For each slice, count the number of SNPs.
-
-This option requires ensembl-variation in perl5lib.
-
-Input data: top level seq regions, variation db
-Output tables: seq_region_attrib (attrib_type code with prefix 'SNPCount')
-
-
-When to run it in the release cycle?
-
-When variation dbs have been handed over
-
-
-Which databases to run it on?
-
-Run on core databases only for new species or if the assembly changed,
-or if the variation positions have changed in the corresponding variation db.
-
-
-How long does it take?
-
-It takes about 20 mins to run for a database in normal queue.
-
-
-
-Usage:
-
- $0 -h host [-port port] -u user -p password \\
- $indent -d database | -pattern pattern \\
- $indent [-mhost ensembl_production host] [-mport ensembl_production host] [-muser ensembl_production host] \\
- $indent -s gene | snp \\
- $indent [-help] \\
-
- -h|host Database host to connect to
-
- -port Database port to connect to (default 3306)
-
- -u|user Database username
-
- -p|pass Password for user
-
- -d|dbname Database name
-
- -pattern Database name regexp
-
- -mhost ensembl_production database host to connect to
-
- -mport ensembl_production database port to connect to
-
- -muser ensembl_production database username
-
- -s|stats 'gene' or 'snp'
-
- -help This message
-
-
-EOF
-
-} ## end sub usage
diff --git a/misc-scripts/density_feature/submit_density_features.pl b/misc-scripts/density_feature/submit_density_features.pl
deleted file mode 100644
index e1ae34ef44..0000000000
--- a/misc-scripts/density_feature/submit_density_features.pl
+++ /dev/null
@@ -1,510 +0,0 @@
-#!/usr/bin/env perl
-
-# The script lists databases/species which should have density features updated at the specified stage in the release cycle. There's an option for submitting a selected script.
-
-use strict;
-use Getopt::Long;
-use DBI qw( :sql_types );
-use Switch;
-
-use FindBin qw($Bin);
-use vars qw($SERVERROOT);
-
-BEGIN {
- $SERVERROOT = "$Bin";
-}
-
-my $gene_gc_path = "$SERVERROOT/..";
-
-my $getdbs;
-my $submit_script;
-my $outdir;
-my @host;
-my @user;
-my @pass;
-my @port;
-
-GetOptions( "getdbs|g", \$getdbs,
- "submit|s=s", \$submit_script,
- "outdir|o=s", \$outdir,
- "host|h=s",\@host,
- "user|u=s",\@user,
- "pass|p=s",\@pass,
- "port=s",\@port,
- "help" , \&usage
- );
-
-my $host_count = @host;
-my $user_count = @user;
-my $pass_count = @pass;
-
-usage() if (!defined @host || !defined @user || !defined @pass );
-
-
-my $port_count = @port;
-my $host_string = join("\|",@host);
-
-# if we have fewer user names specified than hosts copy user name from the first -u parameter
-if ($user_count < $host_count) {
- for (my $i = $user_count; $i < $host_count; $i++) {
- push(@user,$user[0]);
- }
-}
-
-# if we have fewer passwords specified than hosts copy password from the first -p parameter
-if ($pass_count < $host_count) {
- for (my $i = $pass_count; $i < $host_count; $i++) {
- push(@pass,$pass[0]);
- }
-}
-
-if ( (!defined @port) || ($port_count < $host_count) ) {
-
- if (!defined $port[0]) {
- $port[0] = 3306;
- }
-
- for (my $i=1; $i<$host_count;$i++) {
-
- if (!defined $port[$i]) {
- push(@port,$port[0]);
- }
- }
-}
-
-if (!$outdir) {
- $outdir = $ENV{'PWD'};
-} else {
- #strip the final /
- $outdir =~ s/\/$//;
- #test if the directory exists
- if (!-d $outdir) {
- die("Directory $outdir does not exist.");
- }
-
-}
-
-# production/master database location:
-my ( $mhost, $mport ) = ( 'ens-staging1', '3306' );
-my ( $muser, $mpass ) = ( 'ensro', undef );
-my $mdbname = 'ensembl_production';
-
-#connect to the production database
-my $prod_dsn = sprintf( 'DBI:mysql:host=%s;port=%d;database=%s',
- $mhost, $mport, $mdbname );
-my $prod_dbh = DBI->connect( $prod_dsn, $muser, $mpass,
- { 'PrintError' => 1, 'RaiseError' => 1 } );
-
-#get the current release number
-
-my ($current_release ) = $prod_dbh->selectrow_array('select max(db_release) from db where is_current = 1');
-
-
-
-if ( !defined $submit_script ) {
- $getdbs = 1;
-}
-
-if ( defined $getdbs ) {
-# ask the user which stage of the release cycle we're in
- print <<CYCLE;
-
-The output will list density features scripts which can be run at a specified point in the release.
-Select option 3 if variation dbs are handed over. The script will assume that no density features scripts
-have been run yet and list all of them.
-In case you would like to run some scripts earlier in the release, for instance when genebuild and genebuild
-xrefs are complete, select option 1 and the output will only contain density features scripts which can be
-run at this early point in the release. Later, when variation dbs are handed over, run this script again with
-option 3 and bare in mind that it will list all density feature scripts which can be run, inluding
-those that you may have already run before. Make sure that you keep the output from the density feature
-script runs so you don\'t submit those jobs again.
-
-
-Where in the release cycle are we (0,1,2,3)?
-
-0 - None of the below - exit program
-1 - Genebuild and genebuild xrefs complete (excluding projected xrefs), all gene healthchecks cleared
-2 - Compara homologies handed over and core xref projections complete
-3 - Variation dbs handed over
-
-CYCLE
-
-my $response = <>;
-
-if ( !defined $response or $response !~ /[0-4]/ ) {
-
- print <<CYCLE2;
-
-Please specify a valid option: 0, 1, 2 or 3:
-
-Where in the release cycle are we (0,1,2,3)?
-
-0 - None of the below - exit program
-1 - Genebuild and genebuild xrefs complete (excluding projected xrefs), all gene healthchecks cleared
-2 - Compara homologies handed over and core xref projections complete
-3 - Variation dbs hand
-
-CYCLE2
-
- $response = <>;
-}
-
-
-if ($response == 0) {
- $prod_dbh->disconnect;
- exit(0);
-}
-
-
-#dbs for new species or changed assembly
-my @new_sp_assem;
-
-#dbs with changed gene sequence
-
-my @chg_seq;
-
-#dbs with changed repeats
-my @chg_repeats;
-
-#dbs with new variation positions
-my @chg_variation;
-
-#core dbs which have a corresponding variation db
-my @core_with_variation;
-
-if ($response >= 1) {
- #get new dbs, or changed assembly
- @new_sp_assem = map { $_->[0] } @{ $prod_dbh->selectall_arrayref("select distinct concat(full_db_name,'|',db_host) from db_list dl join db d using (db_id) where db_release = $current_release and db_type = 'core' and species_id not in (select distinct species_id from db where db_release <> $current_release and db_type = 'core') union
-select distinct concat(full_db_name,'|',db_host) from db_list dl join db d using (db_id) where db_type = 'core' and is_current = 1 and species_id in (select distinct species_id from changelog_species cs join changelog c using (changelog_id) where release_id = $current_release and status not in ('cancelled','postponed') and assembly = 'Y')") };
- #get dbs with changed sequence
- @chg_seq = map { $_->[0] } @{ $prod_dbh->selectall_arrayref("select distinct concat(full_db_name,'|',db_host) from db_list dl join db d using (db_id) where db_type = 'core' and is_current = 1 and species_id in (select distinct species_id from changelog_species cs join changelog c using (changelog_id) where release_id = $current_release and status not in ('cancelled','postponed') and gene_set = 'Y')") };
-
- #get dbs with changed repeats
- @chg_repeats = map { $_->[0] } @{ $prod_dbh->selectall_arrayref("select distinct concat(full_db_name,'|',db_host) from db_list dl join db d using (db_id) where db_type = 'core' and is_current = 1 and species_id in (select distinct species_id from changelog_species cs join changelog c using (changelog_id) where release_id = $current_release and status not in ('cancelled','postponed') and repeat_masking = 'Y')") };
-
-
- print "1. Density features scripts which can be run when qenebuild and genebuild xrefs (excluding projected xrefs) are complete and all gene healthchecks are cleared: \n\n";
-
- print "gene_gc.pl - run on all core databases (use the commands below or script submit_density_features.pl -submit gene_gc):\n\n";
-
- print "(It takes a total of about 10 hours to run for all core databases in normal queue)\n\n";
-
- for (my $i=0; $i<$host_count;$i++) {
- print "\nbsub -q normal -J genegc_stats -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/core_dbs_$current_release"."_".$host[$i]."_genegc.out -eo $outdir/core_dbs_$current_release" ."_".$host[$i]."_genegc.err perl $gene_gc_path/gene_gc.pl -h ".$host[$i]." -port ".$port[$i]." -u ".$user[$i]." -p ".$pass[$i]." -pattern 'core_$current_release'\n";
- }
-
- print "\npercent_gc_calc.pl – run on core databases for new species, or where sequence or assembly have changed (db names will be stored in file $outdir/percent_gc_data.txt, to submit run submit_density_features.pl -submit percent_gc -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx): \n\n";
-
- print "(It takes about 15 mins to run for a database in normal queue)\n\n";
-
- my %array_union = ();
- foreach my $element (@new_sp_assem, @chg_seq) { $array_union{$element}++ }
- my @dbnames_hosts = sort(keys %array_union);
-
- my $file_path = "$outdir/percent_gc_data.txt";
-
- open(DATAFILE, ">$file_path") or die("Failed to open file $file_path for writing\n");
- foreach my $dbname_host (@dbnames_hosts) {
- my ($db_name, $host) = split(/\|/,$dbname_host);
- if ( ($db_name =~ /core_/) && ( $host_string =~ /$host/) ) {
- print DATAFILE $db_name."\t".$host."\n";
- print $db_name."\t".$host."\n";
- }
- }
- close DATAFILE;
-
- print "\n\nrepeat_coverage_calc.pl – run on core databases for new species, or where sequence, assembly or repeats have changed (db names will be stored in file $outdir/repeat_coverage_data.txt, to submit run submit_density_features.pl -submit repeat_coverage -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx): \n\n";
-
-
- print "(It takes about 1 hour to run for a database in the long queue. The script is slowed down considerably for very fragmented genomes with thousands of toplevel seqs)\n\n";
-
- foreach my $element (@chg_repeats) { $array_union{$element}++ }
- @dbnames_hosts = sort(keys %array_union);
-
- my $file_path = "$outdir/repeat_coverage_data.txt";
- open(DATAFILE, ">$file_path") or die("Failed to open file $file_path for writing\n");
- foreach my $dbname_host (@dbnames_hosts) {
- my ($db_name, $host) = split(/\|/,$dbname_host);
- if ( ($db_name =~ /core_/) && ( $host_string =~ /$host/) ) {
- print DATAFILE $db_name."\t".$host."\n";
- print $db_name."\t".$host."\n";
- }
- }
- close DATAFILE;
-
-}
-
-
-if ($response >= 2) {
- print "\n\n2. Density features scripts which can be run when Compara homologies are handed over and core xref projections are complete:\n\n";
- print "gene_density_calc.pl - run on all core dbs (use the commands below or script submit_density_features.pl -submit gene_density -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx)\n\n";
-
- print "(It takes about 10 mins to run for a database in normal queue)\n\n";
-
- for (my $i=0; $i<$host_count;$i++) {
- print "\nbsub -q normal -J gene_density -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/core_dbs_$current_release"."_".$host[$i]."_gene.out -eo $outdir/core_dbs_$current_release"."_".$host[$i]."_gene.err perl $SERVERROOT/gene_density_calc.pl -h ".$host[$i]." -port ".$port[$i]." -u ".$user[$i]." -p ".$pass[$i]." -pattern 'core_$current_release'\n";
- }
-
- print "\n\nseq_region_stats.pl (gene stats option only) - run on all core databases (use the commands below or script submit_density_features.pl -submit seq_region_stats_gene -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx)\n\n";
-
- print "(It takes about 20 mins to run for a database in normal queue)\n\n";
-
- for (my $i=0; $i<$host_count;$i++) {
- print "\nbsub -q normal -J seqreg_stats_gene -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/core_dbs_$current_release"."_".$host[$i]."_seqreg_gene.out -eo $outdir/core_dbs_$current_release"."_".$host[$i]."_seqreg_gene.err perl $SERVERROOT/seq_region_stats.pl -h ".$host[$i]." -port ".$port[$i]." -u ".$user[$i]." -p ".$pass[$i]." -pattern 'core_$current_release' -s gene\n";
- }
-}
-
-
-if ($response == 3) {
-
- print "\n\n3. Density features scripts which can be run when Variation dbs are handed over:\n";
-
- print "\nvariation_density.pl - run for new species or where the core assembly has changed, or if there are any changes to variation positions in the variation database (species will be stored in file $outdir/variation_density_data.txt, to submit run submit_density_features.pl -submit variation_density -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx):\n\n";
-
- print "(It takes about 25 mins to run for a database in normal queue)\n\n";
-
- #get species for new dbs or changed assembly or where variation positions have changed
- @core_with_variation = map { $_->[0] } @{ $prod_dbh->selectall_arrayref("select distinct concat(full_db_name,'|',db_host) from db_list dl join db d using (db_id) where db_release = $current_release and db_type = 'core' and species_id in (select distinct species_id from db where db_release = $current_release and db_type = 'variation');") };
-
- #this sql relies on variation changed flag in change log
- @chg_variation = map { $_->[0] } @{ $prod_dbh->selectall_arrayref("select distinct concat(full_db_name,'|',db_host) from db_list dl join db d using (db_id) where db_type = 'core' and is_current = 1 and species_id in (select distinct species_id from changelog_species cs join changelog c using (changelog_id) where release_id = $current_release and status not in ('cancelled','postponed') and variation_pos_changed = 'Y')") };
-
- my %core_with_variation = map { $_ => 1 } @core_with_variation;
-
-
- my @new_sp_assem_var;
-
- foreach my $element (@new_sp_assem) {
-
- if ( exists($core_with_variation{$element}) ) {
- push (@new_sp_assem_var, $element);
- }
- }
-
- #make sure dbs flagged having variation positions changed, have a varition db
- my @chg_var;
- foreach my $element (@chg_variation) {
- if (exists($core_with_variation{$element}) ) {
- push (@chg_var, $element);
- }
- }
-
- my %array_union = ();
- foreach my $element (@new_sp_assem_var, @chg_var) { $array_union{$element}++ }
- my @dbnames_hosts = sort(keys %array_union);
-
- my $file_path = "$outdir/variation_density_data.txt";
-
- open(DATAFILE, ">$file_path") or die("Failed to open file $file_path for writing\n");
- foreach my $dbname_host (@dbnames_hosts) {
- my ($db_name, $host) = split(/\|/,$dbname_host);
- if ( ($db_name =~ /([^\s]+)_core/) && ( $host_string =~ /$host/) ) {
- $db_name =~ /([^\s]+)_core/;
- print DATAFILE $1 ."\t".$host."\n";
- print $1 ."\t".$host."\n";
- }
- }
- close DATAFILE;
-
- print "\n\nseq_region_stats.pl (snp stats option only) - run on core databases for new species or if the assembly changed, or if the variation positions have changed in the corresponding variation db (db names will be stored in file $outdir/seq_region_stats_snp_data.txt, to submit run submit_density_features.pl -submit seq_region_stats_snp -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx):\n\n";
-
- print "(It takes about 20 mins to run for a database in normal queue)\n\n";
-
- my $file_path = "$outdir/seq_region_stats_snp_data.txt";
-
- open(DATAFILE, ">$file_path") or die("Failed to open file $file_path for writing\n");
- foreach my $dbname_host (@dbnames_hosts) {
- my ($db_name, $host) = split(/\|/,$dbname_host);
- if ( ($db_name =~ /core_/) && ( $host_string =~ /$host/) ) {
- print DATAFILE $db_name."\t".$host."\n";
- print $db_name."\t".$host."\n";
- }
- }
- close DATAFILE;
-
-}
-
-$prod_dbh->disconnect;
-
-} else {
-
- $prod_dbh->disconnect;
-
-#submit selected script
- my @cmd;
- my $data_file;
- my @print_message;
- my $error;
- my %error_message;
- my $queue;
- my $file_name_end;
- my $script;
- my $script_title;
- my $job_name;
- my $option;
- switch ($submit_script) {
- case 'gene_gc' {
- for (my $i=0; $i<$host_count;$i++) {
- push(@cmd, "bsub -q normal -J genegc_stats -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/core_dbs_$current_release"."_". $host[$i] ."_genegc.out -eo $outdir/core_dbs_$current_release"."_" .$host[$i]. "_genegc.err perl $gene_gc_path/gene_gc.pl -h ".$host[$i]." -port ".$port[$i]." -u ".$user[$i]." -p ".$pass[$i]." -pattern 'core_$current_release'");
- push(@print_message,"Submitting gene GC calculation for host ".$host[$i]." to queue 'normal'. The output from this job goes to the file $outdir/core_dbs_$current_release"."_".$host[$i]."_genegc.out\n");
- }
-
- }
- case 'percent_gc' {
- $data_file = "$outdir/percent_gc_data.txt";
- $queue = "normal";
- $job_name = "gc_calc";
- $file_name_end = "_gc";
- $script = "$SERVERROOT/percent_gc_calc.pl";
- $script_title = "percent GC calculation";
- $option = " -pattern ";
- }
- case 'repeat_coverage' {
- $data_file = "$outdir/repeat_coverage_data.txt";
- $queue = "long";
- $job_name = "repeat_cov";
- $file_name_end = "_repeat";
- $script = "$SERVERROOT/repeat_coverage_calc.pl";
- $script_title = "repeat coverage calculation";
- $option = " -pattern ";
- }
- case 'gene_density' {
- for (my $i=0; $i<$host_count;$i++) {
- push(@cmd, "bsub -q normal -J gene_density -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/core_dbs_$current_release"."_".$host[$i]."_gene.out -eo $outdir/core_dbs_$current_release"."_".$host[$i]."_gene.err perl $SERVERROOT/gene_density_calc.pl -h ".$host[$i]." -port ".$port[$i]." -u ".$user[$i]." -p ".$pass[$i]." -pattern 'core_$current_release'");
- push(@print_message,"Submitting gene density calculation for host ".$host[$i]." to queue 'normal'. The output from this job goes to the file ".$outdir."/core_dbs_".$current_release."_".$host[$i]."_gene.out\n");
- }
- }
- case 'seq_region_stats_gene' {
- for (my $i=0; $i<$host_count;$i++) {
- push(@cmd, "bsub -q normal -J seqreg_stats -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/core_dbs_$current_release"."_".$host[$i]."_seqreg_gene.out -eo $outdir/core_dbs_$current_release"."_".$host[$i]."_seqreg_gene.err perl $SERVERROOT/seq_region_stats.pl -h ".$host[$i]." -port ".$port[$i]." -u ".$user[$i]." -p ".$pass[$i]." -pattern 'core_$current_release' -s gene");
- push(@print_message,"Submitting seq region gene stats for host ".$host[$i]." to queue 'normal'. The output from this job goes to the file ".$outdir."/core_dbs_".$current_release."_".$host[$i]."_seqreg_gene.out\n");
- }
- }
- case 'variation_density' {
- $data_file = "$outdir/variation_density_data.txt";
- $queue = "normal";
- $job_name = "var_density";
- $file_name_end = "_var";
- $script = "$SERVERROOT/variation_density.pl";
- $script_title = "variation density calculation";
- $option = " -s ";
- }
- case 'seq_region_stats_snp' {
- $data_file = "$outdir/seq_region_stats_snp_data.txt";
- $queue = "normal";
- $job_name = "seqreg_stats_snp";
- $file_name_end = "_seqreg_snp";
- $script = "$SERVERROOT/seq_region_stats.pl";
- $script_title = "seq region snp stats";
- $option = " -s snp -pattern ";
- }
- else { usage(); }
- }
-
- if (defined $data_file) {
- open(DATAFILE, "<$data_file") or die("Failed to open file $data_file for reading\n");
- while( my $line = <DATAFILE> ) {
- chomp $line;
- my ($db_name, $host_name) = split(/\t/,$line);
- if ( $host_string =~ /$host_name/) {
- #get user and password for host
- my ( $index )= grep { $host[$_] =~ /$host_name/ } 0..$#host;
- push(@cmd, "bsub -q $queue -J $job_name -M2000000 -R'select[mem>2000] rusage[mem=2000]' -oo $outdir/$db_name$file_name_end.out -eo $outdir/$db_name$file_name_end.err perl $script -h $host_name -port ".$port[$index]." -u ".$user[$index]." -p ".$pass[$index].$option. $db_name);
- push(@print_message,"Submitting ".$script_title." for ".$db_name ." on host ".$host_name." to queue '".$queue."'. The output from this job goes to the file $outdir/$db_name$file_name_end.out\n");
- }
- else {
- $error = 1;
- $error_message{$host_name} = "Host info for $host_name not found in the script's arguments. Please run the script again providing the user and password for host $host_name.\n";
- }
-
- }
- close DATAFILE;
-
- }
-
- #submit jobs to the farm
- if ($error) {
- foreach my $host_in_error (sort(keys %error_message)) {
- print $error_message{$host_in_error};
- }
-
- } else {
- my $cmd_count = 0;
- foreach my $cmd (@cmd) {
- print $print_message[$cmd_count];
- #for testing
- #print "\n\n". $cmd . "\n\n";
- system($cmd);
- $cmd_count++;
- }
- }
-}
-
-sub usage {
- my $indent = ' ' x length($0);
- print <<EOF; exit(0);
-
-
-Options -h -u -p are mandatory and need to be specified for at least one host. When using more than
-one host it\'s possible to leave out the user name and password for the second host and they will
-be copied from the first host: e.g. -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx (user ensadmin
-and password xxxx will be used for both ens-staging1 and ens-staging2).
-The script lists databases/species which should have density features updated at the specified stage in the release cycle when using option -g, e.g.
-
-submit_density_features.pl -g -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx
-
-There\'s an option for submitting a selected script, e.g.
-
-submit_density_features.pl -s gene_gc -h ens-staging1 -h ens-staging2 -u ensadmin -p xxxx
-
-
-Usage:
-
- $0 -h host [-h host] -u user [-u user] -p password [-p password]
- $indent -port port_number [-port port_number]
- $indent [-g] [-s script_name]
- $indent [-o output_directory_path]
- $indent [-help]
-
-
- -h|host Database host (multiple hosts can be specified)
-
- -u|user Database user (each host needs a user specified, if multiple -h|host options
- are given and fewer -u|user options are specified, the first user name will be used
- for the hosts where no user name was given)
-
- -p|pass User password (each host needs a password specified, if multiple -h|host options
- are given and fewer -p|pass options are specified, the first password will be used
- for the hosts where no password was given))
-
- -port Database port (default 3306)
-
- -g|getdbs Use this option to generate input files for the -submit option
-
- -s|submit Use this option to submit a density feature script to the farm:
-
- gene_gc - the script will run on all core databases
-
- percent_gc - the script will run on dbs listed in [outdir]/percent_gc_data.txt
-
- repeat_coverage - the script will run on dbs listed in [outdir]/repeat_coverage_data.txt
-
- gene_density - the script will run on all core databases
-
- seq_region_stats_gene - the script will run on all core databases
-
- variation_density - the script will run for species listed in [outdir]/variation_density_data.txt
-
- seq_region_stats_snp - the script will run on dbs listed in [outdir]/seq_region_stats_snp_data.txt
-
- -o|outdir Path for farm job output and error files as well as data files generated by the script using
- option -g(-getdbs) (current path if not specified)
-
- -help This message
-
-
-EOF
-
-}
-
-
diff --git a/misc-scripts/density_feature/variation_density.pl b/misc-scripts/density_feature/variation_density.pl
deleted file mode 100644
index d7dc86b60c..0000000000
--- a/misc-scripts/density_feature/variation_density.pl
+++ /dev/null
@@ -1,233 +0,0 @@
-#!/usr/local/ensembl/bin/perl -w
-
-# Calculates the variation density from given core database
-
-use strict;
-use Bio::EnsEMBL::Registry;
-use Getopt::Long;
-use Data::Dumper;
-use Bio::EnsEMBL::Utils::ConversionSupport;
-
-$Data::Dumper::Maxdepth = 2;
-
-my $max_slices = 100;
-my $bin_count = 150;
-
-my ($host, $user, $pass, $port, $species);
-# Master database location:
-my ( $mhost, $mport ) = ( 'ens-staging1', '3306' );
-my ( $muser, $mpass ) = ( 'ensro', undef );
-my $mdbname = 'ensembl_production';
-
-my ($block_count, $genome_size, $block_size );
-
-GetOptions( "host|h=s", \$host,
- "user|u=s", \$user,
- "pass|p=s", \$pass,
- "port=i", \$port,
- "mhost=s", \$mhost,
- "mport=i", \$mport,
- "muser=s", \$muser,
- "species|s=s", \$species );
-
-
-usage() if (!$host || !$user || !$pass || !$species );
-
-my $reg = 'Bio::EnsEMBL::Registry';
-$reg->load_registry_from_db(-host => $host, -user => $user, -pass => $pass, -port => $port, -species => $species);
-
-my $density_feature_adaptor = $reg->get_adaptor($species, "core", "DensityFeature") || die "Can't create density feature adaptor";
-my $density_type_adaptor = $reg->get_adaptor($species, "core", "DensityType") || die "Can't create density type adaptor";
-my $analysis_adaptor = $reg->get_adaptor($species, "core", "analysis") || die "Can't create analysis adaptor";
-my $slice_adaptor = $reg->get_adaptor($species, "core", "slice") || die "Can't create slice adaptor";
-
-
-my $variation_feature_adaptor = $reg->get_adaptor($species, "variation", "VariationFeature") || die "Can't create variation feature adaptor";
-
-# TODO - variation from registry
-
-
-# release 63: do not delete analysis, as this is synchronized with production!
-my $sth = $slice_adaptor->dbc->prepare("DELETE df, dt FROM analysis_description ad, density_feature df, density_type dt, analysis a WHERE ad.analysis_id = a.analysis_id AND a.analysis_id=dt.analysis_id AND dt.density_type_id=df.density_type_id AND a.logic_name='snpdensity'");
-$sth->execute();
-
-# Sort slices by coordinate system rank, then by length
-my @sorted_slices = sort( {
- $a->coord_system()->rank() <=> $b->coord_system()->rank()
- || $b->seq_region_length() <=> $a->seq_region_length()
-} @{ $slice_adaptor->fetch_all('toplevel')} );
-# Cannot use this due to PAR regions & basefeatureadaptor's fetch method
-#understanding par region translation
-#} @{ $slice_adaptor->fetch_all('toplevel', '', undef, 1) } );
-
-
-my $analysis = $analysis_adaptor->fetch_by_logic_name('snpdensity');
-
-if ( !defined($analysis) ) {
-
- my $prod_dsn = sprintf( 'DBI:mysql:host=%s;port=%d;database=%s',
- $mhost, $mport, $mdbname );
- my $prod_dbh = DBI->connect( $prod_dsn, $muser, $mpass,
- { 'PrintError' => 1, 'RaiseError' => 1 } );
-
- my ($display_label,$description) = $prod_dbh->selectrow_array("select distinct display_label, description from analysis_description where is_current = 1 and logic_name = 'snpdensity'");
-
- $prod_dbh->disconnect;
-
- $analysis = new Bio::EnsEMBL::Analysis(
- -program => "variation_density.pl",
- -database => "ensembl",
- -gff_source => "variation_density.pl",
- -gff_feature => "density",
- -logic_name => "snpdensity",
- -description => $description,
- -display_label => $display_label,
- -displayable => 1 );
-
- $analysis_adaptor->store($analysis);
-} else {
-
- my $support = 'Bio::EnsEMBL::Utils::ConversionSupport';
- $analysis->created($support->date());
- $analysis_adaptor->update($analysis);
-
-}
-
-# Create and store new density type
-
-my $dt = Bio::EnsEMBL::DensityType->new(-analysis => $analysis,
- -region_features => $bin_count,
- -value_type => 'sum');
-$density_type_adaptor->store($dt);
-
-# Now the actual feature calculation loop
-
-my $slice_count = 0;
-
-my ($current, $current_start, $current_end);
-
-# prepare statement outside of loop
-$sth = $variation_feature_adaptor->prepare("SELECT COUNT(*) FROM variation_feature WHERE seq_region_id = ? AND seq_region_start < ? AND seq_region_end > ?");
-
-my $total_count = 0;
-
-while ( my $slice = shift @sorted_slices){
-
- $block_size = $slice->length() / $bin_count;
-
- my @density_features;
-
- print STDOUT "Calculating SNP densities for ". $slice->seq_region_name() . " with block size $block_size\n";
-
- $current_end = 0;
- $current = 0;
-
- while($current_end < $slice->length) {
-
- $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();
- }
-
- # Get count of SNPs in this region
- $sth->execute($slice->get_seq_region_id(), $current_end, $current_start);
- my $count = ($sth->fetchrow_array())[0];
-
- push @density_features, Bio::EnsEMBL::DensityFeature->new(-seq_region => $slice,
- -start => $current_start,
- -end => $current_end,
- -density_type => $dt,
- -density_value => $count);
- if ($count > 0) {
- #density features with value = 0 are not stored
- $total_count++;
- }
-
- }
-
- $density_feature_adaptor->store(@density_features);
-
- last if ( $slice_count++ > $max_slices );
-
-}
-
-print STDOUT "Written $total_count density features for species $species on server $host\n";
-
-
-sub usage {
- my $indent = ' ' x length($0);
- print <<EOF; exit(0);
-
-What does it do?
-
-Calculates the density of SNP features on top level sequences.
-Deletes out all density_feature and density_type entries having
-analysis logic_name 'snpDensity'. Deletes analysis_description
-where display_label = 'snpDensity'. All toplevel slices are fetched
-and sorted from longest to shortest. Each slice is divided into 150
-bins. For each sub_slice, we count and store the number of
-variation_features (SNPs) on that sub_slice.
-
-
-Deletes out all seq_region_attrib that have attrib_type code of 'SNPCount'.
-Attach variation db if exists. All toplevel slices are fetched.
-For each slice, count the number of SNPs.
-
-Input data: top level seq regions, variation db
-Output tables: updates analysis creation date, density_feature, density_type
-
-The script requires ensembl-variation in perl5lib.
-
-When to run it in the release cycle?
-
-When variation dbs have been handed over
-
-
-Which databases to run it on?
-
-The script updates a core database using data from the corresponding variation database.
-Run it for new species or where the core assembly has changed, or if there are any changes to variation positions in the variation database.
-
-
-How long does it take?
-
-It takes about 25 mins to run for a database in normal queue.
-
-
-
-Usage:
-
- $0 -h host [-port port] -u user -p password \\
- $indent -s species \\
- $indent [-mhost ensembl_production host] [-mport ensembl_production host] [-muser ensembl_production host] \\
- $indent [-help] \\
-
- -h|host Database host to connect to
-
- -port Database port to connect to (default 3306)
-
- -u|user Database username
-
- -p|pass Password for user
-
- -s|species Species name
-
- -mhost ensembl_production database host to connect to
-
- -mport ensembl_production database port to connect to
-
- -muser ensembl_production database username
-
- -help This message
-
-
-EOF
-
-}
--
GitLab