|
| 1 | +#!/usr/bin/perl |
| 2 | + |
| 3 | +##########LICENCE########## |
| 4 | +# PCAP - NGS reference implementations and helper code for the ICGC/TCGA Pan-Cancer Analysis Project |
| 5 | +# Copyright (C) 2014-2016 ICGC PanCancer Project |
| 6 | +# |
| 7 | +# This program is free software; you can redistribute it and/or |
| 8 | +# modify it under the terms of the GNU General Public License |
| 9 | +# as published by the Free Software Foundation; either version 2 |
| 10 | +# of the License, or (at your option) any later version. |
| 11 | +# |
| 12 | +# This program is distributed in the hope that it will be useful, |
| 13 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 14 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 15 | +# GNU General Public License for more details. |
| 16 | +# |
| 17 | +# You should have received a copy of the GNU General Public License |
| 18 | +# along with this program; if not see: |
| 19 | +# http://www.gnu.org/licenses/gpl-2.0.html |
| 20 | +##########LICENCE########## |
| 21 | + |
| 22 | +use Cwd qw(abs_path); |
| 23 | +use strict; |
| 24 | +use English qw( -no_match_vars ); |
| 25 | +use warnings FATAL => 'all'; |
| 26 | + |
| 27 | +use File::Basename; |
| 28 | +use Carp; |
| 29 | +use Getopt::Long; |
| 30 | +use Pod::Usage; |
| 31 | + |
| 32 | +use PCAP; |
| 33 | + |
| 34 | +use Bio::DB::BigWig 'binMean','binStdev'; |
| 35 | + |
| 36 | +my %chr_stats; |
| 37 | +my @chr_order; |
| 38 | + |
| 39 | +{ |
| 40 | + my $options = option_builder(); |
| 41 | + |
| 42 | + my $wig = Bio::DB::BigWig->new(-bigwig=>$options->{'b'}); |
| 43 | + my @chroms = $wig->features(-type=>'summary'); |
| 44 | + |
| 45 | + for my $c (@chroms) { |
| 46 | + my $seqid = $c->seq_id; |
| 47 | + next if(defined $options->{'r'} && $options->{'r'} ne $seqid && 'chr'.$options->{'r'} ne $seqid); |
| 48 | + my $start = $c->start; |
| 49 | + |
| 50 | + my $stats = $c->statistical_summary(1); |
| 51 | + my $bin_width = $c->length/@$stats; |
| 52 | + |
| 53 | + my $s = shift @{$stats}; |
| 54 | + |
| 55 | + my $mean = binMean($s); |
| 56 | + my $stdev = binStdev($s); |
| 57 | + my $end = $start + $bin_width-1; |
| 58 | + |
| 59 | + push @chr_order, $seqid; |
| 60 | + $chr_stats{$seqid} = {'mean' => binMean($s), |
| 61 | + 'stdev' => binStdev($s)}; |
| 62 | + warn sprintf "%s: mean %.2f, stdev %.2f\n", $seqid, $chr_stats{$seqid}{'mean'}, $chr_stats{$seqid}{'stdev'}; |
| 63 | + } |
| 64 | + |
| 65 | + open my $OFH, '>', $options->{'o'} or die "Failed to create $options->{o}: $!\n"; |
| 66 | + for my $chr(@chr_order) { |
| 67 | + my $max_val = $chr_stats{$chr}{'mean'} + ($chr_stats{$chr}{'stdev'} * $options->{'s'}); |
| 68 | + warn sprintf "%s: Max depth permitted = %d\n", $chr, $max_val; |
| 69 | + my $iterator = $wig->get_seq_stream(-seq_id=> $chr); |
| 70 | + while (my $p = $iterator->next_seq) { |
| 71 | + next if($p->score <= $max_val); |
| 72 | + printf $OFH "%s\t%d\t%d\t%d\n", $chr, $p->start-1, $p->end, $p->score; |
| 73 | + } |
| 74 | + } |
| 75 | + close $OFH; |
| 76 | +} |
| 77 | + |
| 78 | +sub option_builder { |
| 79 | + my ($factory) = @_; |
| 80 | + |
| 81 | + my %opts; |
| 82 | + |
| 83 | + &GetOptions ( |
| 84 | + 'h|help' => \$opts{'h'}, |
| 85 | + 'b|bigwig=s' => \$opts{'b'}, |
| 86 | + 'o|output=s' => \$opts{'o'}, |
| 87 | + 'r|ref=s' => \$opts{'r'}, |
| 88 | + 'd|decode=s@' => \$opts{'d'}, |
| 89 | + 's|sd=n' => \$opts{'s'}, |
| 90 | + 'v|version' => \$opts{'v'}, |
| 91 | + ); |
| 92 | + |
| 93 | + if(defined $opts{'v'}) { |
| 94 | + print PCAP->VERSION,"\n"; |
| 95 | + exit 0; |
| 96 | + } |
| 97 | + |
| 98 | + pod2usage(0) if($opts{'h'}); |
| 99 | + |
| 100 | + pod2usage(1) if(!$opts{'b'} || !$opts{'o'}); |
| 101 | + |
| 102 | + croak $opts{'b'}.' was not found or is empty' if(!-e $opts{'b'} || !-s $opts{'b'}); |
| 103 | + |
| 104 | + if($opts{'d'}) { |
| 105 | + if(!$opts{'r'}) { |
| 106 | + croak '-d should not be defined without -r'; |
| 107 | + } |
| 108 | + my %decode; |
| 109 | + foreach my $d_str(@{$opts{'d'}}) { |
| 110 | + if($d_str =~ m/^(\d+)\:(.*)$/) { |
| 111 | + my $num = $1; |
| 112 | + my $chr = $2; |
| 113 | + $decode{$num} = $chr; |
| 114 | + } |
| 115 | + else { |
| 116 | + croak "Decode string of $d_str is invalid see --help"; |
| 117 | + } |
| 118 | + } |
| 119 | + if(defined $decode{$opts{'r'}}) { |
| 120 | + $opts{'r'} = $decode{$opts{'r'}}; |
| 121 | + } |
| 122 | + } |
| 123 | + |
| 124 | + my $fn = fileparse($opts{'b'}); |
| 125 | + $fn =~ s/\.bw$//; |
| 126 | + $opts{'o'} .= '/' if($opts{'o'} !~ m/\/$/); |
| 127 | + $opts{'o'} .= $fn; |
| 128 | + if($opts{'r'}) { |
| 129 | + $opts{'o'} .= '.'.$opts{'r'}; |
| 130 | + } |
| 131 | + $opts{'o'} .= '.bed'; |
| 132 | + |
| 133 | + if(!$opts{'s'}) { |
| 134 | + $opts{'s'} = 12; |
| 135 | + } |
| 136 | + |
| 137 | + return \%opts; |
| 138 | +} |
| 139 | + |
| 140 | +__END__ |
| 141 | +
|
| 142 | +=head1 NAME |
| 143 | +
|
| 144 | +detectExtremeDepth.pl - Generate profile of BigWig file and identify regions outside the normal range |
| 145 | +
|
| 146 | +=head1 SYNOPSIS |
| 147 | +
|
| 148 | + General Options (list OR project must be defined): |
| 149 | +
|
| 150 | + --bigwig (-b) FILE BigWig file path |
| 151 | + --output (-o) DIR Folder to send output to |
| 152 | + - named as input file with '.tab' extension |
| 153 | + - if '-r' defined '.{val}' will prefix '.bed' |
| 154 | +
|
| 155 | + Optional: |
| 156 | + --ref (-r) STR Restrict to this reference (mainly for testing) |
| 157 | + - without 'chr' prefix, will test with and without the 'chr' for you. |
| 158 | + --decode (-d) STR Decode -r to chromosome names (do not include 'chr') |
| 159 | + e.g. -d 23:X -d 24:Y -d 25:MT |
| 160 | + --sd (-s) INT Number of standard deviations above mean for group to be included [12] |
| 161 | + --help (-h) This message |
| 162 | + --version (-v) Version |
| 163 | +
|
| 164 | + Examples: |
| 165 | + perl ~/detectExtremeDepth.pl -o someplace -b sample.bw |
| 166 | +
|
| 167 | +=cut |
0 commit comments