getPileupResultsByPosns.pl

#!/usr/bin/perl -w 

=head1 NAME

perl-template.pl

-head1 SYNOPSIS

=head1 OPTIONS

    -pileupFile <STDIN>     generated by samtools pileup command
    -inFormat|f <string>   {'SNVMix','VCF','SAMtools'}
    -positions|p <string>     filename of file containing single column (chr:posn)
    -out|o <string>

=head1 DESCRIPTION

=head1 CONTACT

Gavin Ha <gha@bccrc.ca>

=cut

use strict;
#use DBI;
use File::Basename;
use Getopt::Long;

sub usage () {
    exec('perldoc', $0);
    exit;
}

my ($snvMixFile, $inFormat, $positions, $out, $help);

GetOptions (
    'positions|p=s' => \$positions,
    'inFormat|f=s' => \$inFormat,
    #'out|o=s' => \$out,
    'help|?' => \$help
    );

if($help) {
    &usage();
}
#print "Parameters:\ninType=$inFormat\npositions=$positions\nout=$out\n";

open POS, $positions || die("Can't open $positions\n");
<POS>; #skip first line
my %pos = ();
while (<POS>) {
    chomp;
    if ($inFormat eq 'SNVMix'){    	
	my ($chr,@rest) = split(/\t/,$_);
	$chr =~ s/chr//g;
	$chr =~ s/23/X/g; $chr =~ s/24/Y/g;
	$pos{$chr} = 1;
    }elsif ($inFormat eq 'VCF'){
	my ($chr,$posn,@rest) = split(/\t/,$_);
	$chr =~ s/chr//g;
	$chr =~ s/23/X/g; $chr =~ s/24/Y/g;
	my $chrStr = $chr . ":" . $posn;
	$pos{$chrStr} = 1;
    }elsif ($inFormat eq 'SAMtools'){
	my ($ps_id,$chr,$posn,@rest) = split(/\t/,$_);
	$chr =~ s/chr//g;
	$chr =~ s/23/X/g; $chr =~ s/24/Y/g;
	my $chrStr = $chr . ":" . $posn;
	$pos{$chrStr} = 1;
    }
    
	
}
close POS;

#open SNV, $snvMixFile || die("Can't open $snvMixFile\n");
#open OUT, ">$out" || die("Can't write to $out\n");
my $linesRead = 0;
while (<STDIN>){
    chomp;
    my ($chr, $posn,$base,$depth,@rest) = split(/\t/,$_);
	$chr =~ s/chr//g;
	my $chrPosn = $chr . ":" . $posn;    
    if (defined $pos{$chrPosn} && $depth>=5) {
	    #print OUT $_ . "\n";
	print $_ . "\n";
    }
    $linesRead +=1;
    #print "Extracting pileup by positions: $linesRead lines read\n" if($linesRead%1000000==0);
}
#close OUT;