美文网首页
把illumina比对后的bam转回双端fastq

把illumina比对后的bam转回双端fastq

作者: 九月_1012 | 来源:发表于2022-01-09 22:57 被阅读0次
#!/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;
}

相关文章

网友评论

      本文标题:把illumina比对后的bam转回双端fastq

      本文链接:https://www.haomeiwen.com/subject/wzmpcrtx.html