RhodeCode

#!/usr/bin/env python
# coding=UTF-8
"""
Run overrepresentation analysis on Panther and return the results.

Usage: panther_bot.py [options] <gene_set>
       panther_bot.py -h | --help

Arguments:
  <gene_set>    file with gene IDs (one per line)

Options:
  -d NAME, --data=NAME   annotation dataset [default: bp]

"""
import time
import requests
import lxml.html

from docopt import docopt


URL = "http://www.pantherdb.org/"
HEADERS = {
	"Accept": "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8",
	"Accept-Encoding": "gzip, deflate",
	"Accept-Language": "en-US,en;q=0.5",
	"User-Agent": "Mozilla/5.0 (X11; Linux x86_64; rv:35.0) Gecko/20100101 Firefox/35.0",
	"Referer": URL
}
DATASETS = {"p":"pathway"}
DATASETS.update((name, "fullgo_{0}_exp".format(name)) for name in ["mf", "bp", "cc"])


def run(file_name, dataset):
	session = requests.session()
	session.get(URL, headers=HEADERS)
	time.sleep(1)

	data = {
		"idField": "",
		"fileType": 10,
		"organism": "Homo sapiens",
		"dataset": "Homo sapiens",
		"resultType": 3
	}
	with open(file_name, "rb") as genes:
		files = {"fileData": (file_name, genes)}
		response = session.post(URL + "geneListAnalysis.do", data=data, files=files,
			headers=HEADERS)
	time.sleep(1)

	params = {
		"numUploaded": 1,
		"numRefList": 1,
		"saveUserChoices": "false",
		"type": DATASETS[dataset],
		"bonferoni": 1
	}
	response = session.post(URL + "servlet/CompareToRefList", params=params,
		headers=HEADERS, allow_redirects=True)

	tree = lxml.html.document_fromstring(response.content)
	table = tree.xpath("//table[@class='scrollTableHead']")[1]

	results = []
	for row in table.iterfind("tr"):
		element = row.xpath("td[1]/a/text()")
		if element:
			term = element[0].strip()
			pvalue = float(row.xpath("td[7]/text()")[0])
			results.append([term, pvalue])

	session.close()
	return results


def main():
	args = docopt(__doc__)

	if args["--data"] not in ["mf", "bp", "cc", "p"]:
		print("Unknown dataset name. Choose one of [mf, bp, cc, p].")
		return

	results = run(args["<gene_set>"], args["--data"])
	results.sort(key=lambda x: x[1])
	for term, pvalue in results:
		print("{0:.5f} {1}".format(pvalue, term))


if __name__ == '__main__':
	main()