tools/linkage/getSeqbyChr.pl

   1 #!/usr/bin/perl -w
   2 use strict;
   3 use warnings;
   4 use lib '/nas/RD_09C/resequencing/soft/lib';
   5 use Time::HiRes qw ( gettimeofday tv_interval );
   6 use Galaxy::ShowHelp;
   7
   8 $main::VERSION=0.0.1;
   9
  10 our $opts='i:o:n:g:v:r:c:f:h:bqd';
  11 our($opt_i, $opt_n, $opt_g, $opt_r, $opt_c, $opt_f, $opt_o, $opt_v, $opt_b, $opt_q, $opt_d, $opt_h);
  12
  13 our $help=<<EOH;
  14 \t-i soaps.lst
  15 \t-c Chromosome NFO file (chr.nfo) in format: /^ChrID\\s+ChrLen\\s?.*\$/
  16 \t-o Output path (./ffq) for [ChrID/Sample.ffq], will mkdir if not exist
  17 \t-v Verbose level (undef=0)
  18 \t-b No pause for batch runs
  19 \t-d Debug Mode, for test only
  20 EOH
  21
  22 ShowHelp();
  23
  24 die "[x]-i $opt_i not exists !\n" unless -f $opt_i;
  25 $opt_c='chr.nfo' if ! $opt_c;
  26 $opt_o='./ffq' if ! $opt_o;
  27
  28 no warnings;
  29 $opt_v=int $opt_v;
  30 use warnings;
  31
  32 $opt_o=~s#/+$##;
  33
  34 print STDERR "From [$opt_i] to [$opt_o] refer to [$opt_c]\n";
  35 print STDERR "DEBUG Mode on !\n" if $opt_d;
  36 print STDERR "Verbose Mode [$opt_v] !\n" if $opt_v;
  37 unless ($opt_b) {print STDERR 'press [Enter] to continue...'; <>;}
  38
  39 my $start_time = [gettimeofday];
  40 #BEGIN
  41 system('mkdir','-p',$opt_o);
  42 my %ChrOutFH;
  43 open ChrNFO,'<',$opt_c or die "[x]Error opening $opt_c: $!\n";
  44 while (<ChrNFO>) {
  45         next if /^#/;
  46         my ($chrid)=split /\t/;
  47         open $ChrOutFH{$chrid},'|-',"gzip -9c - >$opt_o/$chrid.ffq.gz" or die "[x]Error opening $opt_o/$chrid.ffq with gzip: $!\n";
  48 }
  49 close ChrNFO;
  50 open $ChrOutFH{'__UnKnown__'},'|-',"gzip -9c - >$opt_o/__UnKnown__.ffq.gz" or die "[x]Error opening $opt_o/__UnKnown__.ffq with gzip: $!\n";
  51
  52 open LST,'<',$opt_i or die "[x]Error opening $opt_i: $!\n";
  53 while (<LST>) {
  54         chomp;
  55         my ($PESE,$sample,$lib,$FL,$ReadLen,$nfofpath)=split /\t/;
  56         $nfofpath =~ s/\.nfo$//;
  57         my @Files;
  58         if ($PESE eq 'PE') {
  59                 @Files=("$nfofpath.soap","$nfofpath.single");
  60         } else {
  61                 @Files=("$nfofpath.soap");
  62         }
  63         unless (-s "$nfofpath.unmap") {
  64                 open FQ,'-|',"gzip -dc $nfofpath.unmap.gz" or warn "[x]Error opening $PESE [$nfofpath.unmap.gz] with gzip: $!\n" and next;
  65         } else { open FQ,'<',"$nfofpath.unmap" or warn "[x]Error opening $PESE [$nfofpath.unmap]: $!\n" and next; }
  66         my ($name,$seq,$Qstr);
  67         while (<FQ>) {
  68                 chomp;
  69                 s/^[@>]//;
  70                 $name=$_;
  71                 chomp($seq=<FQ>);
  72                 <FQ>;
  73                 chomp($Qstr=<FQ>);
  74                 print ${$ChrOutFH{'__UnKnown__'}} join("\t",$FL,$name,'.','0','0',$seq,$Qstr,$sample,$lib),"\n";
  75         }
  76         close FQ;
  77         for my $file (@Files) {
  78                 unless (-s $file) {
  79                         open SP,'-|',"gzip -dc $file.gz" or warn "[x]Error opening $PESE [$file.gz] with gzip: $!\n" and next;
  80                         print STDERR '[!]',join(', ',$sample,$lib,$FL,"$file.gz");
  81                 } else {
  82                         open SP,'<',$file or warn "[x]Error opening $PESE [$file]: $!\n" and next;
  83                         print STDERR '[!]',join(', ',$sample,$lib,$FL,$file);
  84                 }
  85                 while (<SP>) {
  86                         my ($soapid,$seq,$Qstr,$hit,$strand,$chr,$pos) = (split(/\t/))[0,1,2,3,6,7,8];
  87                         next unless defined $pos;
  88                         $chr='__UnKnown__' unless exists $ChrOutFH{$chr};
  89                         if ($hit == 1) {
  90                                 print ${$ChrOutFH{$chr}} join("\t",$FL,$soapid,$chr,$pos,$hit,$seq,$Qstr,$sample,$lib),"\n";
  91                         } else {
  92                                 print ${$ChrOutFH{'__UnKnown__'}} join("\t",$FL,$soapid,$chr,$pos,$hit,$seq,$Qstr,$sample,$lib),"\n";
  93                         }
  94                 }
  95                 warn ", done.\n";
  96         }
  97 }
  98 close LST;
  99
 100 close $ChrOutFH{$_} for keys %ChrOutFH;
 101
 102 #END
 103 my $stop_time = [gettimeofday];
 104
 105 print STDERR "\nTime Elapsed:\t",tv_interval( $start_time, $stop_time )," second(s).\n";
 106 __END__