perl/popuSNP/conasmP.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4 #use DBI;
   5 use Galaxy::ShowHelp;
   6
   7 $main::VERSION=0.1.3;
   8
   9 our $opts='i:o:c:m:bvne';
  10 our($opt_i, $opt_o, $opt_c, $opt_v, $opt_b, $opt_n, $opt_e, $opt_m);
  11
  12 our $desc='Genome creater for population analyse';
  13 our $help=<<EOH;
  14 \t-c Consensus file with refbase
  15 \t-i filtered Indel file
  16 \t-m merge.list for dividing mixed Chrs if any
  17 \t-o Output file prefix for the sample, directories must exist
  18 \t-n mask noncoverage bases to "n", otherwise little cased
  19 \t-e skip hete indels, otherwise both homo and hete(little cased) used
  20 \t-v show verbose info to STDOUT
  21 \t-b No pause for batch runs
  22 EOH
  23
  24 ShowHelp();
  25
  26 warn "[x]Must specify -c xx.consensus.txt !\n" unless $opt_c;
  27 warn "[x]Must specify -i xx.indel.txt.filter !\n" unless $opt_i;
  28 warn "[x]Must specify -o ./xx !\n" unless $opt_o;
  29 exit 1 unless $opt_i and $opt_o and $opt_c;
  30 if ($opt_m) {
  31     die "[x]Invalid merge.list: [$opt_m] !\n" unless (-s $opt_m);
  32 }
  33
  34 print STDERR "From [$opt_c][$opt_i] to [$opt_o]\n ";
  35 print STDERR "merge.list is [$opt_m]\n " if $opt_m;
  36 print STDERR "noncoverage mask to n.   " if $opt_n;
  37 print STDERR "hete indels will be skiped." if $opt_e;
  38 warn "\n";
  39 unless ($opt_b) {print STDERR 'press [Enter] to continue...'; <>;}
  40
  41 my (%Genome,%Depth,%Merged);
  42 #my %attr = (
  43 #    RaiseError => 0,
  44 #    PrintError => 1,
  45 #    AutoCommit => 0
  46 #);
  47 #my $dbh = DBI->connect('dbi:SQLite:dbname=:memory:','','',\%attr) or die $DBI::errstr;
  48 #my $sql=q/
  49 #CREATE TABLE IF NOT EXISTS merge
  50 #(  chrid TEXT,
  51 #   scaffold TEXT,
  52 #   start INTEGER,
  53 #   end INTEGER  );
  54 #/;
  55 #$dbh->do($sql) or die $dbh->errstr;
  56 #$dbh->commit;
  57 #my $sthi = $dbh->prepare( 'INSERT INTO merge ( chrid,scaffold,start,end ) VALUES ( ?,?,?,? )' );
  58 #my $stho = $dbh->prepare( 'SELECT DISTINCT scaffold,start,end FROM merge WHERE chrid=? AND ? BETWEEN start AND end' );
  59
  60 if ($opt_m) {
  61     open M,'<',$opt_m or die "Error opening $opt_m: $!\n";
  62     while (<M>) {
  63         chomp;
  64         my ($chrid,$scaffold,$start,$end)=split /\t/;
  65         #$sthi->execute($chrid,$scaffold,$start,$end);
  66         $Merged{$chrid}{$start}=[$scaffold,$end];
  67     }
  68     close M;
  69     warn "\n[!]merge.list done.\n";
  70 }
  71
  72 open CNS,'<',$opt_c or die "Error opening $opt_c: $!\n";
  73 my ($chr,$pos,$best,$depth,$ref,$indel,$bases,$homhet,$strand,$t);
  74 while (<CNS>) {
  75     chomp;
  76     ($chr,$pos,undef,$best,undef,undef,$depth,$ref) = split /\t/;
  77     unless ($ref) {   # if file not completed.
  78         print STDERR '.';
  79         next;
  80     }
  81     if ($depth==0) {
  82         unless ($opt_n) {$best=lc $ref;}
  83          else {$best='n';}
  84     } else {$Depth{$chr}->{$pos}=chr $depth;} # string is smaller than double
  85 # WARING: Something might be wrong if $depth > 255 !
  86 # However, people should not have such much money.
  87 #$best = $ref; # DEBUG ONLY !!!
  88     $Genome{$chr}->[$pos]=$best;
  89     #$Depth{$chr}->[$pos]=$depth;
  90 }
  91 close CNS;
  92 warn "[!]CNS done.\n";
  93
  94 open INDEL,'<',$opt_i or die "Error opening $opt_i: $!\n";
  95 while (<INDEL>) {
  96     ($chr,$pos,$indel,$bases,$strand,$homhet,undef,$depth,$t) = split /\t/;
  97     unless ($t) {   # if file not completed.
  98         print STDERR '\'';
  99         next;
 100     }
 101     next unless exists $Genome{$chr};
 102     if ($homhet eq 'hete') {
 103         unless ($opt_e) {$bases = lc $bases;}
 104          else {next;}
 105     }
 106     if ($indel =~ /^([ID])(\d+)$/) {
 107         if ($1 eq 'I') {
 108             $Genome{$chr}->[$pos] .= $bases;  # if length($Genome{$chr}->[$t]) == 1
 109             $Depth{$chr}{$pos} .= chr($depth) x $2;
 110         } else {
 111             my $i=$2-1;
 112             #++$pos;    # deletion including it.
 113             if ($opt_v) {
 114                 print "[D]$indel:${1}_${2}@","$pos $strand\t";
 115                 #print $Genome{$chr}->[$_] for ($pos..$pos+$i);
 116                 print(($t=$Genome{$chr}->[$_])?$t:'.') for ($pos-5..$pos-1);
 117                 print '<';
 118                 print(($t=$Genome{$chr}->[$_])?$t:'.',' ') for ($pos..$pos+$i);
 119                 print '>';
 120                 print(($t=$Genome{$chr}->[$_])?$t:'.') for ($pos+$i+1..$pos+$i+5);
 121                 print ' - ',$bases,"\n";
 122             }
 123             for $t ($pos..$pos+$i) {
 124                 $Genome{$chr}->[$t]='';
 125                 delete $Depth{$chr}{$pos};# if exists $Depth{$chr}{$pos};
 126             }
 127         }
 128     } else {print STDERR '|';next;}
 129 }
 130 close INDEL;
 131 warn "[!]INDEL done.\n";
 132
 133 my ($i,$dep,$scaf,$name);
 134 for $chr (keys %Genome) {
 135         $Genome{$chr}->[0]='';  # no longer undef
 136         $t=0;
 137         $i=-1;
 138         unless (exists $Merged{$chr}) { # exists should faster than defined ?
 139                 $name=$opt_o.'.'.$chr;
 140                 open OUT,'>',$name.'.fa' or die "Error opening ${name}.fa: $!\n";
 141                 print OUT ">$chr\n";
 142                 open DEP,'>',$name.'.dep' or die "Error opening ${name}.dep: $!\n";
 143                 print DEP ">$chr\n";
 144                 #$t=0;
 145                 #$i=-1;
 146                 for (@{$Genome{$chr}}) {
 147                         ++$i;
 148                         next if $_ eq '';
 149                         ++$t;
 150                         print OUT $_;
 151                         $dep=$Depth{$chr}{$i} or $dep="\0";
 152                         $dep = join(' ',map ord,split //,$dep);
 153                         #print DEP '[',$dep,']';
 154                         if ($t > 79) {
 155                                 $t=0;
 156                                 print OUT "\n";
 157                                 print DEP "$dep\n";
 158                         } else {print DEP $dep,' ';}
 159                 }
 160                 print OUT "\n";
 161                 print DEP "\n";
 162                 close OUT;
 163                 close DEP;
 164     } else {
 165                 $t=$opt_o.'.'.$chr;
 166                 system('mkdir','-p',$t);
 167                 my ($end,$scafname,$toPrint)=(-1,'',0); # -1 to fix 'print() on unopened filehandle OUT at ./conasm.pl line 178 & 179 of elsif ($i == $end).'
 168                 $scaf=$Merged{$chr};
 169                 for (@{$Genome{$chr}}) {
 170                         ++$i;
 171                         if (exists $$scaf{$i}) {        # Start point
 172                                 ($scafname,$end)=@{$$scaf{$i}};
 173                                 $name=$opt_o.'.'.$chr.'/'.$scafname;
 174                                 open OUT,'>',$name.'.fa' or die "Error opening ${name}.fa: $!\n";
 175                                 print OUT ">$scafname\n";
 176                                 open DEP,'>',$name.'.dep' or die "Error opening ${name}.dep: $!\n";
 177                                 print DEP ">$scafname\n";
 178                                 $toPrint=1;     # when set to 0, no output for manmade junctions.
 179                         }
 180                         if ($i != $end) {
 181                                 next if ($_ eq '' or $toPrint==0);
 182                                 ++$t;
 183                                 print OUT $_;
 184                                 $dep=$Depth{$chr}{$i} or $dep="\0";
 185                                 $dep = join(' ',map ord,split //,$dep);
 186                                 #print DEP '[',$dep,']';
 187                                 if ($t > 79) {
 188                                         $t=0;
 189                                         print OUT "\n";
 190                                         print DEP "$dep\n";
 191                                 } else {print DEP $dep,' ';}
 192                         } else {        # End point
 193                         if ($_ ne '') {
 194                                 print OUT $_,"\n";
 195                                 $dep=$Depth{$chr}{$i} or $dep="\0";
 196                                 $dep = join(' ',map ord,split //,$dep);
 197                                 print DEP "$dep\n";
 198                         } else {
 199                                 print OUT "\n";
 200                                 print DEP "\n";
 201                         }
 202                         close OUT;
 203                         close DEP;
 204                         $t=$toPrint=0;
 205                         next;
 206                         }
 207                 }
 208         }
 209 }
 210 warn "[!]Output done.\n";
 211
 212 __END__
 213 vf=11.9g for 1.3G Feb 20 07:23 consensus/QRS9.Gm03.txt, 47,781,076
 214
 215 find ./consensus/ -name '*.txt'|perl -lane '$m="";$vf="14g";$a=(split /\//)[-1];$b=(split /\./,$a)[0];open O,">./sh/$a.sh";if ($a=~/SGm/) {$m="-m soybean.merge.list ";$vf="4.7g"};print O "#\$ -N \"C_$a\"\n#\$ -cwd -r y -l vf=$vf,p=1\n#\$ -v PERL5LIB,PATH,PYTHONPATH,LD_LIBRARY_PATH\n#\$ -o /dev/null -e /dev/null\n./conasm.pl ${m}-i ./Indel/${b}.indel.txt.filter -c $_ -bno ./out/$b 2>./log/$a.log";close O'
 216
 217 $ cat sh/QRS21.Gm13.txt.sh
 218 #$ -N "C_QRS21.Gm13.txt"
 219 #$ -cwd -r y -l vf=14g,s_core=1
 220 #$ -v PERL5LIB,PATH,PYTHONPATH,LD_LIBRARY_PATH
 221 #$ -o /dev/null -e /dev/null
 222 ./conasm.pl -i ./Indel/QRS21.indel.txt.filter -c ./consensus/QRS21.Gm13.txt -bo ./out/QRS21 2>./log/QRS21.Gm13.txt.log
 223
 224 ./conasm.pl -m soybean.merge.list -i ./Indel/QRS29.indel.txt.filter -c ./consensus/QRS29.SGm1.txt -bno ./out/QRS29 2>./log/QRS29.SGm1.txt.log