#!/usr/bin/perl 
use strict;
use File::Basename;
use Getopt::Long;
use Switch;

=head1 DESCRIPTION

	var2exomeassisant.pl  -- used to convert the format of your input file (SAMtools pileup, standard VCF or SOAPsnp) to that requered by Exome-assistent.


=head1 VERSION

	Author: Enjian Shen, baby1885@163.com
	Version: 1.0,  Date: 2011-8-30
	Version: 1.1,  Date: 2011-10-5 add minimal sequence depth filter
	

=head1 USAGE

	PARAMETER:
	-type	<str>	The format of input file need to be converted (pileup/soapsnp/vcf). [default:pileup]
	-file	<str>	SNPs/InDels input file.
	-min	<int>	Minimal sequence depth filter. [default:1]
	-name	<str>	Name of project.
	-outdir	<str>	Output directory of the converted result. [default:./]
	-help		This help message.
	
	EXAMPLE:
	perl var2exomeassisant.pl -type pileup -file SNP_InDel.pileup -min 10 -name Cancer -outdir /home/users/
	
=cut

#SAMtools pileip: chr10   61769   g       R       4       4       37      2       A$.     \^ or 
#SOAPsnp: chr21   9928594 A       M       80      A       35      43      44      C       29      12      12      56      0.504044        1.03571 0       18
#VCF: chr10	3165360	.	C	G	2755.12	.	AB=0.50;AC=1;AF=0.50;AN=2;DP=147;Dels=0.00;HRun=1;HaplotypeScore=5.25;MQ=36.73;MQ0=0;QD=18.74;SB=-1095.22	GT:DP:GL:GQ	0/1:141:-321.19,-42.50,-329.14:99.00
#
my ( $type,$file_in,$min_depth,$pro_name,$outdir,$help );
GetOptions(
	"type=s"  => \$type,
	"file=s"  => \$file_in,
	"min:s"   => \$min_depth,
	"name:s"  => \$pro_name,
	"outdir:s"=> \$outdir,
	"help:s"  => \$help,
);

die `pod2text $0` if ( !$type || !$file_in || $help );
die `pod2text $0` if ( $type!~/^pileup$|^soapsnp$|^vcf$/ );

$outdir        ||= "./";
$type          ||= "pileup";
$min_depth     ||=1;
my $snp_indel_bas=basename $file_in;
$pro_name      ||= "$snp_indel_bas";
my $start_time   = times();#Start time
mkdir $outdir unless( -e $outdir );

switch($type){
	case "pileup" {#SAMtools pileup
		pileup( $file_in,$pro_name,$outdir );
	}
	case "soapsnp" {#SOAPsnp for SNPs
		soapsnp( $file_in,$pro_name,$outdir );
	}
	case "vcf" {#VCFv3.3
		vcf( $file_in,$pro_name,$outdir );
	}
}

sub vcf{
	my ( $vcf_file,$vcf_name,$vcf_out ) = @_;
	my $vcf_total = 0;
	my $vcf_snp   = 0;
	input_inf( "vcf",$vcf_file,$vcf_name,$vcf_out );#output information to screen
	
	open HANDPILEUP, "<$vcf_file" or die "$!";
	open OUTSNP,     ">$vcf_out/$vcf_name.snp_list.txt" or die "$!";
	open LOG,        ">$vcf_out/${vcf_name}_2exomeassistant.log" or die "$!";
	open HANDPILEUP, "<$vcf_file" or die "$!";#vcf SNP format for one sample(VCFv3.3)
	while(<HANDPILEUP>){
		chomp;
		next if(/#/);
		$vcf_total++;
		if( $_!~/AF=/ ){
			warn "Unable to distinguish VCF(VCFv3.3) format in line $vcf_total, please check it.";
			next;
		}
		my @vcf_tmp = split/\t/,$_;
		my ( $chr,$pos,$ref,$alt ) = ( $vcf_tmp[0],$vcf_tmp[1],$vcf_tmp[3],$vcf_tmp[4] );
		next if( $ref eq $alt );
		next if( $chr=~/_/ );
		$chr=~s/chr//;
		$_  =~/DP=(\d+);/;
		my $depth = $1;
		next if($depth<$min_depth);
		$vcf_snp++;
		if( /AF=(\d+\.\d+);/ ){
			my $af = $1;
			if( $af eq "1.00" ){#hom
				print OUTSNP "$chr,$pos,$pos,1,$alt/$alt,$depth\n";
			}else{#het
				$_=~/AB=(\S+);/;
				my $fre     = $1;
				my $ref_num = int( $fre*$depth );
				my $alt_num = $depth-$ref_num;
				print OUTSNP "$chr,$pos,$pos,1,$alt/$ref,$alt_num/$ref_num\n";
			}
		}
		
	}
	my $end_time   = times();
	my $speed_time = sprintf "This program took %.2f CPU seconds of user time\n", $end_time - $start_time;
	warn "There are too many SNPs in sample, we'd better suggest you to filter low-quailty or low-coverage SNPs\n" if( $vcf_snp>200000 );
	print LOG "$speed_time";
	print LOG "File:\t$vcf_file\n";
	print LOG "Total VCF of SNPs line:\t$vcf_total\n";
	print LOG "Effective SNPs line:\t$vcf_snp\n\n";
	print LOG "There are too many SNPs in sample, we'd better suggest you to filter low-quailty or low-coverage SNPs\n" if( $vcf_snp>200000 );
	close OUTSNP;
	close LOG;
	
}



sub soapsnp{
	my ( $soap_file,$soap_name,$soap_out ) = @_;
	my ( $alle,$alle_support,$seq_depth );
	my $soap_total = 0;
	my $soap_snp   = 0;
	my $soap_indel = 0;
	input_inf( "soapsnp",$soap_file,$soap_name,$soap_out );#output information to screen 
	open HANDPILEUP, "<$soap_file" or die "$!";
	open OUTSNP,     ">$soap_out/$soap_name.snp_list.txt" or die "$!";
	open LOG,        ">$soap_out/${soap_name}_2exomeassistant.log" or die "$!";
	while(<HANDPILEUP>){
		chomp;
		my @soap_line = split/\t/,$_;
		$soap_total++;
		if(@soap_line != 18){
			warn "Error SOAPsnp format in line $soap_total, please check it.";
			next;
		}
		
#chr21   9928594 A       M       80      A       35      43      44      C       29      12      12      56      0.504044        1.03571 0       18
		my ( $chr,$pos,$ref,$wildtype,$soapbase1,$soapbase1_depth,$soapbase2,$soapbase2_depth ) = ( $soap_line[0],$soap_line[1],$soap_line[2],$soap_line[3],$soap_line[5],$soap_line[8],$soap_line[9],$soap_line[12] );
		next if( $ref eq $wildtype );
		next if( $chr=~/_/ );
		$chr=~s/chr//;
		if( $wildtype=~/[ATCGatcg]/ ){#hom
			$alle		  = "$soapbase1/$soapbase1";
			$alle_support = $soapbase1_depth;
			$seq_depth    = $soapbase1_depth;
		}else{
			$alle		  = "$soapbase1/$soapbase2";
			$alle_support = "$soapbase1_depth/$soapbase2_depth";
			$seq_depth    =$soapbase1_depth+$soapbase2_depth;
		}
		next if($seq_depth<$min_depth);
		$soap_snp++;
		print OUTSNP "$chr,$pos,$pos,1,$alle,$alle_support\n";
	}
	my $end_time   = times();
	my $speed_time = sprintf "This program took %.2f CPU seconds of user time\n", $end_time - $start_time;
	warn "There are too many SNPs in sample, we'd better suggest you to filter low-quailty or low-coverage SNPs\n" if( $soap_snp>200000 );
	print LOG "$speed_time";
	print LOG "File:\t$soap_file\n";
	print LOG "Total SOAPsnp of SNPs line:\t$soap_total\n";
	print LOG "Effective SNPs line:\t$soap_snp\n\n";
	print LOG "There are too many SNPs in sample, we'd better suggest you to filter low-quailty or low-coverage SNPs\n" if( $soap_snp>200000 );
	close OUTSNP;
	close LOG;
	
}


sub pileup{
	my ( $pile_file,$pile_name,$pile_out ) = @_;
	my ( $alle,$alle_support,$seq_depth );
	my $pile_total  = 0;
	my $pile_effect = 0;
	my $pile_snp    = 0;
	my $pile_indel  = 0;
	input_inf( "pileup",$pile_file,$pile_name,$pile_out );#output information to screen 
	open HANDPILEUP, "<$pile_file" or die "$!";
	open OUTSNP,     ">$pile_out/$pile_name.snp_list.txt" or die "$!";
	open OUTINDEL,   ">$pile_out/$pile_name.indel_list.txt" or die "$!";
	open LOG,        ">$pile_out/${pile_name}_2exomeassistant.log" or die "$!";
	
	while(<HANDPILEUP>){
		chomp;
		$pile_total++;
		next if ( /^\s*$/ );
		my @pileup_line = split/\s+/,$_;
#chr10   61769   g       R       4       4       37      2       A$.     \^ or 
#chr12 61334 G 23  ,.....,,.,.,...,,,.,..A <<<;<<<<<<<<<3<=<<<;<<+
		if ( ( scalar(@pileup_line) == 10 ) || (scalar(@pileup_line) == 6 ) ) {#for samtools pileup
			die( "Format error of samtools pileup, please check it at line $pile_total." ) if( $pileup_line[0]!~/chr/ || $pileup_line[1]!~/\d+/ ); 
			my ( $chr,$pos ) = ( $pileup_line[0],$pileup_line[1] );
			next if( $chr=~/_/ );
			$chr=~s/chr//;
			my $ref   = uc( $pileup_line[2] );
			my $depth = $pileup_line[scalar(@pileup_line)-3];
			my $seq   = uc($pileup_line[scalar(@pileup_line)-2]);
			next if( ( $seq!~/[ATCGatcg]/ ) && ( $seq!~/(\+|\-)[0-9]+[ATCGN]+/ ) );
			die( "Format error of samtools pileup, please check it at line $pile_total." ) if( $seq!~/[ATCGatcg,\.]/ );
			$seq =~ s/^\@//g;#remove sequence with '@/$/^' that does not represent a base
			$seq =~ s/\$//g;
			$seq =~ s/\^.//g;
			my %base_nub;
			my %indel_base;
			my %indel_strand;
			my @indel_sort;
			my @base_sort;
			if( $seq!~/(\+|\-)[0-9]+[ATCGN]+/ ){#for SNP line
				next if($seq!~/[ATCGatcg]/);
				foreach my $eva_base( "A","T","C","G","N" ){
					if( $eva_base eq $ref ){
						my $ref_num = ($seq=~s/(,|\.)/$1/g);
						next if( $ref_num eq '' );
						$base_nub{$eva_base} = $ref_num;
					}else{
						my $other_num = ($seq=~s/($eva_base)/$1/g);
						next if( $other_num eq '' );
						$base_nub{$eva_base} = $other_num;
					}
				}
				
				foreach my $eva_base2( sort{$base_nub{$b}<=>$base_nub{$a}} keys %base_nub ){
					push @base_sort,$eva_base2;
				}
				next if( @base_sort == 0 );
				if(@base_sort>=2){
					$alle         = "$base_sort[0]/$base_sort[1]";
					$alle_support = "$base_nub{$base_sort[0]}/$base_nub{$base_sort[1]}";
					$seq_depth    = $base_nub{$base_sort[0]}+$base_nub{$base_sort[1]};
				}else{
					$alle         = "$base_sort[0]/$base_sort[0]";
					$alle_support = "$base_nub{$base_sort[0]}";
					$seq_depth    = $base_nub{$base_sort[0]};
				}
				next if($seq_depth < $min_depth);
				$pile_snp++;
				print OUTSNP "$chr,$pos,$pos,1,$alle,$alle_support\n";
			}else{#for InDel line
				my @indel_str = ($seq=~/([\+-][0-9]+[ATCGN]+)/g);
				foreach my $eva_str(@indel_str){
				
					$eva_str=~/([\+|-])([0-9]+)([ATCGN]+)/;
					my ( $ins_del,$indel_base_num,$indel_bases ) = ( $1,$2,$3 );
					if( $ins_del eq '+' ){
						$ins_del = "\\+";
					}
					my $indel_bases_length_s = length $indel_bases;
					my $indel_support        = ($seq=~s/($ins_del$indel_bases_length_s$indel_bases)/$1/g);
					$ins_del=~s/\\//g;
					$indel_base{"$ins_del$indel_bases"}   = $indel_support;
					$indel_strand{"$ins_del$indel_bases"} = $ins_del;
				}
				my $seq_tmp = $seq;
				$seq_tmp=~s/([,\.])([\+|-])([0-9]+)([ATCGN]+)//g;
				my $indel_nor_num = ($seq_tmp=~s/(,|\.)/$1/g);
				
				foreach my $eva_indel( sort{$indel_base{$b}<=>$indel_base{$a}} keys %indel_base ){
					push @indel_sort,$eva_indel;
				}
				next if( @indel_sort == 0 );
				my $alle_all = $indel_sort[0];
				$alle_all=~/[\+|-]([ATCGN]+)/;
				$alle = $1;
				my ( $strand,$indel_abn_num ) = ( $indel_strand{$alle_all},$indel_base{$alle_all} );
				
				if( $strand eq '+' ){
					my $pos_next = $pos+1;
					if( $indel_nor_num != 0 ){
						$seq_depth=$indel_abn_num+$indel_nor_num;
						next if($seq_depth<$min_depth);
						$pile_indel++;
						print OUTINDEL "$chr,$pos,$pos_next,1,$alle,$indel_abn_num/$indel_nor_num\n";
					}else{
						$seq_depth=$indel_abn_num;
						next if($seq_depth<$min_depth);
						$pile_indel++;
						print OUTINDEL "$chr,$pos,$pos_next,1,$alle,$indel_abn_num\n";
					}
				}else{
					my $pos_start        = $pos;
					my $indel_seq_length = length $alle;
					my $pos_next         = $pos+$indel_seq_length+1;
					if($indel_nor_num != 0){
						$seq_depth=$indel_abn_num+$indel_nor_num;
						next if($seq_depth<$min_depth);
						$pile_indel++;
						print OUTINDEL "$chr,$pos_start,$pos_next,1,-/,$indel_abn_num/$indel_nor_num\n";
					}else{
						$seq_depth=$indel_abn_num;
						next if($seq_depth<$min_depth);
						$pile_indel++;
						print OUTINDEL "$chr,$pos_start,$pos_next,1,-/,$indel_abn_num\n";
					}
					
				}

			}
        }
	}
	
	my $end_time   = times();
	my $speed_time = sprintf "This program took %.2f CPU seconds of user time\n", $end_time - $start_time;
	warn "There are too many SNPs in sample, we'd better suggest you to filter low-quailty or low-coverage SNPs\n" if( $pile_snp>200000 );
	print LOG "$speed_time";
	print LOG "File:\t$pile_file\n";
	print LOG "Total pileup of SNPs/InDel line:\t$pile_total\n";
	print LOG "Effective SNPs line:\t$pile_snp\n";
	print LOG "Effective InDel line:\t$pile_indel\n\n";
	print LOG "There are too many SNPs in sample, we'd better suggest you to filter low-quailty or low-coverage SNPs\n" if( $pile_snp>200000 );
	close OUTINDEL;
	close OUTSNP;
	close LOG;
}

sub input_inf{
	my ( $inp_type,$inp_file,$inp_name,$inp_outdir )=@_;
	print "Input information
File type: $inp_type
File path: $inp_file
Project name: $inp_name
Output document: $inp_outdir
";
}
