comparison TEisotools-1.1.a/commons/core/coord/MatchUtils.py @ 13:feef9a0db09d draft

Uploaded
author urgi-team
date Wed, 20 Jul 2016 09:04:42 -0400
parents
children
comparison
equal deleted inserted replaced
12:22b0494ec883 13:feef9a0db09d
1 # Copyright INRA (Institut National de la Recherche Agronomique)
2 # http://www.inra.fr
3 # http://urgi.versailles.inra.fr
4 #
5 # This software is governed by the CeCILL license under French law and
6 # abiding by the rules of distribution of free software. You can use,
7 # modify and/ or redistribute the software under the terms of the CeCILL
8 # license as circulated by CEA, CNRS and INRIA at the following URL
9 # "http://www.cecill.info".
10 #
11 # As a counterpart to the access to the source code and rights to copy,
12 # modify and redistribute granted by the license, users are provided only
13 # with a limited warranty and the software's author, the holder of the
14 # economic rights, and the successive licensors have only limited
15 # liability.
16 #
17 # In this respect, the user's attention is drawn to the risks associated
18 # with loading, using, modifying and/or developing or reproducing the
19 # software by the user in light of its specific status of free software,
20 # that may mean that it is complicated to manipulate, and that also
21 # therefore means that it is reserved for developers and experienced
22 # professionals having in-depth computer knowledge. Users are therefore
23 # encouraged to load and test the software's suitability as regards their
24 # requirements in conditions enabling the security of their systems and/or
25 # data to be ensured and, more generally, to use and operate it in the
26 # same conditions as regards security.
27 #
28 # The fact that you are presently reading this means that you have had
29 # knowledge of the CeCILL license and that you accept its terms.
30
31 import os
32 import sys
33 import math
34 from commons.core.coord.Match import Match
35 from commons.core.utils.FileUtils import FileUtils
36 from commons.core.checker.RepetException import RepetException
37
38 # # Static methods for the manipulation of Match instances
39 #
40 class MatchUtils (object):
41
42 # # Return a list with Match instances from the given file
43 #
44 # @param inFile name of a file in the Match format
45 # @return a list of Match instances
46 #
47 @staticmethod
48 def getMatchListFromFile(inFile):
49 lMatchInstances = []
50 inFileHandler = open(inFile, "r")
51 while True:
52 line = inFileHandler.readline()
53 if line == "":
54 break
55 if line[0:10] == "query.name":
56 continue
57 m = Match()
58 m.setFromString(line)
59 lMatchInstances.append(m)
60 inFileHandler.close()
61 return lMatchInstances
62
63 # # Split a Match list in several Match lists according to the subject
64 #
65 # @param lMatches a list of Match instances
66 # @return a dictionary which keys are subject names and values Match lists
67 #
68 @staticmethod
69 def getDictOfListsWithSubjectAsKey(lMatches):
70 dSubject2MatchList = {}
71 for iMatch in lMatches:
72 if not dSubject2MatchList.has_key(iMatch.range_subject.seqname):
73 dSubject2MatchList[ iMatch.range_subject.seqname ] = []
74 dSubject2MatchList[ iMatch.range_subject.seqname ].append(iMatch)
75 return dSubject2MatchList
76
77 # # Split a Match list in several Match lists according to the query
78 #
79 # @param lMatches a list of Match instances
80 # @return a dictionary which keys are query names and values Match lists
81 #
82 @staticmethod
83 def getDictOfListsWithQueryAsKey (lMatches):
84 dQuery2MatchList = {}
85 for iMatch in lMatches:
86 if not dQuery2MatchList.has_key(iMatch.range_query.seqname):
87 dQuery2MatchList[ iMatch.range_query.seqname ] = []
88 dQuery2MatchList[ iMatch.range_query.seqname ].append(iMatch)
89 return dQuery2MatchList
90
91 # # Write Match instances contained in the given list
92 #
93 # @param lMatches a list of Match instances
94 # @param fileName name of the file to write the Match instances
95 # @param mode the open mode of the file ""w"" or ""a""
96 #
97 @staticmethod
98 def writeListInFile(lMatches, fileName, mode = "w", header = None):
99 fileHandler = open(fileName, mode)
100 if header:
101 fileHandler.write(header)
102 for iMatch in lMatches:
103 iMatch.write(fileHandler)
104 fileHandler.close()
105
106 # # Give path id list from a list of Match instances
107 #
108 # @param lMatch list of Match instances
109 # @return lId integer list
110 #
111 @staticmethod
112 def getIdListFromMatchList(lMatch):
113 lId = []
114 for iMatch in lMatch:
115 lId.append(iMatch.id)
116 return lId
117
118 # # Remove duplicated matches in a match list
119 # ## replace old PyRepet.MatchDB.rmvDoublons()
120 #
121 # @param lMatch list of Match instances
122 # @return lMatchesUniq match unique list
123 #
124 @staticmethod
125 def rmvDuplicateMatches(lMatch):
126 lMatchesUniq = []
127 for match in lMatch:
128 if len(lMatchesUniq) == 0:
129 lMatchesUniq.append(match)
130 else:
131 nbDoublons = 0
132 for m in lMatchesUniq:
133 if match.isDoublonWith(m):
134 nbDoublons += 1
135 if nbDoublons == 0:
136 lMatchesUniq.append(match)
137
138 for match1 in lMatchesUniq:
139 for match2 in lMatchesUniq:
140 if match1.id != match2.id:
141 if match1.isDoublonWith(match2):
142 raise RepetException ("*** Error: doublon not removed")
143 return lMatchesUniq
144
145 # # Return the list of queries 'included' in subjects when two different databanks are used.
146 # #replace old pyRepet.MatchDB.filterDiffQrySbj()
147 #
148 # @param iBioseqDB bioseqDB databank of queries
149 # @param thresIdentity float identity threshold
150 # @param thresLength float length threshold
151 # @param verbose int verbosity
152 #
153 # @return lMatches match list to keep according to length and identity thresholds
154 # TODO: don't take into account match for sequence against itself. To do ?
155 @staticmethod
156 def filterDiffQrySbj(iBioseqDB, matchFile, thresIdentity = 0.95, thresLength = 0.98, verbose = 0):
157 if verbose > 0:
158 print "filtering matches (id>=%.2f,qlgth>=%.2f)..." % (thresIdentity, thresLength); sys.stdout.flush()
159
160 thresIdentityPerc = math.floor(thresIdentity * 100)
161 lQryToKeep = []
162 dQry2Matches = MatchUtils.getDictOfListsWithQueryAsKey(MatchUtils.getMatchListFromFile(matchFile))
163
164 for seqH in iBioseqDB.idx.keys():
165 # keep it if it has no match
166 if not dQry2Matches.has_key(seqH):
167 if seqH not in lQryToKeep:
168 lQryToKeep.append(seqH)
169 else:
170 isConditionsMet = False
171 for match in dQry2Matches[ seqH ]:
172 # check if they are above the thresholds
173 if match.identity >= thresIdentityPerc and match.query_length_perc >= thresLength:
174 isConditionsMet = True
175 break
176 if not isConditionsMet and seqH not in lQryToKeep:
177 lQryToKeep.append(seqH)
178 return lQryToKeep
179
180 # # Count the number of distinct matches involved in at least one match above the thresholds.
181 # #replace old pyRepet.coord.MatchDB.getNbDistinctSbjWithThres() and pyRepet.coord.MatchDB.getNbDistinctSbjWithThres()
182 #
183 # @param thresIdentity float identity threshold
184 # @param thresLength float length threshold
185 #
186 @staticmethod
187 def getNbDistinctSequencesInsideMatchesWithThresh(lMatches, thresIdentity = 0.95, thresLength = 0.98, whatToCount = "query"):
188 thresIdentityPerc = math.floor(thresIdentity * 100)
189 countSbj = 0
190 if whatToCount.lower() == "query":
191 dMatches = MatchUtils.getDictOfListsWithQueryAsKey(lMatches)
192 else:
193 dMatches = MatchUtils.getDictOfListsWithSubjectAsKey(lMatches)
194
195 for qry in dMatches.keys():
196 countMatch = 0
197 for match in dMatches[ qry ]:
198
199 if match.identity >= thresIdentityPerc and getattr(match, whatToCount.lower() + "_length_perc") >= thresLength:
200 countMatch += 1
201 if countMatch > 0:
202 countSbj += 1
203 return countSbj
204
205 # # Convert a 'match' file (output from Matcher) into an 'align' file
206 # # replace old parser.tab2align
207 #
208 # @param inFileName a string input file name
209 #
210 @staticmethod
211 def convertMatchFileToAlignFile(inFileName):
212 basename = os.path.splitext(inFileName)[0]
213 outFileName = "%s.align" % basename
214 outFile = open(outFileName, "w")
215
216 lMatches = MatchUtils.getMatchListFromFile(inFileName)
217
218 for match in lMatches:
219 string = "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (match.getQueryName(), match.getQueryStart(), match.getQueryEnd(), match.getSubjectName(), match.getSubjectStart(), match.getSubjectEnd(), match.getEvalue(), match.getScore(), match.getIdentity())
220 outFile.write(string)
221
222 outFile.close()
223
224 # # Convert a 'match' file (output from Matcher) into an 'abc' file (MCL input file)
225 # for each match, compute coverage on the smallest seq using this expression
226 # Min (|QMatchEnd - QMatchStart + 1|,|SMatchEnd - SMatchStart + 1|)
227 #coverage = -----------------------------------------------------------------
228 # Min (QLength,SLength)
229 #
230 # Use this coverage for arc value. Matches of a sequence on itself are filtered.
231 #
232 # @param matchFileName string input match file name
233 # @param outFileName string output abc file name
234 # @param coverage float coverage filter threshold
235 #
236 @staticmethod
237 def convertMatchFileIntoABCFileWithCoverageComputeOnSmallestSeq(matchFileName, outFileName, coverageThreshold = 0):
238 sSeqNames = set()
239 sSeqNamesUsed = set()
240 with open(outFileName, "w") as outF:
241 with open(matchFileName) as inF:
242 inF.readline()
243 for inLine in inF:
244 splittedLine = inLine.split("\t")
245 qName = splittedLine[0]
246 sName = splittedLine[6]
247 sSeqNames.add(qName)
248 sSeqNames.add(sName)
249 if qName != sName:
250 matchMin = min(abs(float(splittedLine[2]) - float(splittedLine[1]) + 1), abs(float(splittedLine[8]) - float(splittedLine[7]) + 1))
251 seqMin = min(int(splittedLine[3]) / float(splittedLine[4]), int(splittedLine[9]) / float(splittedLine[10]))
252 coverage = round(matchMin / seqMin, 2)
253 if coverage >= coverageThreshold:
254 outF.write("%s\n" % "\t".join([qName, sName, str(coverage)]))
255 sSeqNamesUsed.add(qName)
256 sSeqNamesUsed.add(sName)
257
258 with open("%s.unused" % outFileName, "w") as outSeqNameOfUnusedMatchesF:
259 for seqName in sorted(sSeqNames - sSeqNamesUsed):
260 outSeqNameOfUnusedMatchesF.write("%s\n" % seqName)
261
262 if FileUtils.isEmpty(outFileName):
263 print "WARNING: '%s' is empty." % outFileName
264
265 # # Convert a 'match' file (output from Matcher) into an 'abc' file (MCL input file)
266 # Use coverage on query for arc value
267 #
268 # @param matchFileName string input match file name
269 # @param outFileName string output abc file name
270 # @param coverage float query coverage filter threshold
271 #
272 @staticmethod
273 def convertMatchFileIntoABCFileOnQueryCoverage(matchFileName, outFileName, coverage = 0):
274 with open(matchFileName) as inF:
275 with open(outFileName, "w") as outF:
276 inF.readline()
277 for inLine in inF:
278 splittedLine = inLine.split("\t")
279 if float(splittedLine[4]) >= coverage:
280 outLine = "\t".join([splittedLine[0], splittedLine[6], splittedLine[4]])
281 outLine += "\n"
282 outF.write(outLine)
283
284 if FileUtils.isEmpty(outFileName):
285 print "WARNING: '%s' is empty." % outFileName
286
287 # # Adapt the path IDs as the input file is the concatenation of several 'Match' files, and remove the extra header lines.
288 # # replace old parser.tabnum2id
289 #
290 # @param fileName a string input file name
291 # @param outputFileName a string output file name (optional)
292 #
293 @staticmethod
294 def generateMatchFileWithNewPathId(fileName, outputFileName = None):
295 if outputFileName is None:
296 outFile = open(fileName, "w")
297 else:
298 outFile = open(outputFileName, "w")
299 outFile.write("query.name\tquery.start\tquery.end\tquery.length\tquery.length.%\tmatch.length.%\tsubject.name\tsubject.start\tsubject.end\tsubject.length\tsubject.length.%\tE.value\tScore\tIdentity\tpath\n")
300
301 lMatches = MatchUtils.getMatchListFromFile(fileName)
302 count = 1
303 dMatchKeyIdcount = {}
304
305 for match in lMatches:
306 key_id = str(match.getIdentifier()) + "-" + match.getQueryName() + "-" + match.getSubjectName()
307 if not key_id in dMatchKeyIdcount.keys():
308 newPath = count
309 count += 1
310 dMatchKeyIdcount[ key_id ] = newPath
311 else:
312 newPath = dMatchKeyIdcount[ key_id ]
313
314 match.id = newPath
315 outFile.write(match.toString() + "\n")
316 outFile.close()