[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