diff SNV/SNVMix2_source/SNVMix2-v0.12.1-rc1/misc/snvmix2summary.pl @ 0:74f5ea818cea

Uploaded
author ryanmorin
date Wed, 12 Oct 2011 19:50:38 -0400
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SNV/SNVMix2_source/SNVMix2-v0.12.1-rc1/misc/snvmix2summary.pl	Wed Oct 12 19:50:38 2011 -0400
@@ -0,0 +1,64 @@
+#!/usr/bin/perl
+
+use strict;
+
+use Getopt::Std;
+my $opt_string = 'hi:c:t:';
+my %opt;
+getopts( "$opt_string", \%opt ) or usage();
+usage() if $opt{h};
+my $SNVMIX_FILE = "-";
+$SNVMIX_FILE = $opt{i} if $opt{i};
+my $TYPE = 2;
+$TYPE = $opt{c} if $opt{c};
+my $THRESHOLD = 0;
+$THRESHOLD = $opt{t} if $opt{t};
+if($TYPE != 2 && $TYPE != 3) { die("ERROR: Unknown class TYPE\n"); }
+
+print STDERR "Reading from ".($SNVMIX_FILE eq "-" ? "STDIN" : $SNVMIX_FILE)."\n";
+print STDERR "Calculating for max between AA".($TYPE == 2 ? " and {AB u BB}" : ", AB and BB")."\n";
+if($THRESHOLD) {
+	print STDERR "Applying threshold of $THRESHOLD, reporting only if ".($TYPE == 2 ? "P{AB u BB}" : "(P{AB} || P{BB})")." >= $THRESHOLD\n";
+}
+
+open(INPUT, "<$SNVMIX_FILE") || die("ERROR: Could not open '$SNVMIX_FILE' for reading\n");
+while(<INPUT>) {
+	chomp;
+	s/
//;
+	my $line = $_;
+	my ($chr_pos, $ref, $nref, $call_str, @extra) = split(/\t/, $line);
+	my ($ref_num, $nref_num, $pAA, $pAB, $pBB, $call) = split(/,/, $call_str);
+	my $snv = 0;
+	if($TYPE == 2) {
+		if($pAA < ($pAB + $pBB)) {
+			if( ($pAB + $pBB) >= $THRESHOLD) {
+				$snv = 1;
+			}
+		}
+	} elsif($TYPE == 3) {
+		if($call == 2 || $call == 3) {
+			if( $pAB >= $THRESHOLD || $pBB >= $THRESHOLD) {
+				$snv = 1;
+			}
+		}
+	} else {
+		die("ERROR, and a weird one, script shouldn't even BE in here...\n");
+	}
+	if($snv) {
+		#print "$chr_pos\t$ref\t".( $snv ? $nref : "-")."\t$snv\n";
+		print "$line\n";
+	}
+}
+close(INPUT);
+
+sub usage() {
+	print "Syntax:\n";
+	print "$0 [-i <file>] -c <TYPE> [-t <THRESHOLD>]\n";
+	print "\tIf file not given, STDIN is read\n";
+	print "\tTYPE is the number of classes to consider\n";
+	print "\t\t'2'\tconsiders only AA and {AB U BB} (default)\n";
+	print "\t\t'3'\tconsiders AA, AB and BB\n";
+	print "\tIf -t THRESHOLD is given, then SNVs will be reported\n";
+	print "\twhen the selected probability exceeds this\n";
+	exit;
+}