bcftools/vcfutils.pl

   1 #!/usr/bin/perl -w
   2
   3 # Author: lh3
   4
   5 use strict;
   6 use warnings;
   7 use Getopt::Std;
   8
   9 &main;
  10 exit;
  11
  12 sub main {
  13   my $version = '0.1.0';
  14   &usage if (@ARGV < 1);
  15   my $command = shift(@ARGV);
  16   my %func = (subsam=>\&subsam);
  17   die("Unknown command \"$command\".\n") if (!defined($func{$command}));
  18   &{$func{$command}};
  19 }
  20
  21 sub subsam {
  22   die(qq/Usage: vcfutils.pl subsam <in.vcf> [samples]\n/) if (@ARGV == 0);
  23   my ($fh, %h);
  24   my $fn = shift(@ARGV);
  25   my @col;
  26   open($fh, ($fn =~ /\.gz$/)? "gzip -dc $fn |" : $fn) || die;
  27   $h{$_} = 1 for (@ARGV);
  28   while (<$fh>) {
  29         if (/^##/) {
  30           print;
  31         } elsif (/^#/) {
  32           my @t = split;
  33           my @s = @t[0..8]; # all fixed fields + FORMAT
  34           for (9 .. $#t) {
  35                 if ($h{$t[$_]}) {
  36                   push(@s, $t[$_]);
  37                   push(@col, $_);
  38                 }
  39           }
  40           pop(@s) if (@s == 9); # no sample selected; remove the FORMAT field
  41           print join("\t", @s), "\n";
  42         } else {
  43           my @t = split;
  44           if (@col == 0) {
  45                 print join("\t", @t[0..7]), "\n";
  46           } else {
  47                 print join("\t", @t[0..8], map {$t[$_]} @col), "\n";
  48           }
  49         }
  50   }
  51   close($fh);
  52 }
  53
  54 sub usage {
  55   die(qq/
  56 Usage:   vcfutils.pl <command> [<arguments>]\n
  57 Command: subsam       get a subset of samples
  58 \n/);
  59 }