Mercurial > repos > ryanmorin > nextgen_variant_identification
comparison 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 |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:74f5ea818cea |
|---|---|
| 1 #!/usr/bin/perl | |
| 2 | |
| 3 use strict; | |
| 4 | |
| 5 use Getopt::Std; | |
| 6 my $opt_string = 'hi:c:t:'; | |
| 7 my %opt; | |
| 8 getopts( "$opt_string", \%opt ) or usage(); | |
| 9 usage() if $opt{h}; | |
| 10 my $SNVMIX_FILE = "-"; | |
| 11 $SNVMIX_FILE = $opt{i} if $opt{i}; | |
| 12 my $TYPE = 2; | |
| 13 $TYPE = $opt{c} if $opt{c}; | |
| 14 my $THRESHOLD = 0; | |
| 15 $THRESHOLD = $opt{t} if $opt{t}; | |
| 16 if($TYPE != 2 && $TYPE != 3) { die("ERROR: Unknown class TYPE\n"); } | |
| 17 | |
| 18 print STDERR "Reading from ".($SNVMIX_FILE eq "-" ? "STDIN" : $SNVMIX_FILE)."\n"; | |
| 19 print STDERR "Calculating for max between AA".($TYPE == 2 ? " and {AB u BB}" : ", AB and BB")."\n"; | |
| 20 if($THRESHOLD) { | |
| 21 print STDERR "Applying threshold of $THRESHOLD, reporting only if ".($TYPE == 2 ? "P{AB u BB}" : "(P{AB} || P{BB})")." >= $THRESHOLD\n"; | |
| 22 } | |
| 23 | |
| 24 open(INPUT, "<$SNVMIX_FILE") || die("ERROR: Could not open '$SNVMIX_FILE' for reading\n"); | |
| 25 while(<INPUT>) { | |
| 26 chomp; | |
| 27 s/ | |
| 28 //; | |
| 29 my $line = $_; | |
| 30 my ($chr_pos, $ref, $nref, $call_str, @extra) = split(/\t/, $line); | |
| 31 my ($ref_num, $nref_num, $pAA, $pAB, $pBB, $call) = split(/,/, $call_str); | |
| 32 my $snv = 0; | |
| 33 if($TYPE == 2) { | |
| 34 if($pAA < ($pAB + $pBB)) { | |
| 35 if( ($pAB + $pBB) >= $THRESHOLD) { | |
| 36 $snv = 1; | |
| 37 } | |
| 38 } | |
| 39 } elsif($TYPE == 3) { | |
| 40 if($call == 2 || $call == 3) { | |
| 41 if( $pAB >= $THRESHOLD || $pBB >= $THRESHOLD) { | |
| 42 $snv = 1; | |
| 43 } | |
| 44 } | |
| 45 } else { | |
| 46 die("ERROR, and a weird one, script shouldn't even BE in here...\n"); | |
| 47 } | |
| 48 if($snv) { | |
| 49 #print "$chr_pos\t$ref\t".( $snv ? $nref : "-")."\t$snv\n"; | |
| 50 print "$line\n"; | |
| 51 } | |
| 52 } | |
| 53 close(INPUT); | |
| 54 | |
| 55 sub usage() { | |
| 56 print "Syntax:\n"; | |
| 57 print "$0 [-i <file>] -c <TYPE> [-t <THRESHOLD>]\n"; | |
| 58 print "\tIf file not given, STDIN is read\n"; | |
| 59 print "\tTYPE is the number of classes to consider\n"; | |
| 60 print "\t\t'2'\tconsiders only AA and {AB U BB} (default)\n"; | |
| 61 print "\t\t'3'\tconsiders AA, AB and BB\n"; | |
| 62 print "\tIf -t THRESHOLD is given, then SNVs will be reported\n"; | |
| 63 print "\twhen the selected probability exceeds this\n"; | |
| 64 exit; | |
| 65 } |