Bio/MapIO/fpc.pm

   1 # fpc.pm,v 1.2.2.1 2005/10/09 15:16:27 jason Exp
   2 #
   3 # BioPerl module for Bio::MapIO::fpc
   4 #
   5 # Cared for by Gaurav Gupta <gaurav@genome.arizona.edu>
   6 #
   7 # Copyright AGCoL
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::MapIO::fpc - A FPC Map reader
  16
  17 =head1 SYNOPSIS
  18
  19 # do not use this object directly it is accessed through the Bio::MapIO system
  20
  21     use Bio::MapIO;
  22
  23      -format  : specifies the format of the file format is "fpc",
  24      -file    : specifies the name of the .fpc file
  25      -readcor : boolean argument, indicating if .cor is to be read
  26                  or not. It looks for the .cor file in the same path
  27                  as .fpc file.
  28                  0 : doesn't read .cor file
  29                  1 : reads the .cor file
  30                  [default 0]
  31      -verbose : indicates the process of loading of fpc file
  32     my $mapio = Bio::MapIO->new(-format  => "fpc",
  33                                -file    => "rice.fpc",
  34                                -readcor => 0,
  35                                -verbose => 0);
  36
  37     my $map = $mapio->next_map();
  38
  39     foreach my $marker ( $map->each_markerid() ) {
  40          # loop through the markers associated with the map
  41          # likewise for contigs, clones, etc.
  42     }
  43
  44
  45 =head1 DESCRIPTION
  46
  47 This object contains code for parsing and processing FPC files and creating
  48 L<Bio::Map::Physical> object from it.
  49
  50 For faster access and better optimization, the data is stored internally in
  51 hashes. The corresponding objects are created on request.
  52
  53 We handle reading of the FPC ourselves, since MapIO module of Bioperl adds
  54 too much overhead.
  55
  56 =cut
  57
  58 # Let the code begin...
  59
  60 package Bio::MapIO::fpc;
  61 use strict;
  62 use POSIX;
  63
  64 use Bio::Map::Physical;
  65 use Bio::Map::Clone;
  66 use Bio::Map::Contig;
  67 use Bio::Map::FPCMarker;
  68 use Bio::Range;
  69
  70 use base qw(Bio::MapIO);
  71
  72 my $_readcor;
  73
  74 =head1 Initializer
  75
  76 =head2 _initialize
  77
  78  Title   : _initialize
  79  Usage   : called implicitly
  80  Function: calls the SUPER::_initialize
  81  Returns : nothing
  82  Args    : species, readcor
  83
  84 =cut
  85
  86 sub _initialize{
  87     my ($self,@args) = @_;
  88     my $species;
  89     $self->SUPER::_initialize(@args);
  90     ($species,$_readcor) = $self->_rearrange([qw(SPECIES READCOR)], @args);
  91     $_readcor = 0 unless (defined($_readcor));
  92 }
  93
  94 =head1 Access Methods
  95
  96 These methods let you get and set the member variables
  97
  98 =head2 next_map
  99
 100  Title   : next_map
 101  Usage   : my $fpcmap = $mapio->next_map();
 102  Function: gets the fpcmap from MapIO
 103  Returns : object of type L<Bio::Map::MapI>
 104  Args    : none
 105
 106 =cut
 107
 108 sub next_map{
 109
 110     my ($self) = @_;
 111
 112     my $line;
 113     my ($name,$fpcver,$moddate,$moduser,$contigcnt,$clonecnt,$markerscnt,
 114         $bandcnt,$marker,$seqclone);
 115     my ($corfile,$corindex,$BUFFER);
 116     my @cordata;
 117     my %fpcmarker;
 118     my ($contig, $contigNumber);
 119     my $curClone  = 0;
 120     my $curMarker = 0;
 121     my $curContig = 0;
 122     my %_clones;
 123     my %_markers;
 124     my %_contigs;
 125     my $ctgzeropos = 1;
 126
 127     my $map = Bio::Map::Physical->new('-units' => 'CB',
 128                                      '-type'  => 'physical');
 129
 130     my $filename = $self->file();
 131     my $fh = $self->{'_filehandle'};
 132
 133     if (defined($_readcor)) {
 134         $map->core_exists($_readcor);
 135     }
 136     else {
 137         $map->core_exists(0);
 138     }
 139
 140     if ($map->core_exists()) {
 141         $corfile = substr($filename,0,length($filename)-3)."cor";
 142         if (open(CORE,$corfile)) {
 143             while(read(CORE,$BUFFER,2)) {
 144                 push(@cordata,unpack('n*', $BUFFER));
 145             }
 146         }
 147         else {
 148             $map->core_exists(0);
 149         }
 150     }
 151
 152     ## Read in the header
 153     while (defined($line = <$fh>)) {
 154         chomp($line);
 155
 156         if ($line =~ m{^//\s+fpc\s+project\s+(.+)}) { $map->name($1); }
 157         if ($line =~ m{^//\s+([\d.]+)}) {
 158             my $version = $1;
 159             $version =~ /((\d+)\.(\d+))(.*)/;
 160             $map->version($1);
 161             if ($line =~ /User:\s+(.+)/) { $map->modification_user($1); }
 162         }
 163
 164         if ($line =~ m{^//\s+Framework\s+(\w+)\s+(\w+)\s+([-\w]+)\s+(\w+)\s+(\w+)\s+(.+)$})
 165         {
 166             $map->group_type($3) if ($2 eq "Label");
 167             $map->group_abbr($5) if ($4 eq "Abbrev");
 168         }
 169
 170         last unless ($line =~ m{^//});
 171     }
 172
 173     if (!defined($map->group_type()) || !defined($map->group_abbr()) ) {
 174         $map->group_type("Chromosome");
 175         $map->group_abbr("Chr");
 176     }
 177
 178     $_contigs{0}{'range'}{'end'}   = 0;
 179     $_contigs{0}{'range'}{'start'} = 0;
 180
 181     ## Read in the clone data
 182     while (defined($line = <$fh>)) {
 183         $marker = 0;
 184         $contig = 0;
 185         $seqclone = 0;
 186         $contigNumber = 0;
 187
 188         my ($type,$name);
 189         my (@amatch,@pmatch,@ematch);
 190
 191         my $bandsread = 0;
 192
 193         last if ($line =~ /^Markerdata/);
 194
 195
 196         $line =~ /^(\w+)\s+:\s+"(.+)"/;
 197
 198         ## these will be set if we did find the clone line
 199         ($type, $name) = ($1, $2);
 200
 201         if ($name =~ /sd1/) {
 202             $seqclone = 1;
 203         }
 204
 205         $_clones{$name}{'type'} = $type;
 206         $_clones{$name}{'contig'} = 0;
 207         $_contigs{'0'}{'clones'}{$name}  = 0;
 208
 209         my $temp;
 210
 211         ## Loop through the following lines, getting attributes for clone
 212         while (defined($line = <$fh>) && $line !~ /^\s*\n$/)  {
 213
 214             if ($line =~ /^Map "ctg(\d+)" Ends (Left|Right) ([-\d]+)/)  {
 215                 $_clones{$name}{'contig'} = $1;
 216                 $_contigs{$1}{'clones'}{$name} = 0;
 217
 218                 delete($_contigs{'0'}{'clones'}{$name});
 219
 220                 $temp = $3;
 221                 $contigNumber = $1;
 222                 $line = <$fh>;
 223                 $line =~ /^Map "ctg(\d+)" Ends (Left|Right) ([\d]+)/;
 224                 $_clones{$name}{'range'}{'start'} = $temp;
 225
 226                 $_contigs{$contigNumber}{'range'}{'start'} = $temp
 227                     if (!exists($_contigs{$contigNumber}{'range'}{'start'})
 228                         || $_contigs{$contigNumber}{'range'}{'start'}
 229                         >  $temp );
 230
 231                 $_clones{$name}{'range'}{'end'} = $3;
 232
 233                 $_contigs{$contigNumber}{'range'}{'end'} = $3
 234                     if (!exists($_contigs{$contigNumber}{'range'}{'end'})
 235                         || $_contigs{$contigNumber}{'range'}{'end'} < $3 );
 236
 237             }
 238             elsif ($line =~ /^([a-zA-Z]+)_match_to_\w+\s+"(.+)"/) {
 239                 my $matchtype = "match" . lc(substr($1, 0, 1));
 240                 $_clones{$name}{$matchtype}{$2} = 0;
 241             }
 242             elsif ($line =~ /^Positive_(\w+)\s+"(.+)"/) {
 243                 $_clones{$name}{'markers'}{$2} = 0;
 244                 $_markers{$2}{'clones'}{$name} = 0;
 245                 $_markers{$2}{'type'} = $1;
 246                 $_markers{$2}{'contigs'}{$contigNumber} = 0;
 247                 $_contigs{$contigNumber}{'markers'}{$2} = 0;
 248             }
 249             elsif ($line =~ /^Bands\s+(\d+)\s+(\d+)/ && !$bandsread) {
 250                 my $i = 0;
 251                 my @numbands;
 252                 $bandsread = 1;
 253
 254                 if ($map->core_exists()) {
 255                     while($i<$2){
 256                         push(@numbands,$cordata[($1-1)+$i]);
 257                         $i++;
 258                     }
 259                     $_clones{$name}{'bands'} = \@numbands;
 260                 }
 261                 else {
 262                     push(@numbands,$1,$2);
 263                     $_clones{$name}{'bands'} = \@numbands;
 264                 }
 265                 if (exists($_contigs{0}{'clones'}{$name})) {
 266                     $_clones{$name}{'range'}{'start'} = $ctgzeropos;
 267                     $_clones{$name}{'range'}{'end'} = $ctgzeropos + $2;
 268                     $_contigs{0}{'range'}{'end'} = $ctgzeropos + $2;
 269                     $ctgzeropos += $2;
 270                 }
 271             }
 272             elsif ($line =~ /^Gel_number\s+(.+)/) {
 273                 $_clones{$name}{'gel'} = $1;
 274             }
 275             elsif ($line =~ /^Remark\s+"(.+)"/)  {
 276                 $_clones{$name}{'remark'} .= $1;
 277                 $_clones{$name}{'remark'} .= "\n";
 278                 if($seqclone == 1 ) {
 279                     if( $1 =~ /\,\s+Chr(\d+)\s+/){
 280                         $_clones{$name}{'group'} = $1;
 281                     }
 282                 }
 283             }
 284             elsif ($line =~ /^Fp_number\s+"(.+)"/) {
 285                 $_clones{$name}{'fp_number'} = $1;
 286             }
 287             elsif ($line =~ /^Shotgun\s+(\w+)\s+(\w+)/) {
 288                 $_clones{$name}{'sequence_type'} = $1;
 289                 $_clones{$name}{'sequence_status'} = $2;
 290             }
 291             elsif ($line =~ /^Fpc_remark\s+"(.+)"/) {
 292                 $_clones{$name}{'fpc_remark'} .= $1;
 293                 $_clones{$name}{'fpc_remark'} .= "\n";
 294             }
 295         }
 296
 297         $curClone++;
 298         print "Adding clone $curClone...\n\r"
 299             if ($self->verbose()  && $curClone % 1000 == 0);
 300     }
 301
 302     $map->_setCloneRef(\%_clones);
 303     $line = <$fh>;
 304
 305     while (defined($line = <$fh>) && $line !~ /Contigdata/) {
 306         my ($type,$name);
 307
 308         last if ($line !~ /^Marker_(\w+)\s+:\s+"(.+)"/);
 309
 310         ($type, $name) = ($1, $2);
 311
 312         $_markers{$name}{'type'}   = $type;
 313         $_markers{$name}{'group'}  = 0;
 314         $_markers{$name}{'global'} = 0;
 315         $_markers{$name}{'anchor'} = 0;
 316
 317         while (defined($line = <$fh>) && $line !~ /^\s*\n$/)  {
 318             if ($line =~ /^Global_position\s+([\d.]+)\s*(Frame)?/)  {
 319                 my $position = $1 - floor($1/1000)*1000;
 320                 $position = sprintf("%.2f",$position);
 321
 322                 $_markers{$name}{'global'} = $position;
 323                 $_markers{$name}{'group'}  = floor($1/1000);
 324                 $_markers{$name}{'anchor'} = 1;
 325
 326                 if(defined($2)) {
 327                     $_markers{$name}{'framework'} = 1;
 328                 }
 329                 else {
 330                     $_markers{$name}{'framework'} = 0;
 331                 }
 332             }
 333             elsif ($line =~ /^Anchor_bin\s+"([\w\d.]+)"/) {
 334                 my $grpmatch = $1;
 335                 my $grptype  = $map->group_type();
 336
 337                 $grpmatch =~ /(\d+|\w)(.*)/;
 338
 339                 my ($group,$subgroup);
 340                 $group    = $1;
 341                 $subgroup = $2;
 342
 343                 $subgroup = substr($subgroup,1) if ($subgroup =~ /^\./);
 344
 345                 $_markers{$name}{'group'}      = $group;
 346                 $_markers{$name}{'subgroup'}   = $subgroup;
 347             }
 348             elsif ($line =~ /^Anchor_pos\s+([\d.]+)\s+(F|P)?/){
 349                 $_markers{$name}{'global'}  = $1;
 350                 $_markers{$name}{'anchor'}  = 1;
 351
 352                 if ($2 eq 'F') {
 353                     $_markers{$name}{'framework'} = 1;
 354                 }
 355                 else {
 356                     $_markers{$name}{'framework'} = 0;
 357                 }
 358             }
 359             elsif ($line =~ /^anchor$/) {
 360                 $_markers{$name}{'anchor'} = 1;
 361             }
 362             elsif ($line =~ /^Remark\s+"(.+)"/)  {
 363                 $_markers{$name}{'remark'} .= $1;
 364                 $_markers{$name}{'remark'} .= "\n";
 365             }
 366         }
 367         $curMarker++;
 368         print "Adding Marker $curMarker...\n"
 369             if ($self->verbose() && $curMarker % 1000 == 0);
 370     }
 371
 372     $map->_setMarkerRef(\%_markers);
 373
 374     my $ctgname;
 375     my $grpabbr = $map->group_abbr();
 376     my $chr_remark;
 377
 378     $_contigs{0}{'group'} = 0;
 379
 380     while (defined($line = <$fh>)) {
 381
 382         if ($line =~ /^Ctg(\d+)/) {
 383             $ctgname = $1;
 384             $_contigs{$ctgname}{'group'}      = 0;
 385             $_contigs{$ctgname}{'anchor'}     = 0;
 386             $_contigs{$ctgname}{'position'}   = 0;
 387
 388             if ($line =~ /#\w*(.*)\w*$/) {
 389                 $_contigs{$ctgname}{'remark'} = $1;
 390                 if ($line =~ /#\s+Chr(\d+)\s+/) {
 391                     $_contigs{$ctgname}{'group'}  = $1;
 392                     $_contigs{$ctgname}{'anchor'} = 1;
 393                 }
 394             }
 395         }
 396         elsif ($line =~ /^Chr_remark\s+"(-|\+|Chr(\d+))\s+(.+)"$/) {
 397
 398             $_contigs{$ctgname}{'anchor'}     = 1;
 399             $_contigs{$ctgname}{'chr_remark'} = $3 if(defined($3));
 400
 401             if (defined($2)) {
 402                 $_contigs{$ctgname}{'group'}  = $2;
 403             }
 404             else {
 405                 $_contigs{$ctgname}{'group'}  = "?";
 406             }
 407         }
 408         elsif ($line =~ /^User_remark\s+"(.+)"/) {
 409             $_contigs{$ctgname}{'usr_remark'} = $1;
 410         }
 411         elsif ($line =~ /^Trace_remark\s+"(.+)"/) {
 412             $_contigs{$ctgname}{'trace_remark'} = $1;
 413         }
 414         elsif ($grpabbr && $line =~ /^Chr_remark\s+"(\W|$grpabbr((\d+)|(\w+)|([.\w\d]+)))\s*(\{(.*)\}|\[(.*)\])?"\s+(Pos\s+((\d.)+|NaN))(NOEDIT)?/)
 415         {
 416             my $grpmatch = $2;
 417             my $pos = $10;
 418             if ($pos eq "NaN") {
 419                 $pos = 0;
 420                 print "Warning: Nan encountered for Contig position \n";
 421             }
 422             $_contigs{$ctgname}{'chr_remark'}   = $6;
 423             $_contigs{$ctgname}{'position'} = $pos;
 424             $_contigs{$ctgname}{'subgroup'} = 0;
 425
 426             if (defined($grpmatch)) {
 427                 $_contigs{$ctgname}{'anchor'} = 1;
 428
 429                 if ($grpmatch =~ /((\d+)((\D\d.\d+)|(.\d+)))|((\w+)(\.\d+))/) {
 430
 431                     my ($group,$subgroup);
 432                     $group    = $2 if($grpabbr eq "Chr");
 433                     $subgroup = $3 if($grpabbr eq "Chr");
 434
 435                     $group    = $7 if($grpabbr eq "Lg");
 436                     $subgroup = $8 if($grpabbr eq "Lg");
 437
 438                     $subgroup = substr($subgroup,1) if ($subgroup =~ /^\./);
 439                     $_contigs{$ctgname}{'group'}     = $group;
 440                     $_contigs{$ctgname}{'subgroup'}  = $subgroup;
 441
 442                 }
 443                 else {
 444                     $_contigs{$ctgname}{'group'} = $grpmatch;
 445                 }
 446             }
 447             else {
 448                 $_contigs{$ctgname}{'anchor'} = 1;
 449                 $_contigs{$ctgname}{'group'}  = "?";
 450             }
 451         }
 452         $curContig++;
 453         print "Adding Contig $curContig...\n"
 454             if ($self->verbose() && $curContig % 100 == 0);
 455     }
 456
 457     $map->_setContigRef(\%_contigs);
 458     $map->_calc_markerposition();
 459     $map->_calc_contigposition() if ($map->version() < 7.0);
 460     $map->_calc_contiggroup() if ($map->version() == 4.6);
 461
 462     return $map;
 463 }
 464
 465
 466 =head2 write_map
 467
 468  Title   : write_map
 469  Usage   : $mapio->write_map($map);
 470  Function: Write a map out
 471  Returns : none
 472  Args    : Bio::Map::MapI
 473
 474 =cut
 475
 476 sub write_map{
 477     my ($self,@args) = @_;
 478     $self->throw_not_implemented();
 479 }
 480
 481 1;
 482
 483 =head1 FEEDBACK
 484
 485 =head2 Mailing Lists
 486
 487 User feedback is an integral part of the evolution of this and other
 488 Bioperl modules. Send your comments and suggestions preferably to
 489 the Bioperl mailing list.  Your participation is much appreciated.
 490
 491   bioperl-l@bioperl.org                  - General discussion
 492   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 493
 494 =head2 Reporting Bugs
 495
 496 Report bugs to the Bioperl bug tracking system to help us keep track
 497 of the bugs and their resolution. Bug reports can be submitted via the
 498 web:
 499
 500   http://bugzilla.open-bio.org/
 501
 502 =head1 AUTHOR - Gaurav Gupta
 503
 504 Email gaurav@genome.arizona.edu
 505
 506 =head1 PROJECT LEADERS
 507
 508 Jamie Hatfield            jamie@genome.arizona.edu
 509
 510 Dr. Cari Soderlund        cari@genome.arizona.edu
 511
 512 =head1 PROJECT DESCRIPTION
 513
 514 The project was done in Arizona Genomics Computational Laboratory
 515 (AGCoL) at University of Arizona.
 516
 517 This work was funded by USDA-IFAFS grant #11180 titled "Web Resources
 518 for the Computation and Display of Physical Mapping Data".
 519
 520 For more information on this project, please refer:
 521   http://www.genome.arizona.edu
 522
 523 =head1 APPENDIX
 524
 525 The rest of the documentation details each of the object methods.
 526 Internal methods are usually preceded with a _
 527
 528 =cut