Newer
Older
use strict;
use DBI;
use IPC::Open3;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use vars '@ISA';
@ISA = qw{ XrefMapper::db };
=head1 NAME
XrefMapper::BasicMapper
=head1 DESCIPTION
This is the basic mapper routine. It will create the necessary fasta files for
both the xref and ensembl sequences. These will then be matched using exonerate
Ian Longden
committed
and the results written to another file. By creating a <species>.pm file and
inheriting from this base class different matching routines, parameters, data
sets etc can be set.
=head1 CONTACT
Post questions to the EnsEMBL development list ensembl-dev@ebi.ac.uk
=cut
=head2 dump_seqs
Arg[1]: xref object which holds info needed for the dump of xref
Description: Dumps out the files for the mapping. Xref object should hold
the value of the databases and source to be used.
Returntype : none
Exceptions : will die if species not known or an error occurs while
: trying to write to files.
Caller : general
=cut
my ($self) = @_;
$self->dump_xref();
$self->dump_ensembl();
Arg[1]: xref object which holds info on method and files.
Description: runs the mapping of the list of files with species methods
Returntype : none
Exceptions : none
Caller : general
=cut
my ($self) = @_;
my @list=();
my $i = 0;
foreach my $method (@{$self->method()}){
my @dna=();
push @dna, $method;
push @dna, $self->xref->dir."/xref_".$i."_dna.fasta";
push @dna, $self->ensembl_dna_file();
push @list, \@dna;
my @pep=();
push @pep, $method;
push @pep, $self->xref->dir."/xref_".$i."_peptide.fasta";
push @pep, $self->ensembl_protein_file();
push @list, \@pep;
$i++;
}
if (!defined($self->use_existing_mappings)) {
$self->run_mapping(\@list);
} else {
print "Using existing mappings";
}
=head2 get_species_id_from_species_name
Arg[1]: species name
Description: get the species_id from the database for the named database.
Example : my $id = get_species_id_from_species_name('homo_sapiens');
Returntype : int (species_id)
Exceptions : will die if species does not exist in given xref database.
Caller : general
=cut
sub get_species_id_from_species_name{
my ($xref,$species) = @_;
my $sql = "select species_id from species where name = '".$species."'";
$sth->execute();
my @row = $sth->fetchrow_array();
my $species_id;
if (defined @row) {
$species_id = $row[0];
} else {
print STDERR "Couldn't get ID for species ".$species."\n";
print STDERR "It must be one of :-\n";
$sql = "select name from species";
$sth = $dbi->prepare($sql);
$sth->execute();
while(my @row = $sth->fetchrow_array()){
print STDERR $row[0]."\n";
}
die("Please try again :-)\n");
return $species_id;
}
=head2 get_set_lists
Description: specifies the list of databases and source to be used in the
: generation of one or more data sets.
Returntype : list of lists
Example : my @lists =@{$self->get_set_lists()};
Exceptions : none
Caller : dump_xref
=cut
sub get_set_lists{
my ($self) = @_;
# return [["ExonerateGappedBest1", ["homo_sapiens","Uniprot/SWISSPROT"]]];
return [["ExonerateGappedBest1", ["homo_sapiens","*"]]];
}
=head2 get_source_id_from_source_name
Arg[1]: source name
Description: get the source_id from the database for the named source.
Example : my $id = get_source_id_from_source_name('RefSeq');
Returntype : int (source_id)
Exceptions : will die if source does not exist in given xref database.
Caller : general
=cut
sub get_source_id_from_source_name{
my ($xref, $source) = @_;
my $source_id;
my $sql = "select source_id from source where name = '".$source."'";
my $dbi = $xref->dbi();
my $sth = $dbi->prepare($sql);
$sth->execute();
my @row = $sth->fetchrow_array();
if (defined $row[0] and $row[0] ne '') {
$source_id = $row[0];
# print $source."\t*".$row[0]."*\n";
} else {
print STDERR "Couldn't get ID for source ".$source."\n";
print STDERR "It must be one of :-\n";
$sql = "select name from source";
$sth = $dbi->prepare($sql);
$sth->execute();
while(my @row = $sth->fetchrow_array()){
print STDERR $row[0]."\n";
}
die("Please try again :-)\n");
}
return $source_id;
}
Arg[1]: xref object which holds info on method and files.
Description: Dumps the Xref data as fasta file(s)
Returntype : none
Exceptions : none
Caller : dump_seqs
=cut
sub dump_xref{
my ($self) = @_;
my $xref =$self->xref();
if(!defined($xref->dir())){
if(defined($self->dir)){
$xref->species($self->dir);
}
else{
$xref->dir(".");
}
}
my @method=();
my @lists =@{$self->get_set_lists()};
my $i=0;
Ian Longden
committed
if(defined($self->dumpcheck())){
my $skip = 1;
foreach my $list (@lists){
Ian Longden
committed
if(!-e $xref->dir()."/xref_".$i."_dna.fasta"){
$skip = 0;
}
if(!-e $xref->dir()."/xref_".$i."_peptide.fasta"){
Ian Longden
committed
$skip = 0;
}
$i++;
}
if($skip){
Ian Longden
committed
return;
}
}
$i=0;
foreach my $list (@lists){
Ian Longden
committed
# print "method->".@$list[0]."\n";
$method[$i] = shift @$list;
my $j = 0;
my @species_id=();
foreach my $element (@$list){
while(my $species = shift(@$element)){
# print $j.")\t".$species."\n";
if($species ne "*"){
$species_id[$j] = get_species_id_from_species_name($xref,$species);
}
else{
$species_id[$j] = -1;
}
my $source = shift(@$element);
if($source ne "*"){
$source_id[$j] = get_source_id_from_source_name($xref,$source);
}
else{
$source_id[$j] = -1;
}
Ian Longden
committed
# print $j."\t".$source. "\t".$source_id[$j] ."\n";
# print $j."\t".$species."\t".$species_id[$j]."\n";
$j++;
}
}
#method data fully defined now
Ian Longden
committed
$self->dump_subset($xref,\@species_id,\@source_id,$i);
$i++;
}
$self->method(\@method);
return;
}
=head2 dump_subset
Arg[1]: xref object which holds info on files.
Arg[2]: list of species to use.
Arg[3]: list of sources to use.
Arg[4]: index to be used in file creation.
Description: Dumps the Xref data for one set of species/databases
Returntype : none
Exceptions : none
Caller : dump_xref
=cut
sub dump_subset{
Ian Longden
committed
my ($self,$xref,$rspecies_id,$rsource_id,$index) = @_;
# generate or condition list for species and sources
my $final_clause;
my $use_all = 0;
my @or_list;
for (my $j = 0; $j < scalar(@$rspecies_id); $j++){
my @condition;
if($$rspecies_id[$j] > 0){
push @condition, "x.species_id=" . $$rspecies_id[$j];
if($$rsource_id[$j] > 0){
push @condition, "x.source_id=" . $$rsource_id[$j];
# note if both source and species are * (-1) there's no need for a final clause
if ( !@condition ) {
$use_all = 1;
last;
}
push @or_list, join (" AND ", @condition);
Ian Longden
committed
}
$final_clause = " AND ((" . join(") OR (", @or_list) . "))" unless ($use_all) ;
for my $sequence_type ('dna', 'peptide') {
my $filename = $xref->dir() . "/xref_" . $index . "_" . $sequence_type . ".fasta";
open(XREF_DUMP,">$filename") || die "Could not open $filename";
my $sql = "SELECT p.xref_id, p.sequence, x.species_id , x.source_id ";
$sql .= " FROM primary_xref p, xref x ";
$sql .= " WHERE p.xref_id = x.xref_id AND ";
$sql .= " p.sequence_type ='$sequence_type' ";
$sql .= $final_clause;
if(defined($self->maxdump())){
$sql .= " LIMIT ".$self->maxdump()." ";
my $sth = $xref->dbi()->prepare($sql);
$sth->execute();
while(my @row = $sth->fetchrow_array()){
$row[1] =~ s/(.{60})/$1\n/g;
print XREF_DUMP ">".$row[0]."\n".$row[1]."\n";
close(XREF_DUMP);
$sth->finish();
=head2 dump_ensembl
Description: Dumps the ensembl data to a file in fasta format.
Returntype : none
Exceptions : none
Caller : dump_seqs
=cut
sub dump_ensembl{
my ($self) = @_;
$self->fetch_and_dump_seq();
=head2 fetch_and_dump_seq
Description: Dumps the ensembl data to a file in fasta format.
Returntype : none
Exceptions : wil die if the are errors in db connection or file creation.
Caller : dump_ensembl
=cut
sub fetch_and_dump_seq{
my $db = new Bio::EnsEMBL::DBSQL::DBAdaptor(-species => $self->species(),
-dbname => $self->dbname(),
-host => $self->host(),
-port => $self->port(),
-password => $self->password(),
-group => 'core');
Ian Longden
committed
# if no directory set then dump in the current directory.
if(!defined($self->dir())){
$self->dir(".");
}
$self->ensembl_dna_file($self->dir."/".$self->species."_dna.fasta");
|| die("Could not open dna file for writing: ".$self->ensembl_dna_file."\n");
#
# store ensembl protein file name and open it
#
$self->ensembl_protein_file($self->dir."/".$self->species."_protein.fasta");
Ian Longden
committed
if(defined($self->dumpcheck()) and -e $self->ensembl_protein_file() and -e $self->ensembl_dna_file()){
return;
}
|| die("Could not open protein file for writing: ".$self->ensembl_protein_file."\n");
my $gene_adap = $db->get_GeneAdaptor();
my @gene_ids = @{$gene_adap->list_dbIDs()};
Ian Longden
committed
my $max = undef;
if(defined($self->maxdump())){
$max = $self->maxdump();
}
my $i =0;
foreach my $gene_id (@gene_ids){
my $gene = $gene_adap->fetch_by_dbID($gene_id);
foreach my $transcript (@{$gene->get_all_Transcripts()}) {
Ian Longden
committed
$i++;
my $seq = $transcript->spliced_seq();
$seq =~ s/(.{60})/$1\n/g;
print DNA ">" . $transcript->dbID() . "\n" .$seq."\n";
my $trans = $transcript->translation();
my $translation = $transcript->translate();
Loading full blame...