diff bin/ApplyCutoff.jie.pl @ 1:adc0f7765d85 draft

planemo upload

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bin/ApplyCutoff.jie.pl	Thu Sep 07 15:06:58 2017 -0400
@@ -0,0 +1,108 @@
+#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)
+}
+