我正在使用 BioPerl 模块Bio::Graphics来绘制一些基因组坐标。我有来自不同支架或轨道的区域，它们具有进化保守的元素，这个想法是将对应的片段与其元素绘制在同一个图中。我怎样才能做到这一点？
这是我的原始数据：
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
其中片段名称为：reftig_235 和 reftig_2。所以当我通过我的代码时：
#!/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;
结果如下：
但它没有绘制名为 reftig_2 的第二个片段。我能做些什么？提前致谢！