Mercurial > repos > devteam > sam_to_bam

comparison sam_to_bam.py @ 1:93f2e3337a33

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update sam_to_bam to use the fasta_indexes data table.
author Dave Bouvier <dave@bx.psu.edu>
date Wed, 11 Dec 2013 12:54:32 -0500
parents 30fdbaccb96b
children ab4c4e07eb3c
comparison
equal deleted inserted replaced
0:30fdbaccb96b 1:93f2e3337a33
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 """ 2 """
3 Converts SAM data to sorted BAM data. 3 Converts SAM data to sorted BAM data.
4 usage: sam_to_bam.py [options] 4 usage: sam_to_bam.py [options]
5 --input1: SAM file to be converted 5 --input1: SAM file to be converted
6 --dbkey: dbkey value 6 --index: path of the indexed reference genome
7 --ref_file: Reference file if choosing from history 7 --ref_file: Reference file if choosing from history
8 --output1: output dataset in bam format 8 --output1: output dataset in bam format
9 --index_dir: GALAXY_DATA_INDEX_DIR
10 """ 9 """
11 10
12 import optparse, os, sys, subprocess, tempfile, shutil, gzip 11 import optparse, os, sys, subprocess, tempfile, shutil
13 from galaxy import eggs
14 import pkg_resources; pkg_resources.require( "bx-python" )
15 from bx.cookbook import doc_optparse
16 from galaxy import util
17 12
18 def stop_err( msg ): 13 def stop_err( msg ):
19 sys.stderr.write( '%s\n' % msg ) 14 sys.stderr.write( '%s\n' % msg )
20 sys.exit() 15 sys.exit()
21 16
22 def 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
35 def __main__(): 17 def __main__():
36 #Parse Command Line 18 #Parse Command Line
37 parser = optparse.OptionParser() 19 parser = optparse.OptionParser()
38 parser.add_option( '', '--input1', dest='input1', help='The input SAM dataset' ) 20 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' ) 21
22 parser.add_option( '', '--index', dest='index', help='The path of the indexed reference genome' )
40 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) 23 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' ) 24 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() 25 ( options, args ) = parser.parse_args()
44 26
45 # output version # of tool 27 # output version # of tool
46 try: 28 try:
47 tmp = tempfile.NamedTemporaryFile().name 29 tmp = tempfile.NamedTemporaryFile().name
59 else: 41 else:
60 raise Exception 42 raise Exception
61 except: 43 except:
62 sys.stdout.write( 'Could not determine Samtools version\n' ) 44 sys.stdout.write( 'Could not determine Samtools version\n' )
63 45
64 cached_seqs_pointer_file = '%s/sam_fa_indices.loc' % options.index_dir
65 if not os.path.exists( cached_seqs_pointer_file ):
66 stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file )
67 # If found for the dbkey, seq_path will look something like /galaxy/data/equCab2/sam_index/equCab2.fa,
68 # and the equCab2.fa file will contain fasta sequences.
69 seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file )
70 tmp_dir = tempfile.mkdtemp( dir='.' ) 46 tmp_dir = tempfile.mkdtemp( dir='.' )
71 if not options.ref_file or options.ref_file == 'None': 47 if not options.ref_file or options.ref_file == 'None':
72 # We're using locally cached reference sequences( e.g., /galaxy/data/equCab2/sam_index/equCab2.fa ). 48 # We're using locally cached reference sequences( e.g., /galaxy/data/equCab2/sam_index/equCab2.fa ).
73 # The indexes for /galaxy/data/equCab2/sam_index/equCab2.fa will be contained in 49 # The indexes for /galaxy/data/equCab2/sam_index/equCab2.fa will be contained in
74 # a file named /galaxy/data/equCab2/sam_index/equCab2.fa.fai 50 # a file named /galaxy/data/equCab2/sam_index/equCab2.fa.fai
75 fai_index_file_base = seq_path 51 fai_index_file_base = seq_path
76 fai_index_file_path = '%s.fai' % seq_path 52 fai_index_file_path = '%s.fai' % options.index
77 if not os.path.exists( fai_index_file_path ): 53 if not os.path.exists( fai_index_file_path ):
78 #clean up temp files 54 #clean up temp files
79 if os.path.exists( tmp_dir ): 55 if os.path.exists( tmp_dir ):
80 shutil.rmtree( tmp_dir ) 56 shutil.rmtree( tmp_dir )
81 stop_err( 'No sequences are available for build (%s), request them by reporting this error.' % options.dbkey ) 57 stop_err( 'Indexed genome %s not present, request it by reporting this error.' % options.index )
82 else: 58 else:
83 try: 59 try:
84 # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will: 60 # Create indexes for history reference ( e.g., ~/database/files/000/dataset_1.dat ) using samtools faidx, which will:
85 # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence 61 # - index reference sequence in the FASTA format or extract subsequence from indexed reference sequence
86 # - if no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk 62 # - if no region is specified, faidx will index the file and create <ref.fasta>.fai on the disk