use warnings;
use Getopt::Long;

# --- command line --- #
my $opthelp=0;
my $infile; #ReadPos_ReadCount FILE 
my $outfile; #Seq and Coverages Vector FILE
my $ref;
my $pos = "all"; 
my $derive=0; #derived RNA由来の形状を作るモード（$maxLengthとペアで使用、また遺伝子の長さ-20のリードも計算にいれない）
my $maxLength=1000000; 
my $minLength=0; 
my $tmpdir = "/tmp";

my $opt=1;  #secret option : 1塩基でもアノテーション上に乗っているリードをCoverageに加えるか(1)、それとも2/3以上アノテーション上にあるリードのみをCoverageに加えるか(2)

GetOptions(
    'help' => \$opthelp,
    'in=s' => \$infile,
    'out=s' => \$outfile,
    'tmp=s' => \$tmpdir,
    'ref=s' => \$ref,
    'opt=s' => \$opt,
    'pos=s' => \$pos,
    'derive' => \$derive,
    'maxlen=i' => \$maxLength,
    'minlen=i' => \$minLength,
    );

if($opthelp){&printhelp; exit(1);}

## --- check --- ##
open(IN,"$infile") || die "input FILE error :\n  Can't open input FILE : $infile $!";
print "input file : $infile\n";
open(OUT,"> $outfile") || die "output FILE error :\n  Can't open output FILE : $outfile $!"; 
print "output file : $outfile\n";
if(!defined($ref))
{
    print "Reference FILE error :\n";
    print "  please input Reference FASTA File : -ref FILE\n"; 
    exit(1);
}
elsif(!-e $ref)
{
    print "Reference FILE error :\n";
    print "  Can't open Reference FILE : $ref $!"; 
    exit(1);
}
print "reference file : $ref\n";
if($pos ne "all" && $pos ne "exon")
{
    print "pos options error :\n";
    print "  please input : -pos <STRING> (possible value=all or exon)\n";
    exit(1);
}
if (!-d "$tmpdir"){mkdir "$tmpdir" || die "can't make directory $tmpdir : please input other directory for tmp directory : -tmp DIR $!";}

if(!defined($opt)){print "input -opt 1 or 2\n"; exit(1);}
## --- --- ##
# --- --- #

# --- main --- #
while(<IN>)
{
    my $geneName=$_; chomp($geneName); #geneName parse
    my @GenePos=&parseGenePos;
    my @Coverages=&makeCoverages(@GenePos);

    print OUT "$geneName\n";
    if($pos eq "all")
    {
	&printExonPos1(@GenePos);
	my @StartEnd = ($GenePos[0],$GenePos[$#GenePos]);
	&printSeq($geneName,\@StartEnd,$pos);
	for(my $i=0; $i<$GenePos[$#GenePos]-$GenePos[0]+1; $i++){print OUT "$Coverages[$i]".",";}
    }
    elsif($pos == "exon")
    {
	&printExonPos2(@GenePos);
	&printSeq($geneName,\@GenePos,$pos);
	for(my $i=0; $i<=$#GenePos; $i+=2)
	{
	    my ($s,$e) = ($GenePos[$i]-$GenePos[0],$GenePos[$i+1]-$GenePos[0]);
	    for(my $j=$s; $j<=$e; $j++){print OUT "$Coverages[$j]".",";}
	}
    }
    print OUT "\n";
    print OUT "end\n";
} 
close(IN);
# --- --- #

# --- sub functions --- #
sub printSeq
{
    my ($geneName, $rGenePos,$pos) = @_;
    my @geneinfos=split(/\//,$geneName);
    my $Seq = "";
    for(my $i=0; $i<=$#{$rGenePos}; $i+=2)
    {
	my $tmpbedfile="$tmpdir/tmp.bed";
	&makeBedInfo($i, $tmpbedfile, \@geneinfos, $rGenePos);
	my $tmpfafile="$tmpdir/tmp.fa";
	system("fastaFromBed -fi $ref -bed $tmpbedfile -fo $tmpfafile -s");
	$Seq.=&returnSeq($tmpfafile);
    }
    print OUT "$Seq\n";
}

sub makeBedInfo
{
    my ($i, $tmpbedfile, $rgeneinfos, $rGenePos) = @_;
    my $start = $$rGenePos[$i]-$$rGenePos[0]+$$rgeneinfos[3]-1;
    my $end = $$rGenePos[$i+1]-$$rGenePos[0]+$$rgeneinfos[3];
    open(TMPBED,"> $tmpbedfile") || die "Can't open tmp bed file : $tmpbedfile : $!"; 
    print TMPBED "$$rgeneinfos[1]\t"; #print strand
    print TMPBED "$start\t"; #print start pos
    print TMPBED "$end\t"; #print end pos
    print TMPBED "$$rgeneinfos[0]\t"; #print name
    print TMPBED "0\t"; #print score
    print TMPBED "$$rgeneinfos[2]\n"; #print strand

    close(TMPBED);
}

sub returnSeq
{
    my $tmpfafile=$_[0];
    open(TMPFA,"$tmpfafile") || die "Can't open tmp fasta file : $tmpfafile $!"; 
    my $head=<TMPFA>;
    my $seq=<TMPFA>; chomp($seq);
    close(TMPFA);
    return $seq;
}

sub parseGenePos
{
    my $GenePosStr=<IN>; chomp($GenePosStr);
    my @temp=split(/;/,$GenePosStr);
    my @GenePos; 
    for(my $i=1; $i<=$#temp; $i++)
    {
	my @temp2=split(/,/,$temp[$i]); #temp2[0]:ExonStart,temp2[1]:ExonEnd
	push(@GenePos,@temp2);
	@temp2=();
    }
    return @GenePos;
}

sub makeCoverages
{
    local @GenePos=@_;
    local @Coverages; &InitCoverages;

    while(<IN>)
    {
	chomp($_);
	if($_ eq "end"){last;}
	else
	{
	    my @ReadPosCount=&returnReadPosCount($_);
	    if($derive)
	    {
		if( &returnReadLength(@ReadPosCount) > $maxLength ){next;} #リード長がmaxLenghを超えるときは、Coverageの計算に入れない。
		if( &returnReadLength(@ReadPosCount) < $minLength ){next;} #リード長がminLenghを超えるときは、Coverageの計算に入れない。
		my $AnnotationLength=$GenePos[$#GenePos]-$GenePos[0]+1; #Annotation上の長さ（イントロンも含めた）
		if( &returnReadLength(@ReadPosCount) > $AnnotationLength - 20){next;} #リード長が Annotation上の長さ-20 を超える場合,Coverageの計算に入れない。
	    }

	    if($opt == "1"){ &plusReadCount_InCoverages1(@ReadPosCount); }
	    elsif($opt == "2"){ &plusReadCount_InCoverages2(@ReadPosCount); }
	}
    }
    return @Coverages;
}

sub InitCoverages
{
    my $GeneStart=$GenePos[0];
    my $GeneEnd=$GenePos[$#GenePos];
    for(my $i=0; $i<$GeneEnd-$GeneStart+1; $i++)
    {
	$Coverages[$i]=0;
    }
}

sub returnReadPosCount
{
    my @temp=split(/\t/,$_[0]);
    my @temp2=split(/;/,$temp[0]);
    my @ReadPosCount;

    for(my $i=0; $i<=$#temp2; $i++)
    {
	my @temp3=split(/,/,$temp2[$i]); #temp3[0]:ReadExonStart,temp3[1]:ReadExonEnd,temp3[2]:ReadExonCount
	push(@ReadPosCount,@temp3);
	@temp3=();
    }

    return @ReadPosCount;
}

sub returnReadLength
{
    my $ReadLength=0;
    for(my $i=0; $i<=$#_; $i+=3)
    {
	$ReadLength+=$_[$i+1]-$_[$i]+1;
    }
    return $ReadLength;
}

sub plusReadCount_InCoverages1
{
    my ($GeneStart,$GeneEnd) = ($GenePos[0],$GenePos[$#GenePos]);
    
    for(my $j=0; $j<=$#_; $j+=3)
    {
	my($ReadExonStart,$ReadExonEnd,$ReadExonCount)=($_[$j],$_[$j+1],$_[$j+2]);
	if($ReadExonStart > $GeneEnd){next;}
	elsif($ReadExonEnd < $GeneStart){next;}
	
	my $s; 
	my $e;
	if($ReadExonStart <= $GeneStart){$s=0;}
	elsif($ReadExonStart > $GeneStart){$s=$ReadExonStart-$GeneStart;}
	if($ReadExonEnd >= $GeneEnd){$e=$GeneEnd-$GeneStart;}
	elsif($ReadExonEnd < $GeneEnd){$e=$ReadExonEnd-$GeneStart;}
	
	for(my $t=$s; $t<=$e; $t++){$Coverages[$t]+=$ReadExonCount;}
    }
}

sub plusReadCount_InCoverages2
{
    my @se; 
    my $count=0;
    my $ReadExonCount=$_[2];
    my $ReadLength=&returnReadPosCount(@_);

    my ($GeneStart,$GeneEnd) = ($GenePos[0],$GenePos[$#GenePos]);

    for(my $j=0; $j<=$#_; $j+=3)
    {
	my($ReadExonStart,$ReadExonEnd)=($_[$j],$_[$j+1]);
	if($ReadExonStart > $GeneEnd){next;}
	elsif($ReadExonEnd < $GeneStart){next;}
	
	if($ReadExonStart <= $GeneStart){$se[$count]=0;}
	elsif($ReadExonStart > $GeneStart){$se[$count]=$ReadExonStart-$GeneStart;}
	if($ReadExonEnd >= $GeneEnd){$se[$count+1]=$GeneEnd-$GeneStart;}
	elsif($ReadExonEnd < $GeneEnd){$se[$count+1]=$ReadExonEnd-$GeneStart;}	    
	$count+=2;
    }
    my $flag=&ReturnFlag($ReadLength,@se);
    for(my $i=0; $i<=$#se; $i+=2)
    {
	my($s,$e)=($se[$i],$se[$i+1]);
	for(my $t=$s; $t<=$e; $t++){$Coverages[$t]+=$ReadExonCount;}
    }
}

sub ReturnFlag
{
    my($ReadLength,@se)=@_;
    my $ReadLengthInGene=0;
    for(my $i=0; $i<=$#se; $i+=2){ $ReadLengthInGene+=$se[$i+1]-$se[$i]+1; }

    my $ratio=$ReadLengthInGene/$ReadLength;
    if($ratio >= 2/3){return 1;}
    else             {return 0;}
}

sub printExonPos1
{
    my @GenePos = @_;
    print OUT "ExonPositions";
    for(my $i=0; $i<=$#GenePos; $i+=2)
    {
	my $s=$GenePos[$i]-$GenePos[0]; my $e=$GenePos[$i+1]-$GenePos[0];
	print OUT ";"."$s".","."$e";
    }
    print OUT "\n";
}

sub printExonPos2
{
    my @GenePos = @_;
    print OUT "ExonPositions";
    my $IntronLength=0;
    for(my $i=2; $i<=$#GenePos; $i+=2)
    {
	$IntronLength+=$GenePos[$i]-$GenePos[$i-1]-1;
    }
    my $revisedGeneEnd=($GenePos[$#GenePos]-$IntronLength)-$GenePos[0];
    print OUT ";"."0".","."$revisedGeneEnd\n";
}
# --- --- #

sub printhelp
{
    print "\nusage : perl makeSeqCoverages.pl -in INFILE -out OUTFILE -ref REFERENCE_FASTA_FILE [other Options]\n\n";
    print "Options : \n";
    print "  -in   FILE       : input ReadPos_ReadCount File ( made by SortReadsByRegions) \n";
    print "  -out  FILE       : output Sequence and Coverages_Vector File ( use this as input file for SHARAKU) \n";
    print "  -ref  FILE       : Reference FASTA file (genome file)\n";
    print "  -tmp  DIR        : tmp directry (defalut = /tmp) \n";
    print "  -pos  <STRING>   : If you want to use only exon region, use this option as -pos exon. (possible option = exon or all, defalut = exon)\n";
    print "  -derive          : Make Coverages Vector consisting of derived-RNA Reads\n";
    print "                     (details as below)\n";
    print "  -maxlen <INT>    : Max length of reference gene which you want to use (default= 1000000) \n";
    print "  -minlen <INT>    : Minimum length of reference gene which you want to use (default= 0) \n";
    print "\n";
}
