#!/usr/bin/perl
use strict;
use Getopt::Long;
if($#ARGV<0){

        my $usage=<<content
usage:

perl $0 \\
	--asIdToNumericId asId.to.numberId.tsv \\
	--jcecA3SS outputDir/JCEC.raw.input.A3SS.txt \\
	--jcecA5SS outputDir/JCEC.raw.input.A5SS.txt\\
	--jcecMXE outputDir/JCEC.raw.input.MXE.txt\\
	--jcecRI outputDir/JCEC.raw.input.RI.txt\\
	--jcecSE outputDir/JCEC.raw.input.SE.txt\\
	--jcA3SS outputDir/JC.raw.input.A3SS.txt\\
	--jcA5SS outputDir/JC.raw.input.A5SS.txt\\
	--jcMXE outputDir/JC.raw.input.MXE.txt\\
	--jcRI outputDir/JCEC.raw.input.RI.txt\\
	--jcSE outputDir/JCEC.raw.input.SE.txt\\
	--outputJCECtsv jcec.tsv \\
	--outputJCtsv jc.tsv

The script is used to calculate PSI values of alternative splicing events, and integrate 5 JCEC files and 5 JC files from rMATS into two single files respectively.

	--asIdToNumericId a mapping file from alternative spicing event to numeric identifer.
	--jcec[asType] the read count files of rMATS JCEC outputs.
	--jc[asType] the read count files of rMATS JC outputs.
	--outputJCECtsv the combined file for JCEC read count files.
	--outputJCtsv the combined file for JC read count files.

content
;
die $usage;
}

my ($asIdToNumericId);
my ($jcecA3SS, $jcecA5SS, $jcecSE, $jcecRI, $jcecMXE);
my ($jcA3SS, $jcA5SS, $jcSE, $jcRI, $jcMXE);
my ($outputJCECtsv, $outputJCtsv);

GetOptions(
        'asIdToNumericId=s'=>\$asIdToNumericId,
        'jcecA3SS=s'=>\$jcecA3SS,
        'jcecA5SS=s'=>\$jcecA5SS,
        'jcecSE=s'=>\$jcecSE,
	'jcecRI=s'=>\$jcecRI,
	'jcecMXE=s'=>\$jcecMXE,
        'jcA3SS=s'=>\$jcA3SS,
        'jcA5SS=s'=>\$jcA5SS,
        'jcSE=s'=>\$jcSE,
	'jcRI=s'=>\$jcRI,
	'jcMXE=s'=>\$jcMXE,
	'outputJCECtsv=s'=>\$outputJCECtsv,
	'outputJCtsv=s'=>\$outputJCtsv,	
);

my ($jeceFile, $jcFile, $asId);
my @type = ("A3SS", "A5SS", "SE", "RI", "MXE");
my @jcecFile = ($jcecA3SS, $jcecA5SS, $jcecSE, $jcecRI, $jcecMXE);
my @jcFile = ($jcA3SS, $jcA5SS, $jcSE, $jcRI, $jcMXE);
my ($line, $num);
my (%a3ssMap, %a5ssMap, %seMap, %mxeMap, %riMap);
open FF, "<$asIdToNumericId";
while($line = <FF>){
	chomp($line);
	($asId, $num) = split(/\t/, $line);
	if(substr($asId, 4, 2) eq "A3"){
		$a3ssMap{$num} = $asId;
	}elsif(substr($asId, 4, 2) eq "A5"){
		$a5ssMap{$num} = $asId;
	}elsif(substr($asId, 4, 2) eq "SE"){
                $seMap{$num} = $asId;
        }elsif(substr($asId, 4, 2) eq "RI"){
                $riMap{$num} = $asId;
        }elsif(substr($asId, 4, 2) eq "MX"){
                $mxeMap{$num} = $asId;
        }
}
close FF;

my (@refA3ssAsId, @refA5ssAsId, @refSeAsId, @refRiAsId, @refMxeAsId);
@refA3ssAsId = keys(%a3ssMap);
@refA5ssAsId = keys(%a5ssMap);
@refSeAsId = keys(%seMap);
@refRiAsId = keys(%riMap);
@refMxeAsId = keys(%mxeMap);
#print join("\t", $#refA3ssAsId, $#refA5ssAsId, $#refSeAsId, $#refRiAsId, $#refMxeAsId) . "\n";
my ($num, $asId, $line);
my ($ID, $IJC_SAMPLE_1, $SJC_SAMPLE_1, $IJC_SAMPLE_2, $SJC_SAMPLE_2, $IncFormLen, $SkipFormLen);
my ($psi);
my ($type, $jcecFile, $jcFile, $asMapHashRef, $refAsNum, $mapAsNum);
open JCEC, ">$outputJCECtsv";
print JCEC join("\t", "ASID", "IJC_SAMPLE", "SJC_SAMPLE", "IncFormLen", "SkipFormLen", "psi") . "\n";
open JC, ">$outputJCtsv";
print JC join("\t", "ASID", "IJC_SAMPLE", "SJC_SAMPLE", "IncFormLen", "SkipFormLen", "psi") . "\n";
for(my $i=0; $i<=4; $i++){
	$type = $type[$i];
	if($type eq "A3SS"){
		$asMapHashRef = \%a3ssMap;
		$refAsNum = $#refA3ssAsId + 1;
	}elsif($type eq "A5SS"){
		$asMapHashRef = \%a5ssMap;
		$refAsNum = $#refA5ssAsId + 1;
	}elsif($type eq "SE"){
		$asMapHashRef = \%seMap;
		$refAsNum = $#refSeAsId + 1;
	}elsif($type eq "RI"){
		$asMapHashRef = \%riMap;
		$refAsNum = $#refRiAsId + 1;
	}elsif($type eq "MXE"){
		$asMapHashRef = \%mxeMap;
		$refAsNum = $#refMxeAsId + 1;
	}

	$jcecFile = $jcecFile[$i];
	open FF, "<$jcecFile";
	$line = <FF>;
	$mapAsNum = 0;
	while($line = <FF>){
		$mapAsNum++;
		chomp($line);
		($ID, $IJC_SAMPLE_1, $SJC_SAMPLE_1, $IJC_SAMPLE_2, $SJC_SAMPLE_2, $IncFormLen, $SkipFormLen) = split(/\t/, $line);
		$asId = $asMapHashRef->{$ID};
		# discard psi in case of no read coverage
		next if($IJC_SAMPLE_1 + $SJC_SAMPLE_1 == 0);
		$psi = ($IJC_SAMPLE_1/$IncFormLen)/(($IJC_SAMPLE_1/$IncFormLen)+($SJC_SAMPLE_1/$SkipFormLen));
		print JCEC join("\t", $asId, $IJC_SAMPLE_1, $SJC_SAMPLE_1, $IncFormLen, $SkipFormLen, $psi) . "\n";
	}
	close FF;
	if($refAsNum != $mapAsNum){
		print "Your reference alternative splicing files are not right.";
		goto OVER;
	}

	$jcFile = $jcFile[$i];
	open FF, "<$jcFile";
	$line = <FF>;
	$mapAsNum = 0;
	while($line = <FF>){
		$mapAsNum++;
		chomp($line);
		($ID, $IJC_SAMPLE_1, $SJC_SAMPLE_1, $IJC_SAMPLE_2, $SJC_SAMPLE_2, $IncFormLen, $SkipFormLen) = split(/\t/, $line);
		$asId = $asMapHashRef->{$ID};
		# discard psi in case of no read coverage
		next if($IJC_SAMPLE_1 + $SJC_SAMPLE_1 == 0);
		$psi = ($IJC_SAMPLE_1/$IncFormLen)/(($IJC_SAMPLE_1/$IncFormLen)+($SJC_SAMPLE_1/$SkipFormLen));
		print JC join("\t", $asId, $IJC_SAMPLE_1, $SJC_SAMPLE_1, $IncFormLen, $SkipFormLen, $psi) . "\n";
	}
	close FF;
	if($refAsNum != $mapAsNum){
		print "Your reference alternative splicing files are not right.";
		goto OVER;
	}

}
close JCEC;
close JC;
OVER: