#!/usr/bin/perl

use strict;
use Getopt::Std;


my $outFile = "result.sam";
my $version = 2.0;

my $inFile = "";
my $GTFfile = "";
my $sortAdd ="";
my $strandness=0;
my $pairOrientation = 0;

my $single_end =1;

my $ex_text = "Input parameters are missing\nUse following:\nperl MapAl.pl -i input_file -g GTF_file [-t temp_path] [-s mem_size] [-S strandness] [-P pair_orientation] [-o output_file]\n";


#### -S 
# 0 - forward
# 1 - both
# 2 - reverse
#### -P
# 0 - opposite
# 1 - both
# 2 - same

$#ARGV > 0 || die ($ex_text);

our($opt_i, $opt_g, $opt_s, $opt_t, $opt_S, $opt_P, $opt_o);

getopt('igstSPo');

if(defined $opt_i && $opt_i ne '1')
{
    $inFile = $opt_i;
}
else
{
    die ($ex_text);
}

if(defined $opt_g && $opt_g ne '1')
{
    $GTFfile = $opt_g;
}
else
{
    die ($ex_text);
}

if(defined $opt_s && $opt_s ne '1')
{
    $sortAdd = $sortAdd." -S ".$opt_s;
}

if(defined $opt_t && $opt_t ne '1')
{
    $sortAdd = $sortAdd." -T ".$opt_t;
}

if(defined $opt_S)
{
    $strandness = $opt_S;
}

if(defined $opt_P)
{
    $pairOrientation = $opt_P;
}

if(defined $opt_o && $opt_o ne '1')
{
    $outFile = $opt_o;
}

#print $strandness."\n"; 

#-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=();

sub exon_sort{
    $a->{start} <=> $b->{start}
}

my @exon_sorted=();

open(GTFIN, $GTFfile) || die ($GTFfile." file missing");

##open(GTFIN, "sort -k3,3 $GTFgile|")

print "Analysing the GTF file\n";
 
while( <GTFIN>) 
{
    chomp;
    my @features = split('\t', $_);

 
    if($features[2] eq "exon")
    {
	$features[8] =~/transcript_id\s\"(.+?)\"/;
	#print $1;
	my $values = $1;

	if($trId eq $values)
	{
	    $exon_count++;
	    #$trans_len+= $feature->length;
	    $trans_len+= $features[4] - $features[3]+1;
	}
	else
	{
	    if($trId ne "")
	    {
		$KGexons{$trId}{count} = $exon_count+1;
		$KGexons{$trId}{length} = $trans_len;
		$KGexons{$trId}{chr} = $chr;
		$KGexons{$trId}{strand} = $strand;

		@exon_sorted = sort exon_sort @temp_arr;


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


		
		
#		if($trId eq "ENST00000353801")
#		{
#		    for(my $jj=0; $jj < $exon_count+1; $jj++)
#		    {
#			print $KGexons{$trId}{exons}[$jj]{start}." ".$KGexons{$trId}{exons}[$jj]{end}." ".($KGexons{$trId}{exons}[$jj]{end}-$KGexons{$trId}{exons}[$jj]{sta#rt}+1)."\n"; 
#		    }
#		}
		
	   
	    }
	    $trId = $values;
	    if($features[6] eq '+')
	    {
		$strand = 1;
	    }
	    else
	    {
		$strand = 0;
	    }
	    $chr = $features[0];
	    $exon_count = 0;
	    $trans_len = $features[4] - $features[3]+1;;
	    @temp_arr=();	
	}
	$temp_arr[$exon_count]{start} = $features[3];
	$temp_arr[$exon_count]{end} = $features[4];
    }

}


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




@exon_sorted = sort exon_sort @temp_arr;

#if($strand == 1)
#{
    $KGexons{$trId}{exons} = [ @exon_sorted ];
#}
#else
#{
#    $KGexons{$trId}{exons} = [ reverse(@exon_sorted) ];
#}
#for(my $jj=0; $jj < $exon_count+1; $jj++)
#{
#    print $exon_sorted[$jj]{start}." ".$exon_sorted[$jj]{end}."\n"; 
#}


#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";
par_exe($inFile, $outFile);

#$t{$inFile} = Thread->new(\&par_exe, $inFile);
#$t{$inFile}->join();

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

sub alignment_sort{
    $a->{chr} cmp $b->{chr}
    ||
	$a->{pos} <=> $b->{pos}
    ||	
	$a->{pos_next} <=> $b->{pos_next}
    ||	
	$a->{nm} <=> $b->{nm}
    ||
	$a->{cigar} cmp $b->{cigar}

}

sub par_exe
{

    my $general_count=0;

    my $out_str = "";
    
    my $in_file = @_[0];
    my $out_file = @_[1].".gen";
    print "Matching alignments to Chromosome locations and storing in SAM format\n";
    print "Analysing file: ".$in_file."\n";
    open(REP, $in_file) || die ($in_file."file missing");

   # my $out_file= $in_file.".gen";
    open(OUT, ">".$out_file) || die ($out_file."file missing");
    
    
    my $read_id="";
    my $old_read_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 $cigar = "";

    my $readLen = 50;
    my $reffLen = 50;
    
    undef @alignments;
    @alignments=();
    

    while(<REP>)
    {
	my @tmpCIGAR=();
	my @CIGAR=();
	
	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)
	    if($fields[1] == 0 || $fields[1] == 16 || ($fields[1] & 2))
	    {
		
		if($fields[1] & 1)
		{
		    my $single_end = 0;
		}
		else
		{
		    my $single_end = 1
		}

		$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];
		my $tmpread_seq = $fields[9];
		my $tmpread_qual = $fields[10];
		my $mate_len = $read_start+$iszie-$mate_start;
		

		$_ =~ /NM:i:(\d+)/;
		$read_mis = $1;


		$readLen=length($tmpread_seq);

		$reffLen = 0;
		my $tmpLen = 0;
		while($read_cigar=~/(\d+)(\w)/g)
		{
		    my $rec={};
		    $rec->{type} = $2;
		    $rec->{length}=$1;
		    push @tmpCIGAR, $rec;
		    
		    if($2 eq "I" or $2 eq "M" or $2 eq "S" or $2 eq "=" or $2 eq "X")
		    {
			$tmpLen+=$1;
		    }
		    if($2 eq "D" or $2 eq "M")
		    {
			$reffLen+=$1;
		    }
		}
		
		if($tmpLen != $readLen)
		{
		    print "read length error\n";
		    #print $trans_id."\t";
		    #print $read_cigar."\t";
		    #print $tmpLen."\t".$readLen."\n";
		    next;		    
		}
		if(not $KGexons{$trans_id})
		{
		    #print "missing transcript\n";
		    next;
		}
		

		#my @CIGAR = reverse(@CIGAR);

		my $read_pos=0;
		my $mate_pos=0;
		
		my $let_in = 0;
		#my $lib_type = 0;

		my $temp_val = $flag & 240;

		if($single_end && (($flag==0 && $strandness < 2) || ($flag==16 && $strandness > 0)))
		{
		    $let_in = 1;
		    #$lib_type = $flag;
		}
		else
		{		   
		    if($pairOrientation <2 && (($strandness < 2 && ($temp_val==96 || $temp_val==144)) || ($strandness > 0 && ($temp_val==80 || $temp_val==160))))
		    {
			$let_in = 1;		    
		    }
		    
		    if($pairOrientation >0 && (($strandness < 2 && ($temp_val==64 || $temp_val==128)) || ($strandness > 0 && ($temp_val==112 || $temp_val==176))))
		    {
			$let_in = 1;
		    }		
		}

## 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 && $strandness < 2) || ($flag==16 && $strandness > 0) )
		if($let_in)
		{

		    
## 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;
		    }
		    

		    if($temp_val == 96 || $temp_val == 160 ||  $flag == 16) ## Setting variables if read was aligned to a forward strand of the transcript (only paired-end)
		    {
			if($KGexons{$trans_id}{strand} == 1)
			{
			    $strand_v = $flag;
			    $read_pos = $read_start-1;
			    if($flag != 16)
			    {
				$mate_pos = $mate_start-1;
			    }
			    @CIGAR = @tmpCIGAR;
			    $read_seq = $tmpread_seq;
			    $read_qual = $tmpread_qual;
			}
			else
			{
			    $strand_v = $flag - 16;
			    $read_pos = $trans_len - $read_start - $reffLen +1;
			    #$mate_pos = $trans_len - $mate_start - $mate_len +1;
			    if($flag != 16)
			    {
				$mate_pos = $trans_len - $mate_start - $mate_len + 1;
			    }
			    @CIGAR = reverse(@tmpCIGAR);
			    $tmpread_seq =~ tr/ATCGatcg/TAGCtagc/;
			    $read_seq = join( "", reverse split( "", $tmpread_seq));
			    $read_qual = join( "", reverse split( "", $tmpread_qual));
			}
		    }
		    elsif($temp_val == 80 || $temp_val == 144 || $flag == 0) ## Setting variables if read was aligned to a rewerse strand of the transcript
		    {						

			if($KGexons{$trans_id}{strand} == 1)
			{
			    $strand_v = $flag;
			    $read_pos = $read_start-1;
			    if($flag != 0)
			    {
				$mate_pos = $mate_start-1;
			    }
			    $read_seq = $tmpread_seq;
			    $read_qual = $tmpread_qual;
			    @CIGAR = @tmpCIGAR;
			}
			else
			{
			    $strand_v = $flag + 16;
			    $read_pos = $trans_len - $read_start - $reffLen + 1;
			    if($flag != 0)
			    {
				$mate_pos = $trans_len - $mate_start + $mate_len + 1;
			    }
			    $tmpread_seq =~ tr/ATCGatcg/TAGCtagc/;
			    $read_seq = join( "", reverse split( "", $tmpread_seq));
			    $read_qual = join( "", reverse split( "", $tmpread_qual));
			    @CIGAR = reverse(@tmpCIGAR);
			}
		    }
		    elsif($temp_val == 64 || $temp_val == 128 || $temp_val == 112 || $temp_val == 176)
		    {						
			
			if($KGexons{$trans_id}{strand} == 1)
			{
			    $strand_v = $flag;
			    $read_pos = $read_start-1;
			    $mate_pos = $mate_start-1;			    
			    $read_seq = $tmpread_seq;
			    $read_qual = $tmpread_qual;
			    @CIGAR = @tmpCIGAR;
			}
			else
			{
			    if($temp_val == 64 || $temp_val == 128)
			    {
				$strand_v = $flag + 48;
			    }
			    else
			    {
				$strand_v = $flag - 48;
			    }
			    $read_pos = $trans_len - $read_start - $reffLen + 1;
			    $mate_pos = $trans_len - $mate_start + $mate_len + 1;			    
			    $tmpread_seq =~ tr/ATCGatcg/TAGCtagc/;
			    $read_seq = join( "", reverse split( "", $tmpread_seq));
			    $read_qual = join( "", reverse split( "", $tmpread_qual));
			    @CIGAR = reverse(@tmpCIGAR);
			}
		    }

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

		    my $oryg_read_pos = $read_pos;

## 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++)
		    {
			#print $KGexons{$trans_id}{exons}[$i]{start}." ".$KGexons{$trans_id}{exons}[$i]{end}."\n";
			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 && $flag != 16)
			{
			    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 > -$reffLen)
			    {

				if(($i+1) < $KGexons{$trans_id}{count})
				{

				    my $temp1 = 0 - $read_pos;
				    my $temp2 = $reffLen + $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($temp2 >0)
				    {
					my $next_exon_len=0;
					my $j=0;
					if($i+$j+1 < $KGexons{$trans_id}{count})
					{
					    $j++;
					    $next_exon_len = $KGexons{$trans_id}{exons}[$i+$j]{end} - $KGexons{$trans_id}{exons}[$i+$j]{start} + 1;
					}

					$cigar = "";
					my $inserted=0;
					for my $href (@CIGAR)
					{
					    if($temp1 >= $href->{length} and $inserted == 0)
					    {
						$temp1=$temp1 - $href->{length};
						$cigar = $cigar.$href->{length}.$href->{type};
					    }
					    elsif($temp1 == 0 and $inserted == 0)
					    {
						$cigar = $cigar.$temp3."N".$href->{length}.$href->{type};
						#$inserted=1;
						####
						$temp1 = $next_exon_len;
						if($i+$j+1 < $KGexons{$trans_id}{count})
						{
						    $j++;
						    $next_exon_len = $KGexons{$trans_id}{exons}[$i+$j]{end} - $KGexons{$trans_id}{exons}[$i+$j]{start} + 1;
						}
						else
						{
						    $inserted=1;
						}
						
					    }
					    elsif($inserted == 0)
					    {
						my $rest = $href->{length} - $temp1;
						$cigar = $cigar.$temp1.$href->{type}.$temp3."N";

						while($rest > $next_exon_len)
						{
						    $temp3 = $KGexons{$trans_id}{exons}[$i+$j+1]{start} - $KGexons{$trans_id}{exons}[$i+$j]{end} - 1;
						 #   print "temp3 ".$temp3."\n";
						    $cigar.=$next_exon_len.$href->{type}.$temp3."N";
						    $rest = $rest - $next_exon_len;

						    if($i+$j+1 < $KGexons{$trans_id}{count})
						    {
							$j++;
							$next_exon_len = $KGexons{$trans_id}{exons}[$i+$j]{end} - $KGexons{$trans_id}{exons}[$i+$j]{start} + 1;
						    }
						    else
						    {
							$inserted=1;
						    }
						}
						$cigar.=$rest.$href->{type};
						$temp1 = $next_exon_len - $rest;

						#print "rest ".$rest."\n";


						if($i+$j+1 < $KGexons{$trans_id}{count})
						{
						    $j++;
						    $next_exon_len = $KGexons{$trans_id}{exons}[$i+$j]{end} - $KGexons{$trans_id}{exons}[$i+$j]{start} + 1;
						}
						else
						{
						    $inserted=1;
						}
					    }
					    else
					    {
						$cigar = $cigar.$href->{length}.$href->{type};
					    }
					    
					}
					
#					$cigar = $temp1."M".$temp3."N".$temp2."M";
					if($KGexons{$trans_id}{strand} == 1)
					{
					    $XS="\tXS:A:+";
					}
					else
					{
					    $XS="\tXS:A:-";
					}
				    }
				    else
				    {
					
					$cigar = "";
					for my $href (@CIGAR)
					{
					    $cigar = $cigar.$href->{length}.$href->{type};
					}
##	$cigar = $readLen."M";
					$XS = "";
				    }
				}
				else
				{
				    $cigar = "*";
				  #print $cigar."\n";  
				}
			    }
			    else
			    {
				
				$cigar = "";
				for my $href (@CIGAR)
				{
				    $cigar = $cigar.$href->{length}.$href->{type};
				}				
				$XS = "";
			    }	
			    $read_stop = 1;
			}
			else	
			{
			    
			    
			}

			if($flag != 0 && $flag != 16)
			{
			    if($mate_pos < 0 && $mate_stop==0)
			    {
				$mate_exon_start = $KGexons{$trans_id}{exons}[$i]{end} + $mate_pos + 1;
				$mate_stop = 1;
			    }
			}					
		    } ### end for loop for exons


		    my $rnext = "=";
		    if($flag == 0 || $flag == 16)
		    {
			$rnext = "*";
		    }

		    
		    if(($flag & 1) && ($mate_exon_start == 0))
		    {
			$cigar="*";
			$general_count++;			
		    }

		    if($cigar ne "*" && $cigar ne "")
		    {
			#$general_count++;
			if($read_id ne $old_read_id)
			{

			    
			    if($old_read_id ne "")
			    {
				my @sorted = ();
				@sorted = sort alignment_sort @alignments;
				my $sort_len = scalar(@sorted);
				my @sorted2_1=();
				my @sorted2_2=();

				for(my $i=0; $i < $sort_len; $i++)			
				{
				    if($i == 0  || ($i > 0 && (($sorted[$i]->{chr} ne $sorted[$i-1]->{chr}) || ($sorted[$i]->{pos} != $sorted[$i-1]->{pos}) || ($sorted[$i]->{pos_next} != $sorted[$i-1]->{pos_next})))) #|| (($sorted[$i]->{cigar} ne $sorted[$i-1]->{cigar}) && ($sorted[$i]->{nm} == $sorted[$i-1]->{nm})))))
				    {
					if($sorted[$i]->{flag} & 64)
					{
					    push @sorted2_1, $sorted[$i];
					}
					elsif($sorted[$i]->{flag} & 128)
					{
					    push @sorted2_2, $sorted[$i];
					}
					elsif($sorted[$i]->{flag} == 0 || $sorted[$i]->{flag} == 16)
					{
					    push @sorted2_1, $sorted[$i];
					}
				    }			 
				}
				my $sort_len2_1 = scalar(@sorted2_1);
				my $sort_len2_2 = scalar(@sorted2_2);
				
				if($sort_len2_1 != $sort_len2_2 && $single_end == 0)
				{
				    print "error\t".$general_count."\t";
				    print $old_read_id." ".$sort_len." ".$sort_len2_1." ".$sort_len2_2."\t";
				    
				    for(my $i=0; $i <= $sort_len; $i++)			
				    {
					print $sorted[$i]->{chr}."\t";
				    }
				    print "\n";
				}
				
				for(my $i=0; $i < $sort_len2_1; $i++)			
				{
				    my $end_1 = "";
				    
				    if($sort_len2_1 >1)		    
				    {
					$end_1.="\tNH:i:".($sort_len2_1);
					if($i>0)
					{
					    $sorted2_1[$i]->{flag}+=256;
					}
				    }

				    if($i < ($sort_len2_1-1))
				    {
					my $cc_1="=";
					if($sorted2_1[$i]->{chr} ne $sorted2_1[$i+1]->{chr})				    
					{
					    $cc_1 = $sorted2_1[$i+1]->{chr};
					}
					$end_1.= "\tCC:Z:".$cc_1."\tCP:i:". $sorted2_1[$i+1]->{pos};
				    }

				    if($sort_len2_1 >1)		    
				    {
					$end_1.="\tHI:i:".$i;
				    }

				    if($version == 1)
				    {
					$end_1="";
				    }
				    
				    print OUT $sorted2_1[$i]->{id}."\t".$sorted2_1[$i]->{flag}."\t".$sorted2_1[$i]->{chr}."\t".$sorted2_1[$i]->{pos}."\t".$sorted2_1[$i]->{score}."\t".$sorted2_1[$i]->{cigar}."\t".$sorted2_1[$i]->{rnext}."\t".$sorted2_1[$i]->{pos_next}."\t".$sorted2_1[$i]->{temp_len}."\t".$sorted2_1[$i]->{seq}."\t".$sorted2_1[$i]->{qual}."\tNM:i:".$sorted2_1[$i]->{nm}.$sorted2_1[$i]->{xs}.$end_1."\n";
				}
				
				for(my $i=0; $i < $sort_len2_2; $i++)			
				{
				    my $end_2 = "";

				    if($sort_len2_2 >1)		    
				    {
					$end_2.="\tNH:i:".($sort_len2_2);
					if($i>0)
					{
					    $sorted2_2[$i]->{flag}+=256;
					}
				    }

				    if($i < ($sort_len2_2-1))
				    {
					my $cc_2="=";
					if($sorted2_2[$i]->{chr} ne $sorted2_2[$i+1]->{chr})				    
					{
					    $cc_2 = $sorted2_2[$i+1]->{chr};
					}
					$end_2.= "\tCC:Z:".$cc_2."\tCP:i:". $sorted2_2[$i+1]->{pos};
				    }

				    if($sort_len2_2 >1)		    
				    {
					$end_2.="\tHI:i:".$i;
				    }

				    if($version == 1)
				    {
					$end_2="";
				    }

				    print OUT $sorted2_2[$i]->{id}."\t".$sorted2_2[$i]->{flag}."\t".$sorted2_2[$i]->{chr}."\t".$sorted2_2[$i]->{pos}."\t".$sorted2_2[$i]->{score}."\t".$sorted2_2[$i]->{cigar}."\t".$sorted2_2[$i]->{rnext}."\t".$sorted2_2[$i]->{pos_next}."\t".$sorted2_2[$i]->{temp_len}."\t".$sorted2_2[$i]->{seq}."\t".$sorted2_2[$i]->{qual}."\tNM:i:".$sorted2_2[$i]->{nm}.$sorted2_2[$i]->{xs}.$end_2."\n";

				    #$general_count++;
				}
			    }			    
			    undef @alignments;
			    @alignments=();

			    #$general_count=0;
			}
			else
			{
			    #print $read_id."\n";
			}


			my $rec= {};
			
			$rec->{id} = $read_id;
			$rec->{flag} = $strand_v;
			$rec->{chr} = $KGexons{$trans_id}{chr};
			$rec->{pos} = $read_exon_start;
			$rec->{score} = $read_score;
			$rec->{cigar} = $cigar;
			$rec->{rnext} = $rnext;
			$rec->{pos_next} = $mate_exon_start;
			$rec->{temp_len} = '0';
#			$rec->{temp_len} = $iszie;
			$rec->{seq} = $read_seq;
			$rec->{qual} = $read_qual;
			$rec->{nm} = $read_mis;
			$rec->{xs} = $XS;
			
			push @alignments, $rec;
			$old_read_id = $read_id;
			
		    }
		    else
		    {
			#$general_count++;
			#print $read_id."\t".$oryg_read_pos."\n";
		    }
		}
		
	    }
	    else
	    {	
	    }
	}
    }    



    my @sorted = ();
    @sorted = sort alignment_sort @alignments;
    my $sort_len = scalar(@sorted);
    my @sorted2_1=();
    my @sorted2_2=();
    
    for(my $i=0; $i < $sort_len; $i++)			
    {
	if($i == 0  || ($i > 0 && (($sorted[$i]->{chr} ne $sorted[$i-1]->{chr}) || ($sorted[$i]->{pos} != $sorted[$i-1]->{pos}) || ($sorted[$i]->{pos_next} != $sorted[$i-1]->{pos_next})))) #|| (($sorted[$i]->{cigar} ne $sorted[$i-1]->{cigar}) && ($sorted[$i]->{nm} == $sorted[$i-1]->{nm})))))
	{
	    if($sorted[$i]->{flag} & 64)
	    {
		push @sorted2_1, $sorted[$i];
	    }
	    elsif($sorted[$i]->{flag} & 128)
	    {
		push @sorted2_2, $sorted[$i];
					}
	    elsif($sorted[$i]->{flag} == 0 || $sorted[$i]->{flag} == 16)
	    {
		push @sorted2_1, $sorted[$i];
	    }
	}			 
    }
				my $sort_len2_1 = scalar(@sorted2_1);
    my $sort_len2_2 = scalar(@sorted2_2);
    
    if($sort_len2_1 != $sort_len2_2 && $single_end == 0)
    {
	print "error\t".$general_count."\t";
	print $old_read_id." ".$sort_len." ".$sort_len2_1." ".$sort_len2_2."\t";
	
	for(my $i=0; $i <= $sort_len; $i++)			
	{
	    print $sorted[$i]->{chr}."\t";
	}
	print "\n";
    }
    
    for(my $i=0; $i < $sort_len2_1; $i++)			
    {
	my $end_1 = "";
	
	if($sort_len2_1 >1)		    
	{
	    $end_1.="\tNH:i:".($sort_len2_1);
	    if($i>0)
	    {
		$sorted2_1[$i]->{flag}+=256;
	    }
	}
	
	if($i < ($sort_len2_1-1))
	{
	    my $cc_1="=";
	    if($sorted2_1[$i]->{chr} ne $sorted2_1[$i+1]->{chr})				    
	    {
		$cc_1 = $sorted2_1[$i+1]->{chr};
	    }
	    $end_1.= "\tCC:Z:".$cc_1."\tCP:i:". $sorted2_1[$i+1]->{pos};
	}
	
	if($sort_len2_1 >1)		    
	{
	    $end_1.="\tHI:i:".$i;
	}
	
	if($version == 1)
	{
	    $end_1="";
	}
	
	print OUT $sorted2_1[$i]->{id}."\t".$sorted2_1[$i]->{flag}."\t".$sorted2_1[$i]->{chr}."\t".$sorted2_1[$i]->{pos}."\t".$sorted2_1[$i]->{score}."\t".$sorted2_1[$i]->{cigar}."\t".$sorted2_1[$i]->{rnext}."\t".$sorted2_1[$i]->{pos_next}."\t".$sorted2_1[$i]->{temp_len}."\t".$sorted2_1[$i]->{seq}."\t".$sorted2_1[$i]->{qual}."\tNM:i:".$sorted2_1[$i]->{nm}.$sorted2_1[$i]->{xs}.$end_1."\n";
    }
    
    for(my $i=0; $i < $sort_len2_2; $i++)			
    {
	my $end_2 = "";

	if($sort_len2_2 >1)		    
	{
	    $end_2.="\tNH:i:".($sort_len2_2);
	    if($i>0)
	    {
		$sorted2_2[$i]->{flag}+=256;
	    }
	}
	
	if($i < ($sort_len2_2-1))
	{
	    my $cc_2="=";
	    if($sorted2_2[$i]->{chr} ne $sorted2_2[$i+1]->{chr})				    
	    {
		$cc_2 = $sorted2_2[$i+1]->{chr};
	    }
	    $end_2.= "\tCC:Z:".$cc_2."\tCP:i:". $sorted2_2[$i+1]->{pos};
	}
	
	if($sort_len2_2 >1)		    
	{
	    $end_2.="\tHI:i:".$i;
	}
	
	if($version == 1)
	{
	    $end_2="";
	}
	
	print OUT $sorted2_2[$i]->{id}."\t".$sorted2_2[$i]->{flag}."\t".$sorted2_2[$i]->{chr}."\t".$sorted2_2[$i]->{pos}."\t".$sorted2_2[$i]->{score}."\t".$sorted2_2[$i]->{cigar}."\t".$sorted2_2[$i]->{rnext}."\t".$sorted2_2[$i]->{pos_next}."\t".$sorted2_2[$i]->{temp_len}."\t".$sorted2_2[$i]->{seq}."\t".$sorted2_2[$i]->{qual}."\tNM:i:".$sorted2_2[$i]->{nm}.$sorted2_2[$i]->{xs}.$end_2."\n";
    }
    close(REP);
    close(OUT);   
}
