perl/etc/radseq/plotphase.pl

   1 #!/usr/bin/env perl
   2 =pod
   3 Author: Hu Xuesong @ BIOPIC <galaxy001@gmail.com>
   4 Version: 1.0.0 @ 20120720
   5 Purpose: Read bcf, get tped for p-link
   6 Notes: rad2marker is deprecated.
   7 =cut
   8 use strict;
   9 use warnings;
  10 use Data::Dump qw(ddx);
  11 use Galaxy::IO;
  12 use Galaxy::SeqTools;
  13 use Data::Dumper;
  14
  15 die "Usage: $0 <tped prefix> [out].chrID.inp\n" if @ARGV<1;
  16 my $prefix=shift;
  17 my $outfs=shift;
  18 unless (defined $outfs) {
  19         warn "Using prefix[$prefix] for both input and output.\n";
  20         $outfs = $prefix;
  21 }
  22 warn "From [$prefix] to [$outfs]\n";
  23
  24 my (@Positions,%HapDat,@IDs,@IDab);
  25 open I,'<',$prefix or die "Error opening $prefix : $!\n";
  26 while (<I>) {
  27         if (/^Positions of loci: /) {
  28                 @Positions = split /\s+/,$_;
  29                 splice @Positions,0,3;
  30         }
  31         if (/^BEGIN BESTPAIRS1/) {
  32                 while (<I>) {
  33                         last if /^END BESTPAIRS1/;
  34                         chomp;
  35                         my $id = (split /\s/,$_)[1];
  36                         $id =~ s/^#//;
  37                         chomp($_ = <I>);
  38                         my @basesA = split /\s+/,$_;
  39                         chomp($_ = <I>);
  40                         my @basesB = split /\s+/,$_;
  41                         for my $base (@basesA,@basesB) {
  42                                 if ($base =~ /^\((\w)\)$/) {
  43                                         $base = lc $1;
  44                                         #last;
  45                                 }
  46                                 $base = '.'.lc($1) if ($base =~ /^\[(\w)\]$/);
  47                         }
  48                         push @IDs,$id;
  49                         push @IDab,"${id}_a";
  50                         push @IDab,"${id}_b";
  51                         $HapDat{$id} = [\@basesA,\@basesB];
  52                 }
  53         }
  54 }
  55 #ddx \%HapDat;
  56 close I;
  57
  58 open O,'>',$prefix.'.tsv' or die "Error opening $prefix.tsv : $!\n";
  59 print O "# From: [$prefix]\n",join("\t",'Pos',@IDab),"\n";
  60 for my $i (0 .. $#Positions) {
  61         my @tmp;
  62         push @tmp,$Positions[$i];
  63         for my $id (@IDs) {
  64                 push @tmp,$HapDat{$id}->[0]->[$i];
  65                 push @tmp,$HapDat{$id}->[1]->[$i];
  66         }
  67         print O join("\t",@tmp),"\n";
  68 }
  69 close O;
  70
  71 __END__