← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to scripts/popgen/bp_composite_LD.pl

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
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") || die("$outfile: $!");
 
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
}
 
105
 
 
106
 
 
107