NAME
GO::TermFinder - identify GO nodes that annotate a group of genes with a significant p-value
DESCRIPTION
This package is intended to provide a method whereby the P-values of a set of GO annotations can be determined for a set of genes, based on the number of genes that exist in the particular genome (or in a selected background distribution from the genome), and their annotation, and the frequency with which the GO nodes are annotated across the provided set of genes. The P-value is simply calculated using either the hypergeometric distribution (the binomial distribution is also provided as an option), as the probability of x or more out of n genes having a given annotation, given that G of N have that annotation in the genome in general. The hypergeometric distribution (sampling without replacement) is more accurate, though slower to calculate than the binomial distibution (sampling with replacement).
In addition, a corrected p-value can be calculated, to correct for multiple hypothesis testing. The correction factor used is the total number of nodes to which the provided list of genes are annotated, excepting any nodes which have only a single annotation in the background, as a priori, we know that these cannot be significantly enriched. The client has access to both the corrected and uncorrected values. It is also possible to correct the p-value using 1000 simulations, which control the Family Wise Error Rate - using this option suggests that the Bonferroni correction is in fact somewhat liberal, rather than conservative, as might be expected. Finally, the False Discovery Rate can also be calculated.
The general idea is that a list of genes may have been identified for some reason, e.g. they are coregulated, and TermFinder can be used to find out if any nodes annotate the set of genes to a level which is extremely improbable if the genes had simply been picked at random.
TODO
1. May want the client to decide the behaviour for ambiguous names, rather than having it hard coded (eg always ignore; use if standard name (current implementation); use all databaseIds for the ambiguous name; decide on a case by case basis (potentially useful if running on command line)).
2. Would probably be a good idea to test using the Math::BigInt module to deal with the factorials and nChooser calculations, so that they are more accurate, rather than the log versions that I currently have. The latest version of Math::BigInt can have a C implementation underneath, if used with the GMP library, that may even speed it up a little.
3. Create new GO::Hypothesis and GO::HypothesisSet objects, so that it is easier to access the information generated about the p-value etc. of any particular GO node that annotates a set of genes.
4. Instead of all the global variables, $k..., replace them with constants, which may improve runtime, as the optimizer should optimize the hash look ups to look like hard-coded strings at runtime, rather than variable lookups.
5. Create a GO::Math class, to extract out all the statistics from this module, so that only the logic of what math is needed to get the p-values, rather than the implementations themselves are left. This will decrease the complexity of this module considerably.
6. Lots of other stuff....
Instance Constructor
new
This is the constructor. It expects to be passed named arguments for an annotationProvider, and an ontologyProvider. In addition, it must be told the aspect of the ontology provider, so that it knows how to query the annotationProvider.
There are also some additional, optional arguments:
population:
This argument allows a client to indicate the population that should used to calculate a background distribution of GO terms. In the absence of population argument, then the background distribution will be drawn from all genes in the annotationProvider. This should be provided as an array reference, and no ambiguous names should be provided (see AnnotationProvider for details of name ambiguity). This option is particularly pertinent in a case where for example you assayed only 2000 genes in a two hybrid experiment, and found 20 interesting ones. To find significant terms, you need to do it in the context of the genes that you assayed, not in the context of all genes with annotation.
totalNumGenes:
This argument allows a client to indicate that the size of the background distribution is in fact larger that the number of genes that exist in the annotation provider, and the extra genes are merely assumed to be entirely unannotated.
NB: This is an API change, as totalNumGenes was previously required.
Thus - if using 'population', the total number of genes considered as the background will be the number of genes in the provided population. If not using 'population', then the number of genes that will be considered as the total population will be the number of genes in the annotationProvider. However, if the totalNumGenes argument is provided, then that number will be used as the size of the population. If it is not larger than the total number of genes in the annotationParser, then the number of genes in the annotationParser will be used. The totalNumGenes and the population arguments are mutually exclusive, and both should not be provided at the same time.
Usage ($num is larger than the number of genes with annotations):
my $termFinder = GO::TermFinder->new(annotationProvider=> $annotationProvider,
ontologyProvider => $ontologyProvider,
totalNumGenes => $num,
aspect => <P|C|F>);Usage (use all annotated genes as population):
my $termFinder = GO::TermFinder->new(annotationProvider=> $annotationProvider,
ontologyProvider => $ontologyProvider,
aspect => <P|C|F>);Usage (use a subset of genes as the background population):
my $termFinder = GO::TermFinder->new(annotationProvider=> $annotationProvider,
ontologyProvider => $ontologyProvider,
population => \@genes,
aspect => <P|C|F>);Instance Methods
findTerms
This method returns an array of hash references, one for each GO::Node that was tested as a hypothesis, that indicates which terms annotate the list of genes with what P-values. The contents of the hashes in the returned array depend on some of the run time options. They are:
key value
-------------------------------------------------------------------------Always Present:
NODE A GO::Node
PVALUE The P-value for having the observed number of
annotations that the provided list of genes
has to that node
NUM_ANNOTATIONS The number of genes within the provided list that
are annotated to the node.
TOTAL_NUM_ANNOTATIONS The number of genes across the genome
annotated to the node
ANNOTATED_GENES A hash reference, whose keys are the
databaseIds that are annotated to the node,
and whose values are the original name
supplied to the findTerms() method.Present if corrected p-values are calculated:
CORRECTED_PVALUE The CORRECTED_PVALUE is the PVALUE, but corrected
for multiple hypothesis testing, due to the
fact that you are more likely to generate
significant looking p-values if you test a
lot of hypotheses. See below for details of
how this pvalue is calculated, and the
options associated with it.Present if p-values were corrected by simulation:
NUM_OBSERVATIONS The number of simulations in which a p-value as
good as this one, or better, was observed.Present if the False Discovery Rate is calculated:
FDR_RATE The False Discovery Rate - this is a fraction
of how many of the nodes with p-values as good or
better than the node with this FDR would be expected
to be false positives.
FDR_OBSERVATIONS The average number of nodes during simulations
that had an uncorrected p-value as good or better
than the p-value of this node.
EXPECTED_FALSE_POSITIVES The expected number of false positives if this node
is chosen as the cut-off.The entries in the returned array are sorted by increasing p-value (i.e. least likely is first). If there is a tie in the p-value, then the sort order is determined by GOID, using a cmp comparison.
findTerm() expects to be passed, by reference, a list of gene names for which terms will be found. If a passed in name is ambiguous (see AnnotationProvider), then the following will occur:
1) If the name can be used as a standard name, it will assume that
it is that.
2) Otherwise it will not use it.Currently a warning will be printed to STDOUT in the case of an ambiguous name being used.
The passed in gene names are converted into a list of databaseIds. If a gene does not map to a databaseId, then an undef is put in the list - however, if the same gene name, which does not map to a databaseId, is used twice then it will produce only one undef in the list. If more than one gene name maps to the same databaseId (either because you used the same name twice, or you used an alias as well), then that databaseId is only put into the list once, and a warning is printed.
If a gene name does not have any information returned from the AnnotationProvider, then it is assumed that the gene is entirely unannotated. For these purposes, TermFinder annotates such genes to the root node (Gene_Ontology), its immediate child (which indicates the aspect of the ontology (such as biological_process), and a dummy go node, corresponding to unannotated. This node will have a goid of 'GO:XXXXXXX', and a term name of 'unannotated'. No other information will be set up for this GO::Node, so you should not count on being able to retrieve it. What it does mean is that you can determine if the predominant feature of a set of genes is that they have no annotation.
If more genes are provided that have been indicated exist in the genome (as provided during object construction), then an error message will be printed out, and an empty list will be returned.
In addition, it is possible that for a small list of genes, that no hypotheses will be tested - in this case, those genes will only have annotated nodes with a count of 1, other than the Gene_Ontology node itself, and the node corresponding to the aspect of the ontology. Neither of these are considered for p-value testing, as a priori they must have a p-value of 1.
MULTIPLE HYPOTHESIS CORRECTION
An optional argument, 'correction' may be used, which indicates what method of multiple hypothesis correction should be used. Multiple hypothesis correction attempts to keep the overall chance of getting any false positives at the same level (e.g. 0.05). Acceptable values are:
bonferroni, none, simulation
: 'bonferroni' will correct the p-values by using as the correction
factor the total number of nodes to which the provided list of
genes are annotated, either directly or indirectly, excepting any
nodes that are annotated only once in the background distribution,
as, a priori, these cannot be overrepresented.
: 'none' will perform no multiple hypothesis correction
: 'simulation' will run 1000 simulations with random lists of genes
(the same size as the originally provided gene list), and determine
a corrected value by how many simulations produced a p-value better
than the p-value associated with one of the real hypotheses.
E.g. if a node from the real data has a p-value of 0.05, but a
p-value that good or better is generated in 500 out of 1000 trials,
the corrected pvalue will be 0.5. In the case that a p-value
generated from a real list of genes is never seen in the
simulations, it will be given a corrected p-value of < 0.001, and
the NUM_OBSERVATIONS attribute of the hypothesis will be 0. Using
this option takes 1000 time as long!The default for this argument, if not provided, is bonferroni.
FALSE DISCOVERY RATE
As a way of pre-empting the potential problems of using p-values corrected for multiple hypothesis testing, the False Discovery Rate can instead be calculated, and you can instead set your cutoff based on an acceptable false discovery rate, such as 0.01 (1%), or 0.05 (5%) etc. Thus, the optional argument 'calculateFDR' can be used. A non-zero value means that the False Discovery Rate will be calculated for each node, such that you can determine, if you chose your p-value cut-off at that node, what the FDR would be. The FDR is calculated by running 50 simulations, and counting the average number of times a p-value as good or better that a p-value generated from the real data is seen. This is used as the numerator. The denominator is the number of p-values in the real data that are as good or better than it.
Usage example - in this example, the default (Bonferroni) correction is used to calculate a corrected p-value, and in addition, the False Discovery Rate is also calculated:
my @pvalueStructures = $termFinder->findTerms(genes => \@genes,
calculateFDR => 1);
my $hypothesis = 1;
foreach my $pvalue (@pvalueStructures){
print "-- $hypothesis of ", scalar @pvalueStructures, "--\n",
"GOID\t", $pvalue->{NODE}->goid, "\n",
"TERM\t", $pvalue->{NODE}->term, "\n",
"P-VALUE\t", $pvalue->{PVALUE}, "\n",
"CORRECTED P-VALUE\t", $pvalue->{CORRECTED_PVALUE}, "\n",
"FALSE DISCOVERY RATE\t", $pvalue->{FDR_RATE}, "\n",
"NUM_ANNOTATIONS\t", $pvalue->{NUM_ANNOTATIONS}, " (of ", $pvalue->{TOTAL_NUM_ANNOTATIONS}, ")\n",
"ANNOTATED_GENES\t", join(", ", values (%{$pvalue->{ANNOTATED_GENES}})), "\n\n";
$hypothesis++;
}aspect
Returns the aspect with the the GO::TermFinder object was constructed.
Usage:
my $aspect = $termFinder->aspect;Authors
Gavin Sherlock; sherlock@genome.stanford.edu
Elizabeth Boyle; ell@mit.edu