prod_reproduction/202307-IDTWES/script/hmm2pcov.pl

#!/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";
}