BioPerl
view release on metacpan or search on metacpan
Bio/Map/Physical.pm view on Meta::CPAN
sub each_markerid {
my ($self) = @_;
return keys (%{$self->{'_markers'}});
}
=head2 get_markerobj
Title : get_markerobj
Usage : my $markerobj = $map->get_markerobj('MARKERA');
Function: returns an object of the marker given in the argument
Returns : object of the marker
Args : scalar representing the marker name
=cut
sub get_markerobj {
my ($self,$marker) = @_;
return 0 if(!defined($marker));
return if($marker eq "");
return if(!exists($self->{'_markers'}{$marker}));
my ($global,$framework,$group,$anchor,$remark,$type,$linkage,$subgroup);
my %mkr = %{$self->{'_markers'}{$marker}};
return $mkr{'marker'} if (ref($mkr{'marker'}) eq 'Bio::Map::FPCMarker');
$type = $mkr{'type'} if(exists($mkr{'type'}));
$global = $mkr{'global'} if(exists($mkr{'global'} ));
$framework = $mkr{'framework'} if(exists($mkr{'framework'}));
$anchor = $mkr{'anchor'} if(exists($mkr{'anchor'}));
$group = $mkr{'group'} if(exists($mkr{'group'}));
$subgroup = $mkr{'subgroup'} if(exists($mkr{'subgroup'}));
$remark = $mkr{'remark'} if(exists($mkr{'remark'}));
my %clones = %{$mkr{'clones'}};
my %contigs = %{$mkr{'contigs'}};
my %markerpos = %{$mkr{'posincontig'}} if(exists($mkr{'posincontig'}));
#*** why doesn't it call Bio::Map::FPCMarker->new ? Seems dangerous...
my $markerobj = bless( {
_name => $marker,
_type => $type,
_global => $global,
_frame => $framework,
_group => $group,
_subgroup => $subgroup,
_anchor => $anchor,
_remark => $remark,
_clones => \%clones,
_contigs => \%contigs,
_position => \%markerpos,
}, 'Bio::Map::FPCMarker');
$self->{'_markers'}{$marker}{'marker'} = $markerobj;
return $markerobj;
}
=head2 each_contigid
Title : each_contigid
Usage : my @contigs = $map->each_contigid();
Function: returns a list of contigs (numbers)
Returns : list of contigs
Args : none
=cut
sub each_contigid {
my ($self) = @_;
return keys (%{$self->{'_contigs'}});
}
=head2 get_contigobj
Title : get_contigobj
Usage : my $contigobj = $map->get_contigobj('CONTIG1');
Function: returns an object of the contig given in the argument
Returns : object of the contig
Args : scalar representing the contig number
=cut
sub get_contigobj {
my ($self,$contig) = @_;
return 0 if(!defined($contig));
return if($contig eq "");
return if(!exists($self->{'_contigs'}{$contig}));
my ($group,$anchor,$uremark,$tremark,$cremark,$startrange,$endrange,
$linkage,$subgroup);
my %ctg = %{$self->{'_contigs'}{$contig}};
my (%position, %pos);
return $ctg{'contig'} if (ref($ctg{'contig'}) eq 'Bio::Map::Contig');
$group = $ctg{'group'} if (exists($ctg{'group'}));
$subgroup = $ctg{'subgroup'} if (exists($ctg{'subgroup'}));
$anchor = $ctg{'anchor'} if (exists($ctg{'anchor'}));
$cremark = $ctg{'chr_remark'} if (exists($ctg{'chr_remark'}));
$uremark = $ctg{'usr_remark'} if (exists($ctg{'usr_remark'}));
$tremark = $ctg{'trace_remark'} if (exists($ctg{'trace_remark'}));
$startrange = $ctg{'range'}{'start'}
if (exists($ctg{'range'}{'start'}));
$endrange = $ctg{'range'}{'end'}
if (exists($ctg{'range'}{'end'}));
my %clones = %{$ctg{'clones'}} if (exists($ctg{'clones'}));
my %markers = %{$ctg{'markers'}} if (exists($ctg{'markers'}));
my $pos = $ctg{'position'};
#*** why doesn't it call Bio::Map::Contig->new ? Seems dangerous...
my $contigobj = bless( {
_group => $group,
_subgroup => $subgroup,
_anchor => $anchor,
_markers => \%markers,
_clones => \%clones,
_name => $contig,
_cremark => $cremark,
_uremark => $uremark,
_tremark => $tremark,
_position => $pos,
_range => Bio::Range->new(-start => $startrange,
-end => $endrange),
}, 'Bio::Map::Contig');
Bio/Map/Physical.pm view on Meta::CPAN
my $score;
$gellen = 3300.0 if (!defined($gellen));
$tol = 7 if (!defined($tol));
if ($numbandsA > $numbandsB) {
$nH = $numbandsA;
$nL = $numbandsB;
}
else {
$nH = $numbandsB;
$nL = $numbandsA;
}
$m = $self->matching_bands($cloneA, $cloneB,$tol);
$logfact[0] = 0.0;
$logfact[1] = 0.0;
for ($i=2; $i<=$nL; $i++) {
$logfact[$i] = $logfact[$i - 1] + log($i);
}
$psmn = 1.0 - ((2*$tol)/$gellen);
$pp = $psmn ** $nH;
$pa = log($pp);
$pb = log(1 - $pp);
$t = 1e-37;
for ($n = $m; $n <= $nL; $n++) {
$c = $logfact[$nL] - $logfact[$nL - $n] - $logfact[$n];
$a = exp($c + ($n * $pb) + (($nL - $n) * $pa));
$t += $a;
}
$score = sprintf("%.e",$t);
return $score;
}
=head2 print_contiglist
Title : print_contiglist
Usage : $map->print_contiglist([showall]); #[Default 0]
Function: prints the list of contigs, markers that hit the contig, the
global position and whether the marker is a placement (P) or
a Framework (F) marker.
Returns : none
Args : [showall] [Default 0], 1 includes all the discrepant markers
=cut
sub print_contiglist{
my ($self,$showall) = @_;
my $pos;
$showall = 0 if (!defined($showall));
my %_contigs = %{$self->{'_contigs'}};
my %_markers = %{$self->{'_markers'}};
my %_clones = %{$self->{'_clones'}};
my @contigs = $self->each_contigid();
my @sortedcontigs = sort {$a <=> $b } @contigs;
print "\n\nContig List\n\n";
foreach my $contig (@sortedcontigs) {
my %list;
my %alist;
my $ctgAnchor = $_contigs{$contig}{'anchor'};
my $ctgGroup = $_contigs{$contig}{'group'};
my @mkr = keys ( %{$_contigs{$contig}{'markers'}} );
foreach my $marker (@mkr) {
my $mrkGroup = $_markers{$marker}{'group'};
my $mrkGlobal = $_markers{$marker}{'global'};
my $mrkFramework = $_markers{$marker}{'framework'};
my $mrkAnchor = $_markers{$marker}{'anchor'};
if($ctgGroup =~ /\d+|\w/ && $ctgGroup != 0) {
if ($mrkGroup eq $ctgGroup) {
if ($mrkFramework == 0) {
$pos = $mrkGlobal."P";
}
else {
$pos = $mrkGlobal."F";
}
$list{$marker} = $pos;
}
elsif ($showall == 1) {
my $chr = $self->group_abbr().$mrkGroup;
$alist{$marker} = $chr;
}
}
elsif ($showall == 1 && $ctgGroup !~ /\d+/) {
my $chr = $self->group_abbr().$mrkGroup;
$alist{$marker} = $chr;
}
}
my $chr = $ctgGroup;
$chr = $self->group_abbr().$ctgGroup if ($ctgGroup =~ /\d+|\w/);
if ($showall == 1 ) {
print " ctg$contig ", $chr, " "
if ($_contigs{$contig}{'group'} !~ /\d+|\w/);
}
elsif ($ctgGroup =~ /\d+|\w/ && $ctgGroup ne 0){
print " ctg",$contig, " ",$chr, " ";
}
while (my ($k,$v) = each %list) {
print "$k/$v ";
}
print "\n" if ($showall == 0 && $ctgGroup =~ /\d+|\w/ &&
$ctgGroup ne 0 );
if ($showall == 1) {
while (my ($k,$v) = each %alist) {
Bio/Map/Physical.pm view on Meta::CPAN
print "Marker List\n\n";
foreach my $marker ($self->each_markerid()) {
print " ",$marker, " ";
my %list;
my %mclones = %{$_markers{$marker}{'clones'}};
foreach my $clone (%mclones) {
if (exists($_clones{$clone}{'contig'}) ) {
my $ctg = $_clones{$clone}{'contig'};
if (exists($list{$ctg})) {
my $clonehits = $list{$ctg};
$clonehits++;
$list{$ctg} = $clonehits;
}
else {
$list{$ctg} = 1;
}
}
}
while (my ($k,$v) = each %list) {
print "$k/$v ";
}
print "\n";
}
}
=head2 print_gffstyle
Title : print_gffstyle
Usage : $map->print_gffstyle([style]);
Function : prints GFF; either Contigwise (default) or Groupwise
Returns : none
Args : [style] default = 0 contigwise, else
1 groupwise (chromosome-wise).
=cut
sub print_gffstyle {
my ($self,$style) = @_;
$style = 0 if(!defined($style));
my %_contigs = %{$self->{'_contigs'}};
my %_markers = %{$self->{'_markers'}};
my %_clones = %{$self->{'_clones'}};
my $i;
my ($depth, $save_depth);
my ($x, $y);
my @stack;
my ($k, $j, $s);
my $pos;
my $contig;
# Calculate the position for the marker in the contig
my @contigs = $self->each_contigid();
my @sortedcontigs = sort {$a <=> $b } @contigs;
my $offset = 0;
my %gffclones;
my %gffcontigs;
my %gffmarkers;
my $basepair = 4096;
foreach my $contig (@sortedcontigs) {
if($_contigs{$contig}{'range'} ) {
$offset = $_contigs{$contig}{'range'}{'start'};
if ($offset <= 0){
$offset = $offset * -1;
$gffcontigs{$contig}{'start'} = 1;
$gffcontigs{$contig}{'end'} =
($_contigs{$contig}{'range'}{'end'} +
$offset ) * $basepair + 1;
}
else {
$offset = 0;
$gffcontigs{$contig}{'start'} =
$_contigs{$contig}{'range'}{'start'} * $basepair;
$gffcontigs{$contig}{'end'} =
$_contigs{$contig}{'range'}{'end'} * $basepair;
}
}
else {
$gffcontigs{$contig}{'start'} = 1;
$gffcontigs{$contig}{'end'} = 1;
}
my @clones = keys %{$_contigs{$contig}{'clones'}};
foreach my $clone (@clones) {
if(exists ($_clones{$clone}{'range'}) ) {
my $gffclone = $clone;
$gffclone =~ s/sd1$//;
$gffclones{$gffclone}{'start'} =
(($_clones{$clone}{'range'}{'start'} + $offset) *
$basepair + 1);
$gffclones{$gffclone}{'end'} =
(($_clones{$clone}{'range'}{'end'}
+ $offset) * $basepair + 1);
}
if(!$contig) {
my %markers = %{$_clones{$clone}{'markers'}}
if (exists($_clones{$clone}{'markers'}));
while (my ($k,$v) = each %markers) {
$gffmarkers{$contig}{$k} =
( ( $_clones{$clone}{'range'}{'start'} +
$_clones{$clone}{'range'}{'end'} ) / 2 ) *
$basepair + 1 ;
}
}
}
Bio/Map/Physical.pm view on Meta::CPAN
$type = "electronicmarker"
if ($_markers{$k}{'type'} eq "eMRK");
if( exists($_markers{$k}{'framework'})) {
$type = "frameworkmarker"
if($_markers{$k}{'framework'} == 1);
$type = "placementmarker"
if($_markers{$k}{'framework'} == 0);
}
print join ("\t","\t$type",$gffmarkers{$contig}{$k}
+ $displacement,$gffmarkers{$contig}{$k}
+ $displacement,".",".",".",
"$type \"$k\"; Name \"$k\"");
my @clonelist;
my @clones = keys %{$_markers{$k}{'clones'}};
foreach my $cl (@clones) {
push (@clonelist, $cl)
if($_clones{$cl}{'contig'} == $contig);
}
$" = " ";
print("; Contig_hit \"ctg$contig - ",
scalar(@clonelist),"\" (@clonelist)\n");
}
}
$displacement += $margin + $gffcontigs{$contig}{'end'};
}
}
}
=head2 _calc_markerposition
Title : _calc_markerposition
Usage : $map->_calc_markerposition();
Function: Calculates the position of the marker in the contig
Returns : none
Args : none
=cut
sub _calc_markerposition {
my ($self) = @_;
my %_contigs = %{$self->{'_contigs'}};
my %_markers = %{$self->{'_markers'}};
my %_clones = %{$self->{'_clones'}};
my $i;
my ($depth, $save_depth);
my ($x, $y);
my @stack;
my ($k, $j, $s);
my $pos;
my $contig;
# Calculate the position for the marker in the contig
my @contigs = $self->each_contigid();
my @sortedcontigs = sort {$a <=> $b } @contigs;
my $offset;
my %gffclones;
my %gffcontigs;
foreach my $marker ($self->each_markerid()) {
my (@ctgmarker, @sortedctgmarker);
my @clones = (keys %{$_markers{$marker}{'clones'}})
if (exists ($_markers{$marker}{'clones'} ));
foreach my $clone (@clones) {
my $record;
$record->{contig} = $_clones{$clone}{'contig'};
$record->{start} = $_clones{$clone}{'range'}{'start'};
$record->{end} = $_clones{$clone}{'range'}{'end'};
push @ctgmarker,$record;
}
# sorting by contig and left position
@sortedctgmarker = sort { $a->{'contig'} <=> $b->{'contig'}
||
$b->{'start'} <=> $a->{'start'}
||
$b->{'end'} <=> $a->{'end'}
} @ctgmarker;
my $ctg = -1;
for ($i=0; $i < scalar(@sortedctgmarker); $i++) {
if ($ctg != $sortedctgmarker[$i]->{'contig'}) {
if ($ctg == -1) {
$ctg = $sortedctgmarker[$i]->{'contig'};
}
else {
if ($depth > $save_depth){
$pos = ($x + $y) >> 1;
$_contigs{$ctg}{'markers'}{$marker} = $pos;
$_markers{$marker}{'posincontig'}{$ctg} = $pos;
}
}
$ctg = $sortedctgmarker[$i]->{'contig'};
$x = $sortedctgmarker[$i]->{'start'};
$y = $sortedctgmarker[$i]->{'end'};
$stack[0] = $y;
$pos = ($x + $y) >> 1;
$_contigs{$ctg}{'markers'}{$marker} = $pos;
$_markers{$marker}{'posincontig'}{$ctg} = $pos;
$depth = $save_depth = 1;
}
elsif ($sortedctgmarker[$i]->{'end'} <= $y) {
$stack[$depth++] = $sortedctgmarker[$i]->{'end'};
# MAX
if ($x < $sortedctgmarker[$i]->{'start'} ) {
$x = $sortedctgmarker[$i]->{'start'};
}
# MIN
Bio/Map/Physical.pm view on Meta::CPAN
}
else {
if ($depth > $save_depth) {
$save_depth = $depth;
$pos = ($x + $y) >> 1;
$_contigs{$ctg}{'markers'}{$marker} = $pos;
$_markers{$marker}{'posincontig'}{$ctg} = $pos;
}
$x = $sortedctgmarker[$i]->{'start'};
$y = $sortedctgmarker[$i]->{'end'};
$stack[$depth++] = $y;
for($j=-1, $k=0, $s=0; $s<$depth; $s++) {
if ($stack[$s] <$x) {
$stack[$s] = -1;
$j = $s if ($j == -1);
}
else {
$k++;
# MIN
$y = $stack[$s] if ($y > $stack[$s]);
if ($stack[$j] == -1) {
$stack[$j] = $stack[$s];
$stack[$s] = -1;
while ($stack[$j] != -1) {$j++;}
}
else {
$j = $s;
}
}
$depth = $k;
}
}
if ($depth > $save_depth) {
$pos = ($x + $y) >> 1;
$_contigs{$ctg}{'markers'}{$marker} = $pos;
$_markers{$marker}{'posincontig'}{$ctg} = $pos;
}
}
}
}
=head2 _calc_contigposition
Title : _calc_contigposition
Usage : $map->_calc_contigposition();
Function: calculates the position of the contig in the group
Returns : none
Args : none
=cut
sub _calc_contigposition{
my ($self) = @_;
my %_contigs = %{$self->{'_contigs'}};
my %_markers = %{$self->{'_markers'}};
my %_clones = %{$self->{'_clones'}};
my @contigs = $self->each_contigid();
my @sortedcontigs = sort {$a <=> $b } @contigs;
foreach my $contig (@sortedcontigs) {
my $position = 0;
my $group;
if (exists($_contigs{$contig}{'group'}) ) {
my %weightedmarkers;
my @mkrs = keys (%{$_contigs{$contig}{'markers'}})
if (exists($_contigs{$contig}{'markers'})) ;
my $chr = $_contigs{$contig}{'group'};
$chr = 0 if ($_contigs{$contig}{'group'} =~ /\?/);
foreach my $mkr (@mkrs) {
if (exists($_markers{$mkr}{'group'})) {
if ( $_markers{$mkr}{'group'} == $chr ) {
my @mkrclones = keys( %{$_markers{$mkr}{'clones'}});
my $clonescount = 0;
foreach my $clone (@mkrclones) {
++$clonescount
if ($_clones{$clone}{'contig'} == $contig);
}
$weightedmarkers{$_markers{$mkr}{'global'}} =
$clonescount;
}
}
}
my $weightedctgsum = 0;
my $totalhits = 0;
while (my ($mpos,$hits) = each %weightedmarkers) {
$weightedctgsum += ($mpos * $hits);
$totalhits += $hits;
}
$position = sprintf("%.2f",$weightedctgsum / $totalhits)
if ($totalhits != 0);
$_contigs{$contig}{'position'} = $position;
}
}
}
=head2 _calc_contiggroup
Title : _calc_contiggroup
Usage : $map->_calc_contiggroup();
Function: calculates the group of the contig
Returns : none
Args : none
=cut
sub _calc_contiggroup {
my ($self) = @_;
my %_contig = %{$self->{'_contigs'}};
my @contigs = $self->each_contigid();
foreach my $ctg (@contigs) {
my $chr = floor($ctg/1000);
$_contig{$ctg}{'group'} = $chr;
}
}
=head2 _setI<E<lt>TypeE<gt>>Ref
Title : _set<Type>Ref
Usage : These are used for initializing the reference of the hash in
Bio::MapIO (fpc.pm) to the corresponding hash in Bio::Map
(physical.pm). Should be used only from Bio::MapIO System.
$map->setCloneRef(\%_clones);
$map->setMarkerRef(\%_markers);
$map->setContigRef(\%_contigs);
Function: sets the hash references to the corresponding hashes
Returns : none
Args : reference of the hash.
=cut
sub _setCloneRef {
my ($self, $ref) = @_;
%{$self->{'_clones'}} = %{$ref};
}
sub _setMarkerRef {
my ($self, $ref) = @_;
%{$self->{'_markers'}} = %{$ref};
}
sub _setContigRef {
my ($self, $ref) = @_;
%{$self->{'_contigs'}} = %{$ref};
}
1;
( run in 0.878 second using v1.01-cache-2.11-cpan-5735350b133 )