#!/usr/bin/perl -w
use strict;
use warnings;
use Getopt::Long;
use Data::Dumper;
use FindBin qw($Bin $Script);
use File::Basename qw(basename dirname);
my $BEGIN_TIME=time();
my $Time_Start = &sub_format_datetime(localtime($BEGIN_TIME));
print "Program Starts Time:$Time_Start\n";
my $version="1.x";
######################################################
# ------------------------------------------------------------------
# GetOptions
# ------------------------------------------------------------------
my ($fIn,$fOut,$read1,$read2);
GetOptions(
"help|?" =>\&USAGE,
"o:s"=>\$fOut,
"bam:s"=>\$fIn,
"fq1:s"=>\$read1,
"fq2:s"=>\$read2,
) or &USAGE;
&USAGE unless ($fIn and $fOut and $read1 and $read2);
my %Bam;
my $name1 =basename($read1);
my $unmap_name1 ="unmap_"."$name1";
my $name2 =basename($read2);
my $unmap_name2 ="unmap_"."$name2";
open (IN,"samtools view $fIn|") or die $!;
while (<IN>) {
chomp;
next if (/^$/);
my $id=(split /\s+/,$_)[0];
$id='@'.$id;
$Bam{$id}=1;
}
close(IN);
#print Dumper %Bam;
open (FQ1,$read1) or die $!;
open (Rice,">$fOut/$name1") or die $!;
open (Flounder,">$fOut/$unmap_name1") or die $!;
while (<FQ1>) {
chomp;
next if (/^$/);
my $id_inf=$_;
my $id=(split /\s+/,$id_inf)[0];
my $seq=<FQ1>;
my $fr=<FQ1>;
my $score=<FQ1>;
if (exists $Bam{$id}) {
print Rice "$id_inf\n$seq$fr$score";
}else{
print Flounder "$id_inf\n$seq$fr$score";
}
}
close(Rice);
close(Flounder);
close(FQ1);
open (FQ2,$read2) or die $!;
open (Rice,">$fOut/$name2") or die $!;
open (Flounder,">$fOut/$unmap_name2") or die $!;
while (<FQ2>) {
chomp;
next if (/^$/);
my $id_inf=$_;
my $id=(split /\s+/,$id_inf)[0];
my $seq=<FQ2>;
my $fr=<FQ2>;
my $score=<FQ2>;
if (exists $Bam{$id}) {
print Rice "$id_inf\n$seq$fr$score";
}else{
print Flounder "$id_inf\n$seq$fr$score";
}
}
close(Rice);
close(Flounder);
close(FQ2);
#######################################################################################
my $Time_End = sub_format_datetime(localtime(time()));
print STDOUT "Program Ends Time:$Time_End\nDone. Total elapsed time : ",time()-$BEGIN_TIME,"s\n";
#######################################################################################
# ------------------------------------------------------------------
# sub function
# ------------------------------------------------------------------
#######################################################################################
sub ABSOLUTE_DIR{ #$pavfile=&ABSOLUTE_DIR($pavfile);
my $cur_dir=`pwd`;chomp($cur_dir);
my ($in)=@_;
my $return="";
if(-f $in){
my $dir=dirname($in);
my $file=basename($in);
chdir $dir;$dir=`pwd`;chomp $dir;
$return="$dir/$file";
}elsif(-d $in){
chdir $in;$return=`pwd`;chomp $return;
}else{
warn "Warning just for file and dir\n";
exit;
}
chdir $cur_dir;
return $return;
}
#######################################################################################
sub max{#&max(lists or arry);
my $max=shift;
my $temp;
while (@_) {
$temp=shift;
$max=$max>$temp?$max:$temp;
}
return $max;
}
#######################################################################################
sub min{#&min(lists or arry);
my $min=shift;
my $temp;
while (@_) {
$temp=shift;
$min=$min<$temp?$min:$temp;
}
return $min;
}
#######################################################################################
sub revcom(){#&revcom($ref_seq);
my $seq=shift;
$seq=~tr/ATCGatcg/TAGCtagc/;
$seq=reverse $seq;
return uc $seq;
}
#######################################################################################
sub GetTime {
my ($sec, $min, $hour, $day, $mon, $year, $wday, $yday, $isdst)=localtime(time());
return sprintf("%4d-%02d-%02d %02d:%02d:%02d", $year+1900, $mon+1, $day, $hour, $min, $sec);
}
#######################################################################################
sub sub_format_datetime {#Time calculation subroutine
my($sec, $min, $hour, $day, $mon, $year, $wday, $yday, $isdst) = @_;
$wday = $yday = $isdst = 0;
sprintf("%4d-%02d-%02d %02d:%02d:%02d", $year+1900, $mon+1, $day, $hour, $min, $sec);
}
sub USAGE {#
my $usage=<<"USAGE";
ProgramName:
Version: $version
Usage:
Options:
-bam <file> input file,bam format,forced
-fq1 <file> input file,rawdata reads1,forced
-fq2 <file> input file,rawdata reads2,forced
-o <dir> output dir,forced
-h Help
USAGE
print $usage;
exit;
}
网友评论