Skip to content

GitLab

Commit f299b420 authored by cvs2git's avatar cvs2git
Browse files

This commit was manufactured by cvs2svn to create branch 'branch- 

vega-49-dev'.

Sprout from master 2008-03-26 14:07:36 UTC Steve Trevanion <st3@sanger.ac.uk> 'merge from vega-48-dev'
Delete:
    misc-scripts/density_feature/vega_repeat_coverage_calc.pl
    misc-scripts/protein_match/process_pmach.pl
    modules/Bio/EnsEMBL/DBSQL/Clone.pm
parent f7ccac98
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 1172 deletions
misc-scripts/density_feature/vega_repeat_coverage_calc.pl deleted 100644 → 0
View file @ f7ccac98
#!/usr/local/bin/perl
=head1 NAME
vega_repeat_coverage_calc.pl - calculate the repeat coverage
=head1 SYNOPSIS
vega_repeat_coverage_calc.pl [options]
General options:
--conffile, --conf=FILE read parameters from FILE
(default: conf/Conversion.ini)
--dbname, db_name=NAME use database NAME
--host, --dbhost, --db_host=HOST use database host HOST
--port, --dbport, --db_port=PORT use database port PORT
--user, --dbuser, --db_user=USER use database username USER
--pass, --dbpass, --db_pass=PASS use database passwort PASS
--logfile, --log=FILE log to FILE (default: *STDOUT)
--logpath=PATH write logfile to PATH (default: .)
--logappend, --log_append append to logfile (default: truncate)
-v, --verbose verbose logging (default: false)
-i, --interactive=0|1 run script interactively (default: true)
-n, --dry_run, --dry=0|1 don't write results to database
-h, --help, -? print help (this message)
--prune undo, i.e. delete from the database changes caused by running the script
=head1 DESCRIPTION
This script calculates the repeat coverage for given database.
The block size is determined so that you have 150 bins for the smallest
chromosome over 5 Mb in length. For chromosomes smaller than 5 Mb, an
additional smaller block size is used to yield 150 bins for the overall
smallest chromosome. This will result in reasonable resolution for small
chromosomes and high performance for big ones.
=head1 LICENCE
This code is distributed under an Apache style licence:
Please see http://www.ensembl.org/code_licence.html for details
=head1 AUTHOR
Patrick Meidl <pm2@sanger.ac.uk>
=head1 CONTACT
Post questions to the EnsEMBL development list ensembl-dev@ebi.ac.uk
=cut
use strict;
use warnings;
no warnings 'uninitialized';
use FindBin qw($Bin);
use vars qw($SERVERROOT);
BEGIN {
$SERVERROOT = "$Bin/../../..";
unshift(@INC, "$SERVERROOT/ensembl/modules");
unshift(@INC, "$SERVERROOT/bioperl-live");
}
use Getopt::Long;
use Pod::Usage;
use Bio::EnsEMBL::Utils::ConversionSupport;
use Bio::EnsEMBL::DensityType;
use Bio::EnsEMBL::DensityFeature;
use Bio::EnsEMBL::Mapper::RangeRegistry;
use POSIX;
$| = 1;
my $support = new Bio::EnsEMBL::Utils::ConversionSupport($SERVERROOT);
# parse options
$support->parse_common_options(@_);
$support->parse_extra_options('prune');
$support->allowed_params($support->get_common_params, 'prune');
if ($support->param('help') or $support->error) {
warn $support->error if $support->error;
pod2usage(1);
}
# ask user to confirm parameters to proceed
$support->confirm_params;
# get log filehandle and print heading and parameters to logfile
$support->init_log;
# connect to database and get adaptors
my $dba = $support->get_database('ensembl');
my $dfa = $dba->get_DensityFeatureAdaptor;
my $dta = $dba->get_DensityTypeAdaptor;
my $aa = $dba->get_AnalysisAdaptor;
my $dbh= $dba->dbc->db_handle;
# check for prune option (undo)
if($support->param('prune')){
#my $query= "delete analysis, density_type, density_feature from analysis, density_type, density_feature where (analysis.program= 'repeat_coverage_calc.pl') and (analysis.analysis_id= density_type.analysis_id) and (density_type.density_type_id= density_feature.density_type_id)";
my $query= "delete analysis, density_type, density_feature from analysis left join density_type on analysis.analysis_id= density_type.analysis_id left join density_feature on density_type.density_type_id= density_feature.density_type_id where analysis.program= 'repeat_coverage_calc.pl'";
if($dbh->do($query)){
$support->log("prune was successfull: any previous entries in the database generated by this script have been deleted\n");
}
else{
$support->log_error("prune failed: any previous entries in the database generated by this script have NOT been deleted\n");
}
}
else{
# Create Analysis object
my $analysis = new Bio::EnsEMBL::Analysis (
-program => "repeat_coverage_calc.pl",
-database => "ensembl",
-gff_source => "repeat_coverage_calc.pl",
-gff_feature => "density",
-logic_name => "PercentageRepeat",
);
$aa->store($analysis) unless ($support->param('dry_run'));
# split chromosomes by size and determine block size
my $chr_slices = $support->split_chromosomes_by_size(5000000);
# loop over block sizes
foreach my $block_size (keys %{ $chr_slices }) {
$support->log("Available chromosomes using block size of $block_size:\n ");
$support->log(join("\n ", map { $_->seq_region_name } @{ $chr_slices->{$block_size} })."\n");
# create DensityType objects
my $density_type = Bio::EnsEMBL::DensityType->new(
-analysis => $analysis,
-block_size => $block_size,
-value_type => 'ratio',
);
$dta->store($density_type) unless ($support->param('dry_run'));
# loop over chromosomes
$support->log_stamped("Looping over chromosomes...\n");
my ($current_start, $current_end);
foreach my $slice (@{ $chr_slices->{$block_size} }) {
$current_start = 1;
my $chr = $slice->seq_region_name;
my $i;
my $bins = POSIX::ceil($slice->end/$block_size);
$support->log_stamped("Chromosome $chr with block size $block_size...\n", 1);
# loop over blocks
while($current_start <= $slice->end) {
$i++;
$current_end = $current_start + $block_size - 1;
if ($current_end > $slice->end) {
$current_end = $slice->end;
}
my $this_block_size = $current_end - $current_start + 1;
my $sub_slice = $slice->sub_Slice($current_start, $current_end);
my $rr = Bio::EnsEMBL::Mapper::RangeRegistry->new or die;
foreach my $repeat (@{ $sub_slice->get_all_RepeatFeatures }) {
$rr->check_and_register("1", $repeat->start, $repeat->end);
}
my $count = 0;
my $non_repeats = $rr->check_and_register("1", 1, $this_block_size);
if (defined $non_repeats) {
foreach my $non_repeat (@{ $non_repeats }) {
$count += ($non_repeat->[1] - $non_repeat->[0]) + 1;
}
}
my $percentage_repeat = (($this_block_size-$count)/$this_block_size)*100;
$support->log_verbose("Chr: $chr | Bin: $i/$bins | ", 2);
$support->log_verbose("\%repeat ".sprintf("%.2f", $percentage_repeat)."\n");
my $df = Bio::EnsEMBL::DensityFeature->new(
-seq_region => $slice,
-start => $current_start,
-end => $current_end,
-density_type => $density_type,
-density_value => $percentage_repeat,
);
$dfa->store($df) unless ($support->param('dry_run'));
$current_start = $current_end + 1;
}
$support->log_stamped("Done.\n", 1);
}
$support->log_stamped("Done.\n");
}
}
# finish logfile
$support->finish_log;
misc-scripts/protein_match/process_pmach.pl deleted 100644 → 0
View file @ f7ccac98
use strict;
=head1 Process pmatch
=head1 Description
=head2 Aims
This script aims to run pmatch and postprocess pmatch to map Ensembl peptides to external databases (currently Swissprot and Refseq but may be extented). The first part of the script runs pmatch, the second part gets the percentage of a match of a unique Ensembl peptide which match to an unique external protein. The third part classify each ensembl match as PRIMARY match (the longest one and the one which will be used for the mapping, PSEUDO, DUPLICATE and REPEAT (pretty arbitrary criterias but may be useful for quality control).
NB: All of the intermediary files are written.
=head2 Options
-ens : Ensembl peptide fasta file
-sp : SP, SPTREMBL fasta file
-refseq: Refseq peptide fasta file
=head2 Contacts
mongin@ebi.ac.uk
birney@ebi.ac.uk
=cut
use Getopt::Long;
my ($ens,$sp,$refseq,$pdb);
&GetOptions(
'ens:s'=>\$ens,
'sp:s'=>\$sp,
'refseq:s'=>\$refseq,
'pdb:s'=>\$pdb
);
&runpmatch();
&postprocesspmatch($sp);
&postprocesspmatch($refseq);
<<<<<<< process_pmach.pl
=======
&postprocesspmatch($pdb);
>>>>>>> 1.8
&finalprocess($sp);
&finalprocess($refseq);
&finalprocess($pdb);
#perl ../../../src/ensembl-live/misc-scripts/protein_match/process_pmach.pl -ens ../primary/SPAN_pepfile -sp ../primary/SPTr.human.expanded -refseq ../primary/hs2.fsa -pdb ../primary/scop_human.fas
sub runpmatch {
print STDERR "Running pmatch\n";
#Run pmatch and store the data in files which will be kept for debugging
my $pmatch1 = "/nfs/griffin2/rd/bin.ALPHA/pmatch -T 14 $sp $ens > ens_sp_rawpmatch";
my $pmatch2 = "/nfs/griffin2/rd/bin.ALPHA/pmatch -T 14 $refseq $ens > ens_refseq_rawpmatch";
#my $pmatch3 = "/nfs/griffin2/rd/bin.ALPHA/pmatch -T 14 $pdb $ens > ens_pdb_rawpmatch";
system($pmatch1); # == 0 or die "$0\Error running '$pmatch1' : $!";
system($pmatch2); #== 0 or die "$0\Error running '$pmatch2' : $!";
#system($pmatch3); #== 0 or die "$0\Error running '$pmatch2' : $!";
}
sub postprocesspmatch {
my ($db) = @_;
my %hash1;
my %hashlength;
#Post process the raw data from pmatch
if ($db eq $sp) {
print STDERR "Postprocessing pmatch for SP mapping\n";
open (OUT, ">ens_sp.processed") || die "Can't open File\n";
open (PROC, "ens_sp_rawpmatch") || die "Can't open File\n";
}
elsif ($db eq $refseq) {
print STDERR "Postprocessing pmatch for REFSEQ mapping\n";
open (OUT, ">ens_refseq.processed") || die "Can't open File\n";;
open (PROC, "ens_refseq_rawpmatch") || die "Can't open file ens_refseq_rawpmatch\n";
}
elsif ($db eq $pdb) {
print STDERR "Postprocessing pmatch for PDB mapping\n";
open (OUT, ">ens_pdb.processed") || die "Can't open File\n";;
open (PROC, "ens_pdb_rawpmatch") || die "Can't open file ens_pdb_rawpmatch\n";
}
while (<PROC>) {
#538 COBP00000033978 1 538 35.3 Q14146 1 538 35.3
my ($len,$id,$start,$end,$tperc,$query,$qst,$qend,$perc) = split;
if ($db eq $refseq) {
#Get only the refseq ac (NP_\d+)
($query) = $query =~ /\w+\|\d+\|\w+\|(\w+)/;
}
my $uniq = "$id:$query";
#Add the percentage of similarity for the Ensembl peptide for a single match
#There is a bug at this step, some similarities can be over 100% !!! This problem may be solved by changing pmatch source code
$hash1{$uniq} += $perc;
$hashlength{$uniq} += $len;
}
#Write out the processed data
foreach my $key ( keys %hash1 ) {
#if (($hashlength{$key} >= 20)) {
if (($hash1{$key} >= 25)) {
($a,$b) = split(/:/,$key);
print OUT "$a\t$b\t$hash1{$key}\n";
}
#else {
# print "$a\t$b\t$hash1{$key}\t$hashlength{$key}\n";
#}
}
close (PROC);
close (OUT);
}
sub finalprocess {
#This final subroutine will use the postprocessed pmatch file and get back the best Ensembl match (labelled as PRIMARY) for a given external known protein.
my ($db) = @_;
if ($db eq $sp) {
print STDERR "Getting final mapping for SP mapping\n";
open (PROC, "ens_sp.processed");
open (OUT, ">ens_sp.final");
}
elsif ($db eq $refseq) {
print STDERR "Getting final mapping for REFSEQ mapping\n";
open (PROC, "ens_refseq.processed") || die "Can' open file ens_refseq.processed\n";
open (OUT, ">ens_refseq.final");
}
elsif ($db eq $pdb) {
print STDERR "Getting final mapping for PDB mapping\n";
open (PROC, "ens_pdb.processed") || die "Can' open file ens_refseq.processed\n";
open (OUT, ">ens_pdb.final");
}
my %hash2;
while (<PROC>) {
my ($ens,$known,$perc) = split;
#if ($perc > 100) {
# print "$ens\t$known\t$perc\n";
#}
if( !defined $hash2{$known} ) {
$hash2{$known} = [];
}
#Each single external protein correspond to an array of objects dealing with the name and the percentage of similarity of the Ensembl peptide matching with the the known external protein.
my $p= NamePerc->new;
$p->name($ens);
$p->perc($perc);
push(@{$hash2{$known}},$p);
}
foreach my $know ( keys %hash2 ) {
my @array = @{$hash2{$know}};
@array = sort { $b->perc <=> $a->perc } @array;
#The Ensembl match to the known protein is labelled as PRIMARY and will be used later for the mapping
my $top = shift @array;
#if ($top->perc >= 20) {
print OUT "$know\t",$top->name,"\t",$top->perc,"\tPRIMARY\n";
foreach $ens ( @array ) {
if( $ens->perc > $top->perc ) {
die "Not good....";
}
}
#If there is more than 20 Ensembl peptides matching a single known protein, these Ensembl peptides are labelled as REPEAT
if (scalar(@array) >= 20) {
foreach my $repeat (@array) {
if( $repeat->perc+1 >= $top->perc ) {
print OUT "$know\t",$repeat->name,"\t",$repeat->perc,"\tDUPLICATE\n";
}
else {
print OUT "$know\t",$repeat->name,"\t",$repeat->perc,"\tREPEAT\n";
}
}
}
#If less than 20, either duplicate if percentage of identity close to the PRIMARY labelled as DUPLICATE or labelled as PSEUDO. DUPLICATEs can also be used for the mapping
if (scalar(@array) < 20) {
foreach my $duplicate (@array) {
if( $duplicate->perc+1 >= $top->perc ) {
print OUT "$know\t",$duplicate->name,"\t",$duplicate->perc,"\tDUPLICATE\n";
}
else {
print OUT "$know\t",$duplicate->name,"\t",$duplicate->perc,"\tPSEUDO\n";
}
}
}
}
#}
close (PROC);
close (OUT);
}
#Set of objects to deal with the script
package NamePerc;
sub new {
my $class= shift;
my $self = {};
bless $self,$class;
return $self;
}
=head2 name
Title : name
Usage : $obj->name($newval)
Function:
Returns : value of name
Args : newvalue (optional)
=cut
sub name{
my $obj = shift;
if( @_ ) {
my $value = shift;
$obj->{'name'} = $value;
}
return $obj->{'name'};
}
=head2 perc
Title : perc
Usage : $obj->perc($newval)
Function:
Returns : value of perc
Args : newvalue (optional)
=cut
sub perc{
my $obj = shift;
if( @_ ) {
my $value = shift;
$obj->{'perc'} = $value;
}
return $obj->{'perc'};
}
modules/Bio/EnsEMBL/DBSQL/Clone.pm deleted 100755 → 0
View file @ f7ccac98
#
# BioPerl module for DB::Clone
#
# Cared for by Ewan Birney <birney@sanger.ac.uk>
#
# Copyright Ewan Birney
#
# You may distribute this module under the same terms as perl itself
# POD documentation - main docs before the code
=head1 NAME
Bio::EnsEMBL::DB::Clone - Object representing one clone
=head1 SYNOPSIS
# $db is Bio::EnsEMBL::DB::Obj
@contig = $db->get_Contigs();
$clone = $db->get_Clone();
@genes = $clone->get_all_Genes();
=head1 DESCRIPTION
Represents information on one Clone
=head1 CONTACT
Describe contact details here
=head1 APPENDIX
The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
=cut
# Let the code begin...
package Bio::EnsEMBL::DBSQL::Clone;
use vars qw(@ISA);
use strict;
# Object preamble - inheriets from Bio::Root::Object
use Bio::Root::Object;
use Bio::EnsEMBL::DBSQL::RawContig;
use Bio::EnsEMBL::DBSQL::Feature_Obj;
use Bio::EnsEMBL::DBSQL::Gene_Obj;
use Bio::EnsEMBL::DB::CloneI;
@ISA = qw(Bio::Root::Object Bio::EnsEMBL::DB::CloneI);
# new() is inherited from Bio::Root::Object
# _initialize is where the heavy stuff will happen when new is called
sub _initialize {
my($self,@args) = @_;
my $make = $self->SUPER::_initialize;
my ($dbobj,$id) = $self->_rearrange([qw(DBOBJ
ID
)],@args);
$id || $self->throw("Cannot make clone db object without id");
$dbobj || $self->throw("Cannot make clone db object without db object");
$dbobj->isa('Bio::EnsEMBL::DBSQL::Obj') || $self->throw("Cannot make clone db object with a $dbobj object");
$self->id($id);
$self->_db_obj($dbobj);
$self->fetch();
# set stuff in self from @args
return $make; # success - we hope!
}
=head2 fetch
Title : fetch
Usage :
Function:
Example :
Returns : nothing
Args :
=cut
sub fetch {
my ($self) = @_;
my $id=$self->id();
my $sth = $self->_db_obj->prepare("select internal_id,id from clone where id = \"$id\";");
my $res = $sth ->execute();
my $rowhash = $sth->fetchrow_hashref;
if( ! $rowhash ) {
# make sure we deallocate sth - keeps DBI happy!
$sth = 0;
$self->throw("Clone $id does not seem to occur in the database!");
}
$self->_internal_id($rowhash->{'internal_id'});
return $self;
}
=head2 delete
Title : delete
Usage : $clone->delete()
Function: Deletes clone (itself), including contigs and features, but not its genes
Example :
Returns : nothing
Args : none
=cut
sub delete {
my ($self) = @_;
#(ref($clone_id)) && $self->throw ("Passing an object reference instead of a variable\n");
my $internal_id = $self->_internal_id;
my @contigs;
my @dnas;
# get a list of contigs to zap
my $sth = $self->_db_obj->prepare("select internal_id,dna from contig where clone = $internal_id");
my $res = $sth->execute;