-
Notifications
You must be signed in to change notification settings - Fork 13
/
Copy pathbs3-align.py
executable file
·556 lines (442 loc) · 28.4 KB
/
bs3-align.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
#!/usr/bin/env python
from itertools import izip
from optparse import OptionParser, OptionGroup
import re
import tempfile
from bs_align import output
from bs_align.bs_rrbs import *
import os
import glob
import subprocess
import marshal
import pickle
if __name__ == '__main__':
parser = OptionParser(usage="Usage: %prog {-i <single> | -1 <mate1> -2 <mate2>} -g <genome.fa> [options]")
# option group 1
opt_group = OptionGroup(parser, "For single end reads")
opt_group.add_option("-i", "--input", type="string", dest="infilename",help="Input read file (FORMAT: fasta, fastq). Ex: read.fa or read.fa.gz", metavar="INFILE")
parser.add_option_group(opt_group)
# option group 2
opt_group = OptionGroup(parser, "For pair end reads")
opt_group.add_option("-1", "--input_1", type="string", dest="infilename_1",help="Input read file, mate 1 (FORMAT: fasta, fastq)", metavar="FILE")
opt_group.add_option("-2", "--input_2", type="string", dest="infilename_2",help="Input read file, mate 2 (FORMAT: fasta, fastq)", metavar="FILE")
opt_group.add_option("-I", "--minins",type = "int",dest = "min_insert_size", help="The minimum insert size for valid paired-end alignments [Default: %default]", default = 100)
opt_group.add_option("-X", "--maxins",type = "int",dest = "max_insert_size", help="The maximum insert size for valid paired-end alignments [Default: %default]", default = 1000)
parser.add_option_group(opt_group)
# option group 3
#opt_group = OptionGroup(parser, "Reduced Representation Bisulfite Sequencing Options")
#opt_group.add_option("-r", "--rrbs", action="store_true", dest="rrbs", default = False, help = 'Map reads to the Reduced Representation genome')
#opt_group.add_option("-c", "--cut-site", type="string",dest="cut_format", help="Cutting sites of restriction enzyme. Ex: MspI(C-CGG), Mael:(C-TAG), double-enzyme MspI&Mael:(C-CGG,C-TAG). [Default: %default]", metavar="pattern", default = "C-CGG")
#opt_group.add_option("-L", "--low", type = "int", dest="rrbs_low_bound", help="Lower bound of fragment length (excluding C-CGG ends) [Default: %default]", default = 20)
#opt_group.add_option("-U", "--up", type = "int", dest="rrbs_up_bound", help="Upper bound of fragment length (excluding C-CGG ends) [Default: %default]", default = 500)
#parser.add_option_group(opt_group)
# option group 4
opt_group = OptionGroup(parser, "General options")
opt_group.add_option("-K", "--alignment_length", dest="K", help="Neglect the alignments with length below this value", type = "int", default = 50)
opt_group.add_option("-t", "--tag", type="string", dest="taginfo",help="[Y]es for undirectional lib, [N]o for directional [Default: %default]", metavar="TAG", default = 'N')
opt_group.add_option("-s","--start_base",type = "int",dest = "cutnumber1", help="The first cycle of the read to be mapped [Default: %default]", default = 1)
opt_group.add_option("-e","--end_base",type = "int",dest = "cutnumber2", help="The last cycle of the read to be mapped [Default: %default]", default = 200)
opt_group.add_option("-a", "--adapter", type="string", dest="adapter_file",help="Input text file of your adaptor sequences (to be trimmed from the 3'end of the reads, ). "
"Input one seq for dir. lib., twon seqs for undir. lib. One line per sequence. "
"Only the first 10bp will be used", metavar="FILE", default = '')
opt_group.add_option("--am",type = "int",dest = "adapter_mismatch", help="Number of mismatches allowed in adapter [Default: %default]", default = 0)
opt_group.add_option("-g", "--genome", type="string", dest="genome",help="Name of the reference genome (should be the same as \"-f\" in bs_seeker2-build.py ) [ex. chr21_hg18.fa]")
opt_group.add_option("-m", "--mismatches",type = "float", dest="no_mismatches",help="Number(>=1)/Percentage([0, 1)) of mismatches in one read. Ex: 4 (allow 8 mismatches) or 0.08 (allow 8% mismatches) [Default: %default]", default = 20)
#opt_group.add_option("--aligner", dest="aligner",help="Aligner program for short reads mapping: " + ', '.join(supported_aligners) + " [Default: %default]", metavar="ALIGNER", default = SNAP)
#opt_group.add_option("-p", "--path", dest="aligner_path", help="Path to the aligner program. Detected: " +' '*70+ '\t'.join(('%s: %s '+' '*70) % (al, aligner_path[al]) for al in sorted(supported_aligners)),
# metavar="PATH"
#)
reference_genome_path = reference_genome_path + '/../../../reference_genome'
opt_group.add_option("-d", "--db", type="string", dest="dbpath",help="Path to the reference genome library (generated in preprocessing genome) [Default: %default]" , metavar="DBPATH", default = reference_genome_path)
opt_group.add_option("-l", "--split_line",type = "int", dest="no_split",help="Number of lines per split (the read file will be split into small files for mapping. The result will be merged. [Default: %default]", default = 12800000, metavar="INT")
opt_group.add_option("-o", "--output", type="string", dest="outfilename",help="The name of output file [INFILE.bs(se|pe|rrbs)]", metavar="OUTFILE")
opt_group.add_option("-f", "--output-format", type="string", dest="output_format",help="Output format: "+', '.join(output.formats)+" [Default: %default]", metavar="FORMAT", default = output.BAM)
opt_group.add_option("--no-header", action="store_true", dest="no_SAM_header",help="Suppress SAM header lines [Default: %default]", default = False)
try:
opt_group.add_option("--temp_dir", type="string", dest="temp_dir",help="The path to your temporary directory [Detected: %default]", metavar="PATH", default = os.environ["TMPDIR"])
except:
opt_group.add_option("--temp_dir", type="string", dest="temp_dir",help="The path to your temporary directory [Detected: %default]", metavar="PATH", default = tempfile.gettempdir())
opt_group.add_option("--XS",type = "string", dest="XS_filter",help="Filter definition for tag XS, format X,Y. X=0.8 and y=5 indicate that for one read, if #(mCH sites)/#(all CH sites)>0.8 and #(mCH sites)>5, then tag XS=1; or else tag XS=0. [Default: %default]", default = "0.5,5") # added by weilong
opt_group.add_option("-M", "--multiple-hit", metavar="FileName", type="string", dest="Output_multiple_hit", default = None, help = 'File to store reads with multiple-hits')
opt_group.add_option("-u", "--unmapped", metavar="FileName", type="string", dest="Output_unmapped_hit", default = None, help = 'File to store unmapped reads')
opt_group.add_option("-v", "--version", action="store_true", dest="version",help="show version of BS-Seeker3", metavar="version", default = False)
parser.add_option_group(opt_group)
# option group 5
opt_group = OptionGroup(parser, "Aligner Options",
"You may specify any additional options for the aligner. You just have to prefix them with " +
', '.join('%s for %s' % (aligner_options_prefixes[aligner], aligner) for aligner in supported_aligners)+
', and BS-Seeker3 will pass them on. Below are some of the options what would alter BSseeker peroformacnce using SNAP:\n\t-d maximum edit distance allowed per read or pair (default: 14)\n-n number of seeds to use per read\n-sc Seed coverage (i.e., readSize/seedSize). Floating point. Exclusive with -n. (default uses -n)\n-h maximum hits to consider per seed (default: 300)')
parser.add_option_group(opt_group)
# option group 6
opt_group = OptionGroup(parser, "Methylation Level Statistics Display Option")
opt_group.add_option("--qcf", type="string", dest="qc_f",help="Supply the length of the reads", default = '3')
parser.add_option_group(opt_group)
#----------------------------------------------------------------
# separate aligner options from BS Seeker options
aligner_options = {}
bs_seeker_options = []
i = 1
while i < len(sys.argv):
arg = sys.argv[i]
m = re.match(r'^%s' % '|'.join('(%s)'% aligner_options_prefixes[al] for al in supported_aligners), arg)
if m:
a_opt = arg.replace(m.group(0),'-',1)
aligner_options[a_opt] = []
while i + 1 < len(sys.argv) and sys.argv[i+1][0] != '-':
aligner_options[a_opt].append(sys.argv[i+1])
i += 1
if len(aligner_options[a_opt]) == 0: # if it is a key-only option
aligner_options[a_opt] = True
else:
bs_seeker_options.append(arg)
i += 1
(options, args) = parser.parse_args(args = bs_seeker_options)
qc_len = int(options.qc_f)
# if no options were given by the user, print help and exit
if len(sys.argv) == 1:
parser.print_help()
exit(0)
if options.version :
show_version()
exit (-1)
else :
show_version()
# check parameters
# input read files
if options.infilename and (options.infilename_1 or options.infilename_2):
error('-i and [-1|-2] options are exclusive. You should use only one of them.')
if not (options.infilename or (options.infilename_1 and options.infilename_2)):
error('You should set either -i or -1 and -2 options.')
# Calculate the length of read
if options.infilename :
read_file = options.infilename
elif options.infilename_1 :
read_file = options.infilename_1
else :
error('You should at least specify -i or -1 options.')
try :
if read_file.endswith(".gz") : # support input file ending with ".gz"
read_inf = gzip.open(read_file, "rb")
else :
read_inf=open(read_file,"r")
except IOError :
print "[Error] Cannot open input file : %s" % read_file
exit(-1)
oneline = read_inf.readline()
oneline = read_inf.readline() # get the second line
read_len = min(len(oneline), (options.cutnumber2-options.cutnumber1))
read_inf.close()
# mismatch allowed: bowtie 1,build-in parameter '-m'; bowtie 2, post-filter paramter
# mismatch should no greater than the read length
no_mismatches = float(options.no_mismatches)
if (no_mismatches < 1) :
int_no_mismatches=int(no_mismatches * read_len)
else :
int_no_mismatches=int(no_mismatches)
str_no_mismatches=str(options.no_mismatches) # pass to specific mode
# -t, directional / un-directional library
asktag=str(options.taginfo).upper()
if asktag not in 'YN':
error('-t option should be either Y or N, not %s' % asktag)
# -a
#if options.aligner not in supported_aligners:
# error('-a option should be: %s' % ' ,'.join(supported_aligners)+'.')
# path for aligner
#aligner_exec = os.path.expanduser( os.path.join(options.aligner_path or aligner_path[options.aligner], options.aligner) )
# -g
if options.genome is None:
error('-g is a required option')
genome = os.path.split(options.genome)[1]
genome_subdir = genome
# try to guess the location of the reference genome for RRBS
'''
if options.rrbs:
if options.rrbs_low_bound and options.rrbs_up_bound:
if options.cut_format == "C-CGG" :
genome_subdir += '_rrbs_%d_%d' % (options.rrbs_low_bound, options.rrbs_up_bound)
else :
genome_subdir += '_rrbs_%s_%d_%d' % ( re.sub(",","-",re.sub("-", "", options.cut_format)), options.rrbs_low_bound, options.rrbs_up_bound)
else:
possible_refs = filter(lambda dir: dir.startswith(genome+'_rrbs_'), os.listdir(options.dbpath))
if len(possible_refs) == 1:
genome_subdir = possible_refs[0]
else:
error('Cannot localize unambiguously the reference genome for RRBS. '
'Please, specify the options \"--low\" and \"--up\" that you used at the index-building step.\n'
'Possible choices are:\n' + '\n'.join([pr.split('_rrbs_')[-1].replace('_',', ') for pr in possible_refs]))
'''
db_path = os.path.expanduser(os.path.join(options.dbpath, genome_subdir + '_' + 'snap'))
if not os.path.isdir(db_path):
error('Index DIR \"' + genome_subdir + '..\" cannot be found in ' + options.dbpath +'.\n\tPlease run the bs_seeker2-build.py '
'to create it with the correct parameters for -g, -r, --low, --up and --aligner.')
# default aligner options
aligner_options_defaults = {
SNAP : { }
}
aligner_options = dict(aligner_options_defaults['snap'], **aligner_options)
aligner_options_string = lambda : ' %s ' % (' '.join(opt_key +
(' ' + ' '.join(map(str,opt_val)) # join all values if the value is an array
if type(opt_val) is list else
('' if type(opt_val) is bool and opt_val # output an empty string if it is a key-only option
else ' ' +str(opt_val)) # output the value if it is a single value
)
for opt_key, opt_val in aligner_options.iteritems() if opt_val not in [None, False]))
options.output_format = options.output_format.lower()
if options.output_format not in output.formats:
error('Output format should be one of: ' + ', '.join(output.formats))
if options.outfilename:
outfilename = options.outfilename
logfilename = outfilename
elif options.infilename is not None:
logfilename = options.infilename+'_'+ ('rr' if False else '') + 'bsse'
outfilename = logfilename + '.' + options.output_format
else:
logfilename = options.infilename_1+'_'+ ('rr' if False else '') + 'bspe'
outfilename = logfilename + '.' + options.output_format
outfilename = os.path.expanduser(outfilename)
logfilename = os.path.expanduser(logfilename)
outfile = output.outfile(outfilename, options.output_format, deserialize(os.path.join(db_path, 'refname')), ' '.join(sys.argv), options.no_SAM_header)
open_log(logfilename+'.bs3_log')
aligner_title = 'snap'
tmp_path = tempfile.mkdtemp(prefix='bs3_%s_-%s-TMP-' % (os.path.split(outfilename)[1], aligner_title ), dir = options.temp_dir)
(XS_x, XS_y) = options.XS_filter.split(",")
XS_pct = float(XS_x)
XS_count = int(XS_y)
if options.infilename is not None:
#if options.aligner == 'snap':
aligner_command = './snap single %(reference_genome)s -fastq %(input_file)s -o -sam %(output_file)s -h 630 -d 20 -t 32 -b' + aligner_options_string()
if 1 == 1: # Normal single end scan
input_fname = os.path.split(options.infilename)[1]
tmp_d = lambda fname: os.path.join(tmp_path, fname)
time1 = time.time()
p = []
j = 0
all_raw_reads = 0
for raw_reads, read_file in isplit_file(options.infilename, tmp_d(input_fname)+'-s-', options.no_split):
all_raw_reads += raw_reads
cmd = ['python bs_single_end3.py ',
'snap',
'temp-' + options.outfilename + '-' + str(j + 1),
read_file,
asktag,
options.adapter_file ,
str(options.cutnumber1) ,
str(options.cutnumber2) ,
str(options.no_split) ,
str(str_no_mismatches),
aligner_command ,
db_path ,
tmp_path ,
str(XS_pct) ,
str(XS_count) ,
str(options.adapter_mismatch) ,
str(options.Output_multiple_hit) ,
str(options.Output_unmapped_hit), str(j), options.qc_f, str(options.K)]
j += 1
with open('command', 'w') as fh:
for line in cmd:
fh.write(line + '\n')
p.append(subprocess.Popen('python bs_align/bs_single_end3.py ', shell=True, close_fds=True))
part = 1
while p:
p.pop().communicate("process " + str(part)+ " running\n")
part += 1
stats = [0 for i in range(13)]
path = 'temp-' + options.outfilename + '-*_stat-*'
for filename in glob.glob(path):
if filename == path :
continue
for i, item in enumerate(open(filename).readline().split()):
stats[i] += float(item)
subprocess.call('rm -r ' + filename, shell = True)
if (qc_len > 5):
qc_bin = [0 for qc_num in range(qc_len)]
path = 'qc-temp-' + options.outfilename + '-*_stat-*'
for filename in glob.glob(path):
if filename == path :
continue
fh = open(filename)
fh.readline()
for i, item in enumerate(fh.readline().split()):
qc_bin[i] += float(item)
subprocess.call('rm -r ' + filename, shell = True)
qc_print = open(options.outfilename + '.qc', 'w')
qc_print.write(str(stats[7]) + '\n')
qc_print.write('\n'.join([str(qc_entry) for qc_entry in qc_bin]))
f = open(options.outfilename + '.stat', 'w')
print "\n\n\tBS-seeker3 Result\n"
f = open(options.outfilename + '.stat', 'w')
f.write("\tFinal Alignment Report\n================================================\n")
print "\tFinal Alignment Report\n================================================\n"
f.write("Number of reads in total: " + str(all_raw_reads) + '\n')
print "Number of reads in total: " + str(all_raw_reads)
f.write("Number of unique-hits reads (before post-filtering): " + str(stats[8]) + '\n')
print "Number of unique-hits reads (before post-filtering): " + str(stats[8])
f.write("Number of reads mapped after post-filtering " + str(stats[7]) + '\n')
print "Number of reads mapped after post-filtering " + str(stats[7])
f.write('Alignment Time: ' + str( time.time() - time1) + 'secs\n')
print '\nAlignment Time: ' + str( time.time() - time1) + 'secs'
print '\n\tFinal Cytosine Report\n================================================\n'
f.write('\n\tFinal Cytosine Report\n================================================\n')
print 'Total Number of Cytosines: ' + str(stats[1] + stats[4] + stats[2] + stats[5] + stats[3] + stats[6] )
f.write('Total Number of Cytosines: ' + str(stats[1] + stats[4] + stats[2] + stats[5] + stats[3] + stats[6] ) + '\n')
print 'Total Number of Cs in CpG context: ' + str(stats[1] + stats[4])
f.write('Total Number of Cs in CpG context: ' + str(stats[1] + stats[4]) + '\n')
print 'Total Number of Cs in CHG context: ' + str(stats[2] + stats[5])
f.write('Total Number of Cs in CHG context: ' + str(stats[2] + stats[5]) + '\n')
print 'Total Number of Cs in CHH context: ' + str(stats[3] + stats[6])
f.write('Total Number of Cs in CHH context: ' + str(stats[3] + stats[6]) + '\n')
n_CG = stats[1] / (stats[1] + stats[4]) * 100 if stats[1] + stats[4] != 0 else 0
n_CHG = stats[2] / (stats[2] + stats[5]) * 100 if stats[2] + stats[5] != 0 else 0
n_CHH = stats[3] / (stats[3] + stats[6]) * 100 if stats[3] + stats[6] != 0 else 0
f.write("Rate of Methylation \n")
print "\nRate of Methylation"
f.write(" mCG %1.3f%%\n" % n_CG)
print " mCG %1.3f%%" % n_CG
f.write(" mCHG %1.3f%%\n" % n_CHG)
print " mCHG %1.3f%%" % n_CHG
f.write(" mCHH %1.3f%%\n" % n_CHH)
print " mCHH %1.3f%%" % n_CHH
f.close()
print '\n\n'
chrom_len = marshal.load(open(os.path.join(db_path, 'refname.data')))
f = open('temp_header', 'w')
f.write('@HD\tVN:1.0\tSO:unsorted'+ '\n')
cmd_line = ' '.join(sys.argv)
chrom_conv = pickle.load(open(os.path.join(db_path, 'chrom_conv.p')))
for c in sorted(chrom_len):
f.write('@SQ\tSN:' + chrom_conv[c] +'\tLN:' + str(chrom_len[c])+'\n')
f.write('@PG\tID:1\tPN:BS Seeker 3\tCL:'+ '\"' + cmd_line+'\"' + '\n')
f.close()
path = 'temp-' + options.outfilename + '-*_sam-*'
sam_f = [filename for filename in glob.glob(path)]
if path in sam_f:
sam_f.remove(path)
subprocess.call('cat ' + ' temp_header ' + ' '.join(sam_f) + ' > ' + options.outfilename, shell = True)
subprocess.call('rm -r ' + tmp_path, shell = True)
subprocess.call('rm -r ' + path, shell = True)
subprocess.call('rm -r temp_header', shell = True)
else:
logm('Pair end')
# pair end specific default options
#if options.aligner == 'snap':
aligner_command = './snap paired %(reference_genome)s -t 32 -b -fastq %(input_file_1)s %(input_file_2)s -o -sam %(output_file)s ' + aligner_options_string() + ' -s ' + str(options.min_insert_size) + ' ' + str(options.max_insert_size)
input_fname = os.path.split(options.infilename_1)[1]
tmp_d = lambda fname: os.path.join(tmp_path, fname)
time1 = time.time()
p = []
j = 0
all_raw_reads = 0
read1 = isplit_file(options.infilename_1, tmp_d(input_fname + '_1')+'-s-', options.no_split) # tmp_d(input_fname + '_1')+'-s-',
read2 = isplit_file(options.infilename_2, tmp_d(input_fname + '_2')+'-s-', options.no_split)
for read_file1, read_file2 in izip(read1, read2):
#raw_reads, read_file1 = read1[0]
#_, read_file2 = read2[0]
all_raw_reads += read_file1[0]
read_file1 = read_file1[1]
read_file2 = read_file2[1]
#print read_file1[0], read_file1[1], read_file2[0], read_file2[1], 'fsa'
cmd = ['python bs_pair_end3.py ',
'snap',
'temp-' + options.outfilename + '-' + str(j + 1),
read_file1,
read_file2,
asktag,
str(options.cutnumber1) ,
str(options.cutnumber2) ,
str(options.no_split) ,
str(str_no_mismatches),
aligner_command ,
db_path ,
tmp_path ,
str(XS_pct) ,
str(XS_count) ,
str(options.Output_multiple_hit) ,
str(options.Output_unmapped_hit),
str(j),
str(options.qc_f),
str(options.K)
]
j += 1
with open('command', 'w') as fh:
for line in cmd:
fh.write(line + '\n')
p.append(subprocess.Popen('python bs_align/bs_pair_end3.py ', shell=True, close_fds=True))
part = 1
while p:
print p.pop().communicate("process " + str(part)+ " running\n")
part += 1
chrom_conv = pickle.load(open(os.path.join(db_path, 'chrom_conv.p')))
stats = [0 for i in range(13)]
path = 'temp-' + options.outfilename + '-*_stat-*'
for filename in glob.glob(path):
if filename == path :
continue
for i, item in enumerate(open(filename).readline().split()):
stats[i] += float(item)
subprocess.call('rm -r ' + filename, shell = True)
if (qc_len > 5):
qc_bin = [0 for qc_num in range(qc_len)]
path = 'qc-temp-' + options.outfilename + '-*_stat-*'
for filename in glob.glob(path):
if filename == path :
continue
fh = open(filename)
fh.readline()
for i, item in enumerate(fh.readline().split()):
qc_bin[i] += float(item)
subprocess.call('rm -r ' + filename, shell = True)
qc_print = open(options.outfilename + '.qc', 'w')
qc_print.write(str(stats[7]) + '\n')
qc_print.write('\n'.join([str(qc_entry) for qc_entry in qc_bin]))
print "\n\n\tBS-seeker3 Result\n"
f = open(options.outfilename + '.stat', 'w')
f.write("\tFinal Alignment Report\n================================================\n")
print "\tFinal Alignment Report\n================================================\n"
f.write("Number of reads in total: " + str(all_raw_reads * 2) + '\n')
print "Number of reads in total: " + str(all_raw_reads * 2)
f.write("Number of unique-hits reads (before post-filtering): " + str(stats[8]) + '\n')
print "Number of unique-hits reads (before post-filtering): " + str(stats[8])
f.write("Number of reads mapped after post-filtering " + str(stats[7]) + '\n')
print "Number of reads mapped after post-filtering " + str(stats[7])
f.write('Alignment Time: ' + str( time.time() - time1) + 'secs\n')
print '\nAlignment Time: ' + str( time.time() - time1) + 'secs'
print '\n\tFinal Cytosine Report\n================================================\n'
f.write('\n\tFinal Cytosine Report\n================================================\n')
print 'Total Number of Cytosines: ' + str(stats[1] + stats[4] + stats[2] + stats[5] + stats[3] + stats[6] )
f.write('Total Number of Cytosines: ' + str(stats[1] + stats[4] + stats[2] + stats[5] + stats[3] + stats[6] ) + '\n')
print 'Total Number of Cs in CpG context: ' + str(stats[1] + stats[4])
f.write('Total Number of Cs in CpG context: ' + str(stats[1] + stats[4]) + '\n')
print 'Total Number of Cs in CHG context: ' + str(stats[2] + stats[5])
f.write('Total Number of Cs in CHG context: ' + str(stats[2] + stats[5]) + '\n')
print 'Total Number of Cs in CHH context: ' + str(stats[3] + stats[6])
f.write('Total Number of Cs in CHH context: ' + str(stats[3] + stats[6]) + '\n')
n_CG = stats[1] / (stats[1] + stats[4]) * 100 if stats[1] + stats[4] != 0 else 0
n_CHG = stats[2] / (stats[2] + stats[5]) * 100 if stats[2] + stats[5] != 0 else 0
n_CHH = stats[3] / (stats[3] + stats[6]) * 100 if stats[3] + stats[6] != 0 else 0
f.write("Rate of Methylation \n")
print "\nRate of Methylation"
f.write(" mCG %1.3f%%\n" % n_CG)
print " mCG %1.3f%%" % n_CG
f.write(" mCHG %1.3f%%\n" % n_CHG)
print " mCHG %1.3f%%" % n_CHG
f.write(" mCHH %1.3f%%\n" % n_CHH)
print " mCHH %1.3f%%" % n_CHH
f.close()
print '\n\n'
chrom_len = marshal.load(open(os.path.join(db_path, 'refname.data')))
f = open('temp_header', 'w')
f.write('@HD\tVN:1.0\tSO:unsorted'+ '\n')
cmd_line = ' '.join(sys.argv)
for c in sorted(chrom_len):
f.write('@SQ\tSN:' + chrom_conv[str(c)] +'\tLN:' + str(chrom_len[c])+'\n')
f.write('@PG\tID:1\tPN:BS Seeker 3\tCL:'+ '\"' + cmd_line+'\"' + '\n')
f.close()
path = 'temp-' + options.outfilename + '-*_sam-*'
sam_f = [filename for filename in glob.glob(path)]
if path in sam_f:
sam_f.remove(path)
subprocess.call('cat ' + ' temp_header ' + ' '.join(sam_f) + ' > ' + options.outfilename, shell = True)
subprocess.call('rm -r ' + tmp_path, shell = True)
subprocess.call('rm -r ' + path, shell = True)
subprocess.call('rm -r temp_header', shell = True)