prod_reproduction/202307-FMBY2/yesuan/prep_FeBY2_plus.pl

#!/usr/bin/perl

use strict;
use warnings;

my $usage="\n\tperl $0 SampleName panel_type Ref_file\n\n";
die $usage if @ARGV<3;
my $sam=shift;
my $panel=shift;
my $ref_file=shift;
my (%call,@cols);
my $flag=0;
my %hash=('FeBY2'      => "Yesuan,Xueshuan",
          'FeBY2_plus' => "Yesuan,Xueshuan,Zigong,Luanchao",
          'IDT_WES'    => "Yesuan,Xueshuan",
          'MBY2'    => "Yesuan"
          );
unless( exists $hash{$panel}){
    print "NO Genotyping for $panel.\n";
    exit;
}

my $type=$hash{$panel};

open my $call_fh, '<' , $sam.'.ye.norm.tab' or die "Cannot open $sam.ye.norm.tab\n$!";
while(<$call_fh>){
    chomp;
    if(/^CHROM/){
        my @cont=split /\t/;
        for my $i(0..$#cont){
            push @cols,$i if $cont[$i] eq 'CHROM' or $cont[$i] eq 'POS' 
                          or $cont[$i] =~ /GT$/ or $cont[$i] =~ /GQ$/ 
                          or $cont[$i] =~ /AD$/;
        }
    next;
    }
    
    my ($chr,$pos,$gt,$gq,$ad)=(split /\t/)[@cols];
    next if $gq eq 'NA' or $gq eq '.';
    my ($nref,$nalt)=split /,/,$ad;
    next if $nref+$nalt <7 ;
    $call{$chr.'_'.$pos}{$gt}=undef if $gq > 20;
}
close($call_fh);

open my $reffh, ' <', $ref_file or die "Cannot open file:$ref_file \n$!";
open my $wfh, ' >', $sam.'.GT.xls' or die "Cannot write file: \n$!";
while(<$reffh>){
    s/\s+$//;
    if(/^##/){
        my ($key)=$_=~/^##(\w+)##\s*$/;
        $flag=1 if  $type=~/$key/ ;
        print $wfh $_,"\n" if $flag;
    }elsif(/^_+\s*$/){        
        print $wfh $_,"\n" if $flag;
        $flag=0;
    }else{
        next unless $flag;
        my ($chr,$pos, $gt, $strand, $info)=split /\t/,$_,5;
        if(/^#/){        
            print  $wfh  '#'.$info,"\tResult\n";
        }else{
            my $res= exists $call{$chr.'_'.$pos}{$gt} ? '+' : '0';
            print $wfh  $info,"\t", $res, "\n" if $res;
        }
    }
            
}
close($reffh);