<html>
  <head>
    <meta content="text/html; charset=ISO-8859-1"
      http-equiv="Content-Type">
  </head>
  <body bgcolor="#FFFFFF" text="#000000">
    Dear Maartje,<br>
    Sorry for the late reply.<br>
    The feature you are printing are the introns, you will need to call
    the method get_all_Exons from a transcript object, and you will need
    to look at the exon boundaries to see which exon correspond to an
    intron.<br>
    <br>
    The extra number you that is printed correspond to the sum of all
    the intron spanning reads.<br>
    <br>
    Regards<br>
    Thibaut<br>
    <br>
    <div class="moz-cite-prefix">On 22/11/12 09:39,
      <a class="moz-txt-link-abbreviated" href="mailto:J.vandeVorst@gen.umcn.nl">J.vandeVorst@gen.umcn.nl</a> wrote:<br>
    </div>
    <blockquote
      cite="mid:CE3C1F4B9384554FB57DA18DA78B7E1F2CA1E9@UMCEXMBX02.umcn.nl"
      type="cite">
      <meta http-equiv="Content-Type" content="text/html;
        charset=ISO-8859-1">
      <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]-->
      <div class="WordSection1">
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">Dear
            Thibaut,<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">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 style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">Regard
            Maartje<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">#!/usr/bin/env
            perl<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">use
            strict;<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">use
            warnings;<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">use
            Bio::EnsEMBL::DBSQL::DBAdaptor;<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">#
            Connection to the rnaseq database<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">my
            $rnaseqdb = new Bio::EnsEMBL::DBSQL::DBAdaptor(<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            -host =>  'ensembldb.ensembl.org',<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            -port =>  5306,<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            -user =>  'anonymous',<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            -dbname =>  'homo_sapiens_rnaseq_69_37');<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">#
            The GeneAdaptor and the DnaFeatureAdaptor are needed<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">my
            $rnaseqsa = $rnaseqdb->get_SliceAdaptor();<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">my
            $rnaseqaa = $rnaseqdb->get_AnalysisAdaptor();<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">my
            $rnaseqdafa = $rnaseqdb->get_DnaAlignFeatureAdaptor();<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">foreach
            my $slice (@{$rnaseqsa->fetch_all('toplevel')}) {<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">   
            foreach my $analysis (@{$rnaseqaa->fetch_all()}) {<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            # To analyse with a logic name<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            my $logic_name = $analysis->logic_name;<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            next if ($logic_name !~ /rnaseq/);<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            print 'Looking at genes from analysis: ', $logic_name, "\n";<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            # To load the gene<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            foreach my $gene (@{$slice->get_all_Genes_by_type(undef,
            $logic_name, 1)}) {<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">           
            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 style="color:#1F497D" lang="EN-US">           
            # Rnaseq models have only one transcript per gene<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">           
            foreach my $transcript (@{$gene->get_all_Transcripts()})
            {<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            print $transcript->display_id, ':',
            $transcript->start, ':', $transcript->end, '  length:
            ', $transcript->length, "\n";<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               my
            $tslice = $transcript->feature_Slice();<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            # To fetch the dna features that overlap the gene<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            my $features = $rnaseqdafa->fetch_all_by_Slice($tslice);<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            my $read_count = 0;<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            foreach my $feature (@{$features}) {<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">                       
            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 style="color:#1F497D" lang="EN-US">                       
            # To count all the reads that span over a splice site<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">                       
            $read_count += $feature->score;<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            }<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">               
            print "\tNumber of spanning reads: ", $read_count, "\n";<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">           
            }<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">           
            print "\n";<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">       
            }<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">   
            }<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US">}<o:p></o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><o:p> </o:p></span></p>
        <p class="MsoNormal"><span style="color:#1F497D" lang="EN-US"><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">
                <a class="moz-txt-link-abbreviated" href="mailto:dev-bounces@ensembl.org">dev-bounces@ensembl.org</a> [<a class="moz-txt-link-freetext" href="mailto:dev-bounces@ensembl.org">mailto:dev-bounces@ensembl.org</a>]
                <b>Namens </b>Thibaut Hourlier<br>
                <b>Verzonden:</b> 19 November 2012 17:18<br>
                <b>Aan:</b> <a class="moz-txt-link-abbreviated" href="mailto:dev@ensembl.org">dev@ensembl.org</a><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
            moz-do-not-send="true"
            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
              moz-do-not-send="true"
              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 moz-do-not-send="true" href="mailto:Dev@ensembl.org">Dev@ensembl.org</a><o:p></o:p></pre>
          <pre>Posting guidelines and subscribe/unsubscribe info: <a moz-do-not-send="true" 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 moz-do-not-send="true" 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">
      <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>
      <br>
      <fieldset class="mimeAttachmentHeader"></fieldset>
      <br>
      <pre wrap="">_______________________________________________
Dev mailing list    <a class="moz-txt-link-abbreviated" href="mailto:Dev@ensembl.org">Dev@ensembl.org</a>
Posting guidelines and subscribe/unsubscribe info: <a class="moz-txt-link-freetext" href="http://lists.ensembl.org/mailman/listinfo/dev">http://lists.ensembl.org/mailman/listinfo/dev</a>
Ensembl Blog: <a class="moz-txt-link-freetext" href="http://www.ensembl.info/">http://www.ensembl.info/</a>
</pre>
    </blockquote>
    <br>
  </body>
</html>