Mercurial > repos > melissacline > ucsc_cancer_browser_stats

changeset 12:fd8529cd1564 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
better t-test
author jingchunzhu
date Mon, 28 Sep 2015 12:36:12 -0700
parents cd4c13ae11ce
children
files ttest/stats.py ttest/stats.xml
diffstat 2 files changed, 34 insertions(+), 59 deletions(-) [+]
line wrap: on
line diff
--- a/ttest/stats.py	Fri Jul 31 13:32:25 2015 -0700
+++ b/ttest/stats.py	Mon Sep 28 12:36:12 2015 -0700
@@ -5,7 +5,7 @@
 import argparse
 import math
 import re
-import sys
+import sys,os
 import numpy
 import scipy
 import scipy.stats
@@ -26,6 +26,7 @@
                                                         args.a, args.b)
     outfile = open(args.outfile, "w")
     outfile.write("%s\tStatistic\tpValue\tMedian1\tMedian2\tDelta\n" % (headerTokens[0]))
+    tempOut = open(args.outfile+"_tmp", "w")
     for row in genomicFp:
         tokens = row.rstrip().split("\t")
         if len(tokens)!= len(headerTokens):
@@ -54,8 +55,11 @@
             outputLine = outputLine + "\tNA"
         else:
             outputLine = "%s\t%f" % (outputLine, medianA - medianB)
-        outfile.write("%s\n" % (outputLine))
+        tempOut.write("%s\n" % (outputLine))
+    tempOut.close()
     outfile.close()
+    os.system("sort -nk 2 "+ args.outfile+"_tmp"+" >> "+args.outfile )
+    os.system("rm "+ args.outfile+"_tmp")
 
 def readClinicalFile(clinicalFile):
     clinicalVector = dict()
--- a/ttest/stats.xml	Fri Jul 31 13:32:25 2015 -0700
+++ b/ttest/stats.xml	Mon Sep 28 12:36:12 2015 -0700
@@ -1,4 +1,4 @@
-<tool id="ucscCancerBrowserStats" description="t-tests of Difference in genomic data" name="Difference between categories (t-test)" version="0.0.1">
+<tool id="ucscCancerBrowserStats" description="t-tests of difference in genomic data" name="Difference between categories (t-test)" version="0.0.1">
   <command interpreter="python">
     stats.py  $genomicMatrix $clinicalFeatures $outFile -a="${category1}" -b="${category2}"
   </command>
@@ -23,76 +23,47 @@
   </tests>
   <help>
 
-This tool performs statistical tests found in the UCSC Cancer Genomics
-Browser.  The input data is a genomic matrix (containing genomic data,
-with rows representing genes or probes and columns representing
-samples or patients), a clinical matrix of two (or more) columns
-assigning categorical values to the samples, and two categorical
-values of interest.  The tool identifies the samples corresponding to
-each categorical value, then identifies the columns in the genomic
-matrix corresponding to those sets of samples, which identifies two
-groups of columns.  For each row in the genomic matrix, it extracts
-the value for those two sets of columns, performs a t-test on the two
-sets of values, and returns the result for the row.  Any values for
-any columns NOT pertaining to one of the categorical values of
-interest are ignored.
+This tool performs t-test on genomic data between two groups of samples, which can be used to identify for example, differentially expressed genes or probes.  The genomic data is in the format of UCSC Xena genomic matrix (a tab-deliminated matrix) with rows representing genes or probes and columns representing samples. The phenotype matrix assigns samples into groups. The tool compares two groups of samples, and computes the t-statistics, p value, and delta of medians for each probe/gene between the two groups. The result can be downloaded to programs such as EXCEL for sorting based on the t-statistics. 
+
+The user runs this tool with the following steps:
 
-The user runs this tool with th following steps:
+1. Specify a genomic matrix.  The expected format is with rows representing genes and columns representing samples, and the first line contains sample names. Matrix can be obtained from UCSC Xena bulk download. See below for an example.
+
+
+2. Specify a phenotype matrix.  Here, rows indicate samples, columns indicate phenotypes or annotations.  Matrix can be obtained from UCSC Xena heatmap download. See below for an example.
 
 
-1. Specify a genomic matrix.  The expected format is with rows representing 
-genes and columns representing samples, and the first line contains sample 
-names.
-
-2. Specify a clinical matrix.  Here, rows indicate samples, columns
-indicate clinical features, and the header row contains feature names.
-The first column MUST indicate the sample names, and MUST correspond
-to the column names of the genomic matrix.  The clinical feature of
-interest MUST be in the second column.  Any other columns will be
-ignored.
+3. Specify the two categorical values that you want to use for defining the two groups.  For example, the two groups could be A and B, 0 and 1, etc.
 
 
-3. Indicate two clinical values that you want to use for defining the
-two groups.  For example, the two groups could be "Red group" and
-"Green group", 0 and 1, or whatever.
+4. The output is, for each probe/gene (in each row), the t-statistics, the p-value, the median value for each group, and the difference between the medians.  If it cannot calculate these values, it returns a vector of NAs.
+
 
-The output indicates, for each row, the t-statistic reporting on the
-difference between the two groups of columns (as specified by the two
-clinical values), the p-value corresponding to that t-statistic, the
-median value for each group, and the difference between the medians.  If it 
-cannot calculate these values, it returns a vector of NAs.
+**Input genomic matrix**::
 
-For example, given the following genomic matrix for (1)::
-
-    Gene  1    2    3    4    5    6    7    8    9    10
+    Gene  s1   s2   s3   s4   s5   s6   s7   s8   s9   s10
     G1    2.0  2.2  3.2  1.1  5.1  8.1  3.2  1.1  8.1  0.2
     G2    0.1  8.2  9.1  4.2  6.1  4.9  3.9  2.3  1.1  0.2
 
-and given the following clinical matrix for (2)::
+**Input phenotyp matrix**::
 
-    sample_id Value
-    1         A
-    2         A
-    3         B
-    4         C
-    5         B
-    6         B
-    7         A
-    8         A
-    9         B
-    10        A
+    sample_id  Value
+    s1         A
+    s2         A
+    s3         B
+    s4         C
+    s5         B
+    s6         B
+    s7         A
+    s8         A
+    s9         B
+    s10        A
     
-and given A for Category 1 and B for Category 2
-
-the tool will assemble the following two groups of values::
+**Category 1 : A**
 
-    G1 A:(2.0, 2.2, 3.2, 1.1, 0.2) B:(3.2, 5.1, 8.1, 8.1)
-    G2 A:(0.1, 8.2, 3.9, 2.3, 0.2) B:(9.1, 6.1, 4.9, 1.1) 
+**Category 2 : B**
 
-Note that the values for sample_id 4 do not appear, because it has a Value
-of C in the second column, which is neither A nor B.
-
-And it will return the output::
+**Output**::
 
     Gene Statistic  pValue    Median1   Median2   Delta
     G1   -4.168999  0.004194  2.000000  6.600000  -4.600000