3 # BioPerl module for Bio::Map::Physical
5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
7 # Cared for by Sendu Bala <bix@sendu.me.uk>
11 # You may distribute this module under the same terms as perl itself
13 # POD documentation - main docs before the code
17 Bio::Map::Physical - A class for handling a Physical Map (such as FPC)
23 # accquire a Bio::Map::Physical using Bio::MapIO::fpc
24 my $mapio = Bio::MapIO->new(-format => "fpc",-file => "rice.fpc",
27 my $physical = $mapio->next_map();
29 # get all the markers ids
30 foreach my $marker ( $physical->each_markerid() ) {
31 print "Marker $marker\n";
33 # acquire the marker object using Bio::Map::FPCMarker
34 my $markerobj = $physical->get_markerobj($marker);
36 # get all the clones hit by this marker
37 foreach my $clone ($markerobj->each_cloneid() ) {
44 This class is basically a continer class for a collection of Contig maps and
45 other physical map information.
47 Bio::Map::Physical has been tailored to work for FPC physical maps, but
48 could probably be used for others as well (with the appropriate MapIO
51 This class also has some methods with specific functionalities:
53 print_gffstyle() : Generates GFF; either Contigwise[Default] or
56 print_contiglist() : Prints the list of Contigs, markers that hit the
57 contig, the global position and whether the marker
58 is a placement (<P>) or a Framework (<F>) marker.
60 print_markerlist() : Prints the markers list; contig and corresponding
63 matching_bands() : Given two clones [and tolerence], this method
64 calculates how many matching bands do they have.
66 coincidence_score() : Given two clones [,tolerence and gellen], this
67 method calculates the Sulston Coincidence score.
69 For faster access and better optimization, the data is stored internally in
70 hashes. The corresponding objects are created on request.
76 User feedback is an integral part of the evolution of this and other
77 Bioperl modules. Send your comments and suggestions preferably to
78 the Bioperl mailing list. Your participation is much appreciated.
80 bioperl-l@bioperl.org - General discussion
81 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
85 Please direct usage questions or support issues to the mailing list:
87 I<bioperl-l@bioperl.org>
89 rather than to the module maintainer directly. Many experienced and
90 reponsive experts will be able look at the problem and quickly
91 address it. Please include a thorough description of the problem
92 with code and data examples if at all possible.
96 Report bugs to the Bioperl bug tracking system to help us keep track
97 of the bugs and their resolution. Bug reports can be submitted via the
100 http://bugzilla.open-bio.org/
102 =head1 AUTHOR - Gaurav Gupta
104 Email gaurav@genome.arizona.edu
108 Sendu Bala bix@sendu.me.uk
110 =head1 PROJECT LEADERS
112 Jamie Hatfield jamie@genome.arizona.edu
113 Dr. Cari Soderlund cari@genome.arizona.edu
115 =head1 PROJECT DESCRIPTION
117 The project was done in Arizona Genomics Computational Laboratory (AGCoL)
118 at University of Arizona.
120 This work was funded by USDA-IFAFS grant #11180 titled "Web Resources for
121 the Computation and Display of Physical Mapping Data".
123 For more information on this project, please refer:
124 http://www.genome.arizona.edu
128 The rest of the documentation details each of the object methods.
129 Internal methods are usually preceded with a _
133 # Let the code begin...
135 package Bio
::Map
::Physical
;
136 use vars
qw($MAPCOUNT);
141 use Bio::Map::Contig;
142 use Bio::Map::FPCMarker;
144 use base qw(Bio::Map::SimpleMap);
145 BEGIN { $MAPCOUNT = 1; }
147 =head1 Access Methods
149 These methods let you get and set the member variables
154 Usage : my $version = $map->version();
155 Function: Get/set the version of the program used to
157 Returns : scalar representing the version
158 Args : none to get, OR string to set
163 my ($self,$value) = @_;
164 if (defined($value)) {
165 $self->{'_version'} = $value;
167 return $self->{'_version'};
170 =head2 modification_user
172 Title : modification_user
173 Usage : my $modification_user = $map->modification_user();
174 Function: Get/set the name of the user who last modified this map
175 Returns : scalar representing the username
176 Args : none to get, OR string to set
180 sub modification_user
{
181 my ($self,$value) = @_;
182 if (defined($value)) {
183 $self->{'_modification_user'} = $value;
185 return $self->{'_modification_user'};
191 Usage : $map->group_type($grptype);
192 my $grptype = $map->group_type();
193 Function: Get/set the group type of this map
194 Returns : scalar representing the group type
195 Args : none to get, OR string to set
200 my ($self,$value) = @_;
201 if (defined($value)) {
202 $self->{'_grouptype'} = $value;
204 return $self->{'_grouptype'};
210 Usage : $map->group_abbr($grpabbr);
211 my $grpabbr = $map->group_abbr();
212 Function: get/set the group abbrev of this map
213 Returns : string representing the group abbrev
214 Args : none to get, OR string to set
219 my ($self,$value) = @_;
220 if (defined($value)) {
221 $self->{'_groupabbr'} = $value;
223 return $self->{'_groupabbr'};
229 Usage : my $core_exists = $map->core_exists();
230 Function: Get/set if the FPC file is accompanied by COR file
232 Args : none to get, OR 1|0 to set
237 my ($self,$value) = @_;
238 if (defined($value)) {
239 $self->{'_corexists'} = $value ?
1 : 0;
241 return $self->{'_corexists'};
247 Usage : my @clones = $map->each_cloneid();
248 Function: returns an array of clone names
249 Returns : list of clone names
256 return keys %{$self->{'_clones'}};
262 Usage : my $cloneobj = $map->get_cloneobj('CLONEA');
263 Function: returns an object of the clone given in the argument
264 Returns : object of the clone
265 Args : scalar representing the clone name
270 my ($self,$clone) = @_;
272 return 0 if(!defined($clone));
273 return if($clone eq "");
274 return if(!exists($self->{'_clones'}{$clone}));
276 my ($type,$contig,$bands,$gel,$group,$remark,$fp_number);
277 my ($sequence_type,$sequence_status,$fpc_remark,@amatch,@pmatch,@ematch,
278 $startrange,$endrange);
279 my %clones = %{$self->{'_clones'}{$clone}};
282 if (ref($clones{'clone'}) eq 'Bio::Map::Clone') {
283 return $clones{'clone'};
286 $type = $clones{'type'} if (exists($clones{'type'}));
287 @markers = (keys %{$clones{'markers'}}) if (exists($clones{'markers'}));
288 $contig = $clones{'contig'} if (exists($clones{'contig'}));
289 $bands = $clones{'bands'} if (exists($clones{'bands'}));
290 $gel = $clones{'gel'} if (exists($clones{'gel'}));
291 $group = $clones{'group'} if (exists($clones{'group'}));
292 $remark = $clones{'remark'} if (exists($clones{'remark'}));
294 $fp_number = $clones{'fp_number'} if (exists($clones{'fp_number'}));
295 $fpc_remark = $clones{'fpc_remark'} if (exists($clones{'fpc_remark'}));
297 $sequence_type = $clones{'sequence_type'}
298 if (exists($clones{'sequence_type'}));
299 $sequence_status = $clones{'sequence_status'}
300 if (exists($clones{'sequence_status'} ));
302 @amatch = (keys %{$clones{'matcha'}}) if (exists($clones{'matcha'}));
303 @ematch = (keys %{$clones{'matche'}}) if (exists($clones{'matche'}));
304 @pmatch = (keys %{$clones{'matchp'}}) if (exists($clones{'matchp'}));
306 $startrange = $clones{'range'}{'start'}
307 if (exists($clones{'range'}{'start'}));
308 $endrange = $clones{'range'}{'end'}
309 if (exists($clones{'range'}{'end'}));
311 #*** why doesn't it call Bio::Map::Clone->new ? Seems dangerous...
312 my $cloneobj = bless( {
314 _markers
=> \
@markers,
321 _fpnumber
=> $fp_number,
322 _sequencetype
=> $sequence_type,
323 _sequencestatus
=> $sequence_status,
324 _fpcremark
=> $fpc_remark,
328 _range
=> Bio
::Range
->new(-start
=> $startrange,
330 }, 'Bio::Map::Clone');
332 $self->{'_clones'}{$clone}{'clone'} = $cloneobj;
338 Title : each_markerid
339 Usage : my @markers = $map->each_markerid();
340 Function: returns list of marker names
341 Returns : list of marker names
348 return keys (%{$self->{'_markers'}});
353 Title : get_markerobj
354 Usage : my $markerobj = $map->get_markerobj('MARKERA');
355 Function: returns an object of the marker given in the argument
356 Returns : object of the marker
357 Args : scalar representing the marker name
362 my ($self,$marker) = @_;
364 return 0 if(!defined($marker));
365 return if($marker eq "");
366 return if(!exists($self->{'_markers'}{$marker}));
368 my ($global,$framework,$group,$anchor,$remark,$type,$linkage,$subgroup);
369 my %mkr = %{$self->{'_markers'}{$marker}};
371 return $mkr{'marker'} if (ref($mkr{'marker'}) eq 'Bio::Map::FPCMarker');
373 $type = $mkr{'type'} if(exists($mkr{'type'}));
374 $global = $mkr{'global'} if(exists($mkr{'global'} ));
375 $framework = $mkr{'framework'} if(exists($mkr{'framework'}));
376 $anchor = $mkr{'anchor'} if(exists($mkr{'anchor'}));
377 $group = $mkr{'group'} if(exists($mkr{'group'}));
378 $subgroup = $mkr{'subgroup'} if(exists($mkr{'subgroup'}));
379 $remark = $mkr{'remark'} if(exists($mkr{'remark'}));
381 my %clones = %{$mkr{'clones'}};
382 my %contigs = %{$mkr{'contigs'}};
384 my %markerpos = %{$mkr{'posincontig'}} if(exists($mkr{'posincontig'}));
386 #*** why doesn't it call Bio::Map::FPCMarker->new ? Seems dangerous...
387 my $markerobj = bless( {
391 _frame
=> $framework,
393 _subgroup
=> $subgroup,
397 _contigs
=> \
%contigs,
398 _position
=> \
%markerpos,
399 }, 'Bio::Map::FPCMarker');
401 $self->{'_markers'}{$marker}{'marker'} = $markerobj;
407 Title : each_contigid
408 Usage : my @contigs = $map->each_contigid();
409 Function: returns a list of contigs (numbers)
410 Returns : list of contigs
417 return keys (%{$self->{'_contigs'}});
422 Title : get_contigobj
423 Usage : my $contigobj = $map->get_contigobj('CONTIG1');
424 Function: returns an object of the contig given in the argument
425 Returns : object of the contig
426 Args : scalar representing the contig number
431 my ($self,$contig) = @_;
433 return 0 if(!defined($contig));
434 return if($contig eq "");
435 return if(!exists($self->{'_contigs'}{$contig}));
437 my ($group,$anchor,$uremark,$tremark,$cremark,$startrange,$endrange,
439 my %ctg = %{$self->{'_contigs'}{$contig}};
440 my (%position, %pos);
442 return $ctg{'contig'} if (ref($ctg{'contig'}) eq 'Bio::Map::Contig');
444 $group = $ctg{'group'} if (exists($ctg{'group'}));
445 $subgroup = $ctg{'subgroup'} if (exists($ctg{'subgroup'}));
446 $anchor = $ctg{'anchor'} if (exists($ctg{'anchor'}));
447 $cremark = $ctg{'chr_remark'} if (exists($ctg{'chr_remark'}));
448 $uremark = $ctg{'usr_remark'} if (exists($ctg{'usr_remark'}));
449 $tremark = $ctg{'trace_remark'} if (exists($ctg{'trace_remark'}));
451 $startrange = $ctg{'range'}{'start'}
452 if (exists($ctg{'range'}{'start'}));
453 $endrange = $ctg{'range'}{'end'}
454 if (exists($ctg{'range'}{'end'}));
456 my %clones = %{$ctg{'clones'}} if (exists($ctg{'clones'}));
457 my %markers = %{$ctg{'markers'}} if (exists($ctg{'markers'}));
459 my $pos = $ctg{'position'};
461 #*** why doesn't it call Bio::Map::Contig->new ? Seems dangerous...
462 my $contigobj = bless( {
464 _subgroup
=> $subgroup,
466 _markers
=> \
%markers,
469 _cremark
=> $cremark,
470 _uremark
=> $uremark,
471 _tremark
=> $tremark,
473 _range
=> Bio
::Range
->new(-start
=> $startrange,
475 }, 'Bio::Map::Contig');
477 $self->{'_contigs'}{$contig}{'contig'} = $contigobj;
481 =head2 matching_bands
483 Title : matching_bands
484 Usage : $self->matching_bands('cloneA','cloneB',[$tol]);
485 Function: given two clones [and tolerence], this method calculates how many
486 matching bands do they have.
487 (this method is ported directly from FPC)
488 Returns : scalar representing the number of matching bands
489 Args : names of the clones ('cloneA', 'cloneB') [Default tolerence=7]
494 my($self,$cloneA,$cloneB,$tol) = @_;
495 my($lstart,$kband,$match,$diff,$i,$j);
497 return 0 if(!defined($cloneA) || !defined($cloneB) ||
498 !($self->core_exists()));
500 $tol = 7 if (!defined($tol));
502 my %_clones = %{$self->{'_clones'}};
504 my @bandsA = @
{$_clones{$cloneA}{'bands'}};
505 my @bandsB = @
{$_clones{$cloneB}{'bands'}};
510 for ($i=0; $i<scalar(@bandsA);$i++) {
511 $kband = $bandsA[$i];
512 for ($j = $lstart; $j<scalar(@bandsB); $j++) {
513 $diff = $kband - $bandsB[$j];
514 if (abs($diff) <= $tol ) {
528 =head2 coincidence_score
530 Title : coincidence_score
531 Usage : $self->coincidence_score('cloneA','cloneB'[,$tol,$gellen]);
532 Function: given two clones [,tolerence and gellen], this method calculates
533 the Sulston Coincidence score.
534 (this method is ported directly from FPC)
535 Returns : scalar representing the Sulston coincidence score.
536 Args : names of the clones ('cloneA', 'cloneB')
537 [Default tol=7 gellen=3300.0]
541 sub coincidence_score
{
542 my($self,$cloneA,$cloneB,$tol,$gellen) = @_;
544 return 0 if(!defined($cloneA) || !defined($cloneB) ||
545 !($self->core_exists()));
547 my %_clones = %{$self->{'_clones'}};
549 my $numbandsA = scalar(@
{$_clones{$cloneA}{'bands'}});
550 my $numbandsB = scalar(@
{$_clones{$cloneB}{'bands'}});
552 my ($nL,$nH,$m,$i,$psmn,$pp,$pa,$pb,$t,$c,$a,$n);
556 $gellen = 3300.0 if (!defined($gellen));
557 $tol = 7 if (!defined($tol));
559 if ($numbandsA > $numbandsB) {
568 $m = $self->matching_bands($cloneA, $cloneB,$tol);
572 for ($i=2; $i<=$nL; $i++) {
573 $logfact[$i] = $logfact[$i - 1] + log($i);
576 $psmn = 1.0 - ((2*$tol)/$gellen);
583 for ($n = $m; $n <= $nL; $n++) {
584 $c = $logfact[$nL] - $logfact[$nL - $n] - $logfact[$n];
585 $a = exp($c + ($n * $pb) + (($nL - $n) * $pa));
589 $score = sprintf("%.e",$t);
593 =head2 print_contiglist
595 Title : print_contiglist
596 Usage : $map->print_contiglist([showall]); #[Default 0]
597 Function: prints the list of contigs, markers that hit the contig, the
598 global position and whether the marker is a placement (P) or
599 a Framework (F) marker.
601 Args : [showall] [Default 0], 1 includes all the discrepant markers
605 sub print_contiglist
{
606 my ($self,$showall) = @_;
609 $showall = 0 if (!defined($showall));
610 my %_contigs = %{$self->{'_contigs'}};
611 my %_markers = %{$self->{'_markers'}};
612 my %_clones = %{$self->{'_clones'}};
614 my @contigs = $self->each_contigid();
615 my @sortedcontigs = sort {$a <=> $b } @contigs;
617 print "\n\nContig List\n\n";
618 foreach my $contig (@sortedcontigs) {
622 my $ctgAnchor = $_contigs{$contig}{'anchor'};
623 my $ctgGroup = $_contigs{$contig}{'group'};
625 my @mkr = keys ( %{$_contigs{$contig}{'markers'}} );
627 foreach my $marker (@mkr) {
628 my $mrkGroup = $_markers{$marker}{'group'};
629 my $mrkGlobal = $_markers{$marker}{'global'};
630 my $mrkFramework = $_markers{$marker}{'framework'};
631 my $mrkAnchor = $_markers{$marker}{'anchor'};
633 if($ctgGroup =~ /\d+|\w/ && $ctgGroup != 0) {
634 if ($mrkGroup eq $ctgGroup) {
635 if ($mrkFramework == 0) {
636 $pos = $mrkGlobal."P";
639 $pos = $mrkGlobal."F";
641 $list{$marker} = $pos;
643 elsif ($showall == 1) {
644 my $chr = $self->group_abbr().$mrkGroup;
645 $alist{$marker} = $chr;
648 elsif ($showall == 1 && $ctgGroup !~ /\d+/) {
649 my $chr = $self->group_abbr().$mrkGroup;
650 $alist{$marker} = $chr;
655 $chr = $self->group_abbr().$ctgGroup if ($ctgGroup =~ /\d+|\w/);
657 if ($showall == 1 ) {
659 print " ctg$contig ", $chr, " "
660 if ($_contigs{$contig}{'group'} !~ /\d+|\w/);
662 elsif ($ctgGroup =~ /\d+|\w/ && $ctgGroup ne 0){
663 print " ctg",$contig, " ",$chr, " ";
666 while (my ($k,$v) = each %list) {
670 print "\n" if ($showall == 0 && $ctgGroup =~ /\d+|\w/ &&
674 while (my ($k,$v) = each %alist) {
682 =head2 print_markerlist
684 Title : print_markerlist
685 Usage : $map->print_markerlist();
686 Function : prints the marker list; contig and corresponding number of
687 clones for each marker.
693 sub print_markerlist
{
696 my %_contigs = %{$self->{'_contigs'}};
697 my %_markers = %{$self->{'_markers'}};
698 my %_clones = %{$self->{'_clones'}};
700 print "Marker List\n\n";
702 foreach my $marker ($self->each_markerid()) {
703 print " ",$marker, " ";
706 my %mclones = %{$_markers{$marker}{'clones'}};
708 foreach my $clone (%mclones) {
709 if (exists($_clones{$clone}{'contig'}) ) {
710 my $ctg = $_clones{$clone}{'contig'};
712 if (exists($list{$ctg})) {
713 my $clonehits = $list{$ctg};
715 $list{$ctg} = $clonehits;
722 while (my ($k,$v) = each %list) {
729 =head2 print_gffstyle
731 Title : print_gffstyle
732 Usage : $map->print_gffstyle([style]);
733 Function : prints GFF; either Contigwise (default) or Groupwise
735 Args : [style] default = 0 contigwise, else
736 1 groupwise (chromosome-wise).
741 my ($self,$style) = @_;
743 $style = 0 if(!defined($style));
745 my %_contigs = %{$self->{'_contigs'}};
746 my %_markers = %{$self->{'_markers'}};
747 my %_clones = %{$self->{'_clones'}};
750 my ($depth, $save_depth);
757 # Calculate the position for the marker in the contig
759 my @contigs = $self->each_contigid();
760 my @sortedcontigs = sort {$a <=> $b } @contigs;
767 foreach my $contig (@sortedcontigs) {
768 if($_contigs{$contig}{'range'} ) {
769 $offset = $_contigs{$contig}{'range'}{'start'};
772 $offset = $offset * -1;
773 $gffcontigs{$contig}{'start'} = 1;
774 $gffcontigs{$contig}{'end'} =
775 ($_contigs{$contig}{'range'}{'end'} +
776 $offset ) * $basepair + 1;
780 $gffcontigs{$contig}{'start'} =
781 $_contigs{$contig}{'range'}{'start'} * $basepair;
782 $gffcontigs{$contig}{'end'} =
783 $_contigs{$contig}{'range'}{'end'} * $basepair;
787 $gffcontigs{$contig}{'start'} = 1;
788 $gffcontigs{$contig}{'end'} = 1;
791 my @clones = keys %{$_contigs{$contig}{'clones'}};
792 foreach my $clone (@clones) {
793 if(exists ($_clones{$clone}{'range'}) ) {
794 my $gffclone = $clone;
796 $gffclone =~ s/sd1$//;
798 $gffclones{$gffclone}{'start'} =
799 (($_clones{$clone}{'range'}{'start'} + $offset) *
802 $gffclones{$gffclone}{'end'} =
803 (($_clones{$clone}{'range'}{'end'}
804 + $offset) * $basepair + 1);
808 my %markers = %{$_clones{$clone}{'markers'}}
809 if (exists($_clones{$clone}{'markers'}));
811 while (my ($k,$v) = each %markers) {
812 $gffmarkers{$contig}{$k} =
813 ( ( $_clones{$clone}{'range'}{'start'} +
814 $_clones{$clone}{'range'}{'end'} ) / 2 ) *
821 my %markers = %{$_contigs{$contig}{'markers'}}
822 if (exists($_contigs{$contig}{'markers'}));
824 while (my ($k,$v) = each %markers) {
825 $gffmarkers{$contig}{$k} = ($v + $offset) * $basepair + 1;
831 foreach my $contig (@sortedcontigs) {
833 if(exists ($_contigs{$contig}{'range'} ) ) {
834 print join("\t","ctg$contig","assembly","contig",
835 $gffcontigs{$contig}{'start'},
836 $gffcontigs{$contig}{'end'},".",".",".",
837 "Sequence \"ctg$contig\"; Name \"ctg$contig\"\n"
841 my @clones = (keys %{$_contigs{$contig}{'clones'}} );
843 foreach my $clone (@clones) {
844 if(exists ($_clones{$clone}{'range'}) ) {
845 print join("\t","ctg$contig","FPC");
847 my $type = $_clones{$clone}{'type'};
849 if($clone =~ /sd1$/) {
853 print join ("\t","\t$type",$gffclones{$clone}{'start'},
854 $gffclones{$clone}{'end'},".",".",".",
855 "$type \"$clone\"; Name \"$clone\"");
857 my @markers = keys %{$_clones{$clone}{'markers'}};
858 print "; Marker_hit" if (scalar(@markers));
860 foreach my $mkr(@markers) {
861 if (exists($_markers{$mkr}{'framework'})) {
862 print " \"$mkr ",$_markers{$mkr}{'group'}," ",
863 $_markers{$mkr}{'global'},"\"";
866 print " \"$mkr 0 0\"";
869 print "; Contig_hit \"",$_clones{$clone}{'contig'},"\" "
870 if (defined($_clones{$clone}{'contig'}));
875 if (exists ($_contigs{$contig}{'markers'}) ) {
876 my %list = %{$_contigs{$contig}{'markers'}};
878 while (my ($k,$v) = each %list) {
879 print "ctg", $contig, "\tFPC\t";
880 my $position = $gffmarkers{$contig}{$k};
884 $type = "electronicmarker"
885 if ($_markers{$k}{'type'} eq "eMRK");
887 if( exists($_markers{$k}{'framework'})) {
888 $type = "frameworkmarker"
889 if($_markers{$k}{'framework'} == 1);
891 $type = "placementmarker"
892 if($_markers{$k}{'framework'} == 0);
895 print join ("\t","$type",$position,$position,".",".",
896 ".","$type \"$k\"; Name \"$k\"");
899 my @clones = keys %{$_markers{$k}{'clones'}};
901 foreach my $cl (@clones) {
902 push (@clonelist, $cl)
903 if($_clones{$cl}{'contig'} == $contig);
907 print("; Contig_hit
\"ctg
$contig - ",scalar(@clonelist),
908 "\" (@clonelist)\n");
915 my $margin = 2 * $basepair;
916 my $displacement = 0;
919 foreach my $contig (@sortedcontigs) {
921 my $chr = $_contigs{$contig}{'group'};
922 $chr = 0 if ($chr !~ /\d+|\w+/);
924 $recordchr->{group} = $chr;
925 $recordchr->{contig} = $contig;
926 $recordchr->{position} = $_contigs{$contig}{'position'};
928 push @grouplist, $recordchr;
931 my @chr = keys (%{$_groups{'group'}});
934 if ($self->group_type eq 'Chromosome') {
935 @sortedchr = sort { $a->{'group'} <=> $b->{'group'}
937 $a->{'contig'} <=> $b->{'contig'}
941 @sortedchr = sort { $a->{'group'} cmp $b->{'group'}
943 $a->{'contig'} cmp $b->{'contig'}
949 foreach my $chr (@sortedchr) {
950 my $chrname = $self->group_abbr().$chr->{'group'};
952 if ($lastchr eq -1 || $chr->{'group'} ne $lastchr ) {
953 $lastchr = $chr->{'group'} if ($lastchr eq -1);
956 # caluclate the end position of the contig
961 if ($chr->{contig} != 0) {
962 foreach my $ch (@sortedchr) {
963 if ($ch->{'group'} eq $chr->{'group'}) {
964 if($ch->{'contig'} != 0) {
965 my $ctg = $ch->{'contig'}
966 if($ch->{'contig'} != 0);
968 $chrend += $gffcontigs{$ctg}->{'end'};
973 $chrend += ($ctgcount-1) * $margin;
976 $chrend = $gffcontigs{'0'}->{'end'};
979 $chrname = $self->group_abbr()."ctg0
"
980 if ($chr->{'contig'} == 0);
982 print join ("\t", $chrname,"assembly
","Chromosome
",1,
983 "$chrend",".",".",".",
984 "Sequence
\"$chrname\"; Name
\"$chrname\"\n");
987 print join ("\t", $chrname,"assembly
","Chromosome
",1,
988 "$chrend",".",".",".",
989 "Sequence
\"$chrname\"; Name
\"$chrname\"\n")
990 if ($chr->{'group'} ne $lastchr && $chr->{'group'} eq 0 );
992 $lastchr = $chr->{'group'};
993 $lastchr = -1 if ($chr->{'contig'} == 0);
995 my $contig = $chr->{'contig'};
997 if(exists ($_contigs{$contig}{'range'} ) ) {
999 print join ("\t",$chrname, "FPC
","contig
",
1000 $gffcontigs{$contig}{'start'}+$displacement,
1001 $gffcontigs{$contig}{'end'}+$displacement,
1003 "contig
\"ctg
$contig\"; Name
\"ctg
$contig\"\n");
1006 my @clones = (keys %{$_contigs{$contig}{'clones'}} );
1007 foreach my $clone (@clones) {
1008 if(exists ($_clones{$clone}{'range'}) ) {
1009 print join ("\t",$chrname,"FPC
");
1010 my $type = $_clones{$clone}{'type'};
1012 if ($clone =~ /sd1$/) {
1014 $type = "sequenced
";
1017 print join ("\t","\t$type",$gffclones{$clone}{'start'}
1018 +$displacement,$gffclones{$clone}{'end'}
1019 +$displacement,".",".",".",
1020 "$type \"$clone\"; Name
\"$clone\"");
1022 my @markers = keys %{$_clones{$clone}{'markers'}};
1023 print "; Marker_hit
" if (scalar(@markers));
1025 foreach my $mkr(@markers) {
1026 if (exists($_markers{$mkr}{'framework'})) {
1027 print " \"$mkr ",$_markers{$mkr}{'group'}," ",
1028 $_markers{$mkr}{'global'},"\"";
1031 print (" \"$mkr 0 0\"");
1034 print "; Contig_hit
\"",$_clones{$clone}{'contig'},"\" "
1035 if (defined($_clones{$clone}{'contig'}));
1040 if (exists ($_contigs{$contig}{'markers'}) ) {
1041 my %list = %{$_contigs{$contig}{'markers'}};
1043 while (my ($k,$v) = each %list) {
1044 print join ("\t",$chrname,"FPC
");
1045 my $type = "marker
";
1047 $type = "electronicmarker
"
1048 if ($_markers{$k}{'type'} eq "eMRK
");
1050 if( exists($_markers{$k}{'framework'})) {
1051 $type = "frameworkmarker
"
1052 if($_markers{$k}{'framework'} == 1);
1054 $type = "placementmarker
"
1055 if($_markers{$k}{'framework'} == 0);
1058 print join ("\t","\t$type",$gffmarkers{$contig}{$k}
1059 + $displacement,$gffmarkers{$contig}{$k}
1060 + $displacement,".",".",".",
1061 "$type \"$k\"; Name
\"$k\"");
1064 my @clones = keys %{$_markers{$k}{'clones'}};
1066 foreach my $cl (@clones) {
1067 push (@clonelist, $cl)
1068 if($_clones{$cl}{'contig'} == $contig);
1072 print("; Contig_hit \"ctg$contig - ",
1073 scalar(@clonelist),"\" (@clonelist)\n");
1076 $displacement += $margin + $gffcontigs{$contig}{'end'};
1081 =head2 _calc_markerposition
1083 Title : _calc_markerposition
1084 Usage : $map->_calc_markerposition();
1085 Function: Calculates the position of the marker in the contig
1091 sub _calc_markerposition
{
1093 my %_contigs = %{$self->{'_contigs'}};
1094 my %_markers = %{$self->{'_markers'}};
1095 my %_clones = %{$self->{'_clones'}};
1098 my ($depth, $save_depth);
1105 # Calculate the position for the marker in the contig
1107 my @contigs = $self->each_contigid();
1108 my @sortedcontigs = sort {$a <=> $b } @contigs;
1113 foreach my $marker ($self->each_markerid()) {
1114 my (@ctgmarker, @sortedctgmarker);
1116 my @clones = (keys %{$_markers{$marker}{'clones'}})
1117 if (exists ($_markers{$marker}{'clones'} ));
1119 foreach my $clone (@clones) {
1121 $record->{contig
} = $_clones{$clone}{'contig'};
1122 $record->{start
} = $_clones{$clone}{'range'}{'start'};
1123 $record->{end
} = $_clones{$clone}{'range'}{'end'};
1124 push @ctgmarker,$record;
1127 # sorting by contig and left position
1128 @sortedctgmarker = sort { $a->{'contig'} <=> $b->{'contig'}
1130 $b->{'start'} <=> $a->{'start'}
1135 for ($i=0; $i < scalar(@sortedctgmarker); $i++) {
1136 if ($ctg != $sortedctgmarker[$i]->{'contig'}) {
1138 $ctg = $sortedctgmarker[$i]->{'contig'};
1141 if ($depth > $save_depth){
1142 $pos = ($x + $y) >> 1;
1143 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1144 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1148 $ctg = $sortedctgmarker[$i]->{'contig'};
1149 $x = $sortedctgmarker[$i]->{'start'};
1150 $y = $sortedctgmarker[$i]->{'end'};
1153 $pos = ($x + $y) >> 1;
1154 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1155 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1157 $depth = $save_depth = 1;
1159 elsif ($sortedctgmarker[$i] <= $y) {
1160 $stack[$depth++] = $sortedctgmarker[$i]->{'end'};
1162 if ($x < $sortedctgmarker[$i]->{'start'} ) {
1163 $x = $sortedctgmarker[$i]->{'start'};
1166 if ($y > $sortedctgmarker[$i]->{'end'}) {
1167 $y = $sortedctgmarker[$i]->{'end'};
1171 if ($depth > $save_depth) {
1172 $save_depth = $depth;
1173 $pos = ($x + $y) >> 1;
1174 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1175 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1178 $x = $sortedctgmarker[$i]->{'start'};
1179 $y = $sortedctgmarker[$i]->{'end'};
1180 $stack[$depth++] = $y;
1182 for($j=-1, $k=0, $s=0; $s<$depth; $s++) {
1183 if ($stack[$s] <$x) {
1185 $j = $s if ($j == -1);
1190 $y = $stack[$s] if ($y > $stack[$s]);
1191 if ($stack[$j] == -1) {
1192 $stack[$j] = $stack[$s];
1194 while ($stack[$j] != -1) {$j++;}
1203 if ($depth > $save_depth) {
1204 $pos = ($x + $y) >> 1;
1205 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1206 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1212 =head2 _calc_contigposition
1214 Title : _calc_contigposition
1215 Usage : $map->_calc_contigposition();
1216 Function: calculates the position of the contig in the group
1222 sub _calc_contigposition
{
1225 my %_contigs = %{$self->{'_contigs'}};
1226 my %_markers = %{$self->{'_markers'}};
1227 my %_clones = %{$self->{'_clones'}};
1229 my @contigs = $self->each_contigid();
1230 my @sortedcontigs = sort {$a <=> $b } @contigs;
1232 foreach my $contig (@sortedcontigs) {
1236 if (exists($_contigs{$contig}{'group'}) ) {
1238 my %weightedmarkers;
1239 my @mkrs = keys (%{$_contigs{$contig}{'markers'}})
1240 if (exists($_contigs{$contig}{'markers'})) ;
1242 my $chr = $_contigs{$contig}{'group'};
1243 $chr = 0 if ($_contigs{$contig}{'group'} =~ /\?/);
1245 foreach my $mkr (@mkrs) {
1246 if (exists($_markers{$mkr}{'group'})) {
1247 if ( $_markers{$mkr}{'group'} == $chr ) {
1248 my @mkrclones = keys( %{$_markers{$mkr}{'clones'}});
1249 my $clonescount = 0;
1250 foreach my $clone (@mkrclones) {
1252 if ($_clones{$clone}{'contig'} == $contig);
1254 $weightedmarkers{$_markers{$mkr}{'global'}} =
1260 my $weightedctgsum = 0;
1263 while (my ($mpos,$hits) = each %weightedmarkers) {
1264 $weightedctgsum += ($mpos * $hits);
1265 $totalhits += $hits;
1268 $position = sprintf("%.2f",$weightedctgsum / $totalhits)
1269 if ($totalhits != 0);
1271 $_contigs{$contig}{'position'} = $position;
1276 =head2 _calc_contiggroup
1278 Title : _calc_contiggroup
1279 Usage : $map->_calc_contiggroup();
1280 Function: calculates the group of the contig
1286 sub _calc_contiggroup
{
1288 my %_contig = %{$self->{'_contigs'}};
1289 my @contigs = $self->each_contigid();
1291 foreach my $ctg (@contigs) {
1292 my $chr = floor
($ctg/1000);
1293 $_contig{$ctg}{'group'} = $chr;
1297 =head2 _setI<E<lt>TypeE<gt>>Ref
1299 Title : _set<Type>Ref
1300 Usage : These are used for initializing the reference of the hash in
1301 Bio::MapIO (fpc.pm) to the corresponding hash in Bio::Map
1302 (physical.pm). Should be used only from Bio::MapIO System.
1303 $map->setCloneRef(\%_clones);
1304 $map->setMarkerRef(\%_markers);
1305 $map->setContigRef(\%_contigs);
1306 Function: sets the hash references to the corresponding hashes
1308 Args : reference of the hash.
1313 my ($self, $ref) = @_;
1314 %{$self->{'_clones'}} = %{$ref};
1318 my ($self, $ref) = @_;
1319 %{$self->{'_markers'}} = %{$ref};
1323 my ($self, $ref) = @_;
1324 %{$self->{'_contigs'}} = %{$ref};