美文网首页
把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