1234567891011121314151617181920212223242526272829303132333435363738394041424344454647 |
- #!/usr/bin/perl
- use strict;
- use warnings;
- my (@cols,$total,$total_read,$nline);
- my $hmmf=shift;
- my $out=shift;
- my $sam=$hmmf;
- $sam=~s/\..+$//;
- open my $hfh, '<' , $hmmf or die "Cannot open file:$hmmf \n$!";
- while(<$hfh>){
- chomp;
- if(/^space\s+/){
- my @conts=split /\t/;
- for my $i (0..$#conts){
- push @cols ,$i if $conts[$i] eq 'space' or $conts[$i] eq 'start'
- or $conts[$i] eq 'end' or $conts[$i] eq 'gene'
- or $conts[$i] eq 'cor.gc' or $conts[$i] eq 'reads';
- }
- next;
- }
- my ($chr,$start,$end,$reads,$gene,$val)=(split /\t/)[@cols];
- $total += $val if $val=~/\d+\.?\d*$/;
- $total_read+=$reads;
- $nline++;
- }
- close($hfh);
- open $hfh, '<' , $hmmf or die "Cannot open file:$hmmf \n$!";
- open my $wfh , '>', $out or die "Cannot write file:$!\n";
- print $wfh join("\t",'contig','start','stop','name',$sam),"\n";
- while(<$hfh>){
- chomp;
- next if /^space\s+/;
- my ($chr,$start,$end,$reads,$gene,$val)=(split /\t/)[@cols];
- my $pval;
- if($val=~/\d+\.?\d*$/){
- $pval=sprintf("%.9f",$val/$total);
- }else{
- $pval=sprintf("%.9f",$reads/$total_read);
- }
- print $wfh join("\t",$chr,$start,$end,$gene,$pval ),"\n";
- }
|