vcf_rs_grep

   1 #!/usr/bin/perl
   2 # vcf_rs_grep greps RS from a VCF file
   3 # and is released under the terms of the GNU GPL version 3, or any
   4 # later version, at your option. See the file README and COPYING for
   5 # more information.
   6 # Copyright 2017 by Don Armstrong <don@donarmstrong.com>.
   7
   8
   9 use warnings;
  10 use strict;
  11
  12 use Getopt::Long;
  13 use Pod::Usage;
  14
  15 =head1 NAME
  16
  17 vcf_rs_grep - greps RS from a VCF file
  18
  19 =head1 SYNOPSIS
  20
  21 vcf_rs_grep [options] vcf_file.gz < rs_list |gzip -c > vcf_greped.gz
  22
  23  Options:
  24    --debug, -d debugging level (Default 0)
  25    --help, -h display this help
  26    --man, -m display manual
  27
  28 =head1 OPTIONS
  29
  30 =over
  31
  32 =item B<--debug, -d>
  33
  34 Debug verbosity. (Default 0)
  35
  36 =item B<--help, -h>
  37
  38 Display brief usage information.
  39
  40 =item B<--man, -m>
  41
  42 Display this manual.
  43
  44 =back
  45
  46 =head1 EXAMPLES
  47
  48 vcf_rs_grep vcf_file.gz < rs_list |gzip -c > vcf_greped.gz
  49
  50 =cut
  51
  52
  53 use vars qw($DEBUG);
  54
  55 my %options = (debug           => 0,
  56                help            => 0,
  57                man             => 0,
  58               );
  59
  60 GetOptions(\%options,
  61            'debug|d+','help|h|?','man|m');
  62
  63 pod2usage() if $options{help};
  64 pod2usage({verbose=>2}) if $options{man};
  65
  66 $DEBUG = $options{debug};
  67
  68 my @USAGE_ERRORS;
  69 if (@ARGV!=1) {
  70     push @USAGE_ERRORS,"You must provide exactly one VCF file to read";
  71 }
  72
  73 pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;
  74
  75
  76 sub open_compressed_file {
  77     my ($file) = @_;
  78     my $fh;
  79     my $mode = '<:encoding(UTF-8)';
  80     my @opts;
  81     if ($file =~ /\.gz$/) {
  82         $mode = '-|:encoding(UTF-8)';
  83         push @opts,'gzip','-dc';
  84     }
  85     if ($file =~ /\.xz$/) {
  86         $mode = '-|:encoding(UTF-8)';
  87         push @opts,'xz','-dc';
  88     }
  89     if ($file =~ /\.bz2$/) {
  90         $mode = '-|:encoding(UTF-8)';
  91         push @opts,'bzip2','-dc';
  92     }
  93     open($fh,$mode,@opts,$file);
  94     return $fh;
  95 }
  96
  97 my $vcf = open_compressed_file($ARGV[0]) or
  98     die "Unable to open file $ARGV[0]";
  99
 100 my %rsids;
 101 while (<STDIN>) {
 102     chomp;
 103     $rsids{$_} = 1;
 104 }
 105
 106 while (<$vcf>) {
 107     if (/^#/o) {
 108         print $_;
 109         next;
 110     }
 111     $_ =~ /^\S+\s+\S+\s+(\S+)/o;
 112     next unless $1;
 113     next unless exists $rsids{$1} and $rsids{$1};
 114     print $_;
 115 }
 116
 117 __END__