perl/etc/radseq/ped2hap.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 Data::Dumper;
  12
  13 =pod
  14 206000000
  15 208386088
  16 211361701
  17 214812903
  18 216154106
  19 220000000
  20 =cut
  21 my ($borderL,$borderR) = (206000000,220000000);
  22
  23 die "Usage: $0 <hap scaffold list> <tped sets prefix> <out>\n" if @ARGV<3;
  24 my $lstfs=shift;
  25 my $infs=shift;
  26 my $outfs=shift;
  27
  28 my (%Stat,$t);
  29 open L,'<',$lstfs or die "Error opening $lstfs : $!\n";
  30 my (%Scaff,@Scaffs);
  31 while (<L>) {
  32         chomp;
  33         push @Scaffs,$_ unless exists $Scaff{$_};
  34         ++$Scaff{$_};
  35 }
  36 close L;
  37 print "Total ",scalar keys %Scaff," Scaffolds: [",join('],[',sort keys %Scaff),"]\n";
  38
  39 open LOG,'>',"$outfs.log" or die "Error opening $outfs.log : $!\n";
  40 print LOG "Total ",scalar keys %Scaff," Scaffolds: [",join('],[',sort keys %Scaff),"]\n";
  41
  42 $t=1;
  43 for (@Scaffs) {
  44         $Scaff{$_} = $t;
  45         print LOG "$_: $t\n";
  46         ++$t;
  47 }
  48 die if $t > 22;
  49
  50 open D,'<',$infs.'.dict' or die "Error opening $infs.dict : $!\n";
  51 my %Rid;
  52 while (<D>) {
  53         chomp;
  54         my ($scaffid,$pos,$rid) = split /\t/;
  55         #if (exists $Scaff{$scaffid}) {
  56         if (exists $Scaff{$scaffid} and $pos>=$borderL and $pos<=$borderR) {
  57                 $Rid{$rid} = $scaffid;
  58         }
  59 }
  60 close D;
  61 print "Total ",scalar keys %Rid," Markers.\n";
  62
  63 for my $type (qw[case control all]) {
  64         my $inprefix;
  65         if ($type eq 'all') {
  66                 $inprefix = $infs;
  67         } else {
  68                 $inprefix = "${infs}.$type";
  69         }
  70         for my $scaffid (keys %Scaff) {
  71                 open C,'<',"$inprefix.tped" or die "Error opening $inprefix.tped : $!\n";
  72                 open O,'>',"$outfs.$type.$scaffid.tped" or die "Error opening $outfs.$type.$scaffid.tped : $!\n";
  73                 #symlink "$inprefix.tfam","$outfs.$type.$scaffid.tfam";
  74                 while (<C>) {
  75                         my @items = split /\t/;
  76                         if (exists $Rid{$items[1]} and $Rid{$items[1]} eq $scaffid) {
  77                                 $items[0] = $Scaff{$Rid{$items[1]}};
  78                                 print O join("\t",@items);
  79                         }
  80                 }
  81                 close C;
  82                 close O;
  83                 my $cmd = "p-link --tped \'$outfs.$type.$scaffid.tped\' --tfam $inprefix.tfam --out \'${outfs}$type.$scaffid\' --recode";
  84                 print "Running [$cmd].\n";
  85                 print LOG "$cmd\n";
  86                 system($cmd);
  87         }
  88 }
  89 close LOG;
  90
  91 __END__
  92
  93 ./ped2hap.pl tighap.lst radseq17 dohap17
  94
  95 p-link --tfile dohap17.case --out dohap17case --recode