Which is the best option to plot two different tracks (Chromosomes, fragments, scaffolds...) with Bio::Graphics?

Question

I am working with the BioPerl module Bio::Graphics to plot some genomic coordinates. I have regions from different scaffolds or tracks that have evolutionary conserved elements, the idea is to plot the correspondent fragment with its elements in the same plot. How can I do this?

This is my raw data:

reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_235  SNORD100    86265   86176   36171   136261  52121   51889   -   RNAzCandidate   217 240680
reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_235  SNORD100    86265   86176   36171   136261  73633   73397   -   RNAzCandidate   224 240680
reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_235  SNORD100    86265   86176   36171   136261  114129  114369  +   RNAzCandidate   232 240680
reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_235  SNORD100    86265   86176   36171   136261  118835  119071  +   RNAzCandidate   224 240680
reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_235  SNORD100    86265   86176   36171   136261  130220  130454  +   RNAzCandidate   221 240680
reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_235  SNORD100    86265   86176   36171   136261  129525  129758  +   RNAzCandidate   219 240680
reftig_235  SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_2    SNORD100    86265   86176   36171   136261  124395  124625  +   RNAzCandidate   213 240680
reftig_2    SNORD100    86265   86176   36171   136261  37251   37449   -   protein_coding.exon .   240680
reftig_2    SNORD100    86265   86176   36171   136261  124687  124923  +   RNAzCandidate   224 240680

Where the fragment name is: reftig_235 and reftig_2. So when I pass my code:

#!/usr/bin/perl

# This is code example 3 in the Graphics-HOWTO
use strict;
use lib '/home/lstein/projects/bioperl-live';
use Bio::Graphics;
use Bio::SeqFeature::Generic;
use Data::Dumper;

my (@fill, @comparison, @db, @temp);
open my $IN, "< $ARGV[0]";

while (<$IN>){
    chomp $_;
    push @db, $_;
}

@temp = split /\t/, $db[0];
my $chrLen = $temp[11];
my $chrS = $temp[4];
my $chrE = $temp[5];
my $chr = $temp[0];
my $ann = $temp[9];
my $typ = $temp[1];

my $panel = Bio::Graphics::Panel->new(
  -length    => $chrLen, #longitud scala
  -width => 1500, #Ancho del gráfico
  -grid => 1,
  -pad_left  => 20,
  -pad_right => 20,
  -truetype => 1
  );

my $full_length = Bio::SeqFeature::Generic->new(
  -start => 0,
  -end   => $chrLen,
  -display_name => $chr,
                        );
my $partial_length = Bio::SeqFeature::Generic->new(
  -start => $chrS,
  -end   => $chrE,
  -display_name => "$chrS\-$chrE",
  );

$panel->add_track($full_length,
                  -key => $chr,
                  -glyph   => 'extending_arrow',
                  -tick    => 3,
                  -fgcolor => 'black',
                  -double  => 1,
                 );

$panel->add_track($partial_length,
                  -key => "$chrS\-$chrE",
                  -glyph   => 'extending_arrow',
                  -tick    => 3,
                  -fgcolor => 'blue',
                  -double  => 1,
                 );

my $track1 = $panel->add_track(
                #-category => 'HMMs'
                -key => 'ncRNAs by HMMs',
                -glyph => 'generic',
                #-glyph       => 'graded_segments',
                -label       => 1,
                -bgcolor     => 'lime',
                -stranded   => 1,
                              #-min_score   => 0,
                              #-max_score   => 240,
                -font2color  => 'black',
                -sort_order => 'left|longest',
                              #-sort_order  => 'high_score',
                              -description => sub {
                                my $feature = shift;
                                my $score   = $feature->score;
                                return "$typ";
                               },
                -bump => 3,
                #d-decorate_introns => 1,
                             );

my $track2 = $panel->add_track(
                #-category => 'Others',
                -key => 'RNAz Predictions',
                              #-glyph       => 'graded_segments',
                 -hbumppad => 1,
                -glyph       => 'transcript',
                #-height       => 7,
                -label       => 1,
                -bgcolor     => 'cyan',
                -min_score   => 0,
                -max_score   => 400,
                -strand_arrow   => 1,
                -label       => 1,
                -font2color  => 'red',
                -sort_order  => 'strand|left|longest',
                #-connector   => 'solid',
                -description => sub {
                                my $feature = shift;
                                my $score   = $feature->score;
                                return "RNAzAnn.; score=$score";
                               },
                             );


my $track3 = $panel->add_track( # Annotations
                -key => 'Genomic Annotations',
                -glyph       => 'transcript2',
                #-glyph       => 'transcript2',
                #-height      => 7,
                -label       => 1,
                -bgcolor     => 'teal',
                -min_score   => 0,
                -max_score   => 400,
                -strand_arrow   => 1,
                -label       => 1,
                -font2color  => 'purple',
                -sort_order  => 'strand|left|longest',
                -description => sub {
                                my $feature = shift;
                               #my $ann1   = $feature->annotation;
                                return "Ann.=$ann";
                               },
                             );


while (<>) { # read blast file
  chomp;
  next if /^\#/;  # ignore comments

#scaffold_401   5_8S_rRNA   3054    3211    0   22710   2620    2831    +   RNAzCandidate   176 22710
#scaffold_401   5_8S_rRNA   3054    3211    0   22710   13649   15877   -   Gaze.mRNA   .   22710
#reftig_235      SNORD100        86265   86176   36171   136261  86850   87250   +       RNAzCandidate   387    240680
#reftig_235      SNORD100        86265   86176   36171   136261  37251   37449   -       protein_coding.exon    .       240680

my($chr,$type,$startncRNA,$endncRNA,$startChr,$endChr, $startAnn, $endAnn, $sense, $annotation, $bitscore, $chrlen) = split /\t/;
    push @comparison, "$chr\.$type\.$startncRNA\.$endncRNA";
    my ($num,$sens, $sen, $se);
    #my $nameAnn = $annotation;
    #my @temp = (split /\;/, $nameAnn);
    #nameAnn =~ s/(^gene\_id\".*\")\;(.*)$/$1/g;
my ($feature1, $feature2, $feature3);
#if($startncRNA < $endncRNA){


if($startncRNA < $endncRNA){
    $sen = "+1";} else {
    $sen = "-1";}
    $feature1 = Bio::SeqFeature::Generic->new(
    -display_name => "\($startncRNA\-$endncRNA\)",
    -start        => $startncRNA,
    -end          => $endncRNA,
    -strand => $sen,
    );

if($annotation eq "RNAzCandidate"){
    if($startAnn < $endAnn){
    $sens = "+1";} else {
    $sens = "-1";}
    $feature2 = Bio::SeqFeature::Generic->new(
    -score        => $bitscore,
    -display_name => "\($startAnn\-$endAnn\)",
    -start        => $startAnn,
    -end          => $endAnn,
    -strand        => $sens,
    );
} elsif ($annotation ne "RNAzCandidate" && $bitscore eq "\."){
    if($sense eq "\+"){
    $se = "+1";} else { $se = "-1";}
    $feature3 = Bio::SeqFeature::Generic->new(
   -display_name => "\($startAnn\-$endAnn\)",
   -start        => $startAnn,
   -end          => $endAnn,
   -strand      => $se,
    );
} else {next;}
    $num = scalar @comparison;
    if ($num == 1){
    $track1->add_feature($feature1); #Anchor ncRNA
    $track2->add_feature($feature2); #RNAzCandidates
    $track3->add_feature($feature3); #Annotations
} else {
    my @split1 = split /\./, $comparison[-2];
    my @split2 = split /\./, $comparison[-1];
    if($split1[0] eq $split2[0] && $split1[1] eq $split2[1] && $split1[2] eq $split2[2] && $split1[3] eq $split2[3]){
        $track2->add_feature($feature2); #RNAzCandidates
        #$track3->add_feature($feature3); #Annotations
        shift @comparison;
        next;
    } else {
    $track1->add_feature($feature1); #Anchor ncRNA
    $track2->add_feature($feature2); #RNAzCandidates
    $track3->add_feature($feature3); #Annotations
    shift @comparison;
}}}
print $panel->gd->gif;

$panel->finished;
#print $panel->gif;

And results in this plot:

Plot generated by BioPerl, represent generic coordinates

But it does not plot the second fragment called reftig_2. What can I do? Thanks in advance!

Please stop adding your code and data as HTML / CSS / JavaScript snippets. That is meant to allow people to post a working web page with which they are having trouble. For bare blocks of Perl and data it is sufficient just to indent the code by four spaces, which will cause it to be highlighted as a code segment. You can select the lines you want to indent and press Ctrl-K to indent them if you like. — Borodin, Apr 14 '15 at 01:05

Which is the best option to plot two different tracks (Chromosomes, fragments, scaffolds...) with Bio::Graphics?

0 Answers0