#!/usr/bin/perl 

# author xcg
# the produced file may lost some information if:
#	1. the ref have heterozygous allels
#       2. the IS or INS have heterozygous alles


use strict;
use Getopt::Std;
use List::Util qw[min max];
use FindBin;
use lib "$FindBin::Bin/..";
use File::Basename;
#use fastapro;

my @hist;
my $n=0;
my %opts =( s=>undef, r=>undef);
getopts('s:r:', \%opts);
my $ref_fas =$opts{r};

my $samplename =$opts{s} ;
if ( !defined $samplename) {
     $samplename=basename($ARGV[0]);	
     #my $zz=$ARGV[0];
     #$zz =~ /(.*)\.[^.]*/;
     #$samplename=$1;
}
if ( !defined $samplename) {
	$samplename="unknown";
}

die "usage:sdi2vcf.pl <file>" if (@ARGV<1 && -t STDIN);
my  %riupac= ( "R"=>"AG", "Y"=>"CT", "S"=>"GC", "W"=>"AT","K"=>"GT", "M"=>"AC", "A" =>"AA", "T"=>"TT","G"=>"GG", "C"=>"CC", "N"=>"NN" ); 

print qq(##fileformat=VCFv4.1
##source=IMR
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	$samplename);
print "\n";


while (<>) {
    my @F=split;
    if ( length($F[2]) ==1)  {  #for snp
    	my $bb=$F[3]; chomp $bb;
	my ($alt, $gt);
        $F[2]="N" if ($F[2] =~ /[RYSWKMryswkm]/); 
	if ($bb =~ /[ATGCatgc]/ )  {
	   $alt=$bb;
	   $gt="1/1";
	} else {
	   my @temp=split(//, $riupac{$bb});
	   if ($temp[1] eq $temp[0] && $temp[0] eq $F[2] ) {     # for [YM..S]->N this kind of SNP. hard to solve since only AGCTN is allowed for .vcf
		$alt=join(",",@temp);
		$gt="0/0";
	   } elsif ($temp[0] eq $F[2]) {
		$alt=$temp[1];
		$gt="0/1";
	   } elsif ( $temp[1] eq $F[2] ) {
		$alt=$temp[0];
		$gt="0/1";
	   } else { 
		$alt=join(",",@temp);
		$gt="1/2";
	   }
	}
    	print join("\t", $F[0],$F[1],'.',$F[2], $alt,'.','.','.','GT',$gt),"\n";
    } else {
	my ($chid,$pos, $len, $oldI, $newI)=@F;
	next if ($pos==1);
	my $newp=$pos-1;
	$oldI="" if ($oldI eq "-"); 
	$newI="" if ($newI eq "-"); 
	$oldI =~ s/[RYSWKMryswkm]/N/g;
	$newI =~ s/[RYSWKMryswkm]/N/g;
             #if $newI have heteozygous allels, ther could be problem	
	my $prefix;
	if (!$ref_fas) {
	    $prefix='N';
	} else {
	    $prefix=getNucEx($ref_fas, $chid, $pos-1, 1);
	}
	$prefix =~ s/[RYSWKMryswkm]/N/g;
    	print join("\t", $F[0],$newp,'.',$prefix.$oldI, $prefix.$newI,'.','.','.','GT',"1/1"),"\n";
	next;
    }

}
