perl/etc/radseq/ped2phase.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 $Target = 'scaffold97,scaffold1457';
  25 my %Targets;
  26 my @tTarget = map { s/\s//g;$Targets{$_}=1;$_ } split /\,/,$Target;
  27 #ddx \%Targets,\@tTarget;
  28 #die;
  29
  30 my %mid2pos;
  31 open ID,'<',$prefix.'.dict' or die "Error opening $prefix.dict : $!\n";
  32 open IP,'<',$prefix.'.tped' or die "Error opening $prefix.tped : $!\n";
  33 open IF,'<',$prefix.'.tfam' or die "Error opening $prefix.tfam : $!\n";
  34
  35 my %GTdata;
  36 while (<ID>) {
  37         chomp;
  38         my ($chr,$pos,$id) = split /\t/;
  39         if (exists $Targets{$chr}) {
  40                 $mid2pos{$id} = [$chr,$pos];
  41                 #ddx $mid2pos{$id};
  42         }
  43 }
  44 close ID;
  45
  46 my ($NumberOfIndividuals,@Individuals);
  47 while (<IF>) {
  48         my @tmp = split /\t/;
  49         push @Individuals,$tmp[1];
  50 }
  51 $NumberOfIndividuals = @Individuals;
  52 close IF;
  53
  54 while (<IP>) {
  55         chomp;
  56         my ($chrNO,$id,undef,$pos,@dat) = split /\t/;
  57         if (exists $mid2pos{$id}) {
  58                 #print join(',',@{$mid2pos{$id}},$chrNO,$id,$pos,@dat),"\n";
  59                 my ($chrid,$chrpos) = @{$mid2pos{$id}};
  60                 #$NumberOfIndividuals = @dat;
  61                 for (@dat) {
  62                         s/ //;
  63                         s/0/\?/g;
  64                         #$_ = '??' if $_ eq '00';
  65                 }
  66                 $GTdata{$chrid}{$chrpos} = \@dat;
  67         }
  68 }
  69 close IP;
  70
  71 #ddx \%GTdata;
  72
  73 for my $chrid (keys %GTdata) {
  74         my @Locus = sort { $a <=> $b } keys %{$GTdata{$chrid}};
  75         my $NumberOfLoci = @Locus;
  76         open O,'>',"$outfs.$chrid.inp" or die "Error opening $outfs.$chrid.inp : $!\n";
  77         print O "$NumberOfIndividuals\n$NumberOfLoci\nP ";
  78         print O join(' ',@Locus),"\n",'S' x $NumberOfLoci,"\n";
  79         for my $i (0 .. $#Individuals) {
  80                 print O "#$Individuals[$i]\n";
  81                 my ($m,$n,@t);
  82                 for (@Locus) {
  83                         @t = split //,$GTdata{$chrid}{$_}->[$i];
  84                         $m .= $t[0];
  85                         $n .= $t[1];
  86                 }
  87                 print O "$m\n$n\n";
  88         }
  89
  90         close O;
  91 }
  92 __END__
  93 perl ped2phase.pl sw000-18
  94 ./hp/phase.2.1.1.source/PHASE sw000-18.scaffold1457.inp p18s1457 >p18s1457.log 2>p18s1457.err &
  95 ./hp/phase.2.1.1.source/PHASE sw000-18.scaffold97.inp p18s97 >p18s97.log 2>p18s97.err &