<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 12 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
color:black;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p
{mso-style-priority:99;
mso-margin-top-alt:auto;
margin-right:0cm;
mso-margin-bottom-alt:auto;
margin-left:0cm;
font-size:12.0pt;
font-family:"Times New Roman","serif";
color:black;}
pre
{mso-style-priority:99;
mso-style-link:"HTML - vooraf opgemaakt Char";
margin:0cm;
margin-bottom:.0001pt;
font-size:10.0pt;
font-family:"Courier New";
color:black;}
span.E-mailStijl17
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.E-mailStijl18
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.HTML-voorafopgemaaktChar
{mso-style-name:"HTML - vooraf opgemaakt Char";
mso-style-priority:99;
mso-style-link:"HTML - vooraf opgemaakt";
font-family:Consolas;
color:black;}
span.E-mailStijl22
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:70.85pt 70.85pt 70.85pt 70.85pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body bgcolor="white" lang="NL" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">Dear Thibaut,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">I appreciate your quick response. I now have an idea how to get the number of spanning reads, but I have another question. Because I would like to have the number of spanning reads per gene, I fetch
the dna features that overlap the gene and count all the reads as you can see in the code below. But as I print the feature IDs, I see only introns. Why are there no exons? And sometimes a big number of reads is shown, like 1728.44444444444. Is this an artefact
or something?<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">Regard Maartje<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">#!/usr/bin/env perl<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">use strict;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">use warnings;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">use Bio::EnsEMBL::DBSQL::DBAdaptor;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"># Connection to the rnaseq database<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">my $rnaseqdb = new Bio::EnsEMBL::DBSQL::DBAdaptor(<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> -host => 'ensembldb.ensembl.org',<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> -port => 5306,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> -user => 'anonymous',<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> -dbname => 'homo_sapiens_rnaseq_69_37');<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"># The GeneAdaptor and the DnaFeatureAdaptor are needed<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">my $rnaseqsa = $rnaseqdb->get_SliceAdaptor();<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">my $rnaseqaa = $rnaseqdb->get_AnalysisAdaptor();<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">my $rnaseqdafa = $rnaseqdb->get_DnaAlignFeatureAdaptor();<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">foreach my $slice (@{$rnaseqsa->fetch_all('toplevel')}) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> foreach my $analysis (@{$rnaseqaa->fetch_all()}) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> # To analyse with a logic name<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> my $logic_name = $analysis->logic_name;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> next if ($logic_name !~ /rnaseq/);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> print 'Looking at genes from analysis: ', $logic_name, "\n";<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> # To load the gene<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> foreach my $gene (@{$slice->get_all_Genes_by_type(undef, $logic_name, 1)}) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> print $gene->display_id, ':', $gene->start, ':', $gene->end, ':', $gene->strand, ' on chromosome ', $gene->seq_region_name, "\n";<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> # Rnaseq models have only one transcript per gene<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> foreach my $transcript (@{$gene->get_all_Transcripts()}) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> print $transcript->display_id, ':', $transcript->start, ':', $transcript->end, ' length: ', $transcript->length, "\n";<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> my $tslice = $transcript->feature_Slice();<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> # To fetch the dna features that overlap the gene<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> my $features = $rnaseqdafa->fetch_all_by_Slice($tslice);<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> my $read_count = 0;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> foreach my $feature (@{$features}) {<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> print "\t", $feature->display_id, ':', $feature->seq_region_start, ':', $feature->seq_region_end, "\tNumber of spanning reads: ", $feature->score, "\n";<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> # To count all the reads that span over a splice site<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> $read_count += $feature->score;<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> print "\tNumber of spanning reads: ", $read_count, "\n";<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> print "\n";<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"> }<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D">}<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="color:#1F497D"><o:p> </o:p></span></p>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm">
<p class="MsoNormal"><b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif";color:windowtext">Van:</span></b><span style="font-size:10.0pt;font-family:"Tahoma","sans-serif";color:windowtext"> dev-bounces@ensembl.org [mailto:dev-bounces@ensembl.org]
<b>Namens </b>Thibaut Hourlier<br>
<b>Verzonden:</b> 19 November 2012 17:18<br>
<b>Aan:</b> dev@ensembl.org<br>
<b>Onderwerp:</b> Re: [ensembl-dev] The number of spanning reads for all genes<o:p></o:p></span></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal" style="margin-bottom:12.0pt">Dear Maartje,<br>
First I have to warn you that for next release which is due in January we will update the models of the RNASeq data for human. These models have been generated with the improved pipeline.<br>
<br>
It is "normal" that you can't get what you want because of the way we put the data in the database. You will need to use a DnaAlignFeatureAdaptor.<br>
Did you have a look at the perl script on this thread: <a href="http://lists.ensembl.org/pipermail/dev/2012-September/008071.html">
http://lists.ensembl.org/pipermail/dev/2012-September/008071.html</a>. It's not exactly what you're trying to do but it is similar and it shows a way of getting the number of spanning reads. It is stored as the score of the DnaAlignFeature linked to the exons.<br>
<br>
Regards<br>
Thibaut<o:p></o:p></p>
<div>
<p class="MsoNormal">On 19/11/12 13:39, <a href="mailto:J.vandeVorst@gen.umcn.nl">
J.vandeVorst@gen.umcn.nl</a> wrote:<o:p></o:p></p>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<p class="MsoNormal">Hello Ensemble team,<o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">I was trying to get the number of spanning reads for one gene by using the API and the database homo_sapiens_rnaseq_69-37. I used the following Perl code, but I didn’t get any results. The API version I use is 69 and
I would like to get the number of spanning reads for all genes. I hope you can help me.</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"><br>
Thank you!</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">Regards Maartje</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">use Bio::EnsEMBL::DBSQL::DBAdaptor;</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> </span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$db = new Bio::EnsEMBL::DBSQL::DBAdaptor(</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -host => 'ensembldb.ensembl.org',</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -port => 5306,</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -user => 'anonymous',</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -dbname => 'homo_sapiens_core_69_37',</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -species => 'Homo_sapiens');</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$ga = $db->get_GeneAdaptor();</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$gene = $ga->fetch_by_stable_id("ENSG00000109061");</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$slice = $gene->slice;</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$rnaseqdb = new Bio::EnsEMBL::DBSQL::DBAdaptor(</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -host => 'ensembldb.ensembl.org',</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -port => 5306,</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -user => 'anonymous',</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> -dbname => 'homo_sapiens_rnaseq_69_37');</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$rnaseqsa = $rnaseqdb->get_SliceAdaptor();</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">$rnaseqslice = $rnaseqsa->fetch_by_name($slice->name);</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">@transcripts = @{$rnaseqslice->get_all_Transcripts('skeletal_rnaseq')};</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">foreach my $transcript (@transcripts) {</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> foreach my $sf (@{$transcript->get_all_supporting_features()}) {</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> #The number of reads that spanned accross the intron</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> print STDOUT $sf->hit_name, ' :', $sf->score, "\n";</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US"> }</span><o:p></o:p></p>
<p class="MsoNormal"><span lang="EN-US">}</span><o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal"><span style="font-size:12.0pt;font-family:"Times New Roman","serif""><br clear="all">
<o:p></o:p></span></p>
<p><span style="font-size:10.0pt;font-family:"Arial","sans-serif"">Het UMC St Radboud staat geregistreerd bij de Kamer van Koophandel in het handelsregister onder nummer 41055629.<br>
The Radboud University Nijmegen Medical Centre is listed in the Commercial Register of the Chamber of Commerce under file number 41055629.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt;font-family:"Times New Roman","serif""><br>
<br>
<br>
<o:p></o:p></span></p>
<pre>_______________________________________________<o:p></o:p></pre>
<pre>Dev mailing list <a href="mailto:Dev@ensembl.org">Dev@ensembl.org</a><o:p></o:p></pre>
<pre>Posting guidelines and subscribe/unsubscribe info: <a href="http://lists.ensembl.org/mailman/listinfo/dev">http://lists.ensembl.org/mailman/listinfo/dev</a><o:p></o:p></pre>
<pre>Ensembl Blog: <a href="http://www.ensembl.info/">http://www.ensembl.info/</a><o:p></o:p></pre>
</blockquote>
<p class="MsoNormal"><span style="font-size:12.0pt;font-family:"Times New Roman","serif""><o:p> </o:p></span></p>
</div>
<br clear=all>
<html>
<body>
<p style="font-size:13px;font-family:arial;">
Het UMC St Radboud staat geregistreerd bij de Kamer van Koophandel in het handelsregister onder nummer 41055629.<BR> The Radboud University Nijmegen Medical Centre is listed in the Commercial Register of the Chamber of Commerce under file number 41055629.<BR>
</p>
</body>
</html>
</body>
</html>