perl/popuSNP/nametrans/transnamed.pl

   1 #!/home/huxuesong/gentoo/usr/bin/perl -w
   2 #use threads;
   3 use strict;
   4 use warnings;
   5 use DBI;
   6 use Time::HiRes qw ( gettimeofday tv_interval );
   7 use Galaxy::ShowHelp;
   8 $main::VERSION=0.1.0;
   9
  10 our $opts='i:o:m:bv';
  11 our($opt_i, $opt_o, $opt_v, $opt_b, $opt_m);
  12
  13 #our $desc='SoapSort library PCR PE Duplication Remover & Merger (Atom Edition)';
  14 our $help=<<EOH;
  15 \t-i Input file, in format: /^ChrID\\tPos\\t.*\\n\$/
  16 \t-m merge list, in format: /^ChrID\\tscaffoldID\\tChrStart\\tChrStop\\n\$/
  17 \t-o Output file
  18 \t-v show verbose info to STDOUT
  19 \t-b No pause for batch runs
  20 EOH
  21
  22 ShowHelp();
  23
  24 warn "[x]Must specify -i !\n" unless $opt_i;
  25 warn "[x]Must specify -m !\n" unless $opt_m;
  26 warn "[x]Must specify -o !\n" unless $opt_o;
  27 exit 1 unless $opt_i and $opt_o and $opt_m;
  28
  29 print STDERR "From [$opt_i] with [$opt_m] to [$opt_o]\n";
  30 unless ($opt_b) {print STDERR 'press [Enter] to continue...'; <>;}
  31
  32 my $start_time = [gettimeofday];
  33 #BEGIN
  34 unless (-s $opt_i) {die "[x]Soaplist [$opt_i] is nothing !\n";}
  35
  36 ### Begin SQL
  37 my %attr = (
  38     RaiseError => 0,
  39     PrintError => 1,
  40     AutoCommit => 0 );
  41 my $dbh = DBI->connect('dbi:SQLite:dbname=:memory:','','',\%attr) or die $DBI::errstr;
  42
  43 # ChrNew2 scaffold90      1       45293
  44 my $sql=q/
  45 CREATE TABLE IF NOT EXISTS merge
  46 (  chrid TEXT,
  47    scaffold TEXT,
  48    start INTEGER,
  49    end INTEGER  );
  50 /;
  51 # The rowid value can be accessed using one of the special names "ROWID", "OID", or "_ROWID_".
  52 for (split /;/,$sql) {
  53         next if /^\s*$/;
  54         $dbh->do($_) or die $dbh->errstr;
  55 }
  56 $dbh->commit;
  57
  58 sub doindex() {
  59         my $sql=q/
  60 CREATE INDEX IF NOT EXISTS se ON merge(start,end);
  61 CREATE INDEX IF NOT EXISTS c ON merge(chrid);
  62 /;
  63         for (split /;/,$sql) {
  64                 next if /^\s*$/;
  65                 $dbh->do($_) or warn $dbh->errstr;
  66         }
  67         $dbh->commit;
  68 }
  69
  70 my $sthi = $dbh->prepare( 'INSERT INTO merge ( chrid,scaffold,start,end ) VALUES ( ?,?,?,? )' );
  71 my $stho = $dbh->prepare( 'SELECT DISTINCT scaffold,start,end FROM merge WHERE chrid=? AND ? BETWEEN start AND end' );
  72 $|=1;
  73
  74 open SAMPLE,'-|',"gzip -dc $opt_m" or die "Error opening $opt_m: $!\n";
  75 while (<SAMPLE>) {
  76         chomp;
  77         my ($chrid,$scaffold,$start,$end)=split /\t/;
  78         $sthi->execute($chrid,$scaffold,$start,$end);
  79 }
  80 close SAMPLE;
  81 &doindex;
  82
  83 open IN,'<',$opt_i or die "Error opening $opt_i: $!\n";
  84 open OUT,'>',$opt_o or die "Error opening $opt_o: $!\n";
  85 my ($chrid,$pos,$pos2,$scaffold,$scaffold2,$start,$end,$qres,$newpos,$newpos2,@last);
  86 while (<IN>) {
  87 # ChrNew1       121900  T T T T W W T T T T T T T T T T W W T T T T T T T - T T T T T T T T T T
  88         ($chrid,$pos,$pos2,@last)=split /\t/;
  89         $stho->execute($chrid,$pos);
  90         $qres = $stho->fetchall_arrayref;
  91         if ($#$qres == -1) {
  92                 warn "No info. for $chrid:$pos, skipped.\n";
  93                 next;
  94         } elsif ($#$qres != 0) { warn "More than 1 hit, using first.\n" }
  95         ($scaffold,$start,$end)=@{$$qres[0]};
  96         $newpos=1+$pos-$start;
  97
  98         $stho->execute($chrid,$pos2);
  99         $qres = $stho->fetchall_arrayref;
 100         if ($#$qres == -1) {
 101                 warn "No info. for $chrid:$pos, skipped.\n";
 102                 next;
 103         } elsif ($#$qres != 0) { warn "More than 1 hit, using first.\n" }
 104         ($scaffold2,$start,$end)=@{$$qres[0]};
 105         $newpos2=1+$pos2-$start;
 106         if ($scaffold ne $scaffold2) {
 107                 print OUT "[$scaffold,$scaffold2] ";
 108                 warn "Diff scaffold happens. ";
 109         }
 110         print OUT join("\t",$scaffold,$newpos,$newpos2,@last);
 111 }
 112 close OUT;
 113 close IN;
 114 #END
 115 my $stop_time = [gettimeofday];
 116
 117 print STDERR "\nTime Elapsed:\t",tv_interval( $start_time, $stop_time )," second(s).\n";
 118 __END__
 119