4 use IO
::Unread
qw(unread);
5 use Data
::Dump
qw(ddx);
7 die "Usage: $0 <fq.gz path> <out path>\n" if @ARGV < 2;
11 my $BwaSai = "bwa aln -l 17 -q 10 $Ref";
12 my $BwaSam = "bwa sampe -a 800 $Ref";
13 my $BwaSort = 'perl pickU_sort.pl';
14 my $BwaRmDup = 'samtools rmdup';
15 my $readsCounter = '/share/users/huxs/git/toGit/c_cpp/faststater/readsCounter';
17 opendir my($dh), $inp or die "Couldn't open dir '$inp': $!";
18 my @files = readdir $dh;
20 # http://perlmeme.org/faqs/file_io/directory_listing.html
22 @files = grep(/\.f(ast|)q\b/i, @files);
26 my (%Pairs,@AllTargets);
28 m
~^([^/]+?
)([\W_
]*?
)([12])?
(\
.f
(ast
|)q
(\
.gz
)?
)$~i
or die $_;
29 my ($m,$rp,$r,$ext) = ($1,$2,$3,$4);
30 $m =~ s/\W//g and die $m;
31 print "$1, $2, $3, $4, $5, $6, $7\t$m,$rp,$r,$ext, $_\n";
32 die unless ($r == 1) or ($r==2); # no SE now
33 push @
{$Pairs{$m}},[$_,$r];
35 @AllTargets = sort keys %Pairs;
40 open M
,'>','Makefile' or die $!;
41 print M
'all: ',join(' ',@AllTargets,'_AdditionalTG_'),"\n";
48 for my $Target ( @AllTargets ) {
49 my @Files = @
{$Pairs{$Target}};
50 die if @Files != 2; # no SE now
51 my (@fqFiles,@newFiles);
53 my ($fqname,$read12) = @
$_;
54 $newFiles[$read12-1] = "${Target}_$read12";
55 $fqFiles[$read12-1] = $fqname;
57 sai/${Target}_$read12.sai: $inp$fqname
58 \t$BwaSai $inp$fqname >sai/${Target}_$read12.sai 2>sai/${Target}_$read12.logsai
59 \tdate >> sai/${Target}_$read12.logsai
61 sai/${Target}_$read12.fqstat: $inp$fqname
62 \t$readsCounter -o sai/${Target}_$read12.fqstat $inp$fqname
64 push @AdditionalTG, "sai/${Target}_$read12.fqstat";
67 sam/${Target}.sam.gz: ",join(' ',map { "sai/$_.sai" } @newFiles),"
68 \t$BwaSam ",join(' ',(map { "sai/$_.sai" } @newFiles),(map { "$inp$_" } @fqFiles) )," 2>sam/${Target}.logsam |gzip -9c >sam/${Target}.sam.gz
69 \tdate >> sam/${Target}.logsam
72 sam/${Target}.sort.bam: sam/${Target}.sam.gz
73 \t$BwaSort sam/${Target}.sam.gz sam/${Target}.sort.bam
76 bam/${Target}.rmdup.bam: sam/${Target}.sort.bam
77 \t$BwaRmDup sam/${Target}.sort.bam bam/${Target}.rmdup.bam
79 print M
"${Target}: bam/${Target}.rmdup.bam\n";
82 print M
join(' ','_AdditionalTG_:',@AdditionalTG),"\n";
86 perl getmakefln
.pl fq
/ .