Are you sick of bugs? Me too!
An easy way to search down bugs is to test your Python code using doctests.
For example, if you have a module like this called 'module.py':
#====================================================================#
# define a function to read GO terms from a line of a file:
def read_go_terms_from_line_of_file(line):
"""read GO terms from a line of a file
>>> read_go_terms_from_line_of_file("791228 GO:0005515")
('791228', ['GO:0005515'])
>>> read_go_terms_from_line_of_file("530242 GO:0000003,GO:0004675,GO:0005024,GO:0016020")
('530242', ['GO:0000003', 'GO:0004675', 'GO:0005024', 'GO:0016020'])
>>> read_go_terms_from_line_of_file("794214 none")
('794214', [])
"""
# eg.
# 791228 GO:0005515
# 5470 GO:0002119,GO:0004652,GO:0006952,GO:0007423,GO:0009792,GO:0010171,GO:0019915,GO:0040002,GO:0040011,GO:0040018,GO:0042302,GO:0042338,GO:0043631,GO:0045138
# 530242 GO:0000003,GO:0004675,GO:0005024,GO:0016020
# 794214 none
temp = line.split()
family = temp[0]
go_terms = [] if temp[1] == "none" else temp[1].split(',')
return(family, go_terms)
#====================================================================#
if __name__ == "__main__":
import doctest
doctest.testmod()
This module has just one function 'read_go_terms_from_line_of_file' that reads the name of a gene family, and GO terms for the genes in that family, from an input file. The function takes as its input the line of a file, and returns the name of the family (a string) and a list of the GO terms. At the top of the function (in red above) are some doctests, which define what the function should return (according to the programmer) for particular example inputs. The end of the module has some code (in pink above) that says that if we run this module, by typing:
% python3 module.py
then the doctests are run. If you get no output when you run 'python3 module.py', then it means that the doctests were run, and all ran fine, and no bugs were found (hurray!) If you do get some output, then it helps you find those nasty bugs, and you can hunt them down and squash them (hurray!)
To get more detailed output regarding which tests were run and which failed/passed, type:
% python3 module.py -v
If you want to include doctests in a program (not a module)
An alternative way is to include doctests in your program, for example, your program could have the same code as above, except that the end of the program would have a 'main' function that calls the functions, and the 'doctest' lines would be removed :
def main():
read_go_terms_from_line_of_file("791228 GO:0005515")
if __name__=="__main__":
main()
Then, to run the doctests we just run the program with the "-m doctest" option:
% python -m doctest program.py
If your function returns a dictionary
Dictionaries are unsorted, so this is a problem. But you can do something like this:
>>> find_midpoint_positions({'scaff1': ['50=100','150=200','250=300'], 'scaff2': ['100=200', '202=500'], 'scaff3': ['300=500']}) == {'scaff1': [125, 225], 'scaff2': [201]}
True
Monday, 27 July 2015
Monday, 13 July 2015
Parsing a phylogenetic tree in Python using ETE2
I've just started learning about ETE2 for parsing and plotting phylogenetic trees in Python, on recommendation from my colleague Diogo Ribeiro (thanks Diogo!)
ETE2 is currently in Python2, but apparently a Python3 version is coming soon (I hope so!)
For example, here's a simple script to read in a tree, and traverse all the nodes in a tree, and for each node print its children, its parent, and all the leaves under (decendants to) the node:
import os
import sys
from ete2 import Tree
# note: this script is in python2, as ete2 is in python2
#====================================================================#
def main():
# check the command-line arguments:
if len(sys.argv) != 2 or os.path.exists(sys.argv[1]) == False:
print("Usage: %s newick_file" % sys.argv[0])
sys.exit(1)
newick_file = sys.argv[1]
# the code below is based on examples on http://pythonhosted.org/ete2/tutorial/tutorial_trees.html
# load a tree from a file. The root node is put into 't'
t = Tree(newick_file, format=1) # format=1 means internal node names are loaded
# iterate over all the nodes in the tree, and get the children of each node:
for node in t.traverse("postorder"):
# get the children of this node:
children = node.get_children()
for child in children:
print "node %s has child %s" % (node.name, child.name)
# iterate over all the nodes in the tree, and get the parent of each node:
for node in t.traverse("postorder"):
# get the parent of this node:
if not node.is_root(): # if it is not the root node
print "node %s has parent %s" % (node.name, node.up.name)
# iterate over all the nodes in the tree, and get the descendant species of each node:
for node in t.traverse("postorder"):
if not node.is_leaf(): # if it is not a leaf
# get the descendant species (leaves):
leaves = node.get_leaves()
leaf_names = [leaf.name for leaf in leaves]
print "node %s has descendants %s" % (node.name, ','.join(leaf_names))
# iterate over all the nodes in the tree, and get the ancestor nodes of each node, back
# to the root:
for node in t.traverse("postorder"):
# get the ancestors of this node:
if not node.is_root():
ancestor = node.up
ancestors = [ancestor]
while not ancestor.is_root():
ancestor = ancestor.up
ancestors.append(ancestor)
ancestor_names = [ancestor.name for ancestor in ancestors]
print "node %s has ancestors %s" % (node.name, ','.join(ancestor_names))
print("FINISHED")
#====================================================================#
if __name__=="__main__":
main()
#====================================================================#
ETE2 is currently in Python2, but apparently a Python3 version is coming soon (I hope so!)
For example, here's a simple script to read in a tree, and traverse all the nodes in a tree, and for each node print its children, its parent, and all the leaves under (decendants to) the node:
import os
import sys
from ete2 import Tree
# note: this script is in python2, as ete2 is in python2
#====================================================================#
def main():
# check the command-line arguments:
if len(sys.argv) != 2 or os.path.exists(sys.argv[1]) == False:
print("Usage: %s newick_file" % sys.argv[0])
sys.exit(1)
newick_file = sys.argv[1]
# the code below is based on examples on http://pythonhosted.org/ete2/tutorial/tutorial_trees.html
# load a tree from a file. The root node is put into 't'
t = Tree(newick_file, format=1) # format=1 means internal node names are loaded
# iterate over all the nodes in the tree, and get the children of each node:
for node in t.traverse("postorder"):
# get the children of this node:
children = node.get_children()
for child in children:
print "node %s has child %s" % (node.name, child.name)
# iterate over all the nodes in the tree, and get the parent of each node:
for node in t.traverse("postorder"):
# get the parent of this node:
if not node.is_root(): # if it is not the root node
print "node %s has parent %s" % (node.name, node.up.name)
# iterate over all the nodes in the tree, and get the descendant species of each node:
for node in t.traverse("postorder"):
if not node.is_leaf(): # if it is not a leaf
# get the descendant species (leaves):
leaves = node.get_leaves()
leaf_names = [leaf.name for leaf in leaves]
print "node %s has descendants %s" % (node.name, ','.join(leaf_names))
# iterate over all the nodes in the tree, and get the ancestor nodes of each node, back
# to the root:
for node in t.traverse("postorder"):
# get the ancestors of this node:
if not node.is_root():
ancestor = node.up
ancestors = [ancestor]
while not ancestor.is_root():
ancestor = ancestor.up
ancestors.append(ancestor)
ancestor_names = [ancestor.name for ancestor in ancestors]
print "node %s has ancestors %s" % (node.name, ','.join(ancestor_names))
print("FINISHED")
#====================================================================#
if __name__=="__main__":
main()
#====================================================================#
Python list comprehensions
List comprehensions are a really nice feature of Python. For example, say you have a list 'leaves' of objects, and each object in the list has an attribute 'name' (a string). To print out a comma-separated string of all the names of all the objects in the list, we just type:
# first get all the names of all the 'leaf' objects in the list 'leaves":
leaf_names = [leaf.name for leaf in leaves]
# now print out a comma-separated string of all the names of all the leaves:
print "The names are %s" % (','.join(leaf_names))
Other examples:
You have a dictionary 'mydict' of lists, get those lists and flatten them into one big list:
The values (not keys) of 'mydict' are lists.
biglist = [x for sublist in mydict.values() for x in sublist]
# mydict.values() gives a list of lists, this flattens it into one big list
You have a dictionary 'fasta_seq' whose keys are protein names, and values are protein sequences, and want the lengths of the proteins:
This gives a list with the lengths of the proteins:
lens = [len(fasta_seq[x]) for x in list(fasta_seq.keys())]
You have a list of words, called 'temp', and want to find the words in that list that do not contain 'status:', 'UniProt:', 'protein_id:', or 'locus:'
matching = [x for x in temp if 'status:' not in x and 'UniProt:' not in x and 'protein_id:' not in x and 'locus:' not in x]
You have a list 'temp' of words, and only want those words that also appear in a second list 'temp2':
seqmatch_words = [x for x in temp if x in temp2]
You have a list 'stack' of tuples, and want the first member of each tuple:
Your list is called 'stack'. This puts the first members of the tuples in a list of their own:
stack_members = [x[0] for x in stack]
You want to create a list containing 9 lists initialised to zero:
Matrix = [[0 for x in xrange(9)] for x in xrange(9)]
You want to create a dictionary with the length of each item in a list of strings:
genes = ['gene1', 'gene2', 'gene32']
mydict = {gene: len(gene) for gene in genes}
You want to create a dictionary, taking the keys from one list and the values from a second list:
genes = ['gene1', 'gene2', 'gene32']
values = [32, 22, 66]
mydict = {genes[i]: values[i] for i in range(len(genes))}
You want to create a set, taking the items in a list of strings that start with 's':
mylist = ['string1', 'hello', 'sausage']
myset = {x for x in mylist if x.startswith('s')}
You have a list of sets, and a dictionary with values for the items in the sets, and want to create a list of sublists, where each sublist have the dictionary values for the objects in a set:
mylist = [{'seq1', 'seq2'}, {'seq3', 'seq4', 'seq5'}, {'seq6'}]
mydict = {'seq1': 5, 'seq2': '4', 'seq3': '4', 'seq4': '3', 'seq5': 2, 'seq6': 2}
mybiglist = [[mydict[x] for x in sublist] for sublist in mylist]
You have two lists of equal lengths and want to merge them into a single list of tuples:
list1 = ['A', 'B', 'C', 'D']
list2 = ['a', 'b', 'c', 'd']
[(list1[i], list2[i]) for i in range(len(list1))]
[('A', 'a'), ('B', 'b'), ('C', 'c'), ('D', 'd')]
# first get all the names of all the 'leaf' objects in the list 'leaves":
leaf_names = [leaf.name for leaf in leaves]
# now print out a comma-separated string of all the names of all the leaves:
print "The names are %s" % (','.join(leaf_names))
Other examples:
You have a dictionary 'mydict' of lists, get those lists and flatten them into one big list:
The values (not keys) of 'mydict' are lists.
biglist = [x for sublist in mydict.values() for x in sublist]
# mydict.values() gives a list of lists, this flattens it into one big list
You have a dictionary 'fasta_seq' whose keys are protein names, and values are protein sequences, and want the lengths of the proteins:
This gives a list with the lengths of the proteins:
lens = [len(fasta_seq[x]) for x in list(fasta_seq.keys())]
You have a list of words, called 'temp', and want to find the words in that list that do not contain 'status:', 'UniProt:', 'protein_id:', or 'locus:'
matching = [x for x in temp if 'status:' not in x and 'UniProt:' not in x and 'protein_id:' not in x and 'locus:' not in x]
You have a list 'temp' of words, and only want those words that also appear in a second list 'temp2':
seqmatch_words = [x for x in temp if x in temp2]
You have a list 'stack' of tuples, and want the first member of each tuple:
Your list is called 'stack'. This puts the first members of the tuples in a list of their own:
stack_members = [x[0] for x in stack]
You want to create a list containing 9 lists initialised to zero:
Matrix = [[0 for x in xrange(9)] for x in xrange(9)]
You want to create a dictionary with the length of each item in a list of strings:
genes = ['gene1', 'gene2', 'gene32']
mydict = {gene: len(gene) for gene in genes}
You want to create a dictionary, taking the keys from one list and the values from a second list:
genes = ['gene1', 'gene2', 'gene32']
values = [32, 22, 66]
mydict = {genes[i]: values[i] for i in range(len(genes))}
You want to create a set, taking the items in a list of strings that start with 's':
mylist = ['string1', 'hello', 'sausage']
myset = {x for x in mylist if x.startswith('s')}
You have a list of sets, and a dictionary with values for the items in the sets, and want to create a list of sublists, where each sublist have the dictionary values for the objects in a set:
mylist = [{'seq1', 'seq2'}, {'seq3', 'seq4', 'seq5'}, {'seq6'}]
mydict = {'seq1': 5, 'seq2': '4', 'seq3': '4', 'seq4': '3', 'seq5': 2, 'seq6': 2}
mybiglist = [[mydict[x] for x in sublist] for sublist in mylist]
You have two lists of equal lengths and want to merge them into a single list of tuples:
list1 = ['A', 'B', 'C', 'D']
list2 = ['a', 'b', 'c', 'd']
[(list1[i], list2[i]) for i in range(len(list1))]
[('A', 'a'), ('B', 'b'), ('C', 'c'), ('D', 'd')]
Making a quick plot of a phylogenetic tree
I've been trying to find nice ways to make a quick plot of a phylogenetic tree.
Here are some nice ones:
- My colleague Bhavana told me about using FigTree [note to self: this is installed on the Sanger farm]
- On the ETE2 webpage, there is a very nice TreeView tool for drawing a tree, where you can just paste in your tree, and it makes a tree picture
- On the Sanger farm, if you logged in using 'ssh -X', you can use the ETE2 command-line tool, by typing 'ete2' with a Newick tree file name, and a GUI pops up
- You can also plot a tree using a Python script that calls ETE2, for example, this script calls ETE2 and pops up a picture of the tree in a GUI (picture below): (again, for this you need to log into the Sanger farm using 'ssh -X')
import os
import sys
from ete2 import Tree
# note: this script is in python2, as ete2 is in python2
#====================================================================#
def main():
# check the command-line arguments:
if len(sys.argv) != 2 or os.path.exists(sys.argv[1]) == False:
print("Usage: %s newick_file" % sys.argv[0])
sys.exit(1)
newick_file = sys.argv[1]
# the code below is based on examples on http://pythonhosted.org/ete2/tutorial/tutorial_drawing.html
# load a tree from a file. The root node is put into 't'
t = Tree(newick_file, format=1) # format=1 means internal node names are loaded
t.show()
#====================================================================#
if __name__=="__main__":
main()
#====================================================================#
Here are some nice ones:
- My colleague Bhavana told me about using FigTree [note to self: this is installed on the Sanger farm]
- On the ETE2 webpage, there is a very nice TreeView tool for drawing a tree, where you can just paste in your tree, and it makes a tree picture
- On the Sanger farm, if you logged in using 'ssh -X', you can use the ETE2 command-line tool, by typing 'ete2' with a Newick tree file name, and a GUI pops up
- You can also plot a tree using a Python script that calls ETE2, for example, this script calls ETE2 and pops up a picture of the tree in a GUI (picture below): (again, for this you need to log into the Sanger farm using 'ssh -X')
import os
import sys
from ete2 import Tree
# note: this script is in python2, as ete2 is in python2
#====================================================================#
def main():
# check the command-line arguments:
if len(sys.argv) != 2 or os.path.exists(sys.argv[1]) == False:
print("Usage: %s newick_file" % sys.argv[0])
sys.exit(1)
newick_file = sys.argv[1]
# the code below is based on examples on http://pythonhosted.org/ete2/tutorial/tutorial_drawing.html
# load a tree from a file. The root node is put into 't'
t = Tree(newick_file, format=1) # format=1 means internal node names are loaded
t.show()
#====================================================================#
if __name__=="__main__":
main()
#====================================================================#
Other tools:
- My colleague Diogo also mentioned a new tool called PhyloCanvas, which he said is very fast.
- My colleage Sam Dougan mentioned Taxonium
Note to self (sorry, this is just for me!):
[To get the Newick version of a 50HG tree, for family 141670:
% get_tree_for_50HG_family_v75.pl 141670 es9_compara_homology_final2_75 newick ]
Thursday, 2 July 2015
Using TopGO to test for GO-term enrichment
I wanted to test whether any GO terms are enriched in my gene list compared to all genes in the genome, so decided to try topGO, since my colleagues recommend it.
How TopGO works
The details of how TopGO are published in a paper by Alexa et al (2006). TopGO takes the GO hierarchy into account when calculating enrichment, and it is thought that this will lead it to report fewer false positive results than methods that do not take the GO hierarchy into account (eg. just using a Fisher's exact test). Grossman et al (2007) et al say that if you don't take the GO hierarchy into account, you get what they call the 'inheritance problem': because the highest (more general) GO terms inherit annotations from more specific descendant terms, and this can lead to false positives.
There are two methods described in the Alexa et al (2006) paper:
(i) 'elim' : this method processes the GO terms by traversing the GO hierarchy from bottom to top, ie. it first assesses the most specific (bottom-most) GO terms, and proceeds later to more general (higher) GO terms. When it assesses a higher (more general) GO term, it discards any genes that are annotated with significantly enriched descendant GO terms (considered significant using a pre-defined P-value threshold). This method does tend to miss some true positives at higher (more general) levels of the GO hierarchy.
(ii) 'weight': this is designed to miss less true positives than the 'elim' method. The significance scores of connected nodes (a parent and its child in the GO hierarchy) are compared in order to detect the locally most significant terms in the GO hierarchy.
The current TopGO actually has four more methods:
(iii) 'classic': each GO term is tested independently, not taking the GO hierarchy into account.
(iv) 'weight01': this is the default method used by TopGO, and is a mixture of the 'elim' and 'weight' methods.
(v) 'lea': this doesn't seem to be described in the TopGO vignette or any paper, so I'm not sure what it does (?)
(vi) 'parentchild': described by Grossman et al (2007) : when assessing a GO term, it takes into accoount the annotation of terms to the current term's parents, and so reduces false positives due to the inheritance problem. You can also use this algorithm in a very nice tool called Ontologizer , which has a lovely GUI (I used it in a collaboration on gene duplicates with Shane Woods, described in his paper Woods et al (2013)).
The 'elim' and 'weight' methods are designed to be more conservative (have less false positives) than the 'classic' method, so should have larger p-values in general. Simulation results reported by Alexa et al (2006) show that the 'weight' algorithm has less false positives than 'classic' and misses few true positives, while 'elim' has even less false positives than 'weight' but misses more true positives. Grossman et al (2007) compare the 'elim', 'weight' and 'parentchild' algorithms, and say that each has advantages in certain scenarios.
How to run TopGO
The vignette described in detail how to use the package, but here's a summary of what I did (so I won't forget!):
R packages
It's important to have the latest version of the packages, to make sure you have bug fixes, but also because the gene ontology is read into TopGO from another R package (GO.db), so you need to make sure you have the latest version of that package.
First you need to make sure that you have the latest version of R.
Then make sure you have the latest version of Bioconductor, and the latest version of the Bioconductor packages topGO:
>
Then check you have the latest versions of everything: (for me it should be R 3.2.0, Bioconductor 3.1 and TopGO 2.2.0 and GO.db (comes with TopGO) 3.1.2:
> sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] topGO_2.20.0 SparseM_1.6 GO.db_3.1.2
[4] RSQLite_1.0.0 DBI_0.3.1 AnnotationDbi_1.30.1
[7] GenomeInfoDb_1.4.1 IRanges_2.2.5 S4Vectors_0.6.1
[10] Biobase_2.28.0 graph_1.46.0 BiocGenerics_0.14.0
[13] BiocInstaller_1.18.3
loaded via a namespace (and not attached):
[1] lattice_0.20-31 grid_3.2.0 tools_3.2.0
Starting topGO
Once you have installed topGO you can start it by typing:
> library("topGO")
Loading your data
Reading in GO annotations for the genes
The file that you provide of gene-to-GO mappings should look be a tab-delimited file looking like this, where the first column is the gene name and the second column contains a list of GO terms (separated by a comma and space):
105778 GO:0016021, GO:0016020
166881 GO:0003674, GO:0016021, GO:0016020, GO:0008150
188691 GO:0016021, GO:0016020, GO:0001584, GO:0007165, GO:0004984, GO:0007186, GO:0004872
184033 GO:0005737, GO:0005515
For example, if your file is called 'annotations.txt', you can read it in using:
> geneID2GO <- readMappings(file = "annotations.txt")
To have the fastest run-time, for a particular gene, it shouldn't be the case that one of the GO terms is an ancestor of one of the other GO terms for that gene. However, TopGO will run fine if this is the case, although it will be slower. In practice, I found that I didn't need to prune my lists of GO terms like this, as TopGO was pretty fast for my data set.
GO hierarchy
The GO hierarchy used by topGO is the one that comes with the GO.db package. Therefore the user does not need to provide this.
Defining your list of genes of interest, and the 'gene universe' to compare it to
You need to tell TopGO what is your list of genes of interest, and what is the list of all genes in the 'gene universe' that you want to compare your list of interest to.
The 'gene universe' can be all the genes that contain GO annotations in your input 'annotations.txt' file (see above):
> geneUniverse <- names(geneID2GO)
We can read in a list of interesting genes from a file, just containing a list of interesting genes, one per line:
> genesOfInterest <- read.table("interestinggenes.txt",header=FALSE)
> genesOfInterest <- as.character(genesOfInterest$V1)
Then we need to tell TopGO where these interesting genes appear in the 'geneUniverse' vector:
> geneList <- factor(as.integer(geneUniverse %in% genesOfInterest))
> names(geneList) <- geneUniverse
The geneList object tells TopGO which genes in the gene universe are your genes of interest.
Putting the data together into an R object
We need to put our data in an object of type 'topGOdata'. This will contain the list of genes of interest, the GO annotations, and the GO hierarchy itself.
> myGOdata <- new("topGOdata", description="My project", ontology="BP", allGenes=geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
The 'ontology' argument can be 'BP' (biological process), 'MF' (molecular function), or 'CC' (cellular component).
The 'description' argument has a short description of your project.
The 'allGenes' argument specifies all the genes in the 'gene universe', and which of those are your genes of interest.
The 'annot' argument tells topGO how to map genes to GO annotations. 'annot'=annFUN.gene2GO means that the user provides gene-to-GO annotations, and we specify here that they are in object 'geneID2GO' (see above for how this was created).
An optional argument is 'nodeSize': this is used to prune the GO hierarchy, eg. nodesize=10 prunes the GO hierarchy, to remove terms which have less than 10 annotated genes.
TopGO generally only seems to use a subset of my 'gene universe' data, which it considers 'feasible'. I think this is because you may have some genes in your input 'gene universe' that don't have any GO terms annotated to them. Also, you might have some terms in the GO annotations for your gene universe that don't exist in the version of the GO hierarchy that TopGO uses (which comes from the GO.db R package). Unfortunately this means that TopGO actually only uses some of your input data, and throws some away. Note that also, you may have chosen to use 'BP' (biological process) terms but many of your input annotations may be molecular function (MF) or cellular component (CC) terms.
If you now type 'myGOdata', you get a summary of the data in it:
> myGOdata
Description:
- My project
Ontology:
- BP
4 available genes (all genes from the array):
- symbol: 105778 166881 188691 184033 ...
- 2 significant genes.
2 feasible genes (genes that can be used in the analysis):
- symbol: 166881 188691 ...
- 1 significant genes.
GO graph (nodes with at least 1 genes):
- a graph with directed edges
- number of nodes = 15
- number of edges = 22
Here 2 of the 4 genes in the 'gene universe' have been mapped to the GO ontology (there are 2 'feasible' genes), so the rest of the analysis will be based on just these genes. TopGO calls the genes of interest 'significant genes'. There are 2 genes of interest in the full gene universe of 4 genes, but just one gene of interest among the 2 feasible genes. The 'myGOdata' object takes some time to build, especially if you have a large gene universe. I found it pretty fast though, less than one minute for a gene universe of ~100,000 genes. I guess the run-time might depend on how much memory (RAM) you have available(?)
The list of genes of interest can be accessed using the method sigGenes():
> sg <- sigGenes(myGOdata)
> str(sg)
> numSigGenes(myGOdata)
Performing enrichment tests
One type of test that topGO does is a Fisher's exact test based on gene counts:
> resultFisher <- runTest(myGOdata, algorithm="classic", statistic="fisher")
Selecting 'algorithm=classic' means that the GO hierarchy isn't taken into account, so each GO term is tested independently. (Usually in fact we want to take the GO hierarchy into account, so would use algorithm='weight01' for example.)
You can type 'resultFisher' to get a summary of the results:
Description: My project
Ontology: BP
'classic' algorithm with the 'fisher' test
15 GO terms scored: 0 terms with p < 0.01
Annotation data:
Annotated genes: 2
Significant genes: 1
Min. no. of genes annotated to a GO: 1
Nontrivial nodes: 1
The p-values have not been corrected for multiple testing.
We can list the top ten significant results found:
> allRes <- GenTable(myGOdata, classicFisher = resultFisher, orderBy = "resultFisher", ranksOf = "classicFisher", topNodes = 10)
GO.ID Term Annotated Significant
1 GO:0005525 GTP binding 3 2
2 GO:0019001 guanyl nucleotide binding 3 2
3 GO:0032561 guanyl ribonucleotide binding 3 2
4 GO:0003746 translation elongation factor activity 1 1
5 GO:0005126 cytokine receptor binding 1 1
6 GO:0005132 interferon-alpha/beta receptor binding 1 1
7 GO:0008135 translation factor activity, nucleic aci... 1 1
8 GO:0008332 low voltage-gated calcium channel activi... 1 1
9 GO:0016667 oxidoreductase activity, acting on a sul... 1 1
10 GO:0016670 oxidoreductase activity, acting on a sul... 1 1
Expected classicFisher
1 0.3 0.026
2 0.3 0.026
3 0.3 0.026
4 0.1 0.101
5 0.1 0.101
6 0.1 0.101
7 0.1 0.101
8 0.1 0.101
9 0.1 0.101
10 0.1 0.101
We can visualise the position of the statistically significant GO terms in the GO hierarchy by using the following functions:
> showSigOfNodes(myGOdata, score(resultFisher), firstSigNodes = 5, useInfo ='all')
> printGraph(myGOdata, resultFisher, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)
The second command makes a pdf file ("tGO_classic_5_all.pdf") with the picture. The significant GO terms are shown as rectangles in the picture. The most significant terms are coloured red and least significant in yellow:
How TopGO works
The details of how TopGO are published in a paper by Alexa et al (2006). TopGO takes the GO hierarchy into account when calculating enrichment, and it is thought that this will lead it to report fewer false positive results than methods that do not take the GO hierarchy into account (eg. just using a Fisher's exact test). Grossman et al (2007) et al say that if you don't take the GO hierarchy into account, you get what they call the 'inheritance problem': because the highest (more general) GO terms inherit annotations from more specific descendant terms, and this can lead to false positives.
There are two methods described in the Alexa et al (2006) paper:
(i) 'elim' : this method processes the GO terms by traversing the GO hierarchy from bottom to top, ie. it first assesses the most specific (bottom-most) GO terms, and proceeds later to more general (higher) GO terms. When it assesses a higher (more general) GO term, it discards any genes that are annotated with significantly enriched descendant GO terms (considered significant using a pre-defined P-value threshold). This method does tend to miss some true positives at higher (more general) levels of the GO hierarchy.
(ii) 'weight': this is designed to miss less true positives than the 'elim' method. The significance scores of connected nodes (a parent and its child in the GO hierarchy) are compared in order to detect the locally most significant terms in the GO hierarchy.
The current TopGO actually has four more methods:
(iii) 'classic': each GO term is tested independently, not taking the GO hierarchy into account.
(iv) 'weight01': this is the default method used by TopGO, and is a mixture of the 'elim' and 'weight' methods.
(v) 'lea': this doesn't seem to be described in the TopGO vignette or any paper, so I'm not sure what it does (?)
(vi) 'parentchild': described by Grossman et al (2007) : when assessing a GO term, it takes into accoount the annotation of terms to the current term's parents, and so reduces false positives due to the inheritance problem. You can also use this algorithm in a very nice tool called Ontologizer , which has a lovely GUI (I used it in a collaboration on gene duplicates with Shane Woods, described in his paper Woods et al (2013)).
The 'elim' and 'weight' methods are designed to be more conservative (have less false positives) than the 'classic' method, so should have larger p-values in general. Simulation results reported by Alexa et al (2006) show that the 'weight' algorithm has less false positives than 'classic' and misses few true positives, while 'elim' has even less false positives than 'weight' but misses more true positives. Grossman et al (2007) compare the 'elim', 'weight' and 'parentchild' algorithms, and say that each has advantages in certain scenarios.
How to run TopGO
The vignette described in detail how to use the package, but here's a summary of what I did (so I won't forget!):
R packages
It's important to have the latest version of the packages, to make sure you have bug fixes, but also because the gene ontology is read into TopGO from another R package (GO.db), so you need to make sure you have the latest version of that package.
First you need to make sure that you have the latest version of R.
Then make sure you have the latest version of Bioconductor, and the latest version of the Bioconductor packages topGO:
>
source("http://bioconductor.org/biocLite.R")> biocLite()> biocLite("topGO")Then check you have the latest versions of everything: (for me it should be R 3.2.0, Bioconductor 3.1 and TopGO 2.2.0 and GO.db (comes with TopGO) 3.1.2:
> sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] topGO_2.20.0 SparseM_1.6 GO.db_3.1.2
[4] RSQLite_1.0.0 DBI_0.3.1 AnnotationDbi_1.30.1
[7] GenomeInfoDb_1.4.1 IRanges_2.2.5 S4Vectors_0.6.1
[10] Biobase_2.28.0 graph_1.46.0 BiocGenerics_0.14.0
[13] BiocInstaller_1.18.3
loaded via a namespace (and not attached):
[1] lattice_0.20-31 grid_3.2.0 tools_3.2.0
Starting topGO
Once you have installed topGO you can start it by typing:
> library("topGO")
Loading your data
Reading in GO annotations for the genes
The file that you provide of gene-to-GO mappings should look be a tab-delimited file looking like this, where the first column is the gene name and the second column contains a list of GO terms (separated by a comma and space):
105778 GO:0016021, GO:0016020
166881 GO:0003674, GO:0016021, GO:0016020, GO:0008150
188691 GO:0016021, GO:0016020, GO:0001584, GO:0007165, GO:0004984, GO:0007186, GO:0004872
184033 GO:0005737, GO:0005515
For example, if your file is called 'annotations.txt', you can read it in using:
> geneID2GO <- readMappings(file = "annotations.txt")
To have the fastest run-time, for a particular gene, it shouldn't be the case that one of the GO terms is an ancestor of one of the other GO terms for that gene. However, TopGO will run fine if this is the case, although it will be slower. In practice, I found that I didn't need to prune my lists of GO terms like this, as TopGO was pretty fast for my data set.
GO hierarchy
The GO hierarchy used by topGO is the one that comes with the GO.db package. Therefore the user does not need to provide this.
Defining your list of genes of interest, and the 'gene universe' to compare it to
You need to tell TopGO what is your list of genes of interest, and what is the list of all genes in the 'gene universe' that you want to compare your list of interest to.
The 'gene universe' can be all the genes that contain GO annotations in your input 'annotations.txt' file (see above):
> geneUniverse <- names(geneID2GO)
We can read in a list of interesting genes from a file, just containing a list of interesting genes, one per line:
> genesOfInterest <- read.table("interestinggenes.txt",header=FALSE)
> genesOfInterest <- as.character(genesOfInterest$V1)
Then we need to tell TopGO where these interesting genes appear in the 'geneUniverse' vector:
> geneList <- factor(as.integer(geneUniverse %in% genesOfInterest))
> names(geneList) <- geneUniverse
The geneList object tells TopGO which genes in the gene universe are your genes of interest.
Putting the data together into an R object
We need to put our data in an object of type 'topGOdata'. This will contain the list of genes of interest, the GO annotations, and the GO hierarchy itself.
> myGOdata <- new("topGOdata", description="My project", ontology="BP", allGenes=geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
The 'ontology' argument can be 'BP' (biological process), 'MF' (molecular function), or 'CC' (cellular component).
The 'description' argument has a short description of your project.
The 'allGenes' argument specifies all the genes in the 'gene universe', and which of those are your genes of interest.
The 'annot' argument tells topGO how to map genes to GO annotations. 'annot'=annFUN.gene2GO means that the user provides gene-to-GO annotations, and we specify here that they are in object 'geneID2GO' (see above for how this was created).
An optional argument is 'nodeSize': this is used to prune the GO hierarchy, eg. nodesize=10 prunes the GO hierarchy, to remove terms which have less than 10 annotated genes.
TopGO generally only seems to use a subset of my 'gene universe' data, which it considers 'feasible'. I think this is because you may have some genes in your input 'gene universe' that don't have any GO terms annotated to them. Also, you might have some terms in the GO annotations for your gene universe that don't exist in the version of the GO hierarchy that TopGO uses (which comes from the GO.db R package). Unfortunately this means that TopGO actually only uses some of your input data, and throws some away. Note that also, you may have chosen to use 'BP' (biological process) terms but many of your input annotations may be molecular function (MF) or cellular component (CC) terms.
If you now type 'myGOdata', you get a summary of the data in it:
> myGOdata
Description:
- My project
Ontology:
- BP
4 available genes (all genes from the array):
- symbol: 105778 166881 188691 184033 ...
- 2 significant genes.
2 feasible genes (genes that can be used in the analysis):
- symbol: 166881 188691 ...
- 1 significant genes.
GO graph (nodes with at least 1 genes):
- a graph with directed edges
- number of nodes = 15
- number of edges = 22
Here 2 of the 4 genes in the 'gene universe' have been mapped to the GO ontology (there are 2 'feasible' genes), so the rest of the analysis will be based on just these genes. TopGO calls the genes of interest 'significant genes'. There are 2 genes of interest in the full gene universe of 4 genes, but just one gene of interest among the 2 feasible genes. The 'myGOdata' object takes some time to build, especially if you have a large gene universe. I found it pretty fast though, less than one minute for a gene universe of ~100,000 genes. I guess the run-time might depend on how much memory (RAM) you have available(?)
The list of genes of interest can be accessed using the method sigGenes():
> sg <- sigGenes(myGOdata)
> str(sg)
> numSigGenes(myGOdata)
Performing enrichment tests
One type of test that topGO does is a Fisher's exact test based on gene counts:
> resultFisher <- runTest(myGOdata, algorithm="classic", statistic="fisher")
Selecting 'algorithm=classic' means that the GO hierarchy isn't taken into account, so each GO term is tested independently. (Usually in fact we want to take the GO hierarchy into account, so would use algorithm='weight01' for example.)
You can type 'resultFisher' to get a summary of the results:
Description: My project
Ontology: BP
'classic' algorithm with the 'fisher' test
15 GO terms scored: 0 terms with p < 0.01
Annotation data:
Annotated genes: 2
Significant genes: 1
Min. no. of genes annotated to a GO: 1
Nontrivial nodes: 1
The p-values have not been corrected for multiple testing.
We can list the top ten significant results found:
> allRes <- GenTable(myGOdata, classicFisher = resultFisher, orderBy = "resultFisher", ranksOf = "classicFisher", topNodes = 10)
GO.ID Term Annotated Significant
1 GO:0005525 GTP binding 3 2
2 GO:0019001 guanyl nucleotide binding 3 2
3 GO:0032561 guanyl ribonucleotide binding 3 2
4 GO:0003746 translation elongation factor activity 1 1
5 GO:0005126 cytokine receptor binding 1 1
6 GO:0005132 interferon-alpha/beta receptor binding 1 1
7 GO:0008135 translation factor activity, nucleic aci... 1 1
8 GO:0008332 low voltage-gated calcium channel activi... 1 1
9 GO:0016667 oxidoreductase activity, acting on a sul... 1 1
10 GO:0016670 oxidoreductase activity, acting on a sul... 1 1
Expected classicFisher
1 0.3 0.026
2 0.3 0.026
3 0.3 0.026
4 0.1 0.101
5 0.1 0.101
6 0.1 0.101
7 0.1 0.101
8 0.1 0.101
9 0.1 0.101
10 0.1 0.101
We can visualise the position of the statistically significant GO terms in the GO hierarchy by using the following functions:
> showSigOfNodes(myGOdata, score(resultFisher), firstSigNodes = 5, useInfo ='all')
> printGraph(myGOdata, resultFisher, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)
The second command makes a pdf file ("tGO_classic_5_all.pdf") with the picture. The significant GO terms are shown as rectangles in the picture. The most significant terms are coloured red and least significant in yellow:
Finding the genes annotated with significant GO terms
If you use the weight01 (default), elim or parentchild algorithm, TopGO may report a GO term as significantly enriched which wasn't actually in your input annotations, but rather is an ancestor (in the GO hierarchy) of a GO term in your input annotations.
The number of GO terms in the TopGO subset of the GO hierarchy (the
GO terms annotated to genes in the gene 'universe' input file, plus the
ancestors of those GO terms in the GO hierarchy) can be found using:
> length(usedGO(myGOdata))
[1] 2398
This should agree with the 'number of nodes' given for the GO graph when you type 'myGOdata', as it is the number of nodes in TopGO's internal version of the GO hierarchy (which is filtered just to contain your terms annotated to your gene universe, plus their ancestors).
You may find that TopGO says that terms "GO:0007610", "GO:0014070", "GO:0045910", but how do we find the genes that TopGO is considering as annotated with those terms, without having a good look at the GO hierarchy ourselves?
After a bit of (slightly painful, as I'm not too good at R!) fiddling with TopGO, I found out that to list the genes annotated to a particular set of GO terms, we can type something like this:
> myterms = c("GO:0007610", "GO:0014070", "GO:0045910")
> mygenes <- genesInTerm(myGOdata, myterms)
> for (i in 1:length(myterms))
{
myterm <- myterms[i]
mygenesforterm <- mygenes[myterm][[1]]
mygenesforterm <- paste(mygenesforterm, collapse=',')
print(paste("Term",myterm,"genes:",mygenesforterm))
}
[1]
"Term GO:0007610 genes:
173785,194776,210700,213042,222387,257151,272506,315440,347686,357493,361031,367238,397820,404238,435988,436908,533191,574720,575998,604672,615430,698226,718388,721238,759288,759468,759918,858210,866451,94187"
[1]
"Term GO:0014070 genes:
1050444,1051974,1058304,1119532,1137071,12391,163450,263383,364148,454323,725603,737096,737144,765031,768366,892110,935204,961708,989693"
[1] "Term GO:0045910 genes: 215658"
Running TopGO using an R script
You can run TopGO using an R script by using a command such as:
Rscript Rscript1.R annotations.txt interestinggenes.txt output_file.txt
where Rscript1.R is your R script, annotations.txt and interestinggenes.txt are the files defined as above, and output_file.txt is the file to save R output in.
Rscript1.R has this text, and executable permissions:
# Rscript1.R
# eg. run:
# Rscript Rscript1.R annotations.txt interestinggenes.txt output_file
args <- commandArgs(trailingOnly = TRUE)
universeFile = args[1]
interestingGenesFile = args[2]
output_file = args[3]
# set the output file
sink(output_file)
# load topGO
library("topGO")
# read in the 'gene universe' file
geneID2GO <- readMappings(file = universeFile)
geneUniverse <- names(geneID2GO)
# read in the genes of interest
genesOfInterest <- read.table(interestingGenesFile,header=FALSE)
genesOfInterest <- as.character(genesOfInterest$V1)
geneList <- factor(as.integer(geneUniverse %in% genesOfInterest))
names(geneList) <- geneUniverse
# build the GOdata object in topGO
myGOdata <- new("topGOdata", description="My project", ontology="BP", allGenes=geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
myGOdata
# run the Fisher's exact tests
resultClassic <- runTest(myGOdata, algorithm="classic", statistic="fisher")
resultElim <- runTest(myGOdata, algorithm="elim", statistic="fisher")
resultTopgo <- runTest(myGOdata, algorithm="weight01", statistic="fisher")
resultParentchild <- runTest(myGOdata, algorithm="parentchild", statistic="fisher")
# see how many results we get where weight01 gives a P-value <= 0.001:
mysummary <- summary(attributes(resultTopgo)$score <= 0.001)
numsignif <- as.integer(mysummary[[3]]) # how many terms is it true that P <= 0.001
# print out the top 'numsignif' results:
allRes <- GenTable(myGOdata, classicFisher = resultClassic, elimFisher = resultElim, topgoFisher = resultTopgo, parentchildFisher = resultParentchild, orderBy = "topgoFisher", ranksOf = "classicFisher", topNodes = numsignif)
allRes
# print a graph (to a pdf file) with the top 'numsignif' results:
output_file2 = paste(output_file,"Topgo", sep="_")
printGraph(myGOdata, resultTopgo, firstSigNodes = numsignif, fn.prefix = output_file2, useInfo = "all", pdfSW = TRUE)
# print out the genes that are annotated with the significantly enriched GO terms:
myterms <- allRes$GO.ID
mygenes <- genesInTerm(myGOdata, myterms)
for (i in 1:length(myterms))
{
myterm <- myterms[i]
mygenesforterm <- mygenes[myterm][[1]]
myfactor <- mygenesforterm %in% genesOfInterest # find the genes that are in the list of genes of interest
mygenesforterm2 <- mygenesforterm[myfactor == TRUE]
mygenesforterm2 <- paste(mygenesforterm2, collapse=',')
print(paste("Term",myterm,"genes:",mygenesforterm2))
}
# close the output file
sink()
Note: 15-Feb-2022: I found a nice blogpost by Zhigang Lu, on how to give TopGO scores for each gene as input (e.g. from Seurat P-values or something like that), and use them in a KS test in TopGO to calculate enrichment P-values: see Zhigang's blogpost here.
Monday, 29 June 2015
Creating a GO-Slim, and mapping GO terms to it
Creating a GO-Slim
I wanted to create my own GO-Slim, and found there is a nice tool for creating a GO-Slim at the EBI. You can start with an existing GO-Slim, eg. the 'generic GO-Slim' (has 149 terms), and add terms.
Mapping GO terms to your own GO-Slim
The next thing I wanted to do was to map GO terms to my GO-Slim. I was able to do this using Map2Slim, which is part of Owltools. To install it, I typed:
% wget http://build.berkeleybop.org/userContent/owltools/owltools
% wget http://build.berkeleybop.org/userContent/owltools/owltools-runner-all.jar
This is a Java program, it ran fine for me locally on my Mac laptop.
To run Map2Slim you need:
(i) a list of the GO terms in your GO-Slim
(ii) the GO annotations for your gene set of interest, with respect to the full ontology
(iii) the gene ontology hierarchy file (.obo file) for the full ontology
You need a list of the GO terms in your GO-Slim. If you are using the generic GO-Slim, you can download the generic GO-Slim in obo format from the geneontology.org website. Then to get a list of the terms, you can type:
% grep "id: GO:" goslim_generic.obo | grep -v alt | cut -d" " -f2 > goslim_terms.txt
The GO annotations for your gene set of interest need to be in GAF-2.0 format. In fact, I found that some of the columns aren't necessary for Map2Slim, so the file could look like this, where the columns marked 'optional' and 'unknown' seem to be ignored by Map2Slim:
!gaf-version: 2.0
WB 482159 482159 optional GO:0005515 pubmed unknown optional unknown optional optional protein unknown 482159 WB optional optional
WB 482159 482159 optional GO:0008270 pubmed unknown optional unknown optional optional protein unknown 482159 WB optional optional
WB 644503 644503 optional GO:0008270 pubmed unknown optional unknown optional optional protein unknown 644503 WB optional optional
WB 644503 644503 optional GO:0005230 pubmed unknown optional unknown optional optional protein unknown 644503 WB optional optional
You can download the latest gene ontology hierarchy (.obo) file from the geneontology.org website.
To run Map2Slim you type for example:
% ./owltools go-basic.obo --gaf my_gaf.txt --map2slim --idfile goslim_terms.txt --write-gaf my_slim.txt
The input GAF file was my_gaf.txt, the input obo file was go-basic.obo and the input list of GO-Slim terms was goslim_terms.txt.
The output file my_slim.txt is in GAF format, but has the GO-Slim terms for your genes. Hurray!
I wanted to create my own GO-Slim, and found there is a nice tool for creating a GO-Slim at the EBI. You can start with an existing GO-Slim, eg. the 'generic GO-Slim' (has 149 terms), and add terms.
Mapping GO terms to your own GO-Slim
The next thing I wanted to do was to map GO terms to my GO-Slim. I was able to do this using Map2Slim, which is part of Owltools. To install it, I typed:
% wget http://build.berkeleybop.org/userContent/owltools/owltools
% wget http://build.berkeleybop.org/userContent/owltools/owltools-runner-all.jar
This is a Java program, it ran fine for me locally on my Mac laptop.
To run Map2Slim you need:
(i) a list of the GO terms in your GO-Slim
(ii) the GO annotations for your gene set of interest, with respect to the full ontology
(iii) the gene ontology hierarchy file (.obo file) for the full ontology
You need a list of the GO terms in your GO-Slim. If you are using the generic GO-Slim, you can download the generic GO-Slim in obo format from the geneontology.org website. Then to get a list of the terms, you can type:
% grep "id: GO:" goslim_generic.obo | grep -v alt | cut -d" " -f2 > goslim_terms.txt
The GO annotations for your gene set of interest need to be in GAF-2.0 format. In fact, I found that some of the columns aren't necessary for Map2Slim, so the file could look like this, where the columns marked 'optional' and 'unknown' seem to be ignored by Map2Slim:
!gaf-version: 2.0
WB 482159 482159 optional GO:0005515 pubmed unknown optional unknown optional optional protein unknown 482159 WB optional optional
WB 482159 482159 optional GO:0008270 pubmed unknown optional unknown optional optional protein unknown 482159 WB optional optional
WB 644503 644503 optional GO:0008270 pubmed unknown optional unknown optional optional protein unknown 644503 WB optional optional
WB 644503 644503 optional GO:0005230 pubmed unknown optional unknown optional optional protein unknown 644503 WB optional optional
You can download the latest gene ontology hierarchy (.obo) file from the geneontology.org website.
To run Map2Slim you type for example:
% ./owltools go-basic.obo --gaf my_gaf.txt --map2slim --idfile goslim_terms.txt --write-gaf my_slim.txt
The input GAF file was my_gaf.txt, the input obo file was go-basic.obo and the input list of GO-Slim terms was goslim_terms.txt.
The output file my_slim.txt is in GAF format, but has the GO-Slim terms for your genes. Hurray!
Friday, 13 March 2015
Testing if two correlations are significantly different
I had two variables var1, var2, var3, and wanted to test whether the correlation between var1 and var3 is significantly different from the correlation between var2 and var3. I found a nice blog post by Jeromy Anglim explaining how to do this in R:
> var1 <- c(14.1,31.5,102.9,26.8,30.6,18.4,11.1,6.3,4.9,3.1,4.1,4.0,4.7)
> var2 <- c(54.9,18.3, 151.7, 53.6, 44.8,30.5,26.9,23.3,18.8,35.3,21.4,21.6,27.6)
> var3 <- c(84.7,63.5,272.8,94.1,100.3,74.6,53.0,47.3,49.4,60.5,43.1,42.8,52.2)
> cor.test(var1,var3,method="spearman")
rho=0.8021978
> cor.test(var2,var3,method="spearman")
rho=0.7692308
> cor.test(var1,var2, method="spearman")
rho=0.3736264
Now test whether a correlation of 0.80 is significantly different from 0.77 given that the two predictor variables (var1, var2) are correlated 0.37, assuming a sample size of 13:
> library(psych)
> paired.r(0.7692308,0.8021978,0.3736264,13)
Call: paired.r(xy = 0.7692308, xz = 0.8021978, yz = 0.3736264, n = 13)
[1] "test of difference between two correlated correlations"
t = -0.22 With probability = 0.83
This tells us that the correlation between var1 and var3 (r=0.80) was not significantly different than the correlation between var2 and var3 (r=77), n=13, p=0.83 (t-test).
Thanks to Jeromy Anglim!
> var1 <- c(14.1,31.5,102.9,26.8,30.6,18.4,11.1,6.3,4.9,3.1,4.1,4.0,4.7)
> var2 <- c(54.9,18.3, 151.7, 53.6, 44.8,30.5,26.9,23.3,18.8,35.3,21.4,21.6,27.6)
> var3 <- c(84.7,63.5,272.8,94.1,100.3,74.6,53.0,47.3,49.4,60.5,43.1,42.8,52.2)
> cor.test(var1,var3,method="spearman")
rho=0.8021978
> cor.test(var2,var3,method="spearman")
rho=0.7692308
> cor.test(var1,var2, method="spearman")
rho=0.3736264
Now test whether a correlation of 0.80 is significantly different from 0.77 given that the two predictor variables (var1, var2) are correlated 0.37, assuming a sample size of 13:
> library(psych)
> paired.r(0.7692308,0.8021978,0.3736264,13)
Call: paired.r(xy = 0.7692308, xz = 0.8021978, yz = 0.3736264, n = 13)
[1] "test of difference between two correlated correlations"
t = -0.22 With probability = 0.83
This tells us that the correlation between var1 and var3 (r=0.80) was not significantly different than the correlation between var2 and var3 (r=77), n=13, p=0.83 (t-test).
Thanks to Jeromy Anglim!
Subscribe to:
Posts (Atom)