#!/usr/bin/perl

use strict;
use Thread;
use Bio::Tools::GFF;


my $outFile = "result.sam";
my $inFile = "";
my $GTFfile = "";
my $sortAdd ="";

$#ARGV > 0 || die ("Input parameters are missing\nUse following:\n./combined.pl input_file GTF_file [-s temp_path mem_size] [output_file]\n");

if($#ARGV == 2)
{
    $outFile = $ARGV[2];
}
if($#ARGV == 5 || $#ARGV == 4)
{   
    if($ARGV[2] eq "-s")
    {
	$sortAdd = "-T ".$ARGV[3]." -S ".$ARGV[4];
    }
    if($#ARGV == 5)
    {
	$outFile = $ARGV[5];
    }
}

my $inFile = $ARGV[0];
my $GTFfile = $ARGV[1];


open(REP, $inFile) || die ($inFile." file missing");
close(REP);

 
#-T /scratch/NFS/bi/ -S 32G
#"GRCh37_58.gtf";
############################################
## Reading the gene models from GTF file  ##
############################################
## hash of transcripts with exons info 
my %KGexons;

#          {transcript_id}{exons}[i]{start} - exon start
#          {transcript_id}{exons}[i]{end} - exon end
#          {transcript_id}{exons}[i]{length} - exon length
#          {transcript_id}{count} - no of exons
#          {transcript_id}{length} - trans length
#          {transcript_id}{chr} - chromosome id
#          {transcript_id}{strand} - strand

## temp variables for reading the GFF/GTF file
my $feature;
my $trId="";
my $exon_count=0;
my $trans_len=0;
my $strand="";
my $chr="";
my @temp_arr=();

open(GTFIN, $GTFfile) || die ($GTFfile." file missing");
my $gffio = Bio::Tools::GFF->new(-fh => \*GTFIN, -gff_version => 2.5);


print "Analysing the GTF file\n";
while($feature = $gffio->next_feature()) 
{

    if($feature->primary_tag() eq "exon")
    {
	my @values = $feature->get_tag_values("transcript_id");

	if($trId eq $values[0])
	{
	    $exon_count++;
	    $trans_len+= $feature->length;
	}
	else
	{
	    if($trId ne "")
	    {
		$KGexons{$trId}{count} = $exon_count+1;
		$KGexons{$trId}{length} = $trans_len;
		$KGexons{$trId}{chr} = $chr;
		$KGexons{$trId}{strand} = $strand;
		
		if($strand == 1)
		{
		    $KGexons{$trId}{exons} = [ @temp_arr ];
		}
		else
		{
		    $KGexons{$trId}{exons} = [ reverse(@temp_arr) ];
		}	   
	    }
	    $trId = $values[0];
	    $strand = $feature->strand;
	    $chr = $feature->seq_id;
	    $exon_count = 0;
	    $trans_len = $feature->length;
	    @temp_arr=();	
	}
	$temp_arr[$exon_count]{start} = $feature->start;
	$temp_arr[$exon_count]{end} = $feature->end;
    }
}

$KGexons{$trId}{count} = $exon_count+1;
$KGexons{$trId}{length} = $trans_len;
$KGexons{$trId}{chr} = $chr;
$KGexons{$trId}{strand} = $strand;

if($strand == 1)
{
    $KGexons{$trId}{exons} = [ @temp_arr ];
}
else
{
    $KGexons{$trId}{exons} = [ reverse(@temp_arr) ];
}

$gffio->close();
close(GTFIN);



############################################
# Processing input SAM file              ###
############################################
my %t;

#print "Mapping the alignments\n";
$t{$inFile} = Thread->new(\&par_exe, $inFile);
$t{$inFile}->join();

print "Sorting the alignments\n";
system("sort  -u -k 3,3 -k 4,4n -k 1,1 -k 6,6 -k 8,8 ".$sortAdd." ".$inFile.".gen > ".$outFile) == 0 or die ("system sort failed: $?");
print "Deleting temporary file\n";
system("rm ".$inFile.".gen");

sub par_exe
{
    
    my $in_file = @_[0];
    print "Matching alignments to Chromosome locations and storing in SAM format\n";
    print $in_file."\n";
    
    my $out_file= $in_file.".gen";
    open(REP, $in_file) || die ($in_file."file missing");
    open(OUT, ">".$out_file) || die ($out_file."file missing");
    
    my $read_id="";
    my $old_id="";
    my $strand="";
    my $trans_id=0;
    my $read_start=0;
    my $read_seq="";
    my $read_qual="";
    my $read_score = 0;
    my $read_cigar = 0;
    my $mate_ref = 0;
    my $mate_start = 0;
    my $iszie = 0;
    my $read_mis=0;
    my $strand_v = 0;
    my $align_count = 0;
    my @alignments = ();
    my $flag=0;

    my $readLen = 50;
    
    while(<REP>)
    {
	chomp;
## Ommit the SAM format header lines
	if($_!~/^@/)
	{
	    my @fields = split(/\t/);

## Ommit lines with values 77 and 141 (for paired-end) 4 and 16 (for single-end) reads
	    if($fields[1] != 77 && $fields[1] != 141 && $fields[1] != 4 && $fields[1] != 16)
	    {

		$read_id = $fields[0];
		$flag = $fields[1];
		$trans_id = $fields[2];
		$read_start = $fields[3];
		$read_score = $fields[4];
		$read_cigar = $fields[5];
		$mate_ref = $fields[6];
		$mate_start = $fields[7];
		$iszie = $fields[8];
		$read_seq = $fields[9];
		$read_qual = $fields[10];
		$read_mis = $fields[13];


		$readLen=length($read_seq);

		my $read_pos=0;
		my $mate_pos=0;

## Lines with values 99, 163, 83, and 147 are accepted for pair-end; 0 for single-end
		if($flag == 99 || $flag == 163 || $flag == 83 || $flag == 147 || $flag==0)
		{
## Calculation of transcript length		    
		    my $trans_len=0;
		    for(my $i=0; $i < $KGexons{$trans_id}{count}; $i++)
		    {
			$trans_len = $trans_len + $KGexons{$trans_id}{exons}[$i]{end} - $KGexons{$trans_id}{exons}[$i]{start} + 1;
		    }
		    
## Setting variables if read was aligned to a forward strand of the transcript
		    if($flag == 99 || $flag == 163)
		    {
			if($KGexons{$trans_id}{strand} == 1)
			{
			    $strand_v = $flag;
			    $read_pos = $read_start-1;
			    $mate_pos = $mate_start-1;
			}
			else
			{
			    $strand_v = $flag - 16;
			    $read_pos = $trans_len - $read_start - $readLen +1;
			    $mate_pos = $trans_len - $mate_start - $readLen +1;
			}
		    }

## Setting variables if read was aligned to a revers strand of the transcript (only paired-end)
		    if($flag == 83 || $flag == 147 || $flag == 0)
		    {
			if($KGexons{$trans_id}{strand} == 1)
			{
			    $strand_v = $flag;
			    $read_pos = $read_start-1;
			    if($flag != 0)
			    {
				$mate_pos = $mate_start-1;
			    }
			}
			else
			{
			    $strand_v = $flag + 16;
			    $read_pos = $trans_len - $read_start - $readLen + 1;
			    if($flag != 0)
			    {
				$mate_pos = $trans_len - $mate_start - $readLen + 1;
			    }
			}
		    }
 

		    my $read_exon_start = 0;
		    my $mate_exon_start = 0;
		    my $cigar = "";
		    my $XS = "";
		    
		    my $read_stop=0;
		    my $mate_stop=0;

## Finding on which exon read maps, read genomic position an eventual location on exon-exon junction	    
		    for(my $i=0; $i < $KGexons{$trans_id}{count}; $i++)
		    {
			my $exon_length = $KGexons{$trans_id}{exons}[$i]{end} - $KGexons{$trans_id}{exons}[$i]{start} + 1;
			if($read_pos > 0)
			{
			    $read_pos = $read_pos - $exon_length;
			}

## If paired-end read calculate also the mate read position
			if($flag != 0)
			{
			    if($mate_pos >0)
			    {
				$mate_pos = $mate_pos - $exon_length;
			    }
			}
			
			if($read_pos < 0 && $read_stop == 0)
			{
			    $read_exon_start = $KGexons{$trans_id}{exons}[$i]{end} + $read_pos + 1;
			    if($read_pos > -$readLen)
			    {
				if(($i+1) < $KGexons{$trans_id}{count})
				{
				    my $temp1 = 0 - $read_pos;
				    my $temp2 = $readLen + $read_pos;
				    my $temp3 = $KGexons{$trans_id}{exons}[$i+1]{start} - $KGexons{$trans_id}{exons}[$i]{end} - 1;
				    
				    ## If read is located on a exon-exon junction create an appropriate CIGAR string
				    if($temp3 >0)
				    {
					$cigar = $temp1."M".$temp3."N".$temp2."M";
					if($KGexons{$trans_id}{strand} == 1)
					{
					    $XS="\tXS:A:+";
					}
					else
					{
					    $XS="\tXS:A:-";
					}
				    }
				    else
				    {
					$cigar = $readLen."M";
					$XS = "";
				    }
				}
				else
				{
				    $cigar = "";
				    
				}
			    }
			    else
			    {
				$cigar = $readLen."M";
				$XS = "";
			    }	
			    $read_stop = 1;
			}

			if($flag != 0)
			{
			    if($mate_pos < 0 && $mate_stop==0)
			    {
				$mate_exon_start = $KGexons{$trans_id}{exons}[$i]{end} + $mate_pos + 1;
				$mate_stop = 1;
			    }
			}
			
		    }
		    my $rnext = "=";
		    if($flag == 0)
		    {
			$rnext = "*";
		    }
		    if($cigar ne "")
		    {
			print OUT $read_id."\t".$strand_v."\t".$KGexons{$trans_id}{chr}."\t".$read_exon_start."\t255\t".$cigar."\t".$rnext."\t".$mate_exon_start."\t0\t".$read_seq."\t".$read_qual."\t".$read_mis.$XS."\n";	
		    }
		}
	    }    
	}
    }    
    close(REP);
    close(OUT);   
}
