[Bioperl-l] PopGen modules

[Jason Stajich]

Whoops.
Forgot to call the proper haploid_pop call on the population first.  
Apply this patch and it should work.

Index: Bio/PopGen/Statistics.pm
===================================================================
RCS file: /home/repository/bioperl/bioperl-live/Bio/PopGen/ 
Statistics.pm,v
retrieving revision 1.28
diff -a -u -r1.28 Statistics.pm
--- Bio/PopGen/Statistics.pm    23 Feb 2005 04:51:59 -0000      1.28
+++ Bio/PopGen/Statistics.pm    10 Nov 2005 14:51:08 -0000
@@ -750,7 +750,7 @@
         # This will handle the case when we pass in a PopulationI  
object
         my $pop = $n;
         $totalsites = $seg_sites; # shift the arguments over by one
-       $n = $pop->get_number_individuals;
+       $n = $pop->haploid_population->get_number_individuals;
         $seg_sites = $self->segregating_sites_count($pop);
      }
      my $a1 = 0;

On Nov 10, 2005, at 1:15 AM, Bingshan Li wrote:

> Hi folks,
>
> I recently started to play with PopGen modules but am confused by  
> the difference between "number of individuals" and "sample size".  
> My understanding is that sample size is the number of haploids or  
> chromosomes, and number of individuals is the number of diploids.  
> For example, 100 humans are genotyped, then sample size should be  
> 200 and number of individuals is 100. Am I right? I could be  
> completely wrong but assume I am right for now.
>
> I constructed a population object (named $pop) using prettybase  
> format. Then
>
> $stats = new Bio::PopGen::Statistics();
> $number_individuals = $pop->get_number_individuals();
> $seg_sites = $stats->segregating_sites_count($pop);
> $theta1 = $stats->theta($pop);
> $theta2 = $stats->theta($number_individuals, $seg_sites);
> $theta3 = $stats->theta($number_individuals*2, $seg_sites);
>
> In the above code, $theta1 == $theta2 != $theta3, and I think  
> $theta3 should be the correct answer.
>
> I used "ms" program of Hudson to simulate 200 chromosomes and I  
> used 200 as sample size which gives correct answers (double  
> confirmed with other programs).
>
> Please let me know if I am too naive about this.
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at portal.open-bio.org
> http://portal.open-bio.org/mailman/listinfo/bioperl-l

--
Jason Stajich
Duke University
http://www.duke.edu/~jes12