view release on metacpan or  search on metacpan

Genetics/API/Analysis.pm  view on Meta::CPAN

                             and/or FrequencySource objects.
                             The source(s) for allele frequencies.
                             Required, for obvious reasons.
              ALLELETYPE => The type of alleles whose frequencies are to be graphed.
                            Optional, the default value is "Code".
  Returns   : N/A
  Scope     : Public
  Comments  : Calls xv to display the graphic.

=cut

sub graphAlleleFreqs {
  my($self, %param) = @_ ;
  my($marker, $sourceListPtr, $alleleType, $markerName, @alleleNames, @data, 
     @sourceNames, $source, $alleleName, $freqsPtr, @freqs, $graph, $fileName) ;
  
  # Get the input parameters and check a few things
  defined($marker = $param{MARKER}) or 
                                 croak "ERROR [graphAlleleFreqs]: No Marker!" ;
  $markerName = $marker->field("name") ;
  defined($sourceListPtr = $param{FREQSOURCES}) or 
    croak "ERROR [graphAlleleFreqs]: No source(s) for allele frequencies!" ;
  defined($alleleType = $param{ALLELETYPE}) or $alleleType = "Code" ;

  # Get the frequencies
  @alleleNames = $self->getAllelesByType($marker, $alleleType) ;
  push @data, [ @alleleNames ] ;
  foreach $source (@$sourceListPtr) {
    if ( ref($source) eq "Genetics::FrequencySource" ) {
      $freqsPtr = $source->getAlleleFreqsByMarkerName($markerName, $alleleType) ;
    } elsif ( ref($source) eq "Genetics::Cluster" ) {
      $freqsPtr = $self->getAlleleFreqs($marker, $alleleType, $source) ;
    } else {
      carp "WARNING [graphAlleleFreqs]: $source is an invalid source for allele frequencies!" ;
      next ;
    }
    @freqs = () ;
    push @sourceNames, $source->field("name") ;
    foreach $alleleName (@alleleNames) {
      push @freqs, $$freqsPtr{$alleleName} ;
    }
    push @data, [ @freqs ] ;
  }
  defined($data[1]) or croak "ERROR [graphAlleleFreqs]: No valid allele frequency data!" ;

  # Create and display the graph
  $graph = GD::Graph::bars->new() ;
  $graph->set( 
	      x_label         => 'Allele',
	      y_label         => 'Frequency',
	      title           => "$markerName Allele Frequencies",
	      long_ticks      => 1,
	      y_max_value     => 1.0,
	      y_tick_number   => 20,
	      y_label_skip    => 2,
	      bar_spacing     => 3,
#	      shadow_depth    => 2,
#	      shadowclr       => 'black',
	      transparent     => 0,
	     ) or carp "WARNING [graphAlleleFreqs]: " . $graph->error ;
  $graph->set_legend( @sourceNames );
  $graph->plot( \@data ) or croak "ERROR [graphAlleleFreqs]: " . $graph->error ;
  $fileName = "/tmp/$$.png" ;
  &_saveGraphToFile($graph, $fileName) ;
  system("xv $fileName &") == 0 or 
         croak "ERROR [graphAlleleFreqs]: system xv $fileName & failed: $?" ;

  return(1) ;
}

=head2 calculateHet

 Function  : Calculate the heterozygosity for a Marker or SNP.
 Arguments : A Marker object, a string containing an allele type, and one of 
              the following defining the Subject group:
                - a Subject Cluster object 
                - an array reference to a list of Subject objects
                - a Kindred Cluster object 
                - an array reference to a list of Kindred objects
 Returns   : A scalar float
 Scope     : Public
 Comments  : Arguments are passed directly to API::DB::Query::getAlleleFreqs()

=cut

sub calculateHet {
  my($self, $marker, $alleleType, $subjGroup) = @_ ;
  my($freqsPtr, $allele, $sumFreqsSquared, $het) ;

  # Parameter checking, to the extent that it happens at all, is taken 
  # care of by API::DB::Query::getAlleleFreqs()
  # Get the allele frequencies
  $freqsPtr = $self->getAlleleFreqs($marker, $alleleType, $subjGroup) ;
  # Calculate heterozygosity
  foreach $allele (sort keys %$freqsPtr) {
    next if ($alleleType eq "Nucleotide" and $allele eq "N") ;
    $sumFreqsSquared += ($$freqsPtr{$allele})**2 ;
  }
  $het = 1 - $sumFreqsSquared ;

  return( _formatFloat($het) ) ;
}

=head2 calculateSnpHW

 Function  : 
 Arguments : 
 Returns   : 
 Example   : calculateHW()
 Scope     : 
 Comments  : 

=cut

sub calculateSnpHW {
  my($self, $po, $sc) = @_ ;
  my($gtCountsPtr, $alleleCountsPtr, @alleles, $i, @allelePairs, 
     $n11, $n12, $n22, $n, $obs, $exp, $chiSquare) ;

  if ( ref($po) !~ /^Genetics::(Marker|SNP)$/ ) {
    croak "ERROR [calculateHW]: Invalid input Marker/SNP $po." ;