root/galaxy-central/tools/samtools/sam_to_bam.py

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

import galaxy-central

行番号 
1#!/usr/bin/env python
2"""
3Converts SAM data to sorted BAM data.
4usage: sam_to_bam.py [options]
5   --input1: SAM file to be converted
6   --dbkey: dbkey value
7   --ref_file: Reference file if choosing from history
8   --output1: output dataset in bam format
9   --index_dir: GALAXY_DATA_INDEX_DIR
10"""
11
12import optparse, os, sys, subprocess, tempfile, shutil, gzip
13from galaxy import eggs
14import pkg_resources; pkg_resources.require( "bx-python" )
15from bx.cookbook import doc_optparse
16from galaxy import util
17
18def stop_err( msg ):
19    sys.stderr.write( '%s\n' % msg )
20    sys.exit()
21
22def check_seq_file( dbkey, cached_seqs_pointer_file ):
23    seq_path = ''
24    for line in open( cached_seqs_pointer_file ):
25        line = line.rstrip( '\r\n' )
26        if line and not line.startswith( '#' ) and line.startswith( 'index' ):
27            fields = line.split( '\t' )
28            if len( fields ) < 3:
29                continue
30            if fields[1] == dbkey:
31                seq_path = fields[2].strip()
32                break
33    return seq_path
34
35def __main__():
36    #Parse Command Line
37    parser = optparse.OptionParser()
38    parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' )
39    parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' )
40    parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' )
41    parser.add_option( '', '--output1', dest='output1', help='The output BAM dataset' )
42    parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' )
43    ( options, args ) = parser.parse_args()
44
45    cached_seqs_pointer_file = '%s/sam_fa_indices.loc' % options.index_dir
46    if not os.path.exists( cached_seqs_pointer_file ):
47        stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file )
48    # If found for the dbkey, seq_path will look something like /galaxy/data/equCab2/sam_index/equCab2.fa,
49    # and the equCab2.fa file will contain fasta sequences.
50    seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file )
51    tmp_dir = tempfile.mkdtemp()
52    if options.ref_file == 'None':
53        # We're using locally cached reference sequences( e.g., /galaxy/data/equCab2/sam_index/equCab2.fa ).
54        # The indexes for /galaxy/data/equCab2/sam_index/equCab2.fa will be contained in
55        # a file named /galaxy/data/equCab2/sam_index/equCab2.fa.fai
56        fai_index_file_base = seq_path
57        fai_index_file_path = '%s.fai' % seq_path
58        if not os.path.exists( fai_index_file_path ):
59            #clean up temp files
60            if os.path.exists( tmp_dir ):
61                shutil.rmtree( tmp_dir )
62            stop_err( 'No sequences are available for build (%s), request them by reporting this error.' % options.dbkey )
63    else:
64        try:
65            # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will:
66            # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence
67            # - if no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk
68            # - if regions are specified, the subsequences will be retrieved and printed to stdout in the FASTA format
69            # - the input file can be compressed in the RAZF format.
70            # IMPORTANT NOTE: a real weakness here is that we are creating indexes for the history dataset
71            # every time we run this tool.  It would be nice if we could somehow keep track of user's specific
72            # index files so they could be re-used.
73            fai_index_file_base = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
74            # At this point, fai_index_file_path will look something like /tmp/dataset_13.dat
75            os.symlink( options.ref_file, fai_index_file_base )
76            fai_index_file_path = '%s.fai' % fai_index_file_base
77            command = 'samtools faidx %s' % fai_index_file_base
78            tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
79            tmp_stderr = open( tmp, 'wb' )
80            proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
81            returncode = proc.wait()
82            tmp_stderr.close()
83            # get stderr, allowing for case where it's very large
84            tmp_stderr = open( tmp, 'rb' )
85            stderr = ''
86            buffsize = 1048576
87            try:
88                while True:
89                    stderr += tmp_stderr.read( buffsize )
90                    if not stderr or len( stderr ) % buffsize != 0:
91                        break
92            except OverflowError:
93                pass
94            tmp_stderr.close()
95            if returncode != 0:
96                raise Exception, stderr
97            if os.path.getsize( fai_index_file_path ) == 0:
98                raise Exception, 'Index file empty, there may be an error with your reference file or settings.'
99        except Exception, e:
100            #clean up temp files
101            if os.path.exists( tmp_dir ):
102                shutil.rmtree( tmp_dir )
103            stop_err( 'Error creating indexes from reference (%s), %s' % ( options.ref_file, str( e ) ) )
104    try:
105        # Extract all alignments from the input SAM file to BAM format ( since no region is specified, all the alignments will be extracted ).
106        tmp_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
107        tmp_aligns_file_name = tmp_aligns_file.name
108        tmp_aligns_file.close()
109        # IMPORTANT NOTE: for some reason the samtools view command gzips the resulting bam file without warning,
110        # and the docs do not currently state that this occurs ( very bad ).
111        command = 'samtools view -bt %s -o %s %s' % ( fai_index_file_path, tmp_aligns_file_name, options.input1 )
112        tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
113        tmp_stderr = open( tmp, 'wb' )
114        proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
115        returncode = proc.wait()
116        tmp_stderr.close()
117        # get stderr, allowing for case where it's very large
118        tmp_stderr = open( tmp, 'rb' )
119        stderr = ''
120        buffsize = 1048576
121        try:
122            while True:
123                stderr += tmp_stderr.read( buffsize )
124                if not stderr or len( stderr ) % buffsize != 0:
125                    break
126        except OverflowError:
127            pass
128        tmp_stderr.close()
129        if returncode != 0:
130            raise Exception, stderr
131        if len( open( tmp_aligns_file_name ).read() ) == 0:
132            raise Exception, 'Initial BAM file empty'
133    except Exception, e:
134        #clean up temp files
135        if os.path.exists( tmp_dir ):
136            shutil.rmtree( tmp_dir )
137        stop_err( 'Error extracting alignments from (%s), %s' % ( options.input1, str( e ) ) )
138    try:
139        # Sort alignments by leftmost coordinates. File <out.prefix>.bam will be created. This command
140        # may also create temporary files <out.prefix>.%d.bam when the whole alignment cannot be fitted
141        # into memory ( controlled by option -m ).
142        tmp_sorted_aligns_file = tempfile.NamedTemporaryFile( dir=tmp_dir )
143        tmp_sorted_aligns_file_name = tmp_sorted_aligns_file.name
144        tmp_sorted_aligns_file.close()
145        command = 'samtools sort %s %s' % ( tmp_aligns_file_name, tmp_sorted_aligns_file_name )
146        tmp = tempfile.NamedTemporaryFile( dir=tmp_dir ).name
147        tmp_stderr = open( tmp, 'wb' )
148        proc = subprocess.Popen( args=command, shell=True, cwd=tmp_dir, stderr=tmp_stderr.fileno() )
149        returncode = proc.wait()
150        tmp_stderr.close()
151        # get stderr, allowing for case where it's very large
152        tmp_stderr = open( tmp, 'rb' )
153        stderr = ''
154        buffsize = 1048576
155        try:
156            while True:
157                stderr += tmp_stderr.read( buffsize )
158                if not stderr or len( stderr ) % buffsize != 0:
159                    break
160        except OverflowError:
161            pass
162        tmp_stderr.close()
163        if returncode != 0:
164            raise Exception, stderr
165    except Exception, e:
166        #clean up temp files
167        if os.path.exists( tmp_dir ):
168            shutil.rmtree( tmp_dir )
169        stop_err( 'Error sorting alignments from (%s), %s' % ( tmp_aligns_file_name, str( e ) ) )
170    # Move tmp_aligns_file_name to our output dataset location
171    sorted_bam_file = '%s.bam' % tmp_sorted_aligns_file_name
172    if os.path.getsize( sorted_bam_file ) == 0:
173        #clean up temp files
174        if os.path.exists( tmp_dir ):
175            shutil.rmtree( tmp_dir )
176        stop_err( 'Error creating sorted version of BAM file' )
177    shutil.move( sorted_bam_file, options.output1 )
178    #clean up temp files
179    if os.path.exists( tmp_dir ):
180        shutil.rmtree( tmp_dir )
181    # check that there are results in the output file
182    if os.path.getsize( options.output1 ) > 0:
183        sys.stdout.write( 'SAM file converted to BAM' )
184    else:
185        stop_err( 'The output file is empty, there may be an error with your input file or settings.' )
186
187if __name__=="__main__": __main__()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。