Commit d66354af authored by cvs2git's avatar cvs2git
Browse files

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

ensembl-50-dev'.

Sprout from master 2008-02-07 11:53:09 UTC Andreas Kusalananda Kähäri <ak4@sanger.ac.uk> 'Fixup POD for method translation().'
Delete:
    misc-scripts/protein_match/process_pmach.pl
    modules/Bio/EnsEMBL/DBSQL/Clone.pm
parent 6cbce1cc
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'};
}
#
# 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;
while( my $rowhash = $sth->fetchrow_hashref) {
push(@contigs,$rowhash->{'internal_id'});
push(@dnas,$rowhash->{'dna'});
}
# Delete from DNA table, Contig table, Clone table
foreach my $contig ( @contigs ) {
my $sth = $self->_db_obj->prepare("delete from contig where internal_id = $contig");
my $res = $sth->execute;
}
foreach my $dna (@dnas) {
$sth = $self->_db_obj->prepare("delete from dna where id = $dna");
$res = $sth->execute;
# Mysql does not optimise or statements in where clauses
$sth = $self->_db_obj->prepare("delete from contigoverlap where dna_a_id = $dna;");
$res = $sth ->execute;
$sth = $self->_db_obj->prepare("delete from contigoverlap where dna_b_id = $dna;");
$res = $sth ->execute;
}
$sth = $self->_db_obj->prepare("delete from clone where internal_id = $internal_id");
$res = $sth->execute;
}
=head2 get_all_Genes
Title : get_all_Genes
Usage :
Function:
Example :
Returns :
Args :
=cut
sub get_all_Genes {
my ($self, $supporting) = @_;
my @out;
my $clone_id = $self->_internal_id();
my %got;
# prepare the SQL statement
my $sth = $self->_db_obj->prepare("
SELECT t.gene
FROM transcript t,
exon_transcript et,
exon e,
contig c
WHERE e.contig = c.internal_id
AND et.exon = e.id
AND t.id = et.transcript
AND c.clone = $clone_id
");
my $res = $sth->execute();
while (my $rowhash = $sth->fetchrow_hashref) {
if( ! exists $got{$rowhash->{'gene'}}) {
my $gene_obj = Bio::EnsEMBL::DBSQL::Gene_Obj->new($self->_db_obj);
my $gene = $gene_obj->get($rowhash->{'gene'}, $supporting);
if ($gene) {
push(@out, $gene);
}
$got{$rowhash->{'gene'}} = 1;
}
}
if (@out) {
return @out;
}
return;
}
=head2 get_Contig
Title : get_Contig
Usage :
Function:
Example :
Returns :
Args :
=cut
sub get_Contig {
my ($self,$contigid) = @_;
# should check this contig is in this clone?
my $contig = $self->_db_obj->get_Contig($contigid);
return $contig->fetch();
}
=head2 get_all_geneid
Title : get_all_geneid
Usage :
Function:
Example :
Returns :
Args :
=cut
sub get_all_my_geneid {
my ($self) = @_;
my $cloneid = $self->_internal_id;
my $sth = $self->_db_obj->prepare("select count(*),cont.clone ,ex.contig,tran.gene " .
"from contig as cont, ".
" transcript as tran, " .
" exon_transcript as et, " .
" exon as ex " .
"where ex.id = et.exon " .
"and tran.id = et.transcript " .
"and cont.clone = $cloneid " .
"and cont.internal_id = ex.contig " .
"group by tran.gene");
my @out;
$sth->execute;
while( my $rowhash = $sth->fetchrow_hashref) {
push(@out,$rowhash->{'gene'});
}
return @out;
}
=head2 get_all_Contigs
Title : get_Contigs
Usage : foreach $contig ( $clone->get_all_Contigs )
Function:
Example :
Returns :
Args :
=cut
sub get_all_Contigs {
my ($self) = @_;
my $sth;
my @res;
my $internal_id = $self->_internal_id();
my $sql = "select id,internal_id from contig where clone = $internal_id";
warn $sql;
$sth= $self->_db_obj->prepare($sql);
my $res = $sth->execute();
my $seen = 0;
my $count = 0;
my $total = 0;
my $version = $self->embl_version();
while( my $rowhash = $sth->fetchrow_hashref) {
my $contig = $self->_db_obj->get_Contig( $rowhash->{'id'});
$contig->internal_id($rowhash->{internal_id});
$contig->seq_version($version);
push(@res,$contig);
$seen = 1;
}
if( $seen == 0 ) {
$self->throw("Clone [$internal_id] has no contigs in the database. Should be impossible, but clearly isn't...");
}
return @res;
}
=head2 get_all_RawContigs
Title : get_rawcontig_by_position
Usage : $obj->get_rawcontig_by_position($position)
Function:
Example :
Returns : returns a raw contig object or undef on error
Args : a position (basepair) in clone
=cut
sub get_rawcontig_by_position {
my ($self, $pos) = @_;
if( !ref $self || ! $self->isa('Bio::EnsEMBL::DB::CloneI') ) {
$self->throw("Must supply a clone to get_all_RawContigs: Bailing out...");
}
if ($pos < 1 ){
$self->throw("get_rawcontig_by_position error: Position must be > 0");
}
my @contigs = $self->get_all_Contigs();
@contigs = sort { $a->embl_offset <=> $b->embl_offset } @contigs;