prod_reproduction/202307-IDTWES/script/plot/get_roh.pl

#!/usr/bin/perl
use strict;
use warnings;

my $sam=shift;
my $vaf=$sam.".ps.tab";
my $roh=$sam."_ROHTableCall.bed";
my (%dat,$num);
open my $rohfh, '<', $roh or die "Cannot open file: $roh $!\n";
while(<$rohfh>){
       chomp;
       my ($chr, $start, $end, $flag, $snp)=split /\t/;
       $chr=~s/chr//;
       if($snp>150 or ($end-$start > 5000000 and $snp>50 ) ){
           $num++;
           $dat{$chr}{$start.'_'.$end}='ROH'.$num;
       }
}

open my $vaffh, '<', $vaf or die "Cannot open file:$vaf $!\n";
open my $wfh, '>' , $vaf.".xls" or die "Cannot write file:$! \n";
$num=0;
while(<$vaffh>){
       chomp;
       next if /^#/;
       if(/^contig/){
              print $wfh $_,"\tVAFadj\tROH\n";
              next;
       }
       my ($chr,$pos, $nref, $nalt, $other, $af)=split /\t/;
       $chr=~s/chr//;
       next if $nalt+$nref<40;
       my $vaf=$nalt/($nalt+$nref);
       $vaf=$vaf >0.5 ? 1-$vaf : $vaf;
       $vaf=sprintf("%.3f", $vaf);
       my $roh;
       for my $reg (keys %{$dat{$chr} }){
            my ($seg_start,$seg_end)=split /_/,$reg;
            if($pos>= $seg_start and $pos <= $seg_end){
                $roh=$dat{$chr}{$reg};
                ($num)=$roh=~/ROH(\d+)$/;
                last;
            }
       }
       if (defined $roh){
              print $wfh $_,"\t", $vaf,"\t", $roh,"\n";
       }else{
              print $wfh $_,"\t", $vaf, "\t", 'HET'.$num,"\n";
       }
       
}