2 # BioPerl module for Bio::Map::Physical
4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
6 # Cared for by Sendu Bala <bix@sendu.me.uk>
10 # You may distribute this module under the same terms as perl itself
12 # POD documentation - main docs before the code
16 Bio::Map::Physical - A class for handling a Physical Map (such as FPC)
22 # accquire a Bio::Map::Physical using Bio::MapIO::fpc
23 my $mapio = Bio::MapIO->new(-format => "fpc",-file => "rice.fpc",
26 my $physical = $mapio->next_map();
28 # get all the markers ids
29 foreach my $marker ( $physical->each_markerid() ) {
30 print "Marker $marker\n";
32 # acquire the marker object using Bio::Map::FPCMarker
33 my $markerobj = $physical->get_markerobj($marker);
35 # get all the clones hit by this marker
36 foreach my $clone ($markerobj->each_cloneid() ) {
43 This class is basically a continer class for a collection of Contig maps and
44 other physical map information.
46 Bio::Map::Physical has been tailored to work for FPC physical maps, but
47 could probably be used for others as well (with the appropriate MapIO
50 This class also has some methods with specific functionalities:
52 print_gffstyle() : Generates GFF; either Contigwise[Default] or
55 print_contiglist() : Prints the list of Contigs, markers that hit the
56 contig, the global position and whether the marker
57 is a placement (<P>) or a Framework (<F>) marker.
59 print_markerlist() : Prints the markers list; contig and corresponding
62 matching_bands() : Given two clones [and tolerence], this method
63 calculates how many matching bands do they have.
65 coincidence_score() : Given two clones [,tolerence and gellen], this
66 method calculates the Sulston Coincidence score.
68 For faster access and better optimization, the data is stored internally in
69 hashes. The corresponding objects are created on request.
75 User feedback is an integral part of the evolution of this and other
76 Bioperl modules. Send your comments and suggestions preferably to
77 the Bioperl mailing list. Your participation is much appreciated.
79 bioperl-l@bioperl.org - General discussion
80 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
84 Please direct usage questions or support issues to the mailing list:
86 I<bioperl-l@bioperl.org>
88 rather than to the module maintainer directly. Many experienced and
89 reponsive experts will be able look at the problem and quickly
90 address it. Please include a thorough description of the problem
91 with code and data examples if at all possible.
95 Report bugs to the Bioperl bug tracking system to help us keep track
96 of the bugs and their resolution. Bug reports can be submitted via the
99 https://redmine.open-bio.org/projects/bioperl/
101 =head1 AUTHOR - Gaurav Gupta
103 Email gaurav@genome.arizona.edu
107 Sendu Bala bix@sendu.me.uk
109 =head1 PROJECT LEADERS
111 Jamie Hatfield jamie@genome.arizona.edu
112 Dr. Cari Soderlund cari@genome.arizona.edu
114 =head1 PROJECT DESCRIPTION
116 The project was done in Arizona Genomics Computational Laboratory (AGCoL)
117 at University of Arizona.
119 This work was funded by USDA-IFAFS grant #11180 titled "Web Resources for
120 the Computation and Display of Physical Mapping Data".
122 For more information on this project, please refer:
123 http://www.genome.arizona.edu
127 The rest of the documentation details each of the object methods.
128 Internal methods are usually preceded with a _
132 # Let the code begin...
134 package Bio
::Map
::Physical
;
135 use vars
qw($MAPCOUNT);
140 use Bio::Map::Contig;
141 use Bio::Map::FPCMarker;
143 use base qw(Bio::Map::SimpleMap);
144 BEGIN { $MAPCOUNT = 1; }
146 =head1 Access Methods
148 These methods let you get and set the member variables
153 Usage : my $version = $map->version();
154 Function: Get/set the version of the program used to
156 Returns : scalar representing the version
157 Args : none to get, OR string to set
162 my ($self,$value) = @_;
163 if (defined($value)) {
164 $self->{'_version'} = $value;
166 return $self->{'_version'};
169 =head2 modification_user
171 Title : modification_user
172 Usage : my $modification_user = $map->modification_user();
173 Function: Get/set the name of the user who last modified this map
174 Returns : scalar representing the username
175 Args : none to get, OR string to set
179 sub modification_user
{
180 my ($self,$value) = @_;
181 if (defined($value)) {
182 $self->{'_modification_user'} = $value;
184 return $self->{'_modification_user'};
190 Usage : $map->group_type($grptype);
191 my $grptype = $map->group_type();
192 Function: Get/set the group type of this map
193 Returns : scalar representing the group type
194 Args : none to get, OR string to set
199 my ($self,$value) = @_;
200 if (defined($value)) {
201 $self->{'_grouptype'} = $value;
203 return $self->{'_grouptype'};
209 Usage : $map->group_abbr($grpabbr);
210 my $grpabbr = $map->group_abbr();
211 Function: get/set the group abbrev of this map
212 Returns : string representing the group abbrev
213 Args : none to get, OR string to set
218 my ($self,$value) = @_;
219 if (defined($value)) {
220 $self->{'_groupabbr'} = $value;
222 return $self->{'_groupabbr'};
228 Usage : my $core_exists = $map->core_exists();
229 Function: Get/set if the FPC file is accompanied by COR file
231 Args : none to get, OR 1|0 to set
236 my ($self,$value) = @_;
237 if (defined($value)) {
238 $self->{'_corexists'} = $value ?
1 : 0;
240 return $self->{'_corexists'};
246 Usage : my @clones = $map->each_cloneid();
247 Function: returns an array of clone names
248 Returns : list of clone names
255 return keys %{$self->{'_clones'}};
261 Usage : my $cloneobj = $map->get_cloneobj('CLONEA');
262 Function: returns an object of the clone given in the argument
263 Returns : object of the clone
264 Args : scalar representing the clone name
269 my ($self,$clone) = @_;
271 return 0 if(!defined($clone));
272 return if($clone eq "");
273 return if(!exists($self->{'_clones'}{$clone}));
275 my ($type,$contig,$bands,$gel,$group,$remark,$fp_number);
276 my ($sequence_type,$sequence_status,$fpc_remark,@amatch,@pmatch,@ematch,
277 $startrange,$endrange);
278 my %clones = %{$self->{'_clones'}{$clone}};
281 if (ref($clones{'clone'}) eq 'Bio::Map::Clone') {
282 return $clones{'clone'};
285 $type = $clones{'type'} if (exists($clones{'type'}));
286 @markers = (keys %{$clones{'markers'}}) if (exists($clones{'markers'}));
287 $contig = $clones{'contig'} if (exists($clones{'contig'}));
288 $bands = $clones{'bands'} if (exists($clones{'bands'}));
289 $gel = $clones{'gel'} if (exists($clones{'gel'}));
290 $group = $clones{'group'} if (exists($clones{'group'}));
291 $remark = $clones{'remark'} if (exists($clones{'remark'}));
293 $fp_number = $clones{'fp_number'} if (exists($clones{'fp_number'}));
294 $fpc_remark = $clones{'fpc_remark'} if (exists($clones{'fpc_remark'}));
296 $sequence_type = $clones{'sequence_type'}
297 if (exists($clones{'sequence_type'}));
298 $sequence_status = $clones{'sequence_status'}
299 if (exists($clones{'sequence_status'} ));
301 @amatch = (keys %{$clones{'matcha'}}) if (exists($clones{'matcha'}));
302 @ematch = (keys %{$clones{'matche'}}) if (exists($clones{'matche'}));
303 @pmatch = (keys %{$clones{'matchp'}}) if (exists($clones{'matchp'}));
305 $startrange = $clones{'range'}{'start'}
306 if (exists($clones{'range'}{'start'}));
307 $endrange = $clones{'range'}{'end'}
308 if (exists($clones{'range'}{'end'}));
310 #*** why doesn't it call Bio::Map::Clone->new ? Seems dangerous...
311 my $cloneobj = bless( {
313 _markers
=> \
@markers,
320 _fpnumber
=> $fp_number,
321 _sequencetype
=> $sequence_type,
322 _sequencestatus
=> $sequence_status,
323 _fpcremark
=> $fpc_remark,
327 _range
=> Bio
::Range
->new(-start
=> $startrange,
329 }, 'Bio::Map::Clone');
331 $self->{'_clones'}{$clone}{'clone'} = $cloneobj;
337 Title : each_markerid
338 Usage : my @markers = $map->each_markerid();
339 Function: returns list of marker names
340 Returns : list of marker names
347 return keys (%{$self->{'_markers'}});
352 Title : get_markerobj
353 Usage : my $markerobj = $map->get_markerobj('MARKERA');
354 Function: returns an object of the marker given in the argument
355 Returns : object of the marker
356 Args : scalar representing the marker name
361 my ($self,$marker) = @_;
363 return 0 if(!defined($marker));
364 return if($marker eq "");
365 return if(!exists($self->{'_markers'}{$marker}));
367 my ($global,$framework,$group,$anchor,$remark,$type,$linkage,$subgroup);
368 my %mkr = %{$self->{'_markers'}{$marker}};
370 return $mkr{'marker'} if (ref($mkr{'marker'}) eq 'Bio::Map::FPCMarker');
372 $type = $mkr{'type'} if(exists($mkr{'type'}));
373 $global = $mkr{'global'} if(exists($mkr{'global'} ));
374 $framework = $mkr{'framework'} if(exists($mkr{'framework'}));
375 $anchor = $mkr{'anchor'} if(exists($mkr{'anchor'}));
376 $group = $mkr{'group'} if(exists($mkr{'group'}));
377 $subgroup = $mkr{'subgroup'} if(exists($mkr{'subgroup'}));
378 $remark = $mkr{'remark'} if(exists($mkr{'remark'}));
380 my %clones = %{$mkr{'clones'}};
381 my %contigs = %{$mkr{'contigs'}};
383 my %markerpos = %{$mkr{'posincontig'}} if(exists($mkr{'posincontig'}));
385 #*** why doesn't it call Bio::Map::FPCMarker->new ? Seems dangerous...
386 my $markerobj = bless( {
390 _frame
=> $framework,
392 _subgroup
=> $subgroup,
396 _contigs
=> \
%contigs,
397 _position
=> \
%markerpos,
398 }, 'Bio::Map::FPCMarker');
400 $self->{'_markers'}{$marker}{'marker'} = $markerobj;
406 Title : each_contigid
407 Usage : my @contigs = $map->each_contigid();
408 Function: returns a list of contigs (numbers)
409 Returns : list of contigs
416 return keys (%{$self->{'_contigs'}});
421 Title : get_contigobj
422 Usage : my $contigobj = $map->get_contigobj('CONTIG1');
423 Function: returns an object of the contig given in the argument
424 Returns : object of the contig
425 Args : scalar representing the contig number
430 my ($self,$contig) = @_;
432 return 0 if(!defined($contig));
433 return if($contig eq "");
434 return if(!exists($self->{'_contigs'}{$contig}));
436 my ($group,$anchor,$uremark,$tremark,$cremark,$startrange,$endrange,
438 my %ctg = %{$self->{'_contigs'}{$contig}};
439 my (%position, %pos);
441 return $ctg{'contig'} if (ref($ctg{'contig'}) eq 'Bio::Map::Contig');
443 $group = $ctg{'group'} if (exists($ctg{'group'}));
444 $subgroup = $ctg{'subgroup'} if (exists($ctg{'subgroup'}));
445 $anchor = $ctg{'anchor'} if (exists($ctg{'anchor'}));
446 $cremark = $ctg{'chr_remark'} if (exists($ctg{'chr_remark'}));
447 $uremark = $ctg{'usr_remark'} if (exists($ctg{'usr_remark'}));
448 $tremark = $ctg{'trace_remark'} if (exists($ctg{'trace_remark'}));
450 $startrange = $ctg{'range'}{'start'}
451 if (exists($ctg{'range'}{'start'}));
452 $endrange = $ctg{'range'}{'end'}
453 if (exists($ctg{'range'}{'end'}));
455 my %clones = %{$ctg{'clones'}} if (exists($ctg{'clones'}));
456 my %markers = %{$ctg{'markers'}} if (exists($ctg{'markers'}));
458 my $pos = $ctg{'position'};
460 #*** why doesn't it call Bio::Map::Contig->new ? Seems dangerous...
461 my $contigobj = bless( {
463 _subgroup
=> $subgroup,
465 _markers
=> \
%markers,
468 _cremark
=> $cremark,
469 _uremark
=> $uremark,
470 _tremark
=> $tremark,
472 _range
=> Bio
::Range
->new(-start
=> $startrange,
474 }, 'Bio::Map::Contig');
476 $self->{'_contigs'}{$contig}{'contig'} = $contigobj;
480 =head2 matching_bands
482 Title : matching_bands
483 Usage : $self->matching_bands('cloneA','cloneB',[$tol]);
484 Function: given two clones [and tolerence], this method calculates how many
485 matching bands do they have.
486 (this method is ported directly from FPC)
487 Returns : scalar representing the number of matching bands
488 Args : names of the clones ('cloneA', 'cloneB') [Default tolerence=7]
493 my($self,$cloneA,$cloneB,$tol) = @_;
494 my($lstart,$kband,$match,$diff,$i,$j);
496 return 0 if(!defined($cloneA) || !defined($cloneB) ||
497 !($self->core_exists()));
499 $tol = 7 if (!defined($tol));
501 my %_clones = %{$self->{'_clones'}};
503 my @bandsA = @
{$_clones{$cloneA}{'bands'}};
504 my @bandsB = @
{$_clones{$cloneB}{'bands'}};
509 for ($i=0; $i<scalar(@bandsA);$i++) {
510 $kband = $bandsA[$i];
511 for ($j = $lstart; $j<scalar(@bandsB); $j++) {
512 $diff = $kband - $bandsB[$j];
513 if (abs($diff) <= $tol ) {
527 =head2 coincidence_score
529 Title : coincidence_score
530 Usage : $self->coincidence_score('cloneA','cloneB'[,$tol,$gellen]);
531 Function: given two clones [,tolerence and gellen], this method calculates
532 the Sulston Coincidence score.
533 (this method is ported directly from FPC)
534 Returns : scalar representing the Sulston coincidence score.
535 Args : names of the clones ('cloneA', 'cloneB')
536 [Default tol=7 gellen=3300.0]
540 sub coincidence_score
{
541 my($self,$cloneA,$cloneB,$tol,$gellen) = @_;
543 return 0 if(!defined($cloneA) || !defined($cloneB) ||
544 !($self->core_exists()));
546 my %_clones = %{$self->{'_clones'}};
548 my $numbandsA = scalar(@
{$_clones{$cloneA}{'bands'}});
549 my $numbandsB = scalar(@
{$_clones{$cloneB}{'bands'}});
551 my ($nL,$nH,$m,$i,$psmn,$pp,$pa,$pb,$t,$c,$a,$n);
555 $gellen = 3300.0 if (!defined($gellen));
556 $tol = 7 if (!defined($tol));
558 if ($numbandsA > $numbandsB) {
567 $m = $self->matching_bands($cloneA, $cloneB,$tol);
571 for ($i=2; $i<=$nL; $i++) {
572 $logfact[$i] = $logfact[$i - 1] + log($i);
575 $psmn = 1.0 - ((2*$tol)/$gellen);
582 for ($n = $m; $n <= $nL; $n++) {
583 $c = $logfact[$nL] - $logfact[$nL - $n] - $logfact[$n];
584 $a = exp($c + ($n * $pb) + (($nL - $n) * $pa));
588 $score = sprintf("%.e",$t);
592 =head2 print_contiglist
594 Title : print_contiglist
595 Usage : $map->print_contiglist([showall]); #[Default 0]
596 Function: prints the list of contigs, markers that hit the contig, the
597 global position and whether the marker is a placement (P) or
598 a Framework (F) marker.
600 Args : [showall] [Default 0], 1 includes all the discrepant markers
604 sub print_contiglist
{
605 my ($self,$showall) = @_;
608 $showall = 0 if (!defined($showall));
609 my %_contigs = %{$self->{'_contigs'}};
610 my %_markers = %{$self->{'_markers'}};
611 my %_clones = %{$self->{'_clones'}};
613 my @contigs = $self->each_contigid();
614 my @sortedcontigs = sort {$a <=> $b } @contigs;
616 print "\n\nContig List\n\n";
617 foreach my $contig (@sortedcontigs) {
621 my $ctgAnchor = $_contigs{$contig}{'anchor'};
622 my $ctgGroup = $_contigs{$contig}{'group'};
624 my @mkr = keys ( %{$_contigs{$contig}{'markers'}} );
626 foreach my $marker (@mkr) {
627 my $mrkGroup = $_markers{$marker}{'group'};
628 my $mrkGlobal = $_markers{$marker}{'global'};
629 my $mrkFramework = $_markers{$marker}{'framework'};
630 my $mrkAnchor = $_markers{$marker}{'anchor'};
632 if($ctgGroup =~ /\d+|\w/ && $ctgGroup != 0) {
633 if ($mrkGroup eq $ctgGroup) {
634 if ($mrkFramework == 0) {
635 $pos = $mrkGlobal."P";
638 $pos = $mrkGlobal."F";
640 $list{$marker} = $pos;
642 elsif ($showall == 1) {
643 my $chr = $self->group_abbr().$mrkGroup;
644 $alist{$marker} = $chr;
647 elsif ($showall == 1 && $ctgGroup !~ /\d+/) {
648 my $chr = $self->group_abbr().$mrkGroup;
649 $alist{$marker} = $chr;
654 $chr = $self->group_abbr().$ctgGroup if ($ctgGroup =~ /\d+|\w/);
656 if ($showall == 1 ) {
658 print " ctg$contig ", $chr, " "
659 if ($_contigs{$contig}{'group'} !~ /\d+|\w/);
661 elsif ($ctgGroup =~ /\d+|\w/ && $ctgGroup ne 0){
662 print " ctg",$contig, " ",$chr, " ";
665 while (my ($k,$v) = each %list) {
669 print "\n" if ($showall == 0 && $ctgGroup =~ /\d+|\w/ &&
673 while (my ($k,$v) = each %alist) {
681 =head2 print_markerlist
683 Title : print_markerlist
684 Usage : $map->print_markerlist();
685 Function : prints the marker list; contig and corresponding number of
686 clones for each marker.
692 sub print_markerlist
{
695 my %_contigs = %{$self->{'_contigs'}};
696 my %_markers = %{$self->{'_markers'}};
697 my %_clones = %{$self->{'_clones'}};
699 print "Marker List\n\n";
701 foreach my $marker ($self->each_markerid()) {
702 print " ",$marker, " ";
705 my %mclones = %{$_markers{$marker}{'clones'}};
707 foreach my $clone (%mclones) {
708 if (exists($_clones{$clone}{'contig'}) ) {
709 my $ctg = $_clones{$clone}{'contig'};
711 if (exists($list{$ctg})) {
712 my $clonehits = $list{$ctg};
714 $list{$ctg} = $clonehits;
721 while (my ($k,$v) = each %list) {
728 =head2 print_gffstyle
730 Title : print_gffstyle
731 Usage : $map->print_gffstyle([style]);
732 Function : prints GFF; either Contigwise (default) or Groupwise
734 Args : [style] default = 0 contigwise, else
735 1 groupwise (chromosome-wise).
740 my ($self,$style) = @_;
742 $style = 0 if(!defined($style));
744 my %_contigs = %{$self->{'_contigs'}};
745 my %_markers = %{$self->{'_markers'}};
746 my %_clones = %{$self->{'_clones'}};
749 my ($depth, $save_depth);
756 # Calculate the position for the marker in the contig
758 my @contigs = $self->each_contigid();
759 my @sortedcontigs = sort {$a <=> $b } @contigs;
766 foreach my $contig (@sortedcontigs) {
767 if($_contigs{$contig}{'range'} ) {
768 $offset = $_contigs{$contig}{'range'}{'start'};
771 $offset = $offset * -1;
772 $gffcontigs{$contig}{'start'} = 1;
773 $gffcontigs{$contig}{'end'} =
774 ($_contigs{$contig}{'range'}{'end'} +
775 $offset ) * $basepair + 1;
779 $gffcontigs{$contig}{'start'} =
780 $_contigs{$contig}{'range'}{'start'} * $basepair;
781 $gffcontigs{$contig}{'end'} =
782 $_contigs{$contig}{'range'}{'end'} * $basepair;
786 $gffcontigs{$contig}{'start'} = 1;
787 $gffcontigs{$contig}{'end'} = 1;
790 my @clones = keys %{$_contigs{$contig}{'clones'}};
791 foreach my $clone (@clones) {
792 if(exists ($_clones{$clone}{'range'}) ) {
793 my $gffclone = $clone;
795 $gffclone =~ s/sd1$//;
797 $gffclones{$gffclone}{'start'} =
798 (($_clones{$clone}{'range'}{'start'} + $offset) *
801 $gffclones{$gffclone}{'end'} =
802 (($_clones{$clone}{'range'}{'end'}
803 + $offset) * $basepair + 1);
807 my %markers = %{$_clones{$clone}{'markers'}}
808 if (exists($_clones{$clone}{'markers'}));
810 while (my ($k,$v) = each %markers) {
811 $gffmarkers{$contig}{$k} =
812 ( ( $_clones{$clone}{'range'}{'start'} +
813 $_clones{$clone}{'range'}{'end'} ) / 2 ) *
820 my %markers = %{$_contigs{$contig}{'markers'}}
821 if (exists($_contigs{$contig}{'markers'}));
823 while (my ($k,$v) = each %markers) {
824 $gffmarkers{$contig}{$k} = ($v + $offset) * $basepair + 1;
830 foreach my $contig (@sortedcontigs) {
832 if(exists ($_contigs{$contig}{'range'} ) ) {
833 print join("\t","ctg$contig","assembly","contig",
834 $gffcontigs{$contig}{'start'},
835 $gffcontigs{$contig}{'end'},".",".",".",
836 "Sequence \"ctg$contig\"; Name \"ctg$contig\"\n"
840 my @clones = (keys %{$_contigs{$contig}{'clones'}} );
842 foreach my $clone (@clones) {
843 if(exists ($_clones{$clone}{'range'}) ) {
844 print join("\t","ctg$contig","FPC");
846 my $type = $_clones{$clone}{'type'};
848 if($clone =~ /sd1$/) {
852 print join ("\t","\t$type",$gffclones{$clone}{'start'},
853 $gffclones{$clone}{'end'},".",".",".",
854 "$type \"$clone\"; Name \"$clone\"");
856 my @markers = keys %{$_clones{$clone}{'markers'}};
857 print "; Marker_hit" if (scalar(@markers));
859 foreach my $mkr(@markers) {
860 if (exists($_markers{$mkr}{'framework'})) {
861 print " \"$mkr ",$_markers{$mkr}{'group'}," ",
862 $_markers{$mkr}{'global'},"\"";
865 print " \"$mkr 0 0\"";
868 print "; Contig_hit \"",$_clones{$clone}{'contig'},"\" "
869 if (defined($_clones{$clone}{'contig'}));
874 if (exists ($_contigs{$contig}{'markers'}) ) {
875 my %list = %{$_contigs{$contig}{'markers'}};
877 while (my ($k,$v) = each %list) {
878 print "ctg", $contig, "\tFPC\t";
879 my $position = $gffmarkers{$contig}{$k};
883 $type = "electronicmarker"
884 if ($_markers{$k}{'type'} eq "eMRK");
886 if( exists($_markers{$k}{'framework'})) {
887 $type = "frameworkmarker"
888 if($_markers{$k}{'framework'} == 1);
890 $type = "placementmarker"
891 if($_markers{$k}{'framework'} == 0);
894 print join ("\t","$type",$position,$position,".",".",
895 ".","$type \"$k\"; Name \"$k\"");
898 my @clones = keys %{$_markers{$k}{'clones'}};
900 foreach my $cl (@clones) {
901 push (@clonelist, $cl)
902 if($_clones{$cl}{'contig'} == $contig);
906 print("; Contig_hit
\"ctg
$contig - ",scalar(@clonelist),
907 "\" (@clonelist)\n");
914 my $margin = 2 * $basepair;
915 my $displacement = 0;
918 foreach my $contig (@sortedcontigs) {
920 my $chr = $_contigs{$contig}{'group'};
921 $chr = 0 if ($chr !~ /\d+|\w+/);
923 $recordchr->{group} = $chr;
924 $recordchr->{contig} = $contig;
925 $recordchr->{position} = $_contigs{$contig}{'position'};
927 push @grouplist, $recordchr;
930 my @chr = keys (%{$_groups{'group'}});
933 if ($self->group_type eq 'Chromosome') {
934 @sortedchr = sort { $a->{'group'} <=> $b->{'group'}
936 $a->{'contig'} <=> $b->{'contig'}
940 @sortedchr = sort { $a->{'group'} cmp $b->{'group'}
942 $a->{'contig'} cmp $b->{'contig'}
948 foreach my $chr (@sortedchr) {
949 my $chrname = $self->group_abbr().$chr->{'group'};
951 if ($lastchr eq -1 || $chr->{'group'} ne $lastchr ) {
952 $lastchr = $chr->{'group'} if ($lastchr eq -1);
955 # caluclate the end position of the contig
960 if ($chr->{contig} != 0) {
961 foreach my $ch (@sortedchr) {
962 if ($ch->{'group'} eq $chr->{'group'}) {
963 if($ch->{'contig'} != 0) {
964 my $ctg = $ch->{'contig'}
965 if($ch->{'contig'} != 0);
967 $chrend += $gffcontigs{$ctg}->{'end'};
972 $chrend += ($ctgcount-1) * $margin;
975 $chrend = $gffcontigs{'0'}->{'end'};
978 $chrname = $self->group_abbr()."ctg0
"
979 if ($chr->{'contig'} == 0);
981 print join ("\t", $chrname,"assembly
","Chromosome
",1,
982 "$chrend",".",".",".",
983 "Sequence
\"$chrname\"; Name
\"$chrname\"\n");
986 print join ("\t", $chrname,"assembly
","Chromosome
",1,
987 "$chrend",".",".",".",
988 "Sequence
\"$chrname\"; Name
\"$chrname\"\n")
989 if ($chr->{'group'} ne $lastchr && $chr->{'group'} eq 0 );
991 $lastchr = $chr->{'group'};
992 $lastchr = -1 if ($chr->{'contig'} == 0);
994 my $contig = $chr->{'contig'};
996 if(exists ($_contigs{$contig}{'range'} ) ) {
998 print join ("\t",$chrname, "FPC
","contig
",
999 $gffcontigs{$contig}{'start'}+$displacement,
1000 $gffcontigs{$contig}{'end'}+$displacement,
1002 "contig
\"ctg
$contig\"; Name
\"ctg
$contig\"\n");
1005 my @clones = (keys %{$_contigs{$contig}{'clones'}} );
1006 foreach my $clone (@clones) {
1007 if(exists ($_clones{$clone}{'range'}) ) {
1008 print join ("\t",$chrname,"FPC
");
1009 my $type = $_clones{$clone}{'type'};
1011 if ($clone =~ /sd1$/) {
1013 $type = "sequenced
";
1016 print join ("\t","\t$type",$gffclones{$clone}{'start'}
1017 +$displacement,$gffclones{$clone}{'end'}
1018 +$displacement,".",".",".",
1019 "$type \"$clone\"; Name
\"$clone\"");
1021 my @markers = keys %{$_clones{$clone}{'markers'}};
1022 print "; Marker_hit
" if (scalar(@markers));
1024 foreach my $mkr(@markers) {
1025 if (exists($_markers{$mkr}{'framework'})) {
1026 print " \"$mkr ",$_markers{$mkr}{'group'}," ",
1027 $_markers{$mkr}{'global'},"\"";
1030 print (" \"$mkr 0 0\"");
1033 print "; Contig_hit
\"",$_clones{$clone}{'contig'},"\" "
1034 if (defined($_clones{$clone}{'contig'}));
1039 if (exists ($_contigs{$contig}{'markers'}) ) {
1040 my %list = %{$_contigs{$contig}{'markers'}};
1042 while (my ($k,$v) = each %list) {
1043 print join ("\t",$chrname,"FPC
");
1044 my $type = "marker
";
1046 $type = "electronicmarker
"
1047 if ($_markers{$k}{'type'} eq "eMRK
");
1049 if( exists($_markers{$k}{'framework'})) {
1050 $type = "frameworkmarker
"
1051 if($_markers{$k}{'framework'} == 1);
1053 $type = "placementmarker
"
1054 if($_markers{$k}{'framework'} == 0);
1057 print join ("\t","\t$type",$gffmarkers{$contig}{$k}
1058 + $displacement,$gffmarkers{$contig}{$k}
1059 + $displacement,".",".",".",
1060 "$type \"$k\"; Name
\"$k\"");
1063 my @clones = keys %{$_markers{$k}{'clones'}};
1065 foreach my $cl (@clones) {
1066 push (@clonelist, $cl)
1067 if($_clones{$cl}{'contig'} == $contig);
1071 print("; Contig_hit \"ctg$contig - ",
1072 scalar(@clonelist),"\" (@clonelist)\n");
1075 $displacement += $margin + $gffcontigs{$contig}{'end'};
1080 =head2 _calc_markerposition
1082 Title : _calc_markerposition
1083 Usage : $map->_calc_markerposition();
1084 Function: Calculates the position of the marker in the contig
1090 sub _calc_markerposition
{
1092 my %_contigs = %{$self->{'_contigs'}};
1093 my %_markers = %{$self->{'_markers'}};
1094 my %_clones = %{$self->{'_clones'}};
1097 my ($depth, $save_depth);
1104 # Calculate the position for the marker in the contig
1106 my @contigs = $self->each_contigid();
1107 my @sortedcontigs = sort {$a <=> $b } @contigs;
1112 foreach my $marker ($self->each_markerid()) {
1113 my (@ctgmarker, @sortedctgmarker);
1115 my @clones = (keys %{$_markers{$marker}{'clones'}})
1116 if (exists ($_markers{$marker}{'clones'} ));
1118 foreach my $clone (@clones) {
1120 $record->{contig
} = $_clones{$clone}{'contig'};
1121 $record->{start
} = $_clones{$clone}{'range'}{'start'};
1122 $record->{end
} = $_clones{$clone}{'range'}{'end'};
1123 push @ctgmarker,$record;
1126 # sorting by contig and left position
1127 @sortedctgmarker = sort { $a->{'contig'} <=> $b->{'contig'}
1129 $b->{'start'} <=> $a->{'start'}
1134 for ($i=0; $i < scalar(@sortedctgmarker); $i++) {
1135 if ($ctg != $sortedctgmarker[$i]->{'contig'}) {
1137 $ctg = $sortedctgmarker[$i]->{'contig'};
1140 if ($depth > $save_depth){
1141 $pos = ($x + $y) >> 1;
1142 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1143 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1147 $ctg = $sortedctgmarker[$i]->{'contig'};
1148 $x = $sortedctgmarker[$i]->{'start'};
1149 $y = $sortedctgmarker[$i]->{'end'};
1152 $pos = ($x + $y) >> 1;
1153 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1154 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1156 $depth = $save_depth = 1;
1158 elsif ($sortedctgmarker[$i] <= $y) {
1159 $stack[$depth++] = $sortedctgmarker[$i]->{'end'};
1161 if ($x < $sortedctgmarker[$i]->{'start'} ) {
1162 $x = $sortedctgmarker[$i]->{'start'};
1165 if ($y > $sortedctgmarker[$i]->{'end'}) {
1166 $y = $sortedctgmarker[$i]->{'end'};
1170 if ($depth > $save_depth) {
1171 $save_depth = $depth;
1172 $pos = ($x + $y) >> 1;
1173 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1174 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1177 $x = $sortedctgmarker[$i]->{'start'};
1178 $y = $sortedctgmarker[$i]->{'end'};
1179 $stack[$depth++] = $y;
1181 for($j=-1, $k=0, $s=0; $s<$depth; $s++) {
1182 if ($stack[$s] <$x) {
1184 $j = $s if ($j == -1);
1189 $y = $stack[$s] if ($y > $stack[$s]);
1190 if ($stack[$j] == -1) {
1191 $stack[$j] = $stack[$s];
1193 while ($stack[$j] != -1) {$j++;}
1202 if ($depth > $save_depth) {
1203 $pos = ($x + $y) >> 1;
1204 $_contigs{$ctg}{'markers'}{$marker} = $pos;
1205 $_markers{$marker}{'posincontig'}{$ctg} = $pos;
1211 =head2 _calc_contigposition
1213 Title : _calc_contigposition
1214 Usage : $map->_calc_contigposition();
1215 Function: calculates the position of the contig in the group
1221 sub _calc_contigposition
{
1224 my %_contigs = %{$self->{'_contigs'}};
1225 my %_markers = %{$self->{'_markers'}};
1226 my %_clones = %{$self->{'_clones'}};
1228 my @contigs = $self->each_contigid();
1229 my @sortedcontigs = sort {$a <=> $b } @contigs;
1231 foreach my $contig (@sortedcontigs) {
1235 if (exists($_contigs{$contig}{'group'}) ) {
1237 my %weightedmarkers;
1238 my @mkrs = keys (%{$_contigs{$contig}{'markers'}})
1239 if (exists($_contigs{$contig}{'markers'})) ;
1241 my $chr = $_contigs{$contig}{'group'};
1242 $chr = 0 if ($_contigs{$contig}{'group'} =~ /\?/);
1244 foreach my $mkr (@mkrs) {
1245 if (exists($_markers{$mkr}{'group'})) {
1246 if ( $_markers{$mkr}{'group'} == $chr ) {
1247 my @mkrclones = keys( %{$_markers{$mkr}{'clones'}});
1248 my $clonescount = 0;
1249 foreach my $clone (@mkrclones) {
1251 if ($_clones{$clone}{'contig'} == $contig);
1253 $weightedmarkers{$_markers{$mkr}{'global'}} =
1259 my $weightedctgsum = 0;
1262 while (my ($mpos,$hits) = each %weightedmarkers) {
1263 $weightedctgsum += ($mpos * $hits);
1264 $totalhits += $hits;
1267 $position = sprintf("%.2f",$weightedctgsum / $totalhits)
1268 if ($totalhits != 0);
1270 $_contigs{$contig}{'position'} = $position;
1275 =head2 _calc_contiggroup
1277 Title : _calc_contiggroup
1278 Usage : $map->_calc_contiggroup();
1279 Function: calculates the group of the contig
1285 sub _calc_contiggroup
{
1287 my %_contig = %{$self->{'_contigs'}};
1288 my @contigs = $self->each_contigid();
1290 foreach my $ctg (@contigs) {
1291 my $chr = floor
($ctg/1000);
1292 $_contig{$ctg}{'group'} = $chr;
1296 =head2 _setI<E<lt>TypeE<gt>>Ref
1298 Title : _set<Type>Ref
1299 Usage : These are used for initializing the reference of the hash in
1300 Bio::MapIO (fpc.pm) to the corresponding hash in Bio::Map
1301 (physical.pm). Should be used only from Bio::MapIO System.
1302 $map->setCloneRef(\%_clones);
1303 $map->setMarkerRef(\%_markers);
1304 $map->setContigRef(\%_contigs);
1305 Function: sets the hash references to the corresponding hashes
1307 Args : reference of the hash.
1312 my ($self, $ref) = @_;
1313 %{$self->{'_clones'}} = %{$ref};
1317 my ($self, $ref) = @_;
1318 %{$self->{'_markers'}} = %{$ref};
1322 my ($self, $ref) = @_;
1323 %{$self->{'_contigs'}} = %{$ref};