annotate mytools/random_interval.py @ 7:f0dc65e7f6c0

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:59:07 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
1 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
2 simulate a random interval set that mimics the size and strand of a reference set
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
3 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
4
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
5 def inferSizeFromRefBed(filename):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
6 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
7 read reference interval set, get chrom size information
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
8 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
9 chrSize = {}
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
10 f = open(filename)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
11 for line in f:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
12 flds = line.strip().split('\t')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
13 if not chrSize.has_key(flds[0]):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
14 chrSize[flds[0]] = int(flds[2])
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
15 elif chrSize[flds[0]] < int(flds[2]):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
16 chrSize[flds[0]] = int(flds[2])
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
17 f.close()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
18 return chrSize
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
19
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
20 def getChrSize(filename):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
21 chrSize = {}
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
22 f = open(filename)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
23 for line in f:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
24 flds = line.strip().split()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
25 if len(flds) >1:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
26 chrSize[flds[0]] = int(flds[1])
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
27 f.close()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
28 return chrSize
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
29
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
30 def makeWeightedChrom(chrSize):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
31 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
32 make a list of chr_id, the freq is proportional to its length
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
33 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
34
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
35 genome_len = 0
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
36
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
37 for chrom in chrSize:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
38 chrom_len = chrSize[chrom]
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
39 genome_len += chrom_len
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
40
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
41 weighted_chrom = []
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
42 for chrom in chrSize:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
43 weight = int(round(1000*float(chrSize[chrom])/genome_len))
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
44 weighted_chrom += [chrom]*weight
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
45
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
46 return weighted_chrom
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
47
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
48 def randomIntervalWithinChrom(infile,outfile,chrSize):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
49 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
50 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
51 fin = open(infile)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
52 fout = open(outfile,'w')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
53 n = 0
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
54 for line in fin:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
55 n = n + 1
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
56 flds = line.strip().split('\t')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
57 interval_size = int(flds[2]) - int(flds[1])
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
58 flds[1] = str(random.randint(0,chrSize[flds[0]]-interval_size))
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
59 flds[2] = str(int(flds[1])+interval_size)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
60 fout.write('\t'.join(flds)+'\n')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
61 fin.close()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
62 fout.close()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
63
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
64 def randomIntervalAcrossChrom(infile,outfile,chrSize,weighted_chrom):
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
65 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
66 '''
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
67 fin = open(infile)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
68 fout = open(outfile,'w')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
69 n = 0
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
70 for line in fin:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
71 n = n + 1
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
72 flds = line.strip().split('\t')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
73 interval_size = int(flds[2]) - int(flds[1])
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
74 # find a random chrom
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
75 flds[0] = weighted_chrom[random.randint(0, len(weighted_chrom) - 1)]
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
76 flds[1] = str(random.randint(0,chrSize[flds[0]]-interval_size))
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
77 flds[2] = str(int(flds[1])+interval_size)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
78 fout.write('\t'.join(flds)+'\n')
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
79 fin.close()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
80 fout.close()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
81
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
82 import sys,random
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
83 def main():
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
84 # python random_interval.py test100.bed testout.bed across human.hg18.genome
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
85
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
86 reference_interval_file = sys.argv[1]
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
87 output_file = sys.argv[2]
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
88 across_or_within_chrom = sys.argv[3] # within or across
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
89 chrom_size_file = sys.argv[4]
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
90 chrSize = getChrSize(chrom_size_file)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
91 print chrSize.keys()
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
92 if across_or_within_chrom == 'within':
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
93 randomIntervalWithinChrom(reference_interval_file,output_file,chrSize)
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
94 else:
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
95 randomIntervalAcrossChrom(reference_interval_file,output_file,chrSize,makeWeightedChrom(chrSize))
f0dc65e7f6c0 Uploaded
xuebing
parents:
diff changeset
96 main()