root/galaxy-central/tools/metag_tools/rmapq_wrapper.py

リビジョン 2, 2.7 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
1#!/usr/bin/env python
2
3import os, sys, tempfile
4
5assert sys.version_info[:2] >= (2.4)
6
7def stop_err( msg ):
8   
9    sys.stderr.write( "%s\n" % msg )
10    sys.exit()
11   
12
13def __main__():
14   
15    # I/O
16    target_path = sys.argv[1]
17    infile = sys.argv[2]
18    scorefile = sys.argv[3]
19    high_score = sys.argv[4]            # -q
20    high_len = sys.argv[5]              # -M
21    read_len = sys.argv[6]              # -w
22    align_len = sys.argv[7]             # -h
23    mismatch = sys.argv[8]              # -m
24    output_file = sys.argv[9]
25   
26    try:
27        float(high_score)
28    except:
29        stop_err('Invalid value for minimal quality score.')
30
31    try:
32        int(high_len)
33    except:
34        stop_err('Invalid value for minimal high quality bases.')
35           
36    # first guess the read length
37    guess_read_len = 0
38    seq = ''
39    for i, line in enumerate(open(infile)):
40        line = line.rstrip('\r\n')
41        if line.startswith('>'):
42            if seq:
43                guess_read_len = len(seq)
44                break
45        else:
46            seq += line
47           
48    try:
49        test = int(read_len)
50        if test == 0:
51            read_len = str(guess_read_len)
52        else:
53            assert test >= 20 and test <= 64
54    except:
55        stop_err('Invalid value for read length. Must be between 20 and 64.')
56
57   
58    try:
59        int(align_len)   
60    except:
61        stop_err('Invalid value for minimal length of a hit.')
62   
63    try:
64        int(mismatch)
65    except:
66        stop_err('Invalid value for mismatch numbers in an alignment.')
67   
68    all_files = []
69    if os.path.isdir(target_path):
70        # check target genome
71        fa_files = os.listdir(target_path)
72           
73        for file in fa_files:
74            file = "%s/%s" % ( target_path, file )
75            file = os.path.normpath(file)
76            all_files.append(file)
77    else:
78        stop_err("No sequences for %s are available for search, please report this error." %(target_path))
79   
80    for detail_file_path in all_files:
81        output_tempfile = tempfile.NamedTemporaryFile().name
82        command = "rmapq -q %s -M %s -h %s -w %s -m %s -Q %s -c %s %s -o %s 2>&1" % ( high_score, high_len, align_len, read_len, mismatch, scorefile, detail_file_path, infile, output_tempfile )
83        #print command
84        try:
85            os.system( command )
86        except Exception, e:
87            stop_err( str( e ) )
88
89        try:
90            assert os.system( 'cat %s >> %s' % ( output_tempfile, output_file ) ) == 0
91        except Exception, e:
92            stop_err( str( e ) )
93       
94        try:
95            os.remove( output_tempfile )
96        except:
97            pass
98
99           
100if __name__ == '__main__': __main__()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。