From b8265695d2e462c08fbbd8c4e26fa42e54e94bd3 Mon Sep 17 00:00:00 2001
From: Wen-Yu Chung <wychung@bx.psu.edu>
Date: Tue, 8 Apr 2008 15:05:43 +0000
Subject: [PATCH] add a new tool (poisson2test) with funcitonal test data to
 taxonomy section.

---
 tool_conf.xml.sample            |   1 +
 tools/taxonomy/poisson2test.py  | 123 ++++++++++++++++++++++++++++++++
 tools/taxonomy/poisson2test.xml | 110 ++++++++++++++++++++++++++++
 3 files changed, 234 insertions(+)
 create mode 100644 tools/taxonomy/poisson2test.py
 create mode 100644 tools/taxonomy/poisson2test.xml

diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample
index 56f212dc2c2..1bc90ed044e 100644
--- a/tool_conf.xml.sample
+++ b/tool_conf.xml.sample
@@ -138,6 +138,7 @@
     <tool file="taxonomy/t2t_report.xml" />
     <tool file="taxonomy/t2ps_wrapper.xml" />
     <tool file="taxonomy/find_diag_hits.xml" />
+    <tool file="taxonomy/poisson2test.xml" />
   </section>
   <!--
   TODO: uncomment the following EMBOSS section whenever 
diff --git a/tools/taxonomy/poisson2test.py b/tools/taxonomy/poisson2test.py
new file mode 100644
index 00000000000..f3e64e9d0c1
--- /dev/null
+++ b/tools/taxonomy/poisson2test.py
@@ -0,0 +1,123 @@
+#!/usr/local/bin/python
+
+import sys
+from math import *
+from rpy import *
+
+
+if ((len(sys.argv)-1) != 5):
+    print 'too few parameters' 
+    print 'usage: inputfile, col1, col2, d-value(not 0), p-val correction method(0 or 1)'
+    sys.exit()
+    
+try:
+    lines_arr = open(sys.argv[1]).readlines()
+except IOError:
+    print'cannot open',sys.argv[1]
+    sys.exit()  
+ 
+try:
+    i = int(sys.argv[2]) #first column to compare
+    j = int(sys.argv[3]) #second colum to compare
+    d = float(sys.argv[4]) #correction factor
+    k = int(sys.argv[5]) #p-val correction method
+    
+    if (i>j):
+        print 'column order not correct col1 < col2'
+        print 'usage: inputfile, col1, col2, d-value, p-val correction method'
+        sys.exit()      
+        
+    try:
+        a = 1 / d
+        assert k in [0,1]
+    except ZeroDivisionError:
+        print 'd cannot be 0'
+        print 'usage: inputfile, col1, col2, d-value, p-val correction method'
+        sys.exit()
+    except:
+        print ' p-val correction should be 0 or 1 (0 = "bonferroni", 1 = "fdr")'
+        print 'usage: inputfile, col1, col2, d-value, p-val correction method'
+        sys.exit()
+except ValueError:
+    print 'parameters are not integers'
+    print 'usage: inputfile, col1, col2, d-value, p-val correction method'
+    sys.exit()
+   
+
+fsize = len(lines_arr)
+
+z1 = []
+z2 = []
+pz1 = []
+pz2 = []
+field = []
+
+if d<1: # Z score calculation
+    for line in lines_arr:
+        line.strip()
+        field = line.split('\t')
+        
+        x = int(field[j-1]) #input column 2
+        y = int(field[i-1]) #input column 1
+        if y>x:
+            z1.append(float((y - ((1/d)*x))/sqrt((1/d)*(x + y))))
+            z2.append(float((2*(sqrt(y+(3/8))-sqrt((1/d)*(x+(3/8)))))/sqrt(1+(1/d))))
+        else:
+            tmp_var1 = x
+            x = y
+            y = tmp_var1
+            z1.append(float((y - (d*x))/sqrt(d*(x + y))))
+            z2.append(float((2*(sqrt(y+(3/8))-sqrt(d*(x+(3/8)))))/sqrt(1+d)))
+            
+else: #d>1 Z score calculation
+    for line in lines_arr:
+        line.strip()
+        field = line.split('\t')
+        x = int(field[i-1]) #input column 1
+        y = int(field[j-1]) #input column 2
+        
+        if y>x:
+            z1.append(float((y - (d*x))/sqrt(d*(x + y))))
+            z2.append(float((2*(sqrt(y+(3/8))-sqrt(d*(x+(3/8)))))/sqrt(1+d)))
+        else:
+            tmp_var2 = x
+            x = y
+            y = tmp_var2
+            z1.append(float((y - ((1/d)*x))/sqrt((1/d)*(x + y))))
+            z2.append(float((2*(sqrt(y+(3/8))-sqrt((1/d)*(x+(3/8)))))/sqrt(1+(1/d))))
+        
+  
+   
+
+
+# P-value caluculation for z1 and z2
+for p in z1:
+    
+    pz1.append(float(r.pnorm(-abs(float(p)))))
+
+for q in z2:
+    
+    pz2.append(float(r.pnorm(-abs(float(q)))))    
+
+# P-value correction for pz1 and pz2
+
+if k == 0:
+    corrz1 = r.p_adjust(pz1,"bonferroni",fsize)
+    corrz2 = r.p_adjust(pz2,"bonferroni",fsize)
+  
+   
+else:
+  
+    corrz1 = r.p_adjust(pz1,"fdr",fsize)
+    corrz2 = r.p_adjust(pz2,"fdr",fsize)
+    
+
+#printing all columns
+for n in range(fsize):
+    print "%s\t%4.3f\t%4.3f\t%8.6f\t%8.6f\t%8.6f\t%8.6f" %(lines_arr[n].strip(),z1[n],z2[n],pz1[n],pz2[n],corrz1[n],corrz2[n])
+
+
+      
+      
+      
+          
diff --git a/tools/taxonomy/poisson2test.xml b/tools/taxonomy/poisson2test.xml
new file mode 100644
index 00000000000..2523e8e36ad
--- /dev/null
+++ b/tools/taxonomy/poisson2test.xml
@@ -0,0 +1,110 @@
+<tool id="poisson2test" name="Poisson 2 test" version="1.0.0">
+  <description>on tab-delimited file</description>
+  <command interpreter="python">poisson2test.py $input1 $input2 $input3 $input4 $input5 > $output1 </command>
+  <inputs>
+    <param name="input1" format="tabular" type="data" label="Input File"/>
+    <param name="input2" type="integer" size="5" value="2" label="First Column"/>
+    <param name="input3" type="integer" size="5" value="3" label="Second Column"/>
+    <param name="input4" type="float" size="5" value="1" label="D value"/>
+    <param name="input5" type="select" label="correction method">
+        <option value="0">Bonferroni</option>
+        <option value="1">FDR</option>
+    </param> 
+  </inputs>
+  <outputs>
+    <data format="tabular" name="output1" />
+  </outputs> 
+  <tests>
+    <test>
+        <param name="input1" value="poisson2test1.txt"/>
+        <param name="input2" value="2" />
+        <param name="input3" value="3" />
+        <param name="input4" value="0.44" />
+        <param name="input5" value="0" />
+        <output name="output1" file="poisson2test1.out" />    
+    </test>
+    <test>
+        <param name="input1" value="poisson2test2.txt"/>
+        <param name="input2" value="2" />
+        <param name="input3" value="3" />
+        <param name="input4" value="0.44" />
+        <param name="input5" value="0" />
+        <output name="output1" file="poisson2test2.out" />    
+    </test>    
+  </tests>
+  <help>
+
+**What it does**
+
+Suppose you have metagenomic samples from two different locations and have classified the reads unique to various taxa. Now you want to test if the number of reads that fall in a particular taxon in location 1 is different from those that fall in the same taxon in location 2. 
+This utility performs this analysis. It assumes that the data comes from a Poisson process and calculates two Z scores (Z1 and Z2) based on the work by Shiue and Bain; 1982 (Z1) and Huffman; 1984 (Z2).
+
+-----
+
+**Z score formula**
+
+Equation 1:
+
+.. image:: ../static/images/poisson2test_eqn1.png 
+
+ 
+Equation 2:
+
+.. image:: ../static/images/poisson2test_eqn2.png
+
+
+X = number of reads falling in a particular taxon in location 1
+ 
+Y = number of reads falling in the same taxon in location 2
+ 
+d = correction factor that accounts for biases in sample collection, DNA concentration, read numbers etc. between the two locations. 
+
+Not only that, this utility also provides corresponding p-values and corrected p-values (using Bonferroni or False Discovery Rate (FDR)). It takes in an input file (a tab delimited file consisting of three or more columns (taxon/category, read counts in location 1, read counts in location 2)), columns to compare, d value and a correction method 0 (Bonferroni) or 1 (FDR).
+
+-----
+
+**Example**
+
+- Input File: phylum, read count in location-1, read count in location-2::
+
+    Annelida            36     2
+    Apicomplexa         17     8
+    Arthropoda        1964   928
+    Ascomycota         436    49
+    Basidiomycota       77    55
+
+- Arguments to be supplied by the user::
+
+    col_i   col_j   d-value    correction-method
+    
+    2       3       0.44       Bonferroni
+
+- Output File: phylum, readcount1, readcount2, z1, z2, p1, p2, corrected p1, corrected p2::
+
+    Annelida            36     2   3.385   4.276  0.000356  0.000010  0.00463  0.00012
+    Apicomplexa         17     8  -0.157  -0.156  0.437707  0.438103  1.00000  1.00000
+    Arthropoda        1964   928  -1.790  -1.777  0.036755  0.037744  0.47782  0.49067
+    Ascomycota         436    49   9.778  11.418  0.000000  0.000000  0.00000  0.00000
+    Basidiomycota       77    55  -2.771  -2.659  0.002792  0.003916  0.03629  0.05091
+
+-----
+
+**Note**
+
+- Input file should be Tab delimited
+- i &lt; j
+- d cannot be 0
+- k = Bonferroni or FDR
+
+-----
+
+**References**
+
+- Shiue, W. and Bain, L. (1982). Experiment Size and Power Comparisons for Two-Sample Poisson Tests. Applied Statistics 31, 130-134.
+
+- Huffman, M. D. (1984). An Improved Approximate Two-Sample Poisson Test. Applied Statistics 33, 224-226.
+
+  </help>
+</tool>
+
+