[ensembl-dev] MXEs labeling the pairs of transcripts (SOS) - With incomplete code.

Julio nunesjulioc at gmail.com
Wed Apr 18 23:33:58 BST 2012


Dear people-at-Ensembl, I hope you can help me with this:

I need to filter the alternative splicing data for mutually exclusive 
exons labeling the pairs of transcripts that are involved in event.

Thanks in advance!

Best regards, Julio.

###


This code is incomplete and I am trying to finish it (I copied from 
http://lists.ensembl.org/pipermail/dev/2011-May/001226.html ).

#!/usr/local/ensembl/bin/perl -w

use lib "/home/julio/src/ensembl/modules";
use strict;
use warnings;
use Getopt::Long;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Feature;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::DBSQL::SplicingEventAdaptor;
use Bio::EnsEMBL::DBSQL::GeneAdaptor;
use Bio::EnsEMBL::DBSQL::AttributeAdaptor;
use Bio::EnsEMBL::DBSQL::SplicingTranscriptPairAdaptor;
use Bio::EnsEMBL::Utils::Exception qw(throw);

my $host   = '';
my $port   = '3306';
my $dbname = '';
my $dbuser = '';
my $dbpass = '';
my $file   = undef;
my $help;
my @coord_system;

&GetOptions(
         'dbhost:s'   => \$host,
         'dbport:n'   => \$port,
         'dbname:s'   => \$dbname,
         'dbuser:s'   => \$dbuser,
         'file:s'      => \$file,
         'dbpass:s'   => \$dbpass,
         'h|help'     => \$help,
         ) or ($help = 1);

if(!$host || !$dbuser || !$dbname) {
         print STDERR "Can't get any information without database 
details\n";
         print STDERR "-dbhost $host -dbuser $dbuser -dbname $dbname ".
                 " -dbpass $dbpass\n";
         $help = 1;
}

if ($help) {
         exec('perldoc', $0);
}

my $db = Bio::EnsEMBL::DBSQL::DBAdaptor->new
         (-dbname => $dbname,
          -host   => $host,
          -user   => $dbuser,
          -port   => $port,
          -pass   => $dbpass);

my $gene_adaptor = $db->get_GeneAdaptor();
my $se_adaptor = $db->get_SplicingEventAdaptor();
my $attrib_adaptor = $db->get_AttributeAdaptor();
my $SplicingTranscriptPairAdaptor = 
$db->get_SplicingTranscriptPairAdaptor();
my @stable_gene_ids = @{$gene_adaptor->list_stable_ids()};
my $size = scalar @stable_gene_ids;
print STDERR "Number of stable ids:\t" . $size . "\n";

my $count = 0;

# example genes
my @genes = ("ENSG00000002016", "ENSG00000002746");

# don't know how to get this information easyly from the core API

my $attrib_names = {
         MXE => "Mutually exclusive exons",
         };

print  " AS name\tType\tShort desc\tlocation\tdetails\n";
print  
"+--------------------------------------------------------------------------+\n";

for my $id (@genes) {

         $count++;
         my $gene = $gene_adaptor->fetch_by_stable_id($id);
         my $gene_id = $gene->display_id();
         my $biotype = $gene->biotype();
         my $chr = $gene->slice->seq_region_name();
         my $strand = $gene->strand();
         my $start = $gene->start();
         my $end = $gene->end();

#And get all the AS events sorted by type and location for a given gene:

         for my $se (  @ { $se_adaptor->fetch_all_by_Gene($gene) } ) {

             # filter on Mutually exclusive events
             next unless ($se->type() eq 'MXE');

             # print location
             my $location = ($se->strand() == 1) ? ("+\t$chr:" . 
$se->start() . "-" . $se->end()) : ("-\t$chr:" . $se->end() . "-" . 
$se->start());
             print  $se->type() .  "\t$location";


             # retrieve the trans
             my $f = $f_map{$splicing_transcript_pair_id};
             my $trans = 
$SplicingTranscriptPairAdaptor->fetch_by_dbID($splicing_transcript_pair_id);

             # concatenate information about this trans


             # depending to which transcript the mutually exclusive exon 
is part of, assign this information to list s1 or s2


             # finally print this information
             print "\tAlternative trans: " . join("; ", @s1) . " / " . 
join(";", @s2) . "\n";
         }
}
exit 0;

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.ensembl.org/pipermail/dev_ensembl.org/attachments/20120418/e14f4225/attachment.html>


More information about the Dev mailing list