annotate tools/human_genome_variation/snpFreq2.pl @ 0:9071e359b9a3

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
1 #!/usr/bin/env perl
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 use strict;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4 use warnings;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6 #expected input: path to file, cols of counts (2 sets of 3), threshold
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7 if (!@ARGV or scalar @ARGV != 9) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 print "usage snpFreq.pl /path/to/snps.txt <6 column numbers(1 based) with counts for alleles, first one group then another> #threshold outfile\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 #get and verify inputs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13 my $file = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 my $a1 = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 if ($a1 =~ /\D/ or $a1 < 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 print "Error the column number, must be an integer greater than or equal to 1. Got $a1\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19 my $a2 = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20 if ($a2 =~ /\D/ or $a2 < 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 print "Error the column number, must be an integer greater than or equal to 1. Got $a2\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24 my $a3 = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25 if ($a3 =~ /\D/ or $a3 < 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 print "Error the column number, must be an integer greater than or equal to 1. Got $a3\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29 my $b1 = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30 if ($b1 =~ /\D/ or $b1 < 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 print "Error the column number, must be an integer greater than or equal to 1. Got $b1\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 my $b2 = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 if ($b2 =~ /\D/ or $b2 < 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36 print "Error the column number, must be an integer greater than or equal to 1. Got $b2\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39 my $b3 = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40 if ($b3 =~ /\D/ or $b3 < 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41 print "Error the column number, must be an integer greater than or equal to 1. Got $b3\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 my $thresh = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 if ($thresh !~ /^\d*\.?\d+$/) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 print "Error the threshold must be a number. Got $thresh\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 }elsif ($thresh > .3) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 print "Error the threshold can not be greater than 0.3 got $thresh\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 exit 1;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 my $outfile = shift @ARGV;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 #run a fishers exact test (using R) on whole table
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 my $cmd = qq|options(warn=-1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 tab <- read.table('$file', sep="\t")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 size <- length(tab[,1])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 width <- length(tab[1,])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 x <- 1:size
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 y <- matrix(data=0, nr=size, nc=6)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 for(i in 1:size) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 m <- matrix(c(tab[i,$a1], tab[i,$b1], tab[i,$a2], tab[i,$b2], tab[i,$a3], tab[i,$b3]), nrow=2)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 t <- fisher.test(m)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 x[i] <- t\$p.value
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65 if (x[i] >= 1) {
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66 x[i] <- .999
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67 }
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68 n <- (tab[i,$a1] + tab[i,$a2] + tab[i,$a3] + tab[i,$b1] + tab[i,$b2] + tab[i,$b3])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 n_a <- (tab[i,$a1] + tab[i,$a2] + tab[i,$a3])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 y[i,1] <- ((tab[i,$a1] + tab[i,$b1])*(n_a))/n
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 y[i,1] <- round(y[i,1],3)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 y[i,2] <- ((tab[i,$a2] + tab[i,$b2])*(n_a))/n
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 y[i,2] <- round(y[i,2],3)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 y[i,3] <- ((tab[i,$a3] + tab[i,$b3])*(n_a))/n
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 y[i,3] <- round(y[i,3],3)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 n_b <- (tab[i,$b1] + tab[i,$b2] + tab[i,$b3])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 y[i,4] <- ((tab[i,$a1] + tab[i,$b1])*(n_b))/n
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 y[i,4] <- round(y[i,4],3)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 y[i,5] <- ((tab[i,$a2] + tab[i,$b2])*(n_b))/n
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 y[i,5] <- round(y[i,5],3)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 y[i,6] <- ((tab[i,$a3] + tab[i,$b3])*(n_b))/n
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 y[i,6] <- round(y[i,6],3)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83 }|;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84 #results <- data.frame(tab[1:size,1:width], x[1:size])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 #write.table(results, file="$outfile", row.names = FALSE ,col.names = FALSE,quote = FALSE, sep="\t")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 #q()|;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 my $cmd2 = qq|suppressPackageStartupMessages(library(qvalue))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 qobj <- qvalue(x[1:size], lambda=seq(0,0.90,$thresh), pi0.method="bootstrap", fdr.level=0.1, robust=FALSE, smooth.log.pi0 = FALSE)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 q <- qobj\$qvalues
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 results <- data.frame(tab[1:size,1:width], y[1:size,1:6], x[1:size], q[1:size])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 write.table(results, file="$outfile", row.names = FALSE ,col.names = FALSE,quote = FALSE, sep="\t")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 q()|;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95 #for TESTING
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96 my $pr = qq|results <- data.frame(tab[1:size,1:width], y[1:size,1:6], x[1:size])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97 write.table(results, file="$outfile", row.names = FALSE ,col.names = FALSE,quote = FALSE, sep="\t")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98 q()|;
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100 open(FT, "| R --slave --vanilla")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 or die "Couldn't call fisher.text, $!\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 print FT $cmd, "\n"; #fisher test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103 print FT $cmd2, "\n"; #qvalues and results
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 #print FT $pr, "\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105 close FT or die "Couldn't finish fisher.test, $!\n";
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107 exit;