perl/etc/mergefa.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <mergedfa> <inputfa>\n" if @ARGV != 2;
   6 my ($out,$in)=@ARGV;
   7 warn "From [$in] to [$out]\n";
   8
   9 sub openfile($) {
  10     my ($filename)=@_;
  11     my $infile;
  12     if ($filename=~/.bz2$/) {
  13             open( $infile,"-|","bzip2 -dc $filename") or die "Error opening $filename: $!\n";
  14     } elsif ($filename=~/.gz$/) {
  15         open( $infile,"-|","gzip -dc $filename") or die "Error opening $filename: $!\n";
  16     } else {open( $infile,"<",$filename) or die "Error opening $filename: $!\n";}
  17     return $infile;
  18 }
  19
  20 my $fh=&openfile($in);
  21 my $SEQ='';
  22 my $Ncnt=0;
  23 while (<$fh>) {
  24     s/^>//;
  25         /^(\S+)/ or next;
  26         my $seqname = $1;
  27     #print STDERR " >$seqname ...";
  28         $/=">";
  29         my $genome=<$fh>;
  30         chomp $genome;
  31         $genome=~s/\s//g;
  32         $/="\n";
  33     $SEQ .= $genome;
  34     my $n=($genome=~s/[^ATCG]/A/ig);
  35     $Ncnt += $n;
  36     #print STDERR "\b\b\b   \t",length $genome,".\n";
  37         $genome='';
  38 }
  39
  40 open OUT,'>',$out or die "Error opening $out: $!\n";
  41 print OUT ">merged $Ncnt\n$SEQ\n";
  42 close OUT;
  43
  44 __END__
  45 ./bwa aln -n 3 -o 0 -I ./ant/ant_meg_free.fa ./ant/ant_error_free_100_500_1.fa >./ant/ant_meg_free.3.sai 2>./ant/ant_meg_free.3.sai.log
  46 ./bwa samse -n 5 -f ./ant/ant_meg_free.3.sam ./ant/ant_meg_free.fa ./ant/ant_meg_free.3.sai ./ant/ant_error_free_100_500_1.fa 2>./ant/ant_meg_free.3.sam.log
  47
  48 XT:A:R
  49
  50 XA:Z:merged,+10484501,100M,0;
  51 XA:Z:merged,-15123200,100M,1;
  52
  53 find . -name '*meg*.fa'|perl -lane '/\/(\w+)\//;next if /ant/;system "./bwa index $_"'
  54 find . -name '*_1.fa'|perl -lane '/\/(\w+)\//;$p=$1;$n=$_;s/_error_([^_]+)_.+/_err_$1/;$err=$1;$m=$_;print STDERR "${m}_?";print "./bwa aln -n $_ -o 0 -N -I ./$p/${p}_meg_$err.fa $n >$m.$_.sai 2>$m.$_.sai.log" for (3,6,10,2,5,1,20)' > sai.cmd
  55
  56 find . -name '*_1.fa'|perl -lane '/\/(\w+)\//;$p=$1;$n=$_;s/_error_([^_]+)_.+/_err_$1/;$err=$1;$m=$_;print STDERR "${m}_?";print "./bwa samse -n 500 -f $m.$_.sam ./$p/${p}_meg_$err.fa $m.$_.sai $n 2>$m.$_.sam.log" for (3,6,10,2,5,1,20)' > sam.cmd
  57
  58 SAMFILE=$(echo $SEED|awk '{print $6}')
  59
  60 grep -P '\tXT:A:R\t' $SAMFILE >$SAMFILE.r
  61
  62 ./bwa aln -n 6 -o 0 -I ./ant/ant_meg_free.fa ./ant/ant_error_free_100_500_1.fa >./ant/ant_meg_free.6.sai 2>./ant/ant_meg_free.6.sai.log
  63 ./bwa aln -n 10 -o 0 -I ./ant/ant_meg_free.fa ./ant/ant_error_free_100_500_1.fa >./ant/ant_meg_free.10.sai 2>./ant/ant_meg_free.10.sai.log
  64
  65 ./bwa samse -n 500 -f ./ant/ant_err_1%.3.sam ./ant/ant_meg_1%.fa ./ant/ant_err_1%.3.sai ./ant/ant_error_1%_100_500_1.fa 2>./ant/ant_err_1%.3.sam.log