Advanced Testcases

Home | Search | FPC | WebAGCoL

Following are some of the more advanced testcases and corresponding outputs showing how to use BioPerl modeles for FPC.

Testcase #1

# Finds all the markers that starts with either "SOG" or whatever prefix
# is given on command line and the number of clones each of them hits.

#!/usr/local/bin/perl -w

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';
my $match = '^'.($ARGV[1] || 'SOG');
my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,-species => "Orysa sativa",
         -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();
my $markercount = 0;

foreach my $marker ($fpcmap->each_markerid()) {
  if($marker =~ /$match/) {
    # Create the marker object of only the required markers
    my $markerobj = $fpcmap->get_markerobj($marker);

    # Get the number of clones in each
    my $clones = scalar ($markerobj->each_cloneid());

    print "Marker: $marker\tClone Hits: $clones\n";
    $markercount++;
  }
}

print "\n\nTotal Markers: $markercount\n\n";

Output #1

Marker: SOG0282 Clone Hits: 7
Marker: SOG0283 Clone Hits: 4
Marker: SOG0284 Clone Hits: 13
Marker: SOG0285 Clone Hits: 4
Marker: SOG0286 Clone Hits: 23
Marker: SOG0287 Clone Hits: 3
Marker: SOG0288 Clone Hits: 1
Marker: SOG0290 Clone Hits: 1
Marker: SOG0292 Clone Hits: 8
Marker: SOG0293 Clone Hits: 50
Marker: SOG0294 Clone Hits: 5
Marker: SOG0295 Clone Hits: 1
Marker: SOG0296 Clone Hits: 1
Marker: SOG0298 Clone Hits: 1
Marker: SOG0301 Clone Hits: 3
Marker: SOG0302 Clone Hits: 124
Marker: SOG0303 Clone Hits: 71
Marker: SOG0305 Clone Hits: 2
.
.
.

Total Markers: 1504

Testcase #2

# Find all the contigs that lie on the Group (chromosome) specified on the command line [default chr1]

#!/usr/local/bin/perl

use Bio::MapIO;
my $file = $ARGV[0] || 'rice.fpc';
my $chrmatch = $ARGV[1] || '1';
my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,-species => "Orysa sativa",
     -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();
my (@ctglist, @sortedctglist);

foreach my $contig ($fpcmap->each_contigid()) {
  next unless ($contig);

  # create the object of the contig
  my $ctgobj = $fpcmap->get_contigobj($contig);

  # check if it lies in the same chromosome
  if ($chrmatch == $ctgobj->group()) {
    push @ctglist, $ctgobj;
  }
}

print "Total Contigs in Group $chrmatch: ", scalar(@ctglist),"\n";
@sortedctglist = sort {$a->name() <=> $b->name()} @ctglist;

if (scalar(@sortedctglist) != 0) {
  print "Contigs in Group $chrmatch:\n";
  foreach $ctg (@sortedctglist) {
    print " +++ctg",$ctg->name(),"\tposition:",$ctg->position($ctg->group()),"\n";
  }
}

Output #2

Total Contigs in Group 1: 28
Contigs in Group 1:
+++ctg1 position:1.588
+++ctg2 position:10.235
+++ctg3 position:25.195
+++ctg4 position:39.362
+++ctg5 position:50.494
+++ctg6 position:53.950
+++ctg7 position:60.525
+++ctg8 position:68.442
+++ctg9 position:72.800
+++ctg10 position:73.250
+++ctg11 position:81.580
+++ctg12 position:102.250
+++ctg13 position:87.463
+++ctg14 position:92.700
+++ctg15 position:98.153
+++ctg16 position:102.745
+++ctg17 position:105.821
.
.
.
.

Testcase #3

# List the contigs based on the number of clones they hit

#!/usr/local/bin/perl -w

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';

my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,-species => "Orysa sativa",
     -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();
my (@reclist,@sortedreclist);

foreach my $contig($fpcmap->each_contigid() ) {
   next unless ($contig);
   my $record;

   # create the object of the contig
   my $ctgobj = $fpcmap->get_contigobj($contig);

   my @clones = $ctgobj->each_cloneid();
   $record->{hits} = scalar(@clones);
   $record->{contig} = $contig;

   push @reclist,$record;
}

@sortedreclist = sort { $b->{hits} <=> $a->{hits} } @reclist;

foreach my $rec (@sortedreclist) {
   print "Ctg",$rec->{contig},"\tNumber of clones: ",$rec->{hits},"\n";
}

Output #3

Ctg183 Number of clones: 1772
Ctg439 Number of clones: 1705
Ctg126 Number of clones: 1518
Ctg81 Number of clones: 1342
Ctg91 Number of clones: 1203
Ctg72 Number of clones: 1097
Ctg99 Number of clones: 1085
Ctg90 Number of clones: 1009
Ctg73 Number of clones: 864
Ctg38 Number of clones: 791
Ctg255 Number of clones: 782
Ctg142 Number of clones: 779
Ctg89 Number of clones: 756
Ctg162 Number of clones: 756
Ctg65 Number of clones: 728
Ctg191 Number of clones: 705
.
.
.
.

Testcase #4

# List the markers based on the number of clones each hit for a given group (chromosome)

#!/usr/local/bin/perl

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';
my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,
    -species => "Orysa sativa",
    -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();
my (@reclist,@sortedreclist);

foreach my $marker ($fpcmap->each_markerid() ) {
  my $record;

  # create the object of the contig
  my $ctgobj = $fpcmap->get_markerobj($marker);

  my @clones = $ctgobj->each_cloneid();
  $record->{hits} = scalar(@clones);
  $record->{marker} = $marker;

  push @reclist,$record;
}

@sortedreclist = sort { $b->{hits} <=> $a->{hits} } @reclist;

foreach my $rec (@sortedreclist) {
  print "Marker: ",$rec->{marker},"\tNumber of clones: ",$rec->{hits},"\n";
}

Output #4

Marker: SOG1574 Number of clones: 769
Marker: SOG1872 Number of clones: 559
Marker: SOG0258 Number of clones: 390
Marker: PRC0214 Number of clones: 389
Marker: PSHR174.1 Number of clones: 318
Marker: SOG1841 Number of clones: 289
Marker: SOG1241 Number of clones: 251
Marker: 10E20fA Number of clones: 250
Marker: SOG1778 Number of clones: 238
Marker: CSU490 Number of clones: 236
Marker: SOG0200 Number of clones: 209
Marker: 35O21f Number of clones: 207
Marker: 10E20fC Number of clones: 200
Marker: SOG0702 Number of clones: 193
Marker: SOG2148 Number of clones: 186
.
.
.
.

Testcase #5

#!/usr/local/bin/perl

# Prints the number of clones and marker for all contigs
# on a given chromosome sorted by the number of clones

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';
my $chrmatch = $ARGV[1] || '1';

my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,
    -species => "Orysa sativa",
    -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();
my (@ctgrec, @sortedctgrec);

foreach my $contig ($fpcmap->each_contigid()) {
  my $record;
  # to skip contig 0   next unless ($contig);

  # create the object of the contig
  my $ctgobj = $fpcmap->get_contigobj($contig);

  # check if it lies in the same chromosome
  if ($ctgobj->group() =~ /\d+/ && $chrmatch == $ctgobj->group()) {
    $record->{clones} = scalar($ctgobj->each_cloneid());
    $record->{markers} = scalar($ctgobj->each_markerid());
    $record->{contig} = $contig;
    push @ctglist, $record;
  }

}

@sortedctgrec = sort {$b->{clones} <=> $a->{clones} } @ctglist;

print "Chromosome #$chrmatch\n";

foreach my $ctgrec (@sortedctgrec) {
  print " +++ctg",$ctgrec->{contig},"\n";
  print " Clones: ",$ctgrec->{clones},"\n";
  print " Markers: ",$ctgrec->{markers},"\n\n";
}

Output #5

Chromosome #1
+++ctg3
  Clones: 623
  Markers: 172

+++ctg2
  Clones: 618
  Markers: 150

+++ctg13
  Clones: 586
  Markers: 81

+++ctg22
  Clones: 574
  Markers: 175

+++ctg21
  Clones: 536
  Markers: 126

+++ctg8
  Clones: 489
  Markers: 84
.
.
.
.

Testcase #6

#!/usr/local/bin/perl -w

# For a given clone if it is a parent or child print all
# markers they share and they dont.

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';
my $clonematch = $ARGV[1] || 'a0060J18';
my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,
    -species => "Orysa sativa",
    -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();

my $cloneobj = $fpcmap->get_cloneobj($clonematch);
my @markerlist = $cloneobj->each_markerid();

my @exactmatch = $cloneobj->each_match('exact');
my @approxmatch = $cloneobj->each_match('approximate');
my @pseudomatch = $cloneobj->each_match('pseudo');

my @childmarkers;

# cheack if its a child
if (scalar(@exactmatch) || scalar(@approxmatch) || scalar(@pseudomatch)) {

  print "$clonematch is a child clone\n";

  foreach $approxclone(@approxmatch) {
    my $approxcloneobj = $fpcmap->get_cloneobj($approxclone);
    @childmarkers = $approxcloneobj->each_markerid();
  }

  foreach $exactclone(@exactmatch) {
    my $exactcloneobj = $fpcmap->get_cloneobj($exactclone);
    @childmarkers = $exactcloneobj->each_markerid();
  }

  foreach $pseudoclone(@pseudomatch) {
    my $pseudocloneobj = $fpcmap->get_cloneobj($pseudoclone);
    @childmarkers = $pseudocloneobj->each_markerid();
  }
}
else {
  my (@childclones,@clones,@child);

  foreach my $clone ($fpcmap->each_cloneid()) {
    my $cobj = $fpcmap->get_cloneobj($clone);

    my @echild = $cobj->each_match('exact');
    @child = grep {$_ eq $clonematch} @echild;
    push(@childclones,$clone) if (scalar(@child));

    my @achild = $cobj->each_match('approximate');
    @child = grep {$_ eq $clonematch} @achild;
    push(@childclones,$clone) if (scalar(@child));

    my @pchild = $cobj->each_match('pseudo');
    @child = grep {$_ eq $clonematch} @pchild;
    push(@childclones,$clone) if (scalar(@child));
  }

  if (scalar(@childclones)) {
    print "$clonematch is a parent clone\n";
    foreach my $clone (@childclones) {
      my $cobj = $fpcmap->get_cloneobj($clone);
      push (@childmarkers,$cobj->each_markerid());
    }
  }
}

# Markers not shared
my %seen;
my @aonly;

@seen{@childmarkers} = ();

foreach $mkr(@markerlist) {
  push(@aonly,$mkr) unless exists ($seen{$mkr});
}

$_ = " ";
print "Markers not shared: @aonly\n";

# Shared markers
@union = @isect = ();
%union = %isect = ();

foreach $mkr (@markerlist, @childmarkers) {$union{$mkr}++ && $isect{$mkr}++}
@isect = keys %isect;

print "Shared markers: @isect\n";

Output #6

a0060J18 is a parent clone
Markers not shared:
Shared markers: SOG1841

Testcase #7

#!/usr/local/bin/perl -w
# For each marker print the contig and its position in each and
# number of clones

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';
my $markermatch = $ARGV[1] || '(.+)';

my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,
    -species => "Orysa sativa",
    -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();

foreach my $marker ($fpcmap->each_markerid()) {

  if($marker =~ /$markermatch/) {
    my $markerobj = $fpcmap->get_markerobj($marker);
    print "Marker $marker\n";

    foreach my $contig ($markerobj->each_contigid()) {
      my $clonecnt;
      print "\tctg$contig\tposition: ",$markerobj->position($contig),"\t#Clones: ";
      foreach my $clone ($markerobj->each_cloneid()){
        my $cloneobj = $fpcmap->get_cloneobj($clone);
        $clonecnt++ if ($cloneobj->contigid() == $contig);
      }
      print "$clonecnt\n";
    }
  }
}

Output #7

Marker E30009S
  ctg253 position: 24 #Clones: 6
Marker RZ284
  ctg331 position: 14 #Clones: 1
Marker SOG0282
  ctg13 position: 513 #Clones: 1
  ctg207 position: 83 #Clones: 2
  ctg135 position: 279 #Clones: 1
  ctg91 position: 305 #Clones: 1
  ctg236 position: 22 #Clones: 1
  ctg68 position: 280 #Clones: 1
.
.
.
.

Testcase #8

#!/usr/local/bin/perl -w

# For a given contig print all marker that hit another contig
# indicate the ones that are framework

use Bio::MapIO;

my $file = $ARGV[0] || 'rice.fpc';
my $ctgmatch = $ARGV[1] || '48';

my $mapio = new Bio::MapIO(-format => "fpc",-file => $file,
    -species => "Orysa sativa",
    -readcor => 0, -verbose => 0);

my $fpcmap = $mapio->next_map();
my $ctgobj = $fpcmap->get_contigobj($ctgmatch);

foreach $marker ($ctgobj->each_markerid()) {
  my $markerobj = $fpcmap->get_markerobj($marker);
  if (scalar($markerobj->each_contigid()) > 1 && $markerobj->framework()) {
    my @ctgs = $markerobj->each_contigid();
    $_ = " ";
    print "$marker\tContigs @ctgs\n";
  }
}

Output #8

R1381 Contigs 48 3
S13984 Contigs 0 48 2
R63 Contigs 0 48
R1737 Contigs 0 48
C1236 Contigs 0 48
E11382SB Contigs 125 255 239 0 118 183 240 48 81 91

Testcase #9

#!/usr/local/bin/perl -w

# Compare the clone matches in two FPC file
# usage: perl compare_fpc [#ofclones default: all clones in fpc1]

use strict;
use Bio::MapIO;

my $file1 = shift;
my $file2 = shift;
my $matches = shift;
my $gellen1 = shift;
my $cutoff1 = shift;
my $gellen2 = shift;
my $cutoff2 = shift;

$cutoff1 = 1e-10 if(!defined($cutoff1));
$cutoff2 = 1e-10 if(!defined($cutoff2));
$gellen1 = 3300.00 if(!defined($gellen1));
$gellen2 = 3300.00 if(!defined($gellen2));

my $clonesmatched;

my $mapio1 = new Bio::MapIO(-format => "fpc", -species => "foo", -readcor => 1,-file => $file1);
my $mapio2 = new Bio::MapIO(-format => "fpc", -species => "foo", -readcor => 1,-file => $file2);

my $fpcmap1 = $mapio1->next_map();
my $fpcmap2 = $mapio2->next_map();

if(!$fpcmap1->core_exists() || !$fpcmap2->core_exists()) {
   print "ERR: make sure COR file is in place\n";
   exit;
}

my @clonesA = $fpcmap1->each_cloneid();
my @clonesB = $fpcmap2->each_cloneid();

my (%clonesA,%clonesB);
my (@aonly,@bonly,@isect,@union);
my (%isect,%union,%seen);
my (@file1only,@file2only,@bothfiles);
my (%match1,%match2);
my ($i,$j,$k);

my $tol = 7;

@clonesA{@clonesA} = ();
@clonesB{@clonesB} = ();

if (!defined($matches)) {
   $clonesmatched = scalar(@clonesA) ;
}
else {
   $clonesmatched = $matches;
}

for ($i=0; $i<$clonesmatched; $i++) {
   if (exists($clonesB{$clonesA[$i]})) {
     for ($j=$i+1; $j < scalar(@clonesA); $j++) {
       my $score = $fpcmap1->coincidence_score($clonesA[$i],$clonesA[$j],$tol,$gellen1);
       if ($score <= $cutoff1) {
         $match1{$clonesA[$j]} = 1;
         $match1{$clonesA[$i]} = 1 if (!exists($clonesA[$i]));
       }
     }
     for ($k=0; $k < scalar(@clonesB); $k++) {
       if ($clonesB[$k] ne $clonesA[$i]) {
         my $score = $fpcmap2->coincidence_score($clonesB[$k],$clonesA[$i],$tol,$gellen2);
         if ($score <= $cutoff2) {
           $match2{$clonesB[$k]} = 1;
           $match2{$clonesA[$i]} = 1 if (!exists($clonesA[$i]));
         }
       }
     }
   }
}

$" = " ";

my @clonematchA = keys ( %match1 );
my @clonematchB = keys ( %match2 );

# clone matches in fpc1 that are not in fpc2
@seen{@clonematchB} = ();
foreach my $item (@clonematchA) {
   push (@aonly,$item) unless exists $seen{$item};
}
print ("Clone matches in $file1 that are not in $file2: \n@aonly");
print "\n\n";

# clone matches in fpc1 that are not in fpc2
@seen{@clonematchA} = ();
foreach my $item (@clonematchB) {
   push (@bonly,$item) unless exists $seen{$item};
}
print ("Clone matches in $file2 that are not in $file1: \n@bonly");
print "\n\n";

# clone matches in both
@union = @isect = ();
%union = %isect = ();

foreach my $match (@clonematchA, @clonematchB) {$union{$match}++ && $isect{$match}++}
@isect = keys %isect;
print "Clone matches in both $file1 and $file2(",scalar(@isect),"): \n @isect\n";

Email Comments To: www@agcol.arizona.edu