|
18
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 # Copyright INRA (Institut National de la Recherche Agronomique)
|
|
|
4 # http://www.inra.fr
|
|
|
5 # http://urgi.versailles.inra.fr
|
|
|
6 #
|
|
|
7 # This software is governed by the CeCILL license under French law and
|
|
|
8 # abiding by the rules of distribution of free software. You can use,
|
|
|
9 # modify and/ or redistribute the software under the terms of the CeCILL
|
|
|
10 # license as circulated by CEA, CNRS and INRIA at the following URL
|
|
|
11 # "http://www.cecill.info".
|
|
|
12 #
|
|
|
13 # As a counterpart to the access to the source code and rights to copy,
|
|
|
14 # modify and redistribute granted by the license, users are provided only
|
|
|
15 # with a limited warranty and the software's author, the holder of the
|
|
|
16 # economic rights, and the successive licensors have only limited
|
|
|
17 # liability.
|
|
|
18 #
|
|
|
19 # In this respect, the user's attention is drawn to the risks associated
|
|
|
20 # with loading, using, modifying and/or developing or reproducing the
|
|
|
21 # software by the user in light of its specific status of free software,
|
|
|
22 # that may mean that it is complicated to manipulate, and that also
|
|
|
23 # therefore means that it is reserved for developers and experienced
|
|
|
24 # professionals having in-depth computer knowledge. Users are therefore
|
|
|
25 # encouraged to load and test the software's suitability as regards their
|
|
|
26 # requirements in conditions enabling the security of their systems and/or
|
|
|
27 # data to be ensured and, more generally, to use and operate it in the
|
|
|
28 # same conditions as regards security.
|
|
|
29 #
|
|
|
30 # The fact that you are presently reading this means that you have had
|
|
|
31 # knowledge of the CeCILL license and that you accept its terms.
|
|
|
32
|
|
|
33 import os
|
|
|
34 import sys
|
|
|
35 import getopt
|
|
|
36
|
|
|
37 ##@file
|
|
|
38 # usage: dbConsensus.py [ options ]
|
|
|
39 # options:
|
|
|
40 # -h: this help
|
|
|
41 # -i: name of the input file (format=aligned fasta)
|
|
|
42 # -n: minimum number of nucleotides in a column to edit a consensus (default=1)
|
|
|
43 # -p: minimum proportion for the major nucleotide to be used, otherwise add 'N' (default=0.0)
|
|
|
44 # -o: name of the output file (default=inFileName+'.cons')
|
|
|
45 # -v: verbose (default=0/1/2)
|
|
|
46
|
|
|
47 if not os.environ.has_key( "REPET_PATH" ):
|
|
|
48 print "ERROR: no environment variable REPET_PATH"
|
|
|
49 sys.exit(1)
|
|
|
50 sys.path.append( os.environ["REPET_PATH"] )
|
|
|
51
|
|
|
52 import commons.core.seq.AlignedBioseqDB
|
|
|
53
|
|
|
54
|
|
|
55 def help():
|
|
|
56 """
|
|
|
57 Give the list of the command-line options.
|
|
|
58 """
|
|
|
59 print
|
|
|
60 print "usage:",sys.argv[0]," [ options ]"
|
|
|
61 print "options:"
|
|
|
62 print " -h: this help"
|
|
|
63 print " -i: name of the input file (format=aligned fasta)"
|
|
|
64 print " -n: minimum number of nucleotides in a column to edit a consensus (default=1)"
|
|
|
65 print " -p: minimum proportion for the major nucleotide to be used, otherwise add 'N' (default=0.0)"
|
|
|
66 print " -o: name of the output file (default=inFileName+'.cons')"
|
|
|
67 print " -H: format the header with pyramid and piles informations (SATannot)"
|
|
|
68 print " -v: verbose (default=0/1/2)"
|
|
|
69 print
|
|
|
70
|
|
|
71
|
|
|
72 def main():
|
|
|
73
|
|
|
74 inFileName = ""
|
|
|
75 minNbNt = 1
|
|
|
76 minPropNt = 0.0
|
|
|
77 outFileName = ""
|
|
|
78 header_SATannot = False
|
|
|
79 verbose = 0
|
|
|
80
|
|
|
81 try:
|
|
|
82 opts, args = getopt.getopt(sys.argv[1:],"hi:n:p:o:v:H")
|
|
|
83 except getopt.GetoptError, err:
|
|
|
84 print str(err); help(); sys.exit(1)
|
|
|
85 for o,a in opts:
|
|
|
86 if o == "-h":
|
|
|
87 help()
|
|
|
88 sys.exit(0)
|
|
|
89 elif o == "-i":
|
|
|
90 inFileName = a
|
|
|
91 elif o == "-n":
|
|
|
92 minNbNt = int(a)
|
|
|
93 elif o == "-p":
|
|
|
94 minPropNt = float(a)
|
|
|
95 elif o == "-o":
|
|
|
96 outFileName = a
|
|
|
97 elif o == "-H":
|
|
|
98 header_SATannot = True
|
|
|
99 elif o == "-v":
|
|
|
100 verbose = int(a)
|
|
|
101
|
|
|
102 if inFileName == "":
|
|
|
103 print "ERROR: missing input file name"
|
|
|
104 help()
|
|
|
105 sys.exit(1)
|
|
|
106
|
|
|
107 if verbose > 0:
|
|
|
108 print "START %s" % (sys.argv[0].split("/")[-1])
|
|
|
109 sys.stdout.flush()
|
|
|
110
|
|
|
111 alnDB = commons.core.seq.AlignedBioseqDB.AlignedBioseqDB( inFileName )
|
|
|
112
|
|
|
113 if alnDB.getSize() < minNbNt:
|
|
|
114 print "WARNING: not enough sequences (<%i)" % ( minNbNt )
|
|
|
115
|
|
|
116 else:
|
|
|
117 consensus = alnDB.getConsensus( minNbNt, minPropNt, verbose, header_SATannot)
|
|
|
118 if consensus != None:
|
|
|
119 consensus.upCase()
|
|
|
120 if outFileName == "":
|
|
|
121 outFileName = "%s.cons" % ( inFileName )
|
|
|
122 outFile = open( outFileName, "w" )
|
|
|
123 consensus.write( outFile )
|
|
|
124 outFile.close()
|
|
|
125
|
|
|
126 if verbose > 0:
|
|
|
127 print "END %s" % (sys.argv[0].split("/")[-1])
|
|
|
128 sys.stdout.flush()
|
|
|
129
|
|
|
130 return 0
|
|
|
131
|
|
|
132 if __name__ == "__main__":
|
|
|
133 main ()
|
