scripts/popgen/bp_composite_LD.pl

   1 #!/usr/bin/perl
   2 # -*-Perl-*-
   3
   4 =head1 NAME
   5
   6 bp_composite_LD -i filename.prettybase.txt --sortbyld E<gt> outfile
   7
   8 =head1 SYNOPSIS
   9
  10   bp_composite_LD -i filename.prettybase [-o out.LD] [-f prettybase/csv] [--sortbyld] [--noconvertindels]
  11
  12 =head2 DESCRIPTION
  13
  14 This a script which allows an easy way to calculate composite LD.
  15
  16 =head2 OPTIONS
  17
  18 -i or --in     filename
  19
  20 -f or --format genotype format (prettybase or CSV)
  21
  22 --sortbyld     To see data sorted by LD instead of just all the
  23                site1/site2 pair LD values.
  24
  25 -o or --out    output filename, otherwise will print to STDOUT
  26
  27 --noconvert    (applicable for prettybase format file only)
  28                if specified will NOT attempt to convert indel
  29                states to 'I' and delete states ('-') to 'D'.
  30
  31 -h or --help   see this documentation
  32
  33 =head2 AUTHOR Jason Stajich, Matthew Hahn
  34
  35 For more information contact:
  36
  37 Matthew Hahn, E<lt>matthew.hahn-at-duke.eduE<gt>
  38 Jason Stajich E<lt>jason-at-bioperl-dot-orgE<gt>
  39
  40
  41 =cut
  42
  43 use strict;
  44 use warnings;
  45
  46 use Bio::PopGen::IO;
  47 use Bio::PopGen::Statistics;
  48 use Bio::PopGen::Population;
  49
  50 use Getopt::Long;
  51
  52 my ($file,$outfile,$sortbyld,$format,$noconvert,$verbose);
  53 $format = 'prettybase'; # default format is prettybase
  54 GetOptions(
  55            'i|in:s'       => \$file, # pass the filename as
  56            'o|out:s'      => \$outfile,
  57            'f|format:s'   => \$format,
  58            'sortbyld'     => \$sortbyld,
  59            'noconvert'    => \$noconvert,
  60            'v|verbose'    => \$verbose,
  61            'h|help'       => sub { system('perldoc', $0);
  62                                    exit; },
  63            );
  64
  65 if( ! $file ) {
  66     $file = shift @ARGV;  # if no -i specified
  67 }
  68
  69 my $io = Bio::PopGen::IO->new(-format => $format,
  70                               -verbose=> $verbose,
  71                               -CONVERT_INDEL_STATES => ! $noconvert,
  72                               -file   => $file);
  73
  74 my $stats = Bio::PopGen::Statistics->new(-verbose => $verbose);
  75 my $pop = $io->next_population;
  76
  77 my %LD_matrix = $stats->composite_LD($pop);
  78
  79 # sites can be ordered by sorting their names
  80
  81 my @sites;
  82 my $out;
  83 if( $outfile ) {
  84     open $out, '>', $outfile or die "Could not write file '$outfile': $!\n";
  85 } else {
  86     $out = \*STDOUT;
  87 }
  88 foreach my $site1 ( sort keys %LD_matrix ) {
  89     foreach my $site2 ( sort keys %{ $LD_matrix{$site1} } ) {
  90         my $LD = $LD_matrix{$site1}->{$site2}; # LD for site1,site2
  91         if( $sortbyld ) {
  92             push @sites, [ $site1,$site2,@$LD];
  93         } else {
  94             printf $out "%s,%s - LD=%.4f chisq=%.4f\n",$site1,$site2,@$LD;
  95         }
  96     }
  97 }
  98
  99 if( $sortbyld ) {
 100     foreach my $s ( sort { $b->[3] <=> $a->[3] } @sites ) {
 101         my ($site1,$site2,$ld,$chisq) = @$s;
 102         print $out "$site1,$site2 - LD=$ld, chisq=$chisq\n";
 103     }
 104 }