view commons/pyRepetUnit/profilesDB/ProfilesDB4Repet.py @ 31:0ab839023fe4

Uploaded

#!/usr/bin/env python

import re
import getopt
import sys
from commons.pyRepetUnit.profilesDB.Profiles import Profiles
from commons.core.LoggerFactory import LoggerFactory

LOG_DEPTH = "commons.pyRepetUnit.profiles"

## Format a profiles DB for pipelines in REPET
#    
class ProfilesDB4Repet( object ):
    
    def __init__(self):
        self.profile = Profiles()
        self._inputFile = "" 
        self._outputFile =  ""
        self._verbosity = 2
        self._log = LoggerFactory.createLogger("%s.%s" % (LOG_DEPTH, self.__class__.__name__), self._verbosity)
        
        
    def _help( self ):
        print
        print "usage: %s.py [ options ]" % ( type(self).__name__ )
        print "options:"
        print "     -h: this help"
        print "     -i: name of the profiles DB to format for Repet"
        print "     -o: name of the output profiles DB for Repet"
        print
        
        
    def _setAttributesFromCmdLine( self ):
        try:
            opts, args = getopt.getopt(sys.argv[1:],"hi:o:")
        except getopt.GetoptError, err:
            print str(err); self._help(); sys.exit(1)
        for o,a in opts:
            if o == "-h":
                self._help(); sys.exit(0)
            elif o == "-i":
                self.setInputFile( a )
            elif o == "-o":
                self.setOutputFile( a )
                
    #TDOD: add nb of each domain in log file, verbose...
    def _searchCurrentDomain(self, profile):
        currentDomain = ""
        #TODO: pattern GAGA and GAGE should be excluded !
        #TODO: add new tag like "ORF1_LTR" for ATHILA as key word in Pfam
        #TODO: add new tags from GypsyDB (MOV etc...)
        if (re.search("[gG][aA][Gg]", profile[1]) or re.search("[gG][aA][Gg]", profile[3])):
            currentDomain = "GAG"
        elif (re.search("Zinc knuckle", profile[1]) or re.search("Zinc knuckle", profile[3])):
            currentDomain = "GAG"
        elif (re.search("PF02813", profile[2]) or re.search("PF01021", profile[2])):
            currentDomain = "GAG"
        elif (re.search("GAG_", profile[1])):
            currentDomain = "GAG"
        elif (re.search("GAGCOAT_", profile[1])):
            currentDomain = "GAG"
            
        elif ((re.search("[aA]spartic", profile[1]) or re.search("[aA]aspartic", profile[3])) and (re.search("[pP]roteinase", profile[1]) or re.search("[pP]roteinase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartic", profile[1]) or re.search("[aA]spartic", profile[3])) and (re.search("[pP]rotease", profile[1]) or re.search("[pP]rotease", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[rR]etrotransposon", profile[1]) or re.search("[rR]etrotransposon", profile[3])) and (re.search("[pP]eptidase", profile[1]) or re.search("[pP]eptidase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartic", profile[1]) or re.search("[aA]spartic", profile[3])) and (re.search("[pP]eptidase", profile[1]) or re.search("[pP]eptidase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartic", profile[1]) or re.search("[aA]spartic", profile[3])) and (re.search("[eE]ndopeptidase", profile[1]) or re.search("[eE]ndopeptidase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartyl", profile[1]) or re.search("[aA]spartyl", profile[3])) and (re.search("[pP]roteinase", profile[1]) or re.search("[pP]roteinase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartyl", profile[1]) or re.search("[aA]spartyl", profile[3])) and (re.search("[pP]rotease", profile[1]) or re.search("[pP]rotease", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartyl", profile[1]) or re.search("[aA]spartyl", profile[3])) and (re.search("[pP]eptidase", profile[1]) or re.search("[pP]eptidase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("[aA]spartyl", profile[1]) or re.search("[aA]spartyl", profile[3])) and (re.search("[eE]ndopeptidase", profile[1]) or re.search("[eE]ndopeptidase", profile[3]))):
            currentDomain = "AP"
        elif ((re.search("AP_", profile[1])) and not (re.search("[eE]ndonuclease", profile[3]))):
            currentDomain = "AP"
            
        elif ((re.search("[iI]ntegrase", profile[1]) or re.search("[iI]ntegrase", profile[3])) and not (re.search(" C ", profile[1])) and not (re.search(" C ", profile[3]))):
            currentDomain = "INT"
        elif (re.search(".*[Cc]hromo.*", profile[1]) or re.search(".*[Cc]hromo.*", profile[3])):
            currentDomain = "INT"
        elif (re.search("INT_", profile[1])):
            currentDomain = "INT"
            
            
        elif ((re.search("[rR]everse", profile[1]) or re.search("[Rr]everse", profile[3])) and (re.search("[tT]ranscriptase", profile[1]) or re.search("[tT]ranscriptase", profile[3]))):
            currentDomain = "RT"
        elif (re.search("RT_", profile[1])):
            currentDomain = "RT"
            
            
        elif ((re.search("R[nN]ase", profile[1]) or re.search("R[nN]ase", profile[3])) and (re.search("H", profile[1]) or re.search("H", profile[3]))):
            currentDomain = "RH"
        elif ((re.search("R[nN]ase", profile[1]) or re.search("R[nN]ase", profile[3])) and (re.search("H ", profile[1]) or re.search("H ", profile[3]))):
            currentDomain = "RH"
        elif (re.search("RNaseH_", profile[1])):
            currentDomain = "RH"
            
            
        elif ((re.search("[eE]nvelope", profile[1]) or re.search("[eE]nvelope", profile[3])) and (re.search("[pP]rotein", profile[1]) or re.search("[pP]rotein", profile[3]))):
            currentDomain = "ENV"
        elif ((re.search("ENV", profile[1]) or re.search("ENV", profile[3])) and (re.search("[pP]olyprotein", profile[1]) or re.search("[pP]olyprotein", profile[3]))):
            currentDomain = "ENV"
        elif ((re.search("[eE]nvelope", profile[1]) or re.search("[eE]nvelope", profile[3])) and (re.search("[pP]olyprotein", profile[1]) or re.search("[pP]olyprotein", profile[3]))):
            currentDomain = "ENV"
        elif ((re.search("[eE]nvelope", profile[1]) or re.search("[eE]nvelope", profile[3])) and (re.search("[gG]lycoprotein", profile[1]) or re.search("[gG]lycoprotein", profile[3]))):
            currentDomain = "ENV"
        elif (re.search("PF07253", profile[2]) or re.search("PF03811", profile[2])):
            currentDomain = "ENV"
        elif (re.search("ENV_", profile[1])):
            currentDomain = "ENV"
            
        elif ((re.search("[tT]yrosine", profile[1]) or re.search("[tT]yrosine", profile[3])) and (re.search("[rR]ecombinase", profile[1]) or re.search("[rR]ecombinase", profile[3]))):
            currentDomain = "YR"
            
        elif ((re.search("[eE]ndonuclease", profile[1]) or re.search("[eE]ndonuclease", profile[3])) and not (re.search("AP ", profile[1]) or re.search("AP ", profile[3])) and not (re.search("[aA]purinic", profile[1]) or re.search("[aA]purinic", profile[3]))):
            currentDomain = "EN"
            
        elif ((re.search("[eE]ndonuclease", profile[1]) or re.search("[eE]ndonuclease", profile[3])) and (re.search("AP ", profile[1]) or re.search("AP ", profile[3]))):
            currentDomain = "APE"
        elif ((re.search("[eE]ndonuclease", profile[1]) or re.search("[eE]ndonuclease", profile[3])) and (re.search("[aA]purinic", profile[1]) or re.search("[aA]purinic", profile[3]))):
            currentDomain = "APE"
            
        elif ((re.search("[tT]ransposase", profile[1]) or re.search("[tT]ransposase", profile[3])) and not (re.search("DDE ", profile[1]) or re.search("DDE ", profile[3]))):
            currentDomain = "Tase"
        elif (re.search("DUF659", profile[1])):
            currentDomain = "Tase"
        elif (re.search("PF01695", profile[2])) or (re.search("PF02316", profile[2])) or (re.search("PF09039", profile[2])) or (re.search("PF04827", profile[2])) or (re.search("PF05699", profile[2])):
            currentDomain = "Tase"
            
        elif ((re.search("[tT]ransposase", profile[1]) or re.search("[tT]ransposase", profile[3])) and (re.search("DDE ", profile[1]) or re.search("DDE ", profile[3]))):
            currentDomain = "Tase*"
            
        elif ((re.search("[rR]eplication", profile[1]) or re.search("[rR]eplication", profile[3])) and (re.search("[pP]rotein", profile[1]) or re.search("[pP]rotein", profile[3])) and (re.search("A ", profile[1]) or re.search("A ", profile[3]))):
            currentDomain = "RPA"
        elif (re.search("[rR]epA ", profile[1]) or re.search("[rR]epA ", profile[3])):
            currentDomain = "RPA"
        elif (re.search("RPA", profile[1]) or re.search("RPA", profile[3])):
            currentDomain = "RPA"
            
        elif (re.search("[cC]-integrase", profile[1]) or re.search("[cC]-integrase", profile[3])):
            currentDomain = "C-INT"
            
        elif ((re.search("[pP]ackaging", profile[1]) or re.search("[pP]ackaging", profile[3])) and (re.search("ATPase", profile[1]) or re.search("ATPase", profile[3]))):
            currentDomain = "ATP"
            
        elif ((re.search("[cC]ysteine", profile[1]) or re.search("[cC]ysteine", profile[3])) and (re.search("[pP]rotease", profile[1]) or re.search("[pP]rotease", profile[3]))):
            currentDomain = "CYP"
        elif ((re.search("[cC]ysteine", profile[1]) or re.search("[cC]ysteine", profile[3])) and (re.search("[pP]eptidase", profile[1]) or re.search("[pP]eptidase", profile[3]))):
            currentDomain = "CYP"
        elif (re.search("[pP]eptidase_C", profile[1]) or re.search("[pP]eptidase_C", profile[3])):
            currentDomain = "CYP"
        elif (re.search("PF00559", profile[2])):
            currentDomain = "CYP"
            
        elif (re.search("[pP]ol\S*_*B", profile[1]) or re.search("[pP]ol\S*_*B", profile[3])):
            currentDomain = "POLB"
        elif ((re.search("[pP]olymerase", profile[1]) or re.search("[pP]olymerase", profile[3])) and (re.search("B ", profile[1]) or re.search("B ", profile[3]))):
            currentDomain = "POLB"
            
        elif (re.search("[hH]elicase", profile[1]) or re.search("[hH]elicase", profile[3])):
            currentDomain = "HEL"
            
        else :
            currentDomain = "OTHER"
        return currentDomain
    
    
    ## Replace the old profile name by accession number, name, domain and gather cut off
    #  
    # @param fout file handle
    # @param profile Profiles instance
    # @param currentDomain string
    #
    def _writeModifiedProfile(self, fout, profile, currentDomain):
        for i in xrange(0, len(profile), 1):
            if i != 1:
                fout.write(profile[i])
            else:
                fout.write("NAME  " + self.profile.accNumber + "_"\
                            + self.profile.name + "_"\
                             + currentDomain + "_"\
                              + str(self.profile.GA_cut_off) + "\n")
                
                
    ## Set input file name
    #  
    # @param inputFileName string
    #
    def setInputFile(self, inputFileName):
        self._inputFile = inputFileName
        
        
    ## Set output file name
    #
    # @param outputFileName string
    #
    def setOutputFile(self, outputFileName):
        self._outputFile = outputFileName   
        
        
    ## Read a profiles DB file, parse it and, write a new profiles DB with TE domain information and GA score cut_off placed side by side of the name
    #    
    def run( self ):
        LoggerFactory.setLevel(self._log, self._verbosity)
        fileIn = open( self._inputFile )
        fout = open( self._outputFile, "w" )
        profile = self.profile.readAndRetrieve( fileIn )
        while profile != None:
            currentDomain = self._searchCurrentDomain(profile)
            if currentDomain == "OTHER":
                self._log.warning(self.profile.accNumber + " " + self.profile.name + " has no associated domain")
            self._writeModifiedProfile(fout, profile, currentDomain)
            profile = self.profile.read( fileIn )
            
            
if __name__ == "__main__":
    i = ProfilesDB4Repet()
    i._setAttributesFromCmdLine()
    i.run()