metastudent 2.0.1-1 / usr / bin / metastudent

#!/usr/bin/python
# -*- coding: utf-8 -*-
#
#   metastudent - predict GO terms from sequence
#   Copyright (C) 2012-2013 Tobias Hamp <hampt@rostlab.org>
#
#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

import tempfile
from getopt import getopt, GetoptError
import sys
import os
from shutil import copy2, rmtree
import random
import math
from ConfigParser import RawConfigParser
import collections
from metastudentPkg.runMethods import runBlast, runMethodA, runMethodB,\
	runMethodC
from metastudentPkg.Logger import Logger
from metastudentPkg.commons import GOConfig, GeneOntology, GOPrediction, encodeFastaHeaders, decodeFastaHeader, p, setSilent, isSilent, getPkgPath, splitBigFastaFile
import commands
import metastudentPkg
from metastudentPkg.version import VERSION
import urllib2

configMap=None
metastudentPath=""
metastudentPkgPath=""
dataPkgPath=""

def usage():
	print "metastudent -i <FASTA_FILE> -o <RESULT_FILE_PREFIX> [--silent] [--debug] [--keep-temp] [--output-blast] [--blast-only] [--all-predictions] [--no-names] [--with-images] [--ontologies=<MFO or BPO or CCO or MFO,BPO or ...>] [--blast-kickstart-databases=<BLAST_RESULT_FILE(S)>] [--temp-dir=<DIR>] [--config=<CONFIG_FILE>]\n!!! Make sure your input fasta file contains at most 500 sequences !!!"

def enrichOutputLinesWithNames(predLinesDict, geneOntology):
	for methodChar, predLines in predLinesDict.iteritems():
		newPredLines = []
		for predLine in predLines:
			targetId, goTerm, rel = predLine.split("\t")
			newPredLine = "%s\t%s\t%s\t%s" % (targetId, goTerm, rel, geneOntology.termToName[goTerm])
			newPredLines.append(newPredLine)
		newPredLinesAll.append(newPredLines)
		predLinesDict[methodChar] = newPredLines

def generateGraphImage(predLines):
	goTerms = []
	for predLine in predLines:
		targetId, goTerm = predLine.split("\t")[:2]
		goTerms.append(goTerm)
	
	chunks = []
	
	if len(goTerms) > 0:
		
		termDataTokens=[]
		for i, goTerm in enumerate(goTerms):
			if i == 0:
				header = "%22" + goTerm.replace(":", "%3A") + "%22%3A"
			else:
				header = "%2C%22" + goTerm.replace(":", "%3A") + "%22%3A"
			openBraket = "{"
			formatString = "%22fill%22%3A%22yellow%22%2C%22font%22%3A++%22black%22%2C%22border%22%3A%22red%22"
			closeBraket = "}"
			fullToken = header + openBraket + formatString + closeBraket
			termDataTokens.append(fullToken)
		termDataOpt = "term_data={" + "".join(termDataTokens) + "}%0D%0A%0D%0A%0D%0A%0D%0A%0D%0A%0D%0A"
		
		baseUrl = "http://amigo1.geneontology.org/cgi-bin/amigo/visualize?"
		modeOpt = "mode=advanced"
		termDataTypeOpt = "term_data_type=json"
		formatOpt = "format=png"
		
		fullURL = baseUrl + "&".join([modeOpt,termDataOpt,termDataTypeOpt,formatOpt])

		req = urllib2.urlopen(fullURL)
		
		CHUNK = 16 * 1024
		while True:
			chunk = req.read(CHUNK)
			if not chunk: 
				break
			chunks.append(chunk)
		
	return chunks

"""
calculate coefficients for linear combinations of results
"""
def metaPredict(preds, go, onto):
	
	targetList = preds.items()[0][1].targetsSorted
	
	protToTermToGroupToScore = collections.defaultdict(dict)
	
	for methodChar, predictioni in preds.iteritems():
		targetToTermToScore = predictioni.targetToTermToScore
		for target, termToScore in targetToTermToScore.iteritems():
			for term, score in termToScore.iteritems():
				protToTermToGroupToScore[target].setdefault((term,None),{})[methodChar] = float(score)	


	weightA = float(configMap["GROUP_A_WEIGHT_%s" % (onto)])
	weightB = float(configMap["GROUP_B_WEIGHT_%s" % (onto)])
	weightC = float(configMap["GROUP_C_WEIGHT_%s" % (onto)])
	intercept = float(configMap["INTERCEPT_%s" % (onto)])

	predLines = []
	for prot in protToTermToGroupToScore.keys():
		for (term, name), groupToScore in protToTermToGroupToScore[prot].iteritems():
			groupAScore = groupToScore.get("A", 0.0)
			groupBScore = groupToScore.get("B", 0.0)
			groupCScore = groupToScore.get("C", 0.0)
			metaScore = (weightA * groupAScore) + (weightB * groupBScore) + (weightC * groupCScore) + intercept
			if metaScore > 0.0:
				if metaScore > 0.95 and onto=="MFO":
					metaScore = 0.95
				elif metaScore > 0.95 and onto=="BPO":
					metaScore = 0.95
				elif metaScore > 0.95 and onto=="CCO":
					metaScore = 0.95
				
				predLines.append( prot + "\t" + term + "\t%.2f" % metaScore )
	
	metaPred = GOPrediction(predLines, go)
	
	return metaPred	
	
def runIt(tempfile, inputFastaFilePath, outputFilePath, outputBlast, blastKickstartDatabasePaths, ontologies, blastOnly, keepTemp, allPreds, debug, noNames, withImages):
	predLinesDict = {}
	predsDict = {}
	
	try:
		p("Creating tmpDir...")
		tmpDirPath=""
		tmpDirPath=tempfile.mkdtemp(prefix=configMap["MYNAME"])
		p("\tUsing " + tmpDirPath)
		
		p("Copying input file to tmpDir...")
		fastaFilePathLocal=""
		try:
			fastaFilePathLocal = os.path.join(tmpDirPath, os.path.basename(inputFastaFilePath)) 
			copy2(inputFastaFilePath, os.path.join(tmpDirPath, fastaFilePathLocal))
			encodeFastaHeaders(os.path.join(tmpDirPath, fastaFilePathLocal))
		except (IOError, os.error), why:
			Logger.log("Error copying files " + str(why))
		
		p("Loading Gene Ontology...")
		try:
			goTreeFilePath = os.path.join(configMap["DATABASE_BASE_PATH"], "goGraph.txt")	
			closureFilePath = os.path.join(configMap["DATABASE_BASE_PATH"], "fullTransitiveClosureGO.txt")
			nameMappingFilePath = os.path.join(configMap["DATABASE_BASE_PATH"], "nameMapping.txt")
			goConfig = GOConfig()
			goConfig.setTreeFilePath(goTreeFilePath)
			goConfig.setClosureFilePath(closureFilePath)
			goConfig.setNameMappingFilePath(nameMappingFilePath)
			geneOntology = GeneOntology(goConfig)
		except (IOError, os.error), why:
			Logger.log("Error loading Gene Ontology " + str(why))
		
		for i, ontology in enumerate(ontologies):
			
			database = os.path.join(os.path.join(configMap["DATABASE_BASE_PATH"], ontology), configMap["BLAST_SRC_%s" % (ontology) ] + ".fasta")
			
			suffix="_eval%s_iters%s_src%s.%s.blast" % (configMap["BLAST_EVAL_%s" % (ontology)], configMap["BLAST_ITERS_%s" % (ontology)], configMap["BLAST_SRC_%s" % (ontology)], ontology.lower())	
			blastOutputFilePathLocal = fastaFilePathLocal+suffix
			blastKickstartDatabasePath = blastKickstartDatabasePaths[i]
			
			if blastKickstartDatabasePath == "": #no kickstart
				p("Running Blast %s" % (ontology))
				runBlast(fastaFilePathLocal, database, blastOutputFilePathLocal, tmpDirPath, configMap["BLAST_EVAL_%s" % (ontology)], configMap["BLAST_ITERS_%s" % (ontology)], configMap)
				blastKickstartDatabasePath = blastOutputFilePathLocal
			if not blastOnly:
				predLinesDict["A"] = runMethodA(blastKickstartDatabasePath, fastaFilePathLocal, tmpDirPath, configMap["GROUP_A_THRESHOLD_%s" % (ontology)], configMap["GROUP_A_K_%s" % (ontology) ], configMap["BLAST_ITERS_%s" % (ontology)], ontology, configMap)
				predLinesDict["B"] = runMethodB(blastKickstartDatabasePath, fastaFilePathLocal, tmpDirPath, configMap["GROUP_B_K_%s" % (ontology) ], ontology, configMap)
				predLinesDict["C"] = runMethodC(blastKickstartDatabasePath, fastaFilePathLocal, tmpDirPath, configMap["GROUP_C_SCORING_%s" % (ontology) ], ontology, configMap, debug)
				

				
				for methodChar, predLines in predLinesDict.iteritems():

					predsDict[methodChar] = GOPrediction(predLines, geneOntology)

					predsDict[methodChar].propagatePrediction()

					if not noNames:
						predLinesDict[methodChar] = predsDict[methodChar].toOutputLinesWithNames()
					else:

						predLinesDict[methodChar] = predsDict[methodChar].toOutputLines()

				if allPreds:
					for methodChar, predLines in predLinesDict.iteritems():

						methodOutputFile = open(outputFilePath+".%s.%s.txt" % (ontology, methodChar),'w')
						methodOutputFile.write("".join(predLinesDict[methodChar]))
						methodOutputFile.close()
						
						if withImages:
							methodImageFileChunks = generateGraphImage(predLinesDict[methodChar])
							if len(methodImageFileChunks) > 0:
								methodImageFile = open(outputFilePath+".%s.%s.png" % (ontology, methodChar),'w')
								for chunki in methodImageFileChunks:
									methodImageFile.write(chunki)
								methodImageFile.close()					

				predsMeta = metaPredict(predsDict, geneOntology, ontology)

				predsMeta.propagatePrediction()

				predLinesMeta=None
				if not noNames:
					predLinesMeta = predsMeta.toOutputLinesWithNames()
				else:
					predLinesMeta = predsMeta.toOutputLines()
					
				methodOutputFile = open(outputFilePath+".%s.txt" % (ontology),'w')
				methodOutputFile.write("".join(predLinesMeta))
				methodOutputFile.close()

				if withImages:
					methodImageFileChunks = generateGraphImage(predLinesMeta)
					if len(methodImageFileChunks) > 0:
						methodImageFile = open(outputFilePath+".%s.png" % (ontology),'w')
						for chunki in methodImageFileChunks:
							methodImageFile.write(chunki)
						methodImageFile.close()					


			if outputBlast:
				p("copying back blast output")
				try:
					copy2(blastOutputFilePathLocal, outputFilePath + suffix )
				except (IOError, os.error), why:
					Logger.log("Error copying file " + str(why))
			
	except:
		print >> sys.stderr,"Error occurred: " + str(sys.exc_info()[0].__name__)
		raise
		sys.exit(1)
	finally:
		if not keepTemp:
			p("Deleting temp directories...")
			rmtree(tmpDirPath)
		p("\n====== Metastudent finished ======\n")


def formatPath(path):
	return path.replace("<ms-path>",metastudentPath) \
		.replace("<pkg-path>", metastudentPkgPath) \
		.replace("<default-data-package-path>", dataPkgPath) \
		.replace("<python-package-path>", metastudentPkgPath)

def readConfigFile(configMap, configPath):
	config = RawConfigParser()
	try:
		config.read(configPath)
	except:
		Logger.log("Error: Could not properly read config file: " + configPath, level=1)
		Logger.log("Error type: " + str(sys.exc_info()[0].__name__),level=1)
		sys.exit(1)
	
	
	if not config.has_section("General"):
		Logger.log("Error: config file " + configPath + " does not have a [General] section")
	else:
		if config.has_option('General', 'MYNAME'):
			configMap["MYNAME"] =  config.get('General', 'MYNAME')	
		if config.has_option('General', 'TMPDIR'):
			if config.get('General', 'TMPDIR') == "":
				configMap["TMPDIR"]=tempfile.gettempdir()
			elif config.get('General', 'TMPDIR') == "base":
				tmpPath = os.path.join(os.path.join(metastudentPath, "tmp"))
				if not os.path.exists(tmpPath):
					os.mkdir(tmpPath)
				configMap["TMPDIR"] = tmpPath
			else:
				if os.path.isdir(config.get('General', 'TMPDIR')):
					configMap["TMPDIR"]=config.get('General', 'TMPDIR')
				else:
					print >> sys.stderr,"Wrong argument for option TMPDIR"
					sys.exit(1)
		if config.has_option('General', "BLASTPGP_EXE_PATH"):
			if os.path.exists(config.get('General', 'BLASTPGP_EXE_PATH')):
				configMap["BLASTPGP_EXE_PATH"]=config.get('General', 'BLASTPGP_EXE_PATH')
			else:
				print >> sys.stderr,"Wrong argument for option BLASTPGP_EXE_PATH"
				sys.exit(1)
                if config.has_option('General', "JAR_INSTALL_FOLDER_PATH"):
                        if os.path.exists(config.get('General', 'JAR_INSTALL_FOLDER_PATH')):
                                configMap["JAR_INSTALL_FOLDER_PATH"]=config.get('General', 'JAR_INSTALL_FOLDER_PATH')
                        else:   
                                print >> sys.stderr,"Wrong argument for option JAR_INSTALL_FOLDER_PATH"
                                sys.exit(1)
		if config.has_option('General', "DATABASE_BASE_PATH"):
			configMap["DATABASE_BASE_PATH"]=formatPath(config.get('General', 'DATABASE_BASE_PATH'))
		if config.has_option('General', "FASTA_SPLIT_SIZE"):
			try:
				configMap["FASTA_SPLIT_SIZE"]=int(config.get('General', 'FASTA_SPLIT_SIZE'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option FASTA_SPLIT_SIZE"
				sys.exit(1)
		if config.has_option('General', 'GROUP_A_PATH'):
			if os.path.isdir(formatPath(config.get('General', 'GROUP_A_PATH'))):
				configMap["GROUP_A_PATH"]=formatPath(config.get('General', 'GROUP_A_PATH'))
			else:
				print >> sys.stderr,"Wrong argument for option GROUP_A_PATH"
				sys.exit(1)	
		if config.has_option('General', 'GROUP_B_PATH'):
			if os.path.isdir(formatPath(config.get('General', 'GROUP_B_PATH'))):
				configMap["GROUP_B_PATH"]=formatPath(config.get('General', 'GROUP_B_PATH'))
			else:
				print >> sys.stderr,"Wrong argument for option GROUP_B_PATH"
				sys.exit(1)	
		if config.has_option('General', 'GROUP_C_PATH'):
			if os.path.isdir(formatPath(config.get('General', 'GROUP_C_PATH'))):
				configMap["GROUP_C_PATH"]=formatPath(config.get('General', 'GROUP_C_PATH'))
			else:
				print >> sys.stderr,"Wrong argument for option GROUP_C_PATH"
				sys.exit(1)	






		if config.has_option('General', "GROUP_A_WEIGHT_MFO"):
			try:
				configMap["GROUP_A_WEIGHT_MFO"]=float(config.get('General', 'GROUP_A_WEIGHT_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_WEIGHT_MFO"
				sys.exit(1)
		if config.has_option('General', "GROUP_B_WEIGHT_MFO"):
			try:
				configMap["GROUP_B_WEIGHT_MFO"]=float(config.get('General', 'GROUP_B_WEIGHT_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_B_WEIGHT_MFO"
				sys.exit(1)
		if config.has_option('General', "GROUP_C_WEIGHT_MFO"):
			try:
				configMap["GROUP_C_WEIGHT_MFO"]=float(config.get('General', 'GROUP_C_WEIGHT_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_C_WEIGHT_MFO"
				sys.exit(1)
		if config.has_option('General', "INTERCEPT_MFO"):
			try:
				configMap["INTERCEPT_MFO"]=float(config.get('General', 'INTERCEPT_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option INTERCEPT_MFO"
				sys.exit(1)
				
		if config.has_option('General', "GROUP_A_THRESHOLD_MFO"):
			try:
				float(config.get('General', 'GROUP_A_THRESHOLD_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_THRESHOLD_MFO"
				sys.exit(1)
			configMap["GROUP_A_THRESHOLD_MFO"]=config.get('General', 'GROUP_A_THRESHOLD_MFO')
		if config.has_option('General', "GROUP_A_K_MFO"):
			try:
				int(config.get('General', 'GROUP_A_K_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_K_MFO"
				sys.exit(1)
			configMap["GROUP_A_K_MFO"]=config.get('General', 'GROUP_A_K_MFO')
		if config.has_option('General', "GROUP_B_K_MFO"):
			try:
				int(config.get('General', 'GROUP_B_K_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_B_K_MFO"
				sys.exit(1)
			configMap["GROUP_B_K_MFO"]=config.get('General', 'GROUP_B_K_MFO')
		if config.has_option('General', "GROUP_C_SCORING_MFO"):
			try:
				int(config.get('General', 'GROUP_C_SCORING_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_C_SCORING_MFO"
				sys.exit(1)
			configMap["GROUP_C_SCORING_MFO"]=config.get('General', 'GROUP_C_SCORING_MFO')
		if config.has_option('General', "BLAST_EVAL_MFO"):
			try:
				float(config.get('General', 'BLAST_EVAL_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option BLAST_EVAL_MFO"
				sys.exit(1)
			configMap["BLAST_EVAL_MFO"]=config.get('General', 'BLAST_EVAL_MFO')
		if config.has_option('General', "BLAST_ITERS_MFO"):
			try:
				int(config.get('General', 'BLAST_ITERS_MFO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option BLAST_ITERS_MFO"
				sys.exit(1)
			configMap["BLAST_ITERS_MFO"]=config.get('General', 'BLAST_ITERS_MFO')
		if config.has_option('General', "BLAST_SRC_MFO"):
			if config.get('General', 'BLAST_SRC_MFO') in ["goaexp", "goasp", "spall", "spexp", "upexp", "all","exp","expup"]:
				configMap["BLAST_SRC_MFO"] = config.get('General', 'BLAST_SRC_MFO')
			else:
				print >> sys.stderr,"Wrong argument for option BLAST_SRC_MFO"
				sys.exit(1)




		if config.has_option('General', "GROUP_A_WEIGHT_BPO"):
			try:
				configMap["GROUP_A_WEIGHT_BPO"]=float(config.get('General', 'GROUP_A_WEIGHT_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_WEIGHT_BPO"
				sys.exit(1)
		if config.has_option('General', "GROUP_B_WEIGHT_BPO"):
			try:
				configMap["GROUP_B_WEIGHT_BPO"]=float(config.get('General', 'GROUP_B_WEIGHT_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_B_WEIGHT_BPO"
				sys.exit(1)
		if config.has_option('General', "GROUP_C_WEIGHT_BPO"):
			try:
				configMap["GROUP_C_WEIGHT_BPO"]=float(config.get('General', 'GROUP_C_WEIGHT_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_C_WEIGHT_BPO"
				sys.exit(1)
		if config.has_option('General', "INTERCEPT_BPO"):
			try:
				configMap["INTERCEPT_BPO"]=float(config.get('General', 'INTERCEPT_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option INTERCEPT_BPO"
				sys.exit(1)
				
		if config.has_option('General', "GROUP_A_THRESHOLD_BPO"):
			try:
				float(config.get('General', 'GROUP_A_THRESHOLD_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_THRESHOLD_BPO"
				sys.exit(1)
			configMap["GROUP_A_THRESHOLD_BPO"]=config.get('General', 'GROUP_A_THRESHOLD_BPO')
		if config.has_option('General', "GROUP_A_K_BPO"):
			try:
				int(config.get('General', 'GROUP_A_K_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_K_BPO"
				sys.exit(1)
			configMap["GROUP_A_K_BPO"]=config.get('General', 'GROUP_A_K_BPO')
		if config.has_option('General', "GROUP_B_K_BPO"):
			try:
				int(config.get('General', 'GROUP_B_K_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_B_K_BPO"
				sys.exit(1)
			configMap["GROUP_B_K_BPO"]=config.get('General', 'GROUP_B_K_BPO')
		if config.has_option('General', "GROUP_C_SCORING_BPO"):
			try:
				int(config.get('General', 'GROUP_C_SCORING_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_C_SCORING_BPO"
				sys.exit(1)
			configMap["GROUP_C_SCORING_BPO"]=config.get('General', 'GROUP_C_SCORING_BPO')
		if config.has_option('General', "BLAST_EVAL_BPO"):
			try:
				float(config.get('General', 'BLAST_EVAL_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option BLAST_EVAL_BPO"
				sys.exit(1)
			configMap["BLAST_EVAL_BPO"]=config.get('General', 'BLAST_EVAL_BPO')
		if config.has_option('General', "BLAST_ITERS_BPO"):
			try:
				int(config.get('General', 'BLAST_ITERS_BPO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option BLAST_ITERS_BPO"
				sys.exit(1)
			configMap["BLAST_ITERS_BPO"]=config.get('General', 'BLAST_ITERS_BPO')
		if config.has_option('General', "BLAST_SRC_BPO"):
			if config.get('General', 'BLAST_SRC_BPO') in ["goaexp", "goasp", "spall", "spexp", "upexp", "all","exp","expup"]:
				configMap["BLAST_SRC_BPO"] = config.get('General', 'BLAST_SRC_BPO')
			else:
				print >> sys.stderr,"Wrong argument for option BLAST_SRC_BPO"
				sys.exit(1)

				


		if config.has_option('General', "GROUP_A_WEIGHT_CCO"):
			try:
				configMap["GROUP_A_WEIGHT_CCO"]=float(config.get('General', 'GROUP_A_WEIGHT_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_WEIGHT_CCO"
				sys.exit(1)
		if config.has_option('General', "GROUP_B_WEIGHT_CCO"):
			try:
				configMap["GROUP_B_WEIGHT_CCO"]=float(config.get('General', 'GROUP_B_WEIGHT_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_B_WEIGHT_CCO"
				sys.exit(1)
		if config.has_option('General', "GROUP_C_WEIGHT_CCO"):
			try:
				configMap["GROUP_C_WEIGHT_CCO"]=float(config.get('General', 'GROUP_C_WEIGHT_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_C_WEIGHT_CCO"
				sys.exit(1)
		if config.has_option('General', "INTERCEPT_CCO"):
			try:
				configMap["INTERCEPT_CCO"]=float(config.get('General', 'INTERCEPT_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option INTERCEPT_CCO"
				sys.exit(1)
				
		if config.has_option('General', "GROUP_A_THRESHOLD_CCO"):
			try:
				float(config.get('General', 'GROUP_A_THRESHOLD_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_THRESHOLD_CCO"
				sys.exit(1)
			configMap["GROUP_A_THRESHOLD_CCO"]=config.get('General', 'GROUP_A_THRESHOLD_CCO')
		if config.has_option('General', "GROUP_A_K_CCO"):
			try:
				int(config.get('General', 'GROUP_A_K_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_A_K_CCO"
				sys.exit(1)
			configMap["GROUP_A_K_CCO"]=config.get('General', 'GROUP_A_K_CCO')
		if config.has_option('General', "GROUP_B_K_CCO"):
			try:
				int(config.get('General', 'GROUP_B_K_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_B_K_CCO"
				sys.exit(1)
			configMap["GROUP_B_K_CCO"]=config.get('General', 'GROUP_B_K_CCO')
		if config.has_option('General', "GROUP_C_SCORING_CCO"):
			try:
				int(config.get('General', 'GROUP_C_SCORING_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option GROUP_C_SCORING_CCO"
				sys.exit(1)
			configMap["GROUP_C_SCORING_CCO"]=config.get('General', 'GROUP_C_SCORING_CCO')
		if config.has_option('General', "BLAST_EVAL_CCO"):
			try:
				float(config.get('General', 'BLAST_EVAL_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option BLAST_EVAL_CCO"
				sys.exit(1)
			configMap["BLAST_EVAL_CCO"]=config.get('General', 'BLAST_EVAL_CCO')
		if config.has_option('General', "BLAST_ITERS_CCO"):
			try:
				int(config.get('General', 'BLAST_ITERS_CCO'))
			except ValueError:
				print >> sys.stderr,"Wrong argument for option BLAST_ITERS_CCO"
				sys.exit(1)
			configMap["BLAST_ITERS_CCO"]=config.get('General', 'BLAST_ITERS_CCO')
		if config.has_option('General', "BLAST_SRC_CCO"):
			if config.get('General', 'BLAST_SRC_CCO') in ["goaexp", "goasp", "spall", "spexp", "upexp","all","exp","expup"]:
				configMap["BLAST_SRC_CCO"] = config.get('General', 'BLAST_SRC_CCO')
			else:
				print >> sys.stderr,"Wrong argument for option BLAST_SRC_CCO"
				sys.exit(1)



if __name__ == "__main__":
	currentPath = os.getcwd()
	if hasattr(sys, "frozen"):
		metastudentPath=os.path.dirname(unicode(sys.executable, sys.getfilesystemencoding( )))
	else:
		metastudentPath = os.path.dirname(os.path.abspath(__file__))
	metastudentPkgPath=getPkgPath()
	dataPkgPath=os.path.join(sys.prefix, "share", "metastudent-data")

	if commands.getstatusoutput('pp_popcon_cnt -p metastudent')[0] != 0:
		print >> sys.stderr, "The Rost Lab recommends you install the pp-popularity-contest package that provides pp_popcon_cnt:\n\nsudo apt-get install pp-popularity-contest\n"


	try:
		opts, args=getopt(sys.argv[1:], "hi:o:", ["debug", "help", "silent", "keep-temp", "temp-dir=", "version", "config=", "output-blast", "no-names", "with-images", "blast-kickstart-databases=", "ontologies=","blast-only","all-predictions"])
	except GetoptError, err:
		print >> sys.stderr, str(err)
		usage()
		sys.exit(1)
	
	debug=False
	inputFastaFilePath=""
	outputFilePath=""
	tempDirPath=""
	customConfigFilePath=""
	blastKickstartDatabasePaths=["","",""]
	outputBlast=False
	ontologies=["MFO","BPO","CCO"]
	blastOnly=False
	keepTemp=False
	allPreds = False
	noNames = False
	withImages = False

	Logger(os.path.join(metastudentPath, "metastudent.log"))
	
	optArgs = [arg[0] for arg in opts]
	if "--version" in optArgs:
		print VERSION
		sys.exit(0)
	if ("-i" not in optArgs and "--blast-kickstart-databases" not in optArgs) or "-o" not in optArgs:
		usage()
		sys.exit(1)
	
	for o, a in opts:
		if o == "-i":
			if os.path.exists(a):
				inputFastaFilePath=os.path.abspath(a)
			else:
				print >> sys.stderr,"Error: Input fasta file does not exist"
				sys.exit(1)
		elif o == "-o":
			if os.path.exists(os.path.dirname(os.path.abspath(a))):
				try:
					tmpFile=open(os.path.abspath(a), 'w')
					tmpFile.close()
					os.remove(os.path.abspath(a))
					outputFilePath=os.path.abspath(a)
				except:
					print >> sys.stderr,"Error: Cannot open output file for writing"
					sys.exit(1)
			else:
				print >> sys.stderr,"Error: Output file can not be created"
				sys.exit(1)
		elif o == "--temp-dir":
			if os.path.exists(a):
				tempDirPath=os.path.abspath(a)
			else:
				print >> sys.stderr,"Error: Temp directory does not exist."
				sys.exit(1)
		elif o == "--config":
			if os.path.exists(a):
				customConfigFilePath = a
			else:
				print >> sys.stderr,"Error: Config directory does not exist."
				sys.exit(1)
		elif o == "--blast-kickstart-databases":
			blastKickstartDatabasePathsTmp = a.split(",")
			blastKickstartDatabasePaths=[]
			for db in blastKickstartDatabasePathsTmp:
				dbAbs = db
				if db!="":
					dbAbs = os.path.abspath(db)
					if not os.path.exists(dbAbs):
						print >> sys.stderr,"Error: blast-kickstart-database does not exist."
						sys.exit(1)
				blastKickstartDatabasePaths.append(dbAbs)
		elif o == "--debug":
			debug = True;
		elif o in ("-h", "--help"):
			usage()
			sys.exit()
		elif o == "--silent":
			setSilent(True)
		elif o == "--blast-only":
			blastOnly=True
		elif o == "--output-blast":
			outputBlast = True
		elif o == "--keep-temp":
			keepTemp = True
		elif o == "--all-predictions":
			allPreds = True
		elif o == "--no-names":
			noNames = True
		elif o == "--with-images":
			withImages = True
		elif o == "--ontologies":
			ontologies = a.split(",")
			for ontology in ontologies:
				if ontology not in ["MFO", "BPO","CCO"]:
					print >> sys.stderr,"Error: unknown ontology: " % (ontology)
					sys.exit(1)
		else:
			assert False, "Error: Unhandled option: " + o
			sys.exit(1)

	if blastKickstartDatabasePaths != ["","",""]:
		if len(ontologies) != len(blastKickstartDatabasePaths):
			print >> sys.stderr,"Error: number of kickstart databases not equal to number of ontologies."
			sys.exit(1)
		else:
			p("BLAST results pre-calculted")
			for blastKickstartDatabasePath in blastKickstartDatabasePaths:
				if blastKickstartDatabasePath != "":
					p("\t" + blastKickstartDatabasePath)

	configPaths = [os.path.join(sys.prefix, "share", "metastudent", "metastudentrc.default"),
				"/etc/metastudentrc",
				os.path.expanduser("~/.metastudentrc"),
				customConfigFilePath]
	configMap={}
	if len([readConfigFile(configMap, configPath) for configPath in configPaths if os.path.exists(configPath)]) == 0:
		print >> sys.stderr,"Error: No config file found"
		sys.exit(1)
	
	if tempDirPath != "":
		tempfile.tempdir=tempDirPath
	else:
		if not os.path.exists(configMap["TMPDIR"]):
			os.mkdir(configMap["TMPDIR"])
		tempfile.tempdir=configMap["TMPDIR"]
	
	if not os.path.isdir(configMap.get('DATABASE_BASE_PATH')):
		print >> sys.stderr,"Error: DATABASE_BASE_PATH (%s) does not exist." % configMap.get('DATABASE_BASE_PATH')
		sys.exit(1);

	
	numSeqs = len([None for line in open(inputFastaFilePath).read().split("\n") if line.startswith(">")])
	if numSeqs != 1 and withImages:
		print >> sys.stderr,"Error: option --with-images can only be used with exactly one input sequence."
		sys.exit(1);
		
	
	runIt(tempfile, inputFastaFilePath, outputFilePath, outputBlast, blastKickstartDatabasePaths, ontologies, blastOnly, keepTemp, allPreds, debug, noNames, withImages)
	

"""
=pod

=head1 NAME

metastudent - predictor of gene ontology terms from protein sequence

=head1 SYNOPSIS

metastudent -i I<FASTA_FILE> -o I<RESULT_FILE_PREFIX> [--debug] [--keep-temp] [--silent] [--output-blast] [--blast-only] [--all-predictions] [--ontologies=I<MFO or BPO or CCO or MFO,BPO or ...>] [--with-images] [--blast-kickstart-databases=I<BLAST_RESULT_FILE(S)>] [--temp-dir=I<DIR>] [--config=I<CONFIG_FILE>]
!!! Make sure your fasta file contains at most 500 sequences !!!

=head1 DESCRIPTION

Metastudent predicts Gene Ontology (GO) terms from the Molecular Function Ontology (MFO), Biological Process Ontology (BPO) and Cellular Component Ontology (CCO)
for input protein sequences by homology-based inference from already annotated proteins.

Large (1 GB in total) data files necessary for the operation of metastudent are downloaded automatically on the first use of the program.  The download is restartable.  You can also make an explicit call to
       F<metastudentdata> (by default /usr/bin/metastudentdata) to download the data files.  In case the data directory (by default /usr/share/metastudent-data) is not writable and you are not root, the operation is
       reattempted with sudo(8).

=head2 Output format

For each selected ontology (see I<--ontologies>), one output file is produced (see I<-o>).
Each line in each file associates a protein with a GO term and a reliability for the association (0.0 to 1.0). The following format is used:
<PROTEIN ID><TAB><GO_TERM><TAB><RELIABILITY>

=head1 REFERENCES

=over

=item Hamp, T., Kassner, R., Seemayer, S., Vicedo, E., Schaefer, C., Achten, D., ... & Rost, B. (2013). Homology-based inference sets the bar high for protein function prediction. BMC Bioinformatics, 14(Suppl 3), S7.

=back

=head1 OPTIONS

=over

=item -i FASTA_FILE

The input fasta file. Please try to remove any special formattings (e.g. whitespaces) in the sequences before using them as input. Due to high memory usage, make sure your fasta file contains at most 500 sequences.

=item -o RESULT_FILE_PREFIX

The file name prefix of the output files. GO terms are organized in ontologies. Metatstudent treats each ontology differently and outputs one result file for each. For example, if <RESULT_FILE>=./myresult and MFO (Molecular Function Ontology) and BPO (Biological Process Ontology) ontologies are selected (see option B<--ontologies>), then metastudent creates two output files: F<./myresult.MFO.txt> and F<./myresult.BPO.txt>.

=item --debug

Print extra debugging messages.

=item --keep-temp

Whether to keep the temp directories after metastudent has finished (they can be useful when errors occur or in combination with B<--blast-kickstart-databases>).

=item --silent

No progress messages (stdout), only errors (stderr).

=item --output-blast

Whether to output the result of the BLAST runs. Useful in combination with B<--blast-kickstart-databases>. Output file name format is I<RESULT_FILE_PREFIX>.<BLAST_OPTIONS>.blast.

=item --blast-only

Whether to only output the result of the BLAST runs, and nothing else. See options I<--output-blast> and I<--blast-kickstart-databases>.

=item --all-predictions

Whether to output the prediction results of the individual predictors. File name format of the output file is <RESULT_FILE_PREFIX>.<ONTOLOGY>.<METHOD>.txt.

=item --no-names

Whether to exclude the actual GO term names as a column in the output prediction file(s).

=item --with-images

Whether to also generate an image that displays the predicted GO terms as a GO graph. This option can only be used with exactly one input sequence and only when connected to the internet.

=item --ontologies=I<MFO or BPO or CCO or MFO,BPO or ...>

A comma separated list of ontologies to create predictions for. Default is MFO,BPO,CCO. If used in combination with B<--blast-kickstart-databases>, the number and order of the ontologies must correspond to the kickstart files.

=item --blast-kickstart-databases=<BLAST_RESULT_FILES>

Since running BLAST is usually the part that takes the longest in metastudent, this option allows you to re-use the output of a previous run. This is useful to test, for example, different parameters or when you have lost a prediction. The number of kickstart files must correspond to the number of ontologies (see option B<--ontologies>). Separate the file paths with commas. For example: --blast-kickstart-databases=<RESULT_FILE_MFO>,<RESULT_FILE_BPO> (kickstart for both ontologies) or --blast-kickstart-databases=,<RESULT_FILE_BPO> (only kickstart BPO; note the comma).

=item --temp-dir=DIR

The parent temp directory to use instead of the one specified with tmpDir in the metastudent configuration file.

=item --config=FILE

The path to a custom metastudent configuration file; overrides all settings of the configuration files found in the FILES section of this man page.

=back

=head1 FILES

=over

=item F<<< <package_data_dir>/metastudentrc.default >>>

The metastudent configuration file.

=item F<<< <sysconfdir>/metastudentrc >>>

The metastudent configuration file, overrides F<<< <package_data_dir>/metastudentrc.default >>>.

=item F<<< <homedir>/.metastudentrc >>>

The metastudent configuration file, overrides F<<< <sysconfdir>/metastudentrc >>>.

=back

=head1 EXAMPLES

The example F<test.fasta> file can be found in F<<< <package_doc_dir>/examples >>> (usually F</usr/share/doc/metastudent/examples>).

=over

=item Predict the GO terms for the sequences in test.fasta for both the MFO and the BPO ontology:

 metastudent -i test.fasta -o test.result

=item Create the BLAST output to predict the MFO terms for sequences in test.fasta (not the actual predictions, yet; see next example).

 metastudent -i test.fasta -o test.result --blast-only --output-blast --ontologies=MFO

=item Predict the MFO and BPO terms for sequences in test.fasta with a precomputed MFO BLAST output (see previous example; note the comma at the end).

 metastudent -i test.fasta -o test.result --ontologies=MFO,BPO --blast-kickstart-databases=test.result_eval0.001_iters3_srcexp.mfo.blast,

=back

=head1 BUGS

Please use this link to report bugs:

L<https://rostlab.org/bugzilla3/enter_bug.cgi?product=metastudent>

=head1 AUTHOR

Tobias Hamp <hampt@rostlab.org>

=cut
"""