Raw content of Bio::EnsEMBL::Analysis::RunnableDB::WGA2GenesDirect
<![CDATA[# Cared for by Ensembl
#
# Copyright GRL & EBI
#
# You may distribute this module under the same terms as perl itself
#
# POD documentation - main docs before the code
=pod
=head1 NAME
Bio::EnsEMBL::Analysis::RunnableDB::WGA2GenesDirect
=head1 SYNOPSIS
my $db = Bio::EnsEMBL::DBAdaptor->new($locator);
my $genscan = Bio::EnsEMBL::Analysis::RunnableDB::WGA2GenesDirect
->new (-db => $pipelinedb,
-input_id => $input_id
-analysis => $analysis );
$genscan->fetch_input();
$genscan->run();
$genscan->write_output(); #writes to stdout
=head1 DESCRIPTION
=cut
package Bio::EnsEMBL::Analysis::RunnableDB::WGA2GenesDirect;
use vars qw(@ISA);
use strict;
use Bio::SeqIO;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Compara::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Analysis::RunnableDB;
use Bio::EnsEMBL::Analysis::Config::General;
use Bio::EnsEMBL::Analysis::Config::WGA2GenesDirect;
#use Bio::EnsEMBL::Analysis::Tools::Logger;
use Bio::EnsEMBL::Analysis::Tools::WGA2Genes::GeneScaffold;
use Bio::EnsEMBL::Analysis::Tools::WGA2Genes::GeneScaffoldDirect;
use Bio::EnsEMBL::Analysis::Tools::ClusterFilter;
use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptUtils
qw(replace_stops_with_introns);
@ISA = qw(Bio::EnsEMBL::Analysis::RunnableDB);
############################################################
=head2 new
Title : fetch_input
=cut
sub new {
my ($class,@args) = @_;
my $self = $class->SUPER::new(@args);
$self->read_and_check_config($WGA2GENES_CONFIG_BY_LOGIC);
return $self;
}
############################################################
=head2 fetch_input
Title : fetch_input
=cut
sub fetch_input {
my( $self) = @_;
my $input_id = $self->input_id;
throw("No input id") unless defined($input_id);
print "YOUR INPUT ID:",$input_id,"\n";
my $q_dbh = Bio::EnsEMBL::DBSQL::DBAdaptor->
new(%{$self->QUERY_CORE_DB});
my $t_dbh = Bio::EnsEMBL::DBSQL::DBAdaptor->
new(%{$self->TARGET_CORE_DB});
my $compara_dbh = Bio::EnsEMBL::Compara::DBSQL::DBAdaptor->
new(%{$self->COMPARA_DB});
my $query_species =
$q_dbh->get_MetaContainerAdaptor->get_Species->binomial;
my $target_species =
$t_dbh->get_MetaContainerAdaptor->get_Species->binomial;
my $gdb_adap = $compara_dbh->get_GenomeDBAdaptor;
my $q_gdb = $gdb_adap->fetch_by_name_assembly($query_species);
my $t_gdb = $gdb_adap->fetch_by_name_assembly($target_species);
########
# check that the default assembly for the query and target agrees
# with that for the method_link_species_set GenomeDBs
########
my ($q_assembly_version, $t_assembly_version);
eval {
$q_assembly_version = $q_dbh->get_CoordSystemAdaptor->
fetch_by_name('toplevel',
$q_gdb->assembly);
$t_assembly_version = $t_dbh->get_CoordSystemAdaptor->
fetch_by_name('toplevel',
$t_gdb->assembly);
};
$@ and do {
throw("Had trouble fetching coord systems for ".
$q_gdb->assembly . " and " .
$t_gdb->assembly . " from core dbs: $@");
};
########
# fetch the genes; need to work in the coordinate space of the
# top-level slice to be consistent with compara
########
my $gene = $q_dbh->get_GeneAdaptor->fetch_by_stable_id($input_id);
$self->_check_gene($gene);
$self->gene($gene);
#########
# get the compara data: MethodLinkSpeciesSet, reference DnaFrag,
# and all GenomicAlignBlocks
#########
my $mlss = $compara_dbh->get_MethodLinkSpeciesSetAdaptor
->fetch_by_method_link_type_GenomeDBs($self->INPUT_METHOD_LINK_TYPE,
[$q_gdb, $t_gdb]);
throw("No MethodLinkSpeciesSet for :\n" .
$self->INPUT_METHOD_LINK_TYPE . "\n" .
$query_species . "\n" .
$target_species)
if not $mlss;
my $dnafrag = $compara_dbh->get_DnaFragAdaptor->
fetch_by_GenomeDB_and_name($q_gdb,
$self->gene->slice->seq_region_name);
my $gaba = $compara_dbh->get_GenomicAlignBlockAdaptor;
my $gen_al_blocks =
$gaba->fetch_all_by_MethodLinkSpeciesSet_DnaFrag($mlss,
$dnafrag,
$gene->start,
$gene->end);
my (%chains, @chains);
foreach my $block (@$gen_al_blocks) {
my $qga = $block->reference_genomic_align;
my ($tga) = @{$block->get_all_non_reference_genomic_aligns};
###########################################################
###### INVESTIGATE: do we want to just use level 1 chains?
###########################################################
#next if $qga->level_id != 1 or $tga->level_id != 1;
# fetch the target slice for later reference
if (not exists $self->target_slices->{$tga->dnafrag->name}) {
$self->target_slices->{$tga->dnafrag->name} =
$t_dbh->get_SliceAdaptor->fetch_by_region('toplevel',
$tga->dnafrag->name);
}
if ($block->reference_genomic_align->dnafrag_strand < 0) {
$block->reverse_complement;
}
push @{$chains{$block->group_id}}, $block;
}
foreach my $chain_id (keys %chains) {
push @chains, [
sort {
$a->reference_genomic_align->dnafrag_start <=>
$b->reference_genomic_align->dnafrag_start;
} @{$chains{$chain_id}}
];
}
$self->genomic_align_block_chains(\@chains);
}
############################################################
=head2 run
Title : run
=cut
sub run {
my ($self) = @_;
my @res_tran;
my $tran_stable_id;
my @final_tran;
print scalar(@{$self->genomic_align_block_chains}),"\n";
foreach my $chain (@{$self->genomic_align_block_chains}) {
my $gene_scaffold = Bio::EnsEMBL::Analysis::Tools::WGA2Genes::GeneScaffold->new(
# my $gene_scaffold = Bio::EnsEMBL::Analysis::Tools::WGA2Genes::GeneScaffoldDirect->new(
-genomic_align_blocks => $chain,
-from_slice => $self->gene->slice,
-to_slices => $self->target_slices,
-transcripts => $self->good_transcripts,
-max_readthrough_dist => $self->MAX_EXON_READTHROUGH_DIST,
-direct_target_coords => 1,
);
foreach my $tran (@{$self->good_transcripts}) {
$tran_stable_id = $tran->stable_id;
my $proj_trans =
# Remember to remove the 1 in place transcripts as right now this is only used for testing purposes
#$gene_scaffold->place_transcript($tran,1);
$gene_scaffold->place_transcript($tran);
if ($proj_trans) {
push @res_tran, $proj_trans;
}
}
}
if ($self->TRANSCRIPT_FILTER){
@res_tran = @{$self->filter->filter_results(\@res_tran)};
}
foreach my $res_tran (@res_tran){
$res_tran = $self->process_transcript($res_tran,
$tran_stable_id);
push @final_tran, $res_tran;
}
# create new gene object for each transcript
print "At the end of RUN, we had ", scalar(@final_tran), " transcripts\n";
$self->output(\@final_tran);
}
############################################################
=head2 write_output
Title : write_output
=cut
sub write_output {
my ($self) = @_;
my $trans_count = 0;
my $t_dbh = Bio::EnsEMBL::DBSQL::DBAdaptor->
new(%{$self->TARGET_CORE_DB});
my $t_gene_adaptor = $t_dbh->get_GeneAdaptor();
foreach my $t (@{$self->output}) {
$t->analysis($self->analysis);
my $gene = Bio::EnsEMBL::Gene->new( -analysis => $self->analysis,
-biotype => 'protein_coding',);
foreach my $tsf ( @{ $t->get_all_supporting_features }){
$tsf->analysis($self->analysis);
}
foreach my $exon (@{$t->get_all_Exons()}){
$exon->analysis($self->analysis);
foreach my $esf (@{$exon->get_all_supporting_features()}){
$esf->analysis($self->analysis);
}
}
$gene->add_Transcript($t);
$t_gene_adaptor->store($gene);
$trans_count++;
print "TRANSCRIPT:\n";
foreach my $e (@{$t->get_all_Exons}) {
printf("%s\tEnsembl\tExon\t%d\t%d\t%d\t%d\t%d\n", $e->slice->seq_region_name, $e->start, $e->end, $e->strand, $e->phase, $e->end_phase);
}
my $seqio = Bio::SeqIO->new(-format => 'fasta',
-fh => \*STDOUT);
$seqio->write_seq($t->translate);
}
print "For gene " . $self->input_id . " you stored ", $trans_count, " transcripts\n";
# to do: write gene back to core target database
}
######################################
# internal methods
#####################################
sub _check_gene {
my ($self, $gene) = @_;
my @good_transcripts;
foreach my $t (@{$gene->get_all_Transcripts}){
#print "Transcript Start: ",$t->start, " END: ", $t->end,"\n";
#print "translation: ",$t->translateable_seq,"\n";
if ((length($t->translateable_seq) % 3) == 0){
push @good_transcripts, $t;
}
else{
warn ("Gene ", $gene->display_id(), " contains no valid transcripts");
#push @good_transcripts, $t;
#print "Done\n";
#exit(1);
}
}
warn ("Gene ", $gene->display_id(), " contains no valid transcripts") if (scalar(@good_transcripts) == 0);
$self->good_transcripts(\@good_transcripts);
# throw an exception here if:
# 1. gene is not protein-coding (translations in all transcripts)
# 2. gene contains at least one transcript with a coding region
# that is not multiple-of-three in length
}
#################################################################
# FUNCTION : process_transcript
#
# Description:
# Subjects the given transcript to a few tests, returning
# the transcript if they succeed, undef if not. If the
# transcript contains less than $max_stops stops, these
# are "spliced out"; otherwise the transcripts is rejected
#################################################################
sub process_transcript {
my ($self,
$tran,
$source_id
) = @_;
return 0 if not $tran;
my (@processed_transcripts);
my @exons = @{$tran->get_all_Exons};
my ($tsf) = @{$tran->get_all_supporting_features};
my $pep = $tran->translate->seq;
if (CORE::length($pep) == 0) {
logger_info("Rejecting proj of $source_id because was just a stop codon");
return 0;
}
##################
# number of stops is non-zero but acceptable. Need to
# operate on the transcript to jump over the stops
$tran = replace_stops_with_introns($tran);
return $tran;
}
###########################
# gets/sets
###########################
sub genomic_align_block_chains {
my ($self, $val) = @_;
if (defined $val) {
$self->{_gen_al_chains} = $val;
}
return $self->{_gen_al_chains};
}
sub gene {
my ($self, $val) = @_;
if (defined $val) {
$self->{_gene_recs} = $val;
}
return $self->{_gene_recs};
}
sub good_transcripts {
my ($self, $val) = @_;
if (defined $val) {
$self->{_good_transcripts} = $val;
}
return $self->{_good_transcripts};
}
sub target_slices {
my ($self, $val) = @_;
if (defined $val) {
$self->{_target_slices} = $val;
}
if (not exists $self->{_target_slices}) {
$self->{_target_slices} = {};
}
return $self->{_target_slices};
}
####################################
# config variable holders
####################################
sub read_and_check_config {
my ($self, $hash) = @_;
$self->SUPER::read_and_check_config($hash);
my $logic = $self->analysis->logic_name;
foreach my $var (qw(INPUT_METHOD_LINK_TYPE
QUERY_CORE_DB
TARGET_CORE_DB
COMPARA_DB)) {
throw("You must define $var in config for logic '$logic'" .
" or in the DEFAULT entry")
if not $self->$var;
}
# filter does not have to be defined, but if it is, it should
# give details of an object and its parameters
if ($self->TRANSCRIPT_FILTER) {
if (not ref($self->TRANSCRIPT_FILTER) eq "HASH" or
not exists($self->TRANSCRIPT_FILTER->{OBJECT}) or
not exists($self->TRANSCRIPT_FILTER->{PARAMETERS})) {
throw("FILTER in config foR '$logic' must be a hash ref with elements:\n" .
" OBJECT : qualified name of the filter module;\n" .
" PARAMETERS : anonymous hash of parameters to pass to the filter");
} else {
my $module = $self->TRANSCRIPT_FILTER->{OBJECT};
my $pars = $self->TRANSCRIPT_FILTER->{PARAMETERS};
(my $class = $module) =~ s/::/\//g;
eval{
require "$class.pm";
};
throw("Couldn't require ".$class." Exonerate2Genes:require_module $@") if($@);
$self->filter($module->new(%{$pars}));
}
}
}
#
# core options
#
sub INPUT_METHOD_LINK_TYPE {
my ($self, $type) = @_;
if (defined $type) {
$self->{_input_method_link_type} = $type;
}
return $self->{_input_method_link_type};
}
sub COMPARA_DB {
my ($self, $db) = @_;
if (defined $db) {
$self->{_compara_db} = $db;
}
return $self->{_compara_db};
}
sub QUERY_CORE_DB {
my ($self, $db) = @_;
if (defined $db) {
$self->{_query_core_db} = $db;
}
return $self->{_query_core_db};
}
sub TARGET_CORE_DB {
my ($self, $db) = @_;
if (defined $db) {
$self->{_target_core_db} = $db;
}
return $self->{_target_core_db};
}
#
# transcript editing and filtering
#
sub TRANSCRIPT_FILTER {
my ($self, $val) = @_;
if (defined $val) {
$self->{_transcript_filter} = $val;
}
return $self->{_transcript_filter};
}
sub filter {
my ($self, $val) = @_;
if ($val) {
$self->{_filter} = $val;
}
return $self->{_filter};
}
sub MAX_EXON_READTHROUGH_DIST {
my ($self, $val) = @_;
if (defined $val) {
$self->{_max_ex_rt_dist} = $val;
}
return $self->{_max_ex_rt_dist};
}
sub MIN_COVERAGE {
my ($self, $val) = @_;
if (defined $val) {
$self->{_min_coverage} = $val;
}
return $self->{_min_coverage};
}
1;
]]>