perl/popuSNP/snplots.pl

   1 #!/bin/env perl
   2 use lib '/share/raid010/resequencing/resequencing/tmp/bin/annotation/glfsqlite';
   3 use strict;
   4 use warnings;
   5 use DBI;
   6 use Time::HiRes qw ( gettimeofday tv_interval );
   7 use Galaxy::ShowHelp;
   8
   9 $main::VERSION=0.0.2;
  10
  11 our $opts='i:o:s:l:m:bv';
  12 our ($opt_i, $opt_o, $opt_s, $opt_v, $opt_b, $opt_l, $opt_m);
  13
  14 our $help=<<EOH;
  15 \t-i PSNP list (./npsnp.lst) for chrid.individual.finalSNPs
  16 \t-l scaffold.list (./scaffold.N90.id)
  17 \t-m merge.list for scaffold positions (./watermelon.merge.list)
  18 \t-s GLF list (./glf.list), will use \$1 of (([^/]+)/[^/]+\$) for sample names
  19 \t-o Output Prefix (./plots/p_)
  20 \t-v show verbose info to STDOUT
  21 \t-b No pause for batch runs
  22 EOH
  23
  24 ShowHelp();
  25
  26 $opt_o='./plots/p_' unless defined $opt_o;
  27 $opt_i='./npsnp.lst' unless $opt_i;
  28 $opt_s='./glf.list' unless $opt_s;
  29 $opt_l='./scaffold.N90.id' unless $opt_l;
  30 $opt_m='./watermelon.merge.list' unless $opt_m;
  31
  32 $opt_i =~ s/\/$//;
  33 print STDERR "From [$opt_i] to [$opt_o], with [$opt_s][$opt_l][$opt_m]\n";
  34 if (! $opt_b) {print STDERR 'press [Enter] to continue...'; <>;}
  35
  36 my $start_time = [gettimeofday];
  37
  38 system('mkdir','-p',$opt_o);
  39 system('rmdir',$opt_o) if $opt_o =~ m#/[\s.]*[^/\s.]+[^/]*$#;
  40
  41 my @Samples;
  42 open L,'<',$opt_s or die "[x]Error opening $opt_s: $!\n";
  43 print STDERR "[!]Sample Order: ";
  44 while (<L>) {
  45         m#([^/]+)/[^/]+$#;
  46         push @Samples,$1;
  47         print STDERR (scalar @Samples),':[',$1,"] ";
  48 }
  49 close L;
  50 print STDERR "\n";
  51
  52 my ($LEN,%Scaffolds,@Scaffoldo)=(0);
  53 open L,'<',$opt_l or die "[x]Error opening $opt_l: $!\n";
  54 print STDERR "[!]Scaffolds: ";
  55 while (<L>) {
  56         chomp;
  57         --$Scaffolds{$_};
  58         #print STDERR (scalar (keys %Scaffolds)),':[',$_,"] ";
  59 }
  60 close L;
  61 print STDERR scalar (keys %Scaffolds),"\n";
  62
  63 open L,'<',$opt_m or die "[x]Error opening $opt_m: $!\n";
  64 print STDERR "[!]Scaffolds Len: ";
  65 my ($i,$len)=(0);
  66 while (<L>) {
  67         chomp;
  68         my ($chrid,$scaffold,$start,$end)=split /\t/;
  69         next unless exists $Scaffolds{$scaffold};
  70         $len=$end-$start+1;
  71         $Scaffolds{$scaffold}=[$len,0];
  72         $LEN += $len;
  73 warn "$chrid,$scaffold,$start,$end\t$len\t$LEN\n" if $opt_v;
  74         ++$i;
  75 }
  76 close L;
  77 print STDERR "$LEN\n";
  78 @Scaffoldo = sort {${$Scaffolds{$b}}[0] <=> ${$Scaffolds{$a}}[0]} keys %Scaffolds;
  79 my $i=0;
  80 for (@Scaffoldo) {
  81         ${$Scaffolds{$_}}[1]=$i;
  82         $i += ${$Scaffolds{$_}}[0];
  83 warn "$_\t${$Scaffolds{$_}}[0]\t${$Scaffolds{$_}}[1]\n"  if $opt_v;
  84 }
  85
  86 my (%Plot,$POS);
  87 print STDERR "[!]Parsing SNP ";
  88 open P,'<',$opt_i or die "[x]Error opening $opt_i: $!\n";
  89 while (my $file=<P>) {
  90         chomp $file;
  91         open SNP,'<',$file or (warn "\n[!]Error opening $file: $!\n" and next);
  92         print STDERR ".\b";
  93         while (<SNP>) {
  94                 my ($scaffold,$pos,$ref,$tail)=split /\t/;
  95                 next unless exists $Scaffolds{$scaffold};
  96                 $POS = ${$Scaffolds{$_}}[1] + $pos;
  97                 my @Values;
  98                 for (split / /,$tail) {
  99                         next unless /[ACGTRYMKSWHBVDNX-]/;
 100                         if (/[ACGTRYMKSWHBVDNX]/) {
 101                                 push @Values,;
 102                         } elsif ($_ eq '-') {
 103                                 ;
 104                         } else {next;}
 105                 }
 106         }
 107         close SNP;
 108 }
 109 close P;
 110
 111 my $stop_time = [gettimeofday];
 112
 113 print STDERR "\nTime Elapsed:\t",tv_interval( $start_time, $stop_time )," second(s).\n";
 114