annotate bin/ApplyCutoff.jie.pl @ 5:2ebca9da5e42 draft default tip

planemo upload

#apply our cuttoff hash table on the IR calculated by Jie
#Modified from Martin's code
use strict;

use Cwd;
my $PROG = $0;
my $CUR_DIR = Cwd::abs_path(Cwd::cwd());
my $PROG_ABS_PATH = Cwd::abs_path($PROG);
my $SrcFolder=`dirname $PROG_ABS_PATH`;
chomp($SrcFolder);

my %cutoff;
my @Exlen;

my $cutoff_level=$ARGV[1];
my $JunctionCut = $ARGV[2];
my $noirm = $ARGV[3];

my $cutoff_level_index=7;

 $cutoff_level_index=8 if $cutoff_level eq "H";
$cutoff_level_index=6 if $cutoff_level eq "L";

open(CUT,"$SrcFolder/../cutoffs/cutoff.pair.0".$cutoff_level_index.".txt") || die "cutoff file not found $!\n";

while(<CUT>){
	chomp;
	my @a=split(/\t/,$_);
	push @Exlen,$a[0];
	$cutoff{$a[0]}=$a[1];
}
close(CUT);

open(IN,$ARGV[0]);

while(<IN>){
	chomp;
	my @a=split(/\t/,$_);
	my $Ez='Ez=yes';
	my $print=$_;
	if($a[0]=~m/#/g){next}
	my $eventid=substr($a[0],0,2);
	my $bir =$a[2];
	$bir =$a[1] if($noirm eq "noirm");
	my $j12 = $a[8];
	my $j23 = $a[9];
	my $j13 = $a[10];
	my $cov1=$a[11];
	my $cov2=$a[12];
	my $cov3=$a[13];
	my $siz1=$a[15];
	my $siz2=$a[16];
	my $siz3=$a[17];


	my $stat1='exon1='.cutoff($siz1,$cov1,\@Exlen,%cutoff);
	my $stat2='exon2='.cutoff($siz2,$cov2,\@Exlen,%cutoff);
	my $stat3='exon3='.cutoff($siz3,$cov3,\@Exlen,%cutoff);
	if($stat1 eq "exon1=yes" and $stat3 eq "exon3=yes") 
	{
		#$Ez="passed";
		$Ez=$eventid if $eventid eq "AA";
		$Ez=$eventid if $eventid eq "AD";
		$Ez=$eventid if $eventid eq "IR";
		if ($eventid eq "CS" or $eventid eq "CA" or $eventid eq "ME") 
		{
			if($bir >0.9)
			{
				$Ez = "CS";
			}
			else
			{
				$Ez = "CA";
			}
		}
		
	}
	else
	{
		#$Ez="declined";
		$Ez = "na";
	}
	if( ($j12<$JunctionCut or $j23<$JunctionCut) and $j13 <$JunctionCut)
	{
		$Ez = "na";
	}
	print $print,"\t",$stat1,"\t",$stat2,"\t",$stat3,"\t",$Ez,"\n";
}
close(IN);
####################################################################

sub cutoff{
	my($s,$c,$E,%cutoff)=@_;
	my @Exlen=@$E;
	if($c eq 'NA'){return('NA')}
	my $range=$Exlen[$#Exlen];
	foreach my $l(@Exlen){if($s<$l){$range=$l;last}}
        if($c<$cutoff{$range}){return('no')}
	return('yes')
}


sub contain{
	my ($a,@a)=@_;
	foreach(@a){if($a eq $_){return(1)}}
	return(0)
}