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            'merge=s',
  62            'debug|d+','help|h|?','man|m');
  63
  64 pod2usage() if $options{help};
  65 pod2usage({verbose=>2}) if $options{man};
  66
  67 $DEBUG = $options{debug};
  68
  69 my @USAGE_ERRORS;
  70 if (@ARGV!=1) {
  71     push @USAGE_ERRORS,"You must provide exactly one VCF file to read";
  72 }
  73
  74 pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;
  75
  76
  77 sub open_compressed_file {
  78     my ($file) = @_;
  79     my $fh;
  80     my $mode = '<:encoding(UTF-8)';
  81     my @opts;
  82     if ($file =~ /\.gz$/) {
  83         $mode = '-|:encoding(UTF-8)';
  84         push @opts,'gzip','-dc';
  85     }
  86     if ($file =~ /\.xz$/) {
  87         $mode = '-|:encoding(UTF-8)';
  88         push @opts,'xz','-dc';
  89     }
  90     if ($file =~ /\.bz2$/) {
  91         $mode = '-|:encoding(UTF-8)';
  92         push @opts,'bzip2','-dc';
  93     }
  94     open($fh,$mode,@opts,$file);
  95     return $fh;
  96 }
  97
  98 my $vcf = open_compressed_file($ARGV[0]) or
  99     die "Unable to open file $ARGV[0]";
 100
 101 my %rsids;
 102 while (<STDIN>) {
 103     chomp;
 104     $rsids{$_} = 1;
 105 }
 106
 107 my %merge_rsids;
 108 if (defined $options{merge}) {
 109     my $merge = open_compressed_file($options{merge})
 110         or die "Unable to open file $options{merge}: $!";
 111     while (<$merge>) {
 112         chomp;
 113         my ($old,$new,undef) = split /\t/;
 114         next unless exists $rsids{'rs'.$old};
 115         $merge_rsids{'rs'.$old} = 'rs'.$new;
 116         $rsids{'rs'.$new} = 1;
 117     }
 118     close ($merge);
 119 }
 120
 121 while (<$vcf>) {
 122     if (/^#/o) {
 123         print $_;
 124         next;
 125     }
 126     $_ =~ /^\S+\s+\S+\s+(\S+)/o;
 127     next unless $1;
 128     next unless exists $rsids{$1} and $rsids{$1};
 129     $rsids{$1}++;
 130     print $_;
 131 }
 132
 133 my @unused_rsids;
 134 for my $rsid (keys %rsids) {
 135     if ($rsids{$rsid} == 1) {
 136         push @unused_rsids,$rsid;
 137     }
 138 }
 139 if (@unused_rsids) {
 140     print STDERR "The following rsids were not found\n";
 141     print STDERR map {$_."\n"} @unused_rsids;
 142 }
 143
 144
 145 __END__