annotate sampline.py @ 12:2f4ea569f048

Uploaded
author xuebing
date Sat, 10 Mar 2012 08:10:44 -0500
parents b7f1d9f8f3bc
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
1 #!/usr/bin/env python
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
2
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
3 """
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
4 Sampling random records from a file. Each record is defined by a fixed number of lines.
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
5
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
6 Usage: sampline.py [options]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
7
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
8 Options:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
9 -h, --help show this help message and exit
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
10 -r, --replacement Sampling with replacement
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
11 -i INPUT, --input=INPUT
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
12 Input file
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
13 -o OUTPUT, --output=OUTPUT
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
14 Output file
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
15 -k NSAMPLE, --nSample=NSAMPLE
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
16 (required) number of records to be sampled/output
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
17 -m RECSIZE, --recSize=RECSIZE
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
18 (default=1) number of lines spanned by each record
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
19 -n NSKIP, --nSkip=NSKIP
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
20 (default=0) number of comment lines to skip at the
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
21 beginning
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
22
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
23 example:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
24 python sampline.py -i test10000.fastq -o out.txt --nSample=5 --recSize=4 --nSkip=0 --replacement
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
25 """
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
26
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
27 import optparse, string, random,sys,math,itertools
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
28
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
29 assert sys.version_info[:2] >= ( 2, 4 )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
30
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
31 def main():
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
32
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
33 # Parse command line
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
34 parser = optparse.OptionParser( usage="%prog [options] " )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
35 parser.add_option( "-r", "--replacement", action="store_true", dest="replacement",default=False,
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
36 help="Sampling with replacement" )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
37 parser.add_option( "-i", "--input", dest="input", default=None,
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
38 help="Input file" )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
39 parser.add_option( "-o", "--output", dest="output", default=None,
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
40 help="Output file" )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
41 parser.add_option("-k","--nSample", type='int',dest="nSample",default=None,
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
42 help="(required) number of records to be sampled/output" )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
43 parser.add_option("-m","--recSize", type='int',dest="recSize",default=1,
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
44 help="(default=1) number of lines spanned by each record" )
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
45 parser.add_option("-n","--nSkip", type='int',dest="nSkip",default=0,
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
46 help="(default=0) number of comment lines to skip at the beginning")
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
47 options, args = parser.parse_args()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
48 #assert options.region in ( 'coding', 'utr3', 'utr5', 'transcribed' ), "Invalid region argument"
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
49
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
50 sampline(options.input,options.output,options.nSample,options.recSize,options.nSkip,options.replacement)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
51
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
52 def sample_wr(population, k):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
53 "Chooses k random elements (with replacement) from a population"
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
54 n = len(population)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
55 _random, _int = random.random, int # speed hack
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
56 return [_int(_random() * n) for i in itertools.repeat(None, k)]
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
57
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
58 # num of lines
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
59 def readinput(filename):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
60 try:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
61 f = open (filename)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
62 except:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
63 print >> sys.stderr, "can't open file "+str(filename)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
64 sys.exit(0)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
65
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
66 nline = 0
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
67 for line in f:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
68 nline = nline + 1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
69 f.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
70 return nline
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
71
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
72 def sampline(infile,outfile,nSample,recSize,nSkip,replacement):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
73 # sample nSample records from file
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
74 # each record contains recSize lines
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
75 # skip the top nSkip lines
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
76
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
77 nLine = readinput(infile)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
78 print 'num of lines in input: ',nLine
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
79 print 'avoid sampling the first ',nSkip,' lines'
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
80 print 'lines per record: ',recSize
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
81
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
82 if (nLine-nSkip) % recSize:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
83 print >> sys.stderr, "the number of lines is not dividable by record size!"
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
84 sys.exit(0)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
85
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
86 nTotalRecords = (nLine-nSkip) / recSize
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
87 print "total number of records: ",nTotalRecords
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
88
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
89 if replacement or nTotalRecords < nSample:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
90 sel = sample_wr(range(nTotalRecords),nSample)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
91 else:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
92 sel = random.sample(range(nTotalRecords),nSample)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
93
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
94 #print len(sel), sorted(sel)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
95
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
96 # output
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
97 try:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
98 fout = open (outfile,'w')
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
99 except:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
100 print >> sys.stderr, "can't open file "+str(outfile)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
101 sys.exit(0)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
102 fin = open(infile)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
103 n = 0 # index of line
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
104 rec = "" # to store all content of a record
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
105 nrepeat = 0 # number of times a record is sampled
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
106 curr_rec = -1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
107 for line in fin:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
108 if n < nSkip:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
109 n = n + 1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
110 fout.write(line)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
111 continue
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
112
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
113 if not (n-nSkip) % recSize:# a new record
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
114 # print the previous sampled record
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
115 for i in range(nrepeat):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
116 fout.write(rec)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
117 curr_rec = (n-nSkip)/recSize
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
118 nrepeat = sel.count(curr_rec)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
119 if nrepeat: # sampled
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
120 rec = line
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
121 #print curr_rec,nrepeat
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
122 elif (n-nSkip)/recSize == curr_rec:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
123 rec = rec + line
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
124 n = n + 1
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
125 # if the last record is selected
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
126 if curr_rec == nTotalRecords-1:
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
127 for i in range(nrepeat):
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
128 fout.write(rec)
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
129 fin.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
130 fout.close()
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
131
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
132
b7f1d9f8f3bc Uploaded
xuebing
parents:
diff changeset
133 if __name__ == "__main__": main()