use strict;
use Getopt::Long;
use DBI;
if($#ARGV<0){
        my $usage=<<content
usage:
perl $0 \\
	--python path-to-python \\
	--prepDE path-to-prepDE.py \\
	--readLen 101 \\
	--transcriptomeFile transcriptome.gtf \\
	--geneAbundanceFile gene.abundance.tab \\
	--outputGeneExprFile gene.expr.tsv \\
	--outputTrsptExprFile transcript.expr.tsv


The script is used to extract gene and transcript expression data from gene abundance and transcriptome assembly files.
The expression data include FPKM/TPM expression levels, read counts and coverage.

About inputs:

	You can get the path of python through the command "which python", for example "/usr/bin/python". You can also get the version of python through the command "python --version".

	If your python is Python3 version, only prepDE.py3 is available. Both prepDE.py and prepDE.py3 can be downloaded from the official website of stringTie2.

	readLen is the RNA-seq read length, for example 101.

	transcriptomeFile and geneAbundanceFile is the two output files of stringTie2 in the previous step.

content
;
die $usage;
}

my ($python, $prepDE, $readLen, $transcriptomeFile, $geneAbundanceFile, $outputGeneExprFile, $outputTrsptExprFile);

GetOptions(
	'python=s'=>\$python,
	'prepDE=s'=>\$prepDE,
        'readLen=s'=>\$readLen,
	'transcriptomeFile=s'=>\$transcriptomeFile,
        'geneAbundanceFile=s'=>\$geneAbundanceFile,
	'outputGeneExprFile=s'=>\$outputGeneExprFile,
	'outputTrsptExprFile=s'=>\$outputTrsptExprFile,
);

# detect dos2unix
my $cmd = `which dos2unix 2>/dev/null`;
if($cmd eq ""){
	print "dos2unix is not available. please specify dos2unix.\n";
	exit;
}
# detect input files
if(not -e $transcriptomeFile ){
	print "$transcriptomeFile is not available.\n";
	exit;
}
if(not -e $geneAbundanceFile){
	print "$geneAbundanceFile is not available.\n";
	exit;
}
# detect python version
$cmd = `$python --version`;
if($cmd=~/Python 3.*/ and not $prepDE=~/prepDE\.py3/){
	print "Only prepDE.py3 is available for Python 3\n";
	exit;
}


my (%geneExpr, $geneExprHref, %trsptExpr, $trsptExprHref);
my ($line, $geneId, $geneName, $ref, $strand, $start, $end, $cov, $fpkm, $tpm, @field, $trsptId, @trsptId, $cmd, @geneId);
my ($geneSymbol, $readCount, $symbol);

# hash for gene and transcript expression
%geneExpr = ();
$geneExprHref = \%geneExpr;
%trsptExpr = ();
$trsptExprHref = \%trsptExpr;

# read gene abundance file
open FF, "<$geneAbundanceFile";
# Gene ID 	Gene Name       Reference       Strand  Start   End     Coverage        FPKM    	TPM
# Os01g0100100  -       	1       	+       2983    10815   63.793186       6.685385        10.758520
# Os01g0100200  -       	1       	+       11218   12435   4.926498        0.515535        0.829630
<FF>;
while($line=<FF>){
	chomp($line);
	($geneId, $geneName, $ref, $strand, $start, $end, $cov, $fpkm, $tpm) = split(/\t/, $line);
	$geneExprHref->{$geneId}->{"cov"} = $cov;
	$geneExprHref->{$geneId}->{"fpkm"} = $fpkm;
	$geneExprHref->{$geneId}->{"tpm"} = $tpm;
	$geneExprHref->{$geneId}->{"readCount"} = 0;
}
close FF;

# read transcriptome gtf file to obtain transcript coverage, fpkm, tpm
open FF, "<$transcriptomeFile";
while($line=<FF>){
	chomp($line);
	next if($line=~/^#/);
	
	@field = split(/\t/, $line);
	next if($field[2] ne "transcript");
	($geneId, $trsptId, $cov, $fpkm, $tpm) = ("", "", 0, 0, 0);
	&getAttr($field[8], \$geneId, \$trsptId, \$cov, \$fpkm, \$tpm);
	if($trsptId ne ""){
		$trsptExprHref->{$trsptId}->{"cov"} = $cov;
		$trsptExprHref->{$trsptId}->{"fpkm"} = $fpkm;
		$trsptExprHref->{$trsptId}->{"tpm"} = $tpm;
		$trsptExprHref->{$trsptId}->{"readCount"} = 0;
	}
}
close FF;

my $id = time();
open WW, ">/tmp/$id.prepDE.input.txt";
print WW "exptId " . $transcriptomeFile;
close WW;

# call prepDE[3].py to obtain read count on gene and transcript
$cmd = "$python $prepDE -i /tmp/$id.prepDE.input.txt -g /tmp/$id.gene.count.tsv -t /tmp/$id.trspt.count.tsv -l $readLen";
system($cmd);

system("rm -rf /tmp/$id.prepDE.input.txt");

# extract gene read count into gene hash
system("dos2unix /tmp/$id.gene.count.tsv &>/dev/null");
open FF, "</tmp/$id.gene.count.tsv";
<FF>;
while($line = <FF>){
	chomp($line);
	# Os09g0127800|OsWD40-165,1810
	# Os03g0760800|GASR1,190
	($geneSymbol, $readCount) = split(/,/, $line);
	if($geneSymbol=~/(.*)\|(.*)/){
		$geneId = $1;
	}else{
		$geneId = $geneSymbol;
	}
	$geneExprHref->{$geneId}->{"readCount"} += $readCount;
}
close FF;
system("rm -rf /tmp/$id.gene.count.tsv");

# extract transcript read count into transcript hash
system("dos2unix /tmp/$id.trspt.count.tsv &>/dev/null");
open FF, "</tmp/$id.trspt.count.tsv";
<FF>;
while($line = <FF>){
	chomp($line);
	# transcript_id,SRX5808211
	# SRX5173465.28905.3,45
	($trsptId, $readCount) = split(/,/, $line);
	$trsptExprHref->{$trsptId}->{"readCount"} =$readCount;
}
close FF;
system("rm -rf /tmp/$id.trspt.count.tsv");


# write final gene expression file
@geneId = ();
@geneId = keys(%geneExpr);
open WW, ">$outputGeneExprFile";
print WW join("\t", "geneId", "readCount", "Coverage", "FPKM", "TPM") . "\n";
foreach $geneId(@geneId){
	print WW join("\t", $geneId, $geneExprHref->{$geneId}->{"readCount"}, $geneExprHref->{$geneId}->{"cov"}, $geneExprHref->{$geneId}->{"fpkm"}, $geneExprHref->{$geneId}->{"tpm"}) . "\n";
}
close WW;

# write final transcript expression file 
@trsptId = ();
@trsptId = keys(%trsptExpr);
open WW, ">$outputTrsptExprFile";
print WW join("\t", "trsptId", "readCount", "Coverage", "FPKM", "TPM") . "\n";
foreach $trsptId(@trsptId){
	print WW join("\t", $trsptId, $trsptExprHref->{$trsptId}->{"readCount"}, $trsptExprHref->{$trsptId}->{"cov"}, $trsptExprHref->{$trsptId}->{"fpkm"}, $trsptExprHref->{$trsptId}->{"tpm"}) . "\n";
}
close WW;


sub getAttr{
	my ($attrString, $geneId, $trsptId, $cov, $fpkm, $tpm) = @_;
	my (@attr, $attr);
	@attr = split(/;/, $attrString);
	foreach $attr(@attr){
		if($attr=~/gene_id "(.*)"/){
			$$geneId = $1;
		}elsif($attr=~/transcript_id "(.*)"/){
			$$trsptId = $1;
		}elsif($attr=~/cov "(.*)"/){
			$$cov = $1;
		}elsif($attr=~/FPKM "(.*)"/){
			$$fpkm = $1;
		}elsif($attr=~/TPM "(.*)"/){
			$$tpm = $1;
		}
	}
}