Mercurial > repos > yufei-luo > s_mart

comparison SMART/Java/Python/getRandomRegions.py @ 46:169d364ddd91

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Mon, 30 Sep 2013 03:19:26 -0400
parents 769e306b7933
children
comparison
equal deleted inserted replaced
45:e454402ba9d9 46:169d364ddd91
42 repetitions = 100 42 repetitions = 100
43 43
44 44
45 class RandomRegionsGenerator(object): 45 class RandomRegionsGenerator(object):
46 46
47 def __init__(self, verbosity): 47 def __init__(self, verbosity):
48 self.verbosity = verbosity 48 self.verbosity = verbosity
49 self.strands = False 49 self.strands = False
50 self.distribution = "uniform" 50 self.distribution = "uniform"
51 self.transcripts = None 51 self.transcripts = None
52 self.sequenceParser = None 52 self.sequenceParser = None
53 random.seed() 53 random.seed()
54 54
55 55
56 def setInput(self, fileName): 56 def setInput(self, fileName):
57 self.sequenceParser = FastaParser(fileName, self.verbosity) 57 self.sequenceParser = FastaParser(fileName, self.verbosity)
58 58
59 59
60 def setGenomeSize(self, size): 60 def setGenomeSize(self, size):
61 self.genomeSize = size 61 self.genomeSize = size
62 62
63 63
64 def setChromosomeName(self, name): 64 def setChromosomeName(self, name):
65 self.chromosomeName = name 65 self.chromosomeName = name
66 66
67 67
68 def setAnnotation(self, fileName, format): 68 def setAnnotation(self, fileName, format):
69 parser = TranscriptContainer(fileName, format, self.verbosity) 69 parser = TranscriptContainer(fileName, format, self.verbosity)
70 self.transcripts = [] 70 self.transcripts = []
71 for transcript in parser.getIterator(): 71 for transcript in parser.getIterator():
72 self.transcripts.append(transcript) 72 self.transcripts.append(transcript)
73 self.setNumber(len(self.transcripts)) 73 self.setNumber(len(self.transcripts))
74 self.setSize(0) 74 self.setSize(0)
75 75
76 76
77 def setOutputFile(self, fileName): 77 def setOutputFile(self, fileName):
78 self.outputFileName = fileName 78 self.outputFileName = fileName
79 79
80 80
81 def setSize(self, size): 81 def setSize(self, size):
82 self.minSize = size 82 self.minSize = size
83 self.maxSize = size 83 self.maxSize = size
84 84
85 85
86 def setMinSize(self, size): 86 def setMinSize(self, size):
87 self.minSize = size 87 self.minSize = size
88 88
89 89
90 def setMaxSize(self, size): 90 def setMaxSize(self, size):
91 self.maxSize = size 91 self.maxSize = size
92 92
93 93
94 def setNumber(self, number): 94 def setNumber(self, number):
95 self.number = number 95 self.number = number
96 96
97 97
98 def setStrands(self, strands): 98 def setStrands(self, strands):
99 self.strands = strands 99 self.strands = strands
100 100
101 101
102 def setMaxDistribution(self, maxElements): 102 def setMaxDistribution(self, maxElements):
103 if maxElements == None: 103 if maxElements == None:
104 return 104 return
105 self.maxElements = maxElements 105 self.maxElements = maxElements
106 self.distribution = "gaussian" 106 self.distribution = "gaussian"
107 107
108 108
109 def setDeviationDistribution(self, deviation): 109 def setDeviationDistribution(self, deviation):
110 if deviation == None: 110 if deviation == None:
111 return 111 return
112 self.deviation = deviation 112 self.deviation = deviation
113 self.distribution = "gaussian" 113 self.distribution = "gaussian"
114 114
115 115
116 def getSizes(self): 116 def getSizes(self):
117 if self.sequenceParser == None: 117 if self.sequenceParser == None:
118 self.chromosomes = [self.chromosomeName] 118 self.chromosomes = [self.chromosomeName]
119 self.sizes = {self.chromosomeName: self.genomeSize} 119 self.sizes = {self.chromosomeName: self.genomeSize}
120 self.cumulatedSize = self.genomeSize 120 self.cumulatedSize = self.genomeSize
121 self.cumulatedSizes = {self.chromosomeName: self.genomeSize} 121 self.cumulatedSizes = {self.chromosomeName: self.genomeSize}
122 return 122 return
123 self.chromosomes = self.sequenceParser.getRegions() 123 self.chromosomes = self.sequenceParser.getRegions()
124 self.sizes = {} 124 self.sizes = {}
125 self.cumulatedSize = 0 125 self.cumulatedSize = 0
126 self.cumulatedSizes = {} 126 self.cumulatedSizes = {}
127 for chromosome in self.chromosomes: 127 for chromosome in self.chromosomes:
128 self.sizes[chromosome] = self.sequenceParser.getSizeOfRegion(chromosome) 128 self.sizes[chromosome] = self.sequenceParser.getSizeOfRegion(chromosome)
129 self.cumulatedSize += self.sizes[chromosome] 129 self.cumulatedSize += self.sizes[chromosome]
130 self.cumulatedSizes[chromosome] = self.cumulatedSize 130 self.cumulatedSizes[chromosome] = self.cumulatedSize
131 131
132 132
133 def findPosition(self, size = None): 133 def findPosition(self, size = None):
134 if size == None: 134 if size == None:
135 size = random.randint(self.minSize, self.maxSize) 135 size = random.randint(self.minSize, self.maxSize)
136 integer = random.randint(0, self.cumulatedSize) 136 integer = random.randint(0, self.cumulatedSize)
137 for chromosome in self.chromosomes: 137 for chromosome in self.chromosomes:
138 if self.cumulatedSizes[chromosome] > integer: 138 if self.cumulatedSizes[chromosome] > integer:
139 break 139 break
140 start = random.randint(1, self.sizes[chromosome] - size) 140 start = random.randint(1, self.sizes[chromosome] - size)
141 return (chromosome, start, size) 141 return (chromosome, start, size)
142 142
143 143
144 def createTranscript(self, chromosome, start, size, strand, cpt): 144 def createTranscript(self, chromosome, start, size, strand, cpt):
145 transcript = Transcript() 145 transcript = Transcript()
146 transcript.setChromosome(chromosome) 146 transcript.setChromosome(chromosome)
147 transcript.setStart(start) 147 transcript.setEnd(start + size-1)
148 transcript.setEnd(start + size-1) 148 transcript.setStart(start)
149 transcript.setDirection(strand) 149 transcript.setDirection(strand)
150 transcript.setName("rand_%d" % (cpt)) 150 transcript.setName("rand_%d" % (cpt))
151 return transcript 151 return transcript
152 152
153 153
154 def moveTranscript(self, chromosome, start, transcript): 154 def moveTranscript(self, chromosome, start, transcript):
155 while transcript.getEnd() + start - transcript.getStart() > self.cumulatedSizes[chromosome]: 155 while transcript.getEnd() + start - transcript.getStart() > self.cumulatedSizes[chromosome]:
156 chromosome, start, size = self.findPosition(transcript.getEnd() - transcript.getStart()) 156 chromosome, start, size = self.findPosition(transcript.getEnd() - transcript.getStart())
157 transcript.setChromosome(chromosome) 157 newTranscript = Transcript()
158 oldStart, oldEnd = transcript.getStart(), transcript.getEnd() 158 newTranscript.setChromosome(chromosome)
159 if transcript.getNbExons() > 1: 159 newTranscript.tags = transcript.tags
160 for exon in transcript.getNbExons(): 160 if transcript.getNbExons() > 1:
161 oldExonStart, oldExonEnd = exon.getStart(), exon.getEnd() 161 for exon in transcript.getNbExons():
162 exon.setStart(oldExonStart + start - oldStart) 162 newExon = Interval()
163 exon.setEnd(oldExonEnd + start - oldStart) 163 newExon.setChromosome(chromosome)
164 transcript.setStart(start) 164 newExon.setEnd(exon.getEnd() + start - transcript.getStart())
165 transcript.setEnd(oldEnd + start - oldStart) 165 newExon.setStart(exon.getStart() + start - transcript.getStart())
166 return [transcript] 166 newTranscript.addExon(newExon)
167 167 newTranscript.setEnd(transcript.getEnd() + start - transcript.getStart())
168 168 newTranscript.setStart(start)
169 def createUniformCluster(self, chromosome, start, size, strand, cpt): 169 newTranscript.setDirection(transcript.getDirection())
170 transcript = self.createTranscript(chromosome, start, size, strand, cpt) 170 return [newTranscript]
171 return [transcript] 171
172 172
173 173 def createUniformCluster(self, chromosome, start, size, strand, cpt):
174 def findNbTranscripts(self, cpt): 174 transcript = self.createTranscript(chromosome, start, size, strand, cpt)
175 return min(int(round(math.exp(random.random() * math.log(self.maxElements)))), self.number - cpt + 1) 175 return [transcript]
176 176
177 177
178 def getDev(self): 178 def findNbTranscripts(self, cpt):
179 deviation = 0.0 179 return min(int(round(math.exp(random.random() * math.log(self.maxElements)))), self.number - cpt + 1)
180 for j in range(repetitions): 180
181 deviation += random.randint(-self.deviation, self.deviation) 181
182 deviation /= repetitions 182 def getDev(self):
183 deviation = int(round(deviation)) 183 deviation = 0.0
184 return deviation 184 for j in range(repetitions):
185 185 deviation += random.randint(-self.deviation, self.deviation)
186 186 deviation /= repetitions
187 def createGaussianCluster(self, chromosome, start, size, strand, cpt): 187 deviation = int(round(deviation))
188 transcripts = [] 188 return deviation
189 nbTranscripts = self.findNbTranscripts(cpt) 189
190 for i in range(nbTranscripts): 190
191 transcript = self.createTranscript(chromosome, start + self.getDev(), size + self.getDev(), strand, cpt + i) 191 def createGaussianCluster(self, chromosome, start, size, strand, cpt):
192 transcripts.append(transcript) 192 transcripts = []
193 return transcripts 193 nbTranscripts = self.findNbTranscripts(cpt)
194 194 for i in range(nbTranscripts):
195 195 transcript = self.createTranscript(chromosome, start + self.getDev(), size + self.getDev(), strand, cpt + i)
196 def writeRegions(self): 196 transcripts.append(transcript)
197 writer = Gff3Writer(self.outputFileName, self.verbosity) 197 return transcripts
198 outputFile = open(self.outputFileName, "w") 198
199 progress = Progress(self.number, "Writing to %s" % (self.outputFileName), self.verbosity) 199
200 i = 0 200 def writeRegions(self):
201 while i < self.number: 201 writer = Gff3Writer(self.outputFileName, self.verbosity)
202 chromosome, start, size = self.findPosition() 202 outputFile = open(self.outputFileName, "w")
203 strand = random.choice([-1, 1]) if self.strands else 1 203 progress = Progress(self.number, "Writing to %s" % (self.outputFileName), self.verbosity)
204 if self.transcripts != None: 204 i = 0
205 transcripts = self.moveTranscript(chromosome, start, self.transcripts[i]) 205 while i < self.number:
206 elif self.distribution == "uniform": 206 chromosome, start, size = self.findPosition()
207 transcripts = self.createUniformCluster(chromosome, start, size, strand, i+1) 207 strand = random.choice([-1, 1]) if self.strands else 1
208 else: 208 if self.transcripts != None:
209 transcripts = self.createGaussianCluster(chromosome, start, size, strand, i+1) 209 transcripts = self.moveTranscript(chromosome, start, self.transcripts[i])
210 for transcript in transcripts: 210 elif self.distribution == "uniform":
211 writer.addTranscript(transcript) 211 transcripts = self.createUniformCluster(chromosome, start, size, strand, i+1)
212 i += 1 212 else:
213 progress.inc() 213 transcripts = self.createGaussianCluster(chromosome, start, size, strand, i+1)
214 progress.done() 214 for transcript in transcripts:
215 outputFile.close() 215 writer.addTranscript(transcript)
216 writer.write() 216 i += 1
217 writer.close() 217 progress.inc()
218 218 progress.done()
219 219 outputFile.close()
220 def run(self): 220 writer.write()
221 self.getSizes() 221 writer.close()
222 self.writeRegions() 222
223
224 def run(self):
225 self.getSizes()
226 self.writeRegions()
223 227
224 228
225 if __name__ == "__main__": 229 if __name__ == "__main__":
226 230
227 # parse command line 231 # parse command line
228 description = "Get Random Regions v1.0.2: Get some random coordinates on a genome. May use uniform or gaussian distribution (in gaussion distribution, # of element per cluster follows a power law). [Category: Other]" 232 description = "Get Random Regions v1.0.2: Get some random coordinates on a genome. May use uniform or gaussian distribution (in gaussion distribution, # of element per cluster follows a power law). [Category: Other]"
229 233
230 parser = OptionParser(description = description) 234 parser = OptionParser(description = description)
231 parser.add_option("-r", "--reference", dest="reference", action="store", default=None, type="string", help="file that contains the sequences [format: file in FASTA format]") 235 parser.add_option("-r", "--reference", dest="reference", action="store", default=None, type="string", help="file that contains the sequences [format: file in FASTA format]")
232 parser.add_option("-S", "--referenceSize", dest="referenceSize", action="store", default=None, type="int", help="size of the chromosome (when no reference is given) [format: int]") 236 parser.add_option("-S", "--referenceSize", dest="referenceSize", action="store", default=None, type="int", help="size of the chromosome (when no reference is given) [format: int]")
233 parser.add_option("-c", "--chromosome", dest="chromosome", action="store", default=None, type="string", help="name of the chromosome (when no reference is given) [format: string]") 237 parser.add_option("-c", "--chromosome", dest="chromosome", action="store", default=None, type="string", help="name of the chromosome (when no reference is given) [format: string]")
234 parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [compulsory] [format: output file in FASTA format]") 238 parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [compulsory] [format: output file in FASTA format]")
235 parser.add_option("-i", "--input", dest="inputFileName", action="store", default=None, type="string", help="optional file containing regions to shuffle [format: file in transcript format given by -f]") 239 parser.add_option("-i", "--input", dest="inputFileName", action="store", default=None, type="string", help="optional file containing regions to shuffle [format: file in transcript format given by -f]")
236 parser.add_option("-f", "--format", dest="format", action="store", default=None, type="string", help="format of the previous file [format: transcript file format]") 240 parser.add_option("-f", "--format", dest="format", action="store", default=None, type="string", help="format of the previous file [format: transcript file format]")
237 parser.add_option("-s", "--size", dest="size", action="store", default=None, type="int", help="size of the regions (if no region set is provided) [format: int]") 241 parser.add_option("-s", "--size", dest="size", action="store", default=None, type="int", help="size of the regions (if no region set is provided) [format: int]")
238 parser.add_option("-z", "--minSize", dest="minSize", action="store", default=None, type="int", help="minimum size of the regions (if no region set nor a fixed size are provided) [format: int]") 242 parser.add_option("-z", "--minSize", dest="minSize", action="store", default=None, type="int", help="minimum size of the regions (if no region set nor a fixed size are provided) [format: int]")
239 parser.add_option("-Z", "--maxSize", dest="maxSize", action="store", default=None, type="int", help="maximum size of the regions (if no region set nor a fixed size are provided) [format: int]") 243 parser.add_option("-Z", "--maxSize", dest="maxSize", action="store", default=None, type="int", help="maximum size of the regions (if no region set nor a fixed size are provided) [format: int]")
240 parser.add_option("-n", "--number", dest="number", action="store", default=None, type="int", help="number of regions (if no region set is provided) [format: int]") 244 parser.add_option("-n", "--number", dest="number", action="store", default=None, type="int", help="number of regions (if no region set is provided) [format: int]")
241 parser.add_option("-t", "--strands", dest="strands", action="store_true", default=False, help="use both strands (if no region set is provided) [format: boolean]") 245 parser.add_option("-t", "--strands", dest="strands", action="store_true", default=False, help="use both strands (if no region set is provided) [format: boolean]")
242 parser.add_option("-m", "--max", dest="max", action="store", default=None, type="int", help="max. # reads in a cluster (for Gaussian dist.) [format: int]") 246 parser.add_option("-m", "--max", dest="max", action="store", default=None, type="int", help="max. # reads in a cluster (for Gaussian dist.) [format: int]")
243 parser.add_option("-d", "--deviation", dest="deviation", action="store", default=None, type="int", help="deviation around the center of the cluster (for Gaussian dist.) [format: int]") 247 parser.add_option("-d", "--deviation", dest="deviation", action="store", default=None, type="int", help="deviation around the center of the cluster (for Gaussian dist.) [format: int]")
244 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]") 248 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]")
245 (options, args) = parser.parse_args() 249 (options, args) = parser.parse_args()
246 250
247 rrg = RandomRegionsGenerator(options.verbosity) 251 rrg = RandomRegionsGenerator(options.verbosity)
248 if options.reference == None: 252 if options.reference == None:
249 rrg.setGenomeSize(options.referenceSize) 253 rrg.setGenomeSize(options.referenceSize)
250 rrg.setChromosomeName(options.chromosome) 254 rrg.setChromosomeName(options.chromosome)
251 else: 255 else:
252 rrg.setInput(options.reference) 256 rrg.setInput(options.reference)
253 rrg.setOutputFile(options.outputFileName) 257 rrg.setOutputFile(options.outputFileName)
254 if options.inputFileName == None: 258 if options.inputFileName == None:
255 if options.size != None: 259 if options.size != None:
256 rrg.setSize(options.size) 260 rrg.setSize(options.size)
257 else: 261 else:
258 rrg.setMinSize(options.minSize) 262 rrg.setMinSize(options.minSize)
259 rrg.setMaxSize(options.maxSize) 263 rrg.setMaxSize(options.maxSize)
260 rrg.setNumber(options.number) 264 rrg.setNumber(options.number)
261 rrg.setStrands(options.strands) 265 rrg.setStrands(options.strands)
262 else: 266 else:
263 rrg.setAnnotation(options.inputFileName, options.format) 267 rrg.setAnnotation(options.inputFileName, options.format)
264 rrg.setMaxDistribution(options.max) 268 rrg.setMaxDistribution(options.max)
265 rrg.setDeviationDistribution(options.deviation) 269 rrg.setDeviationDistribution(options.deviation)
266 rrg.run() 270 rrg.run()
267 271