Mercurial > repos > devteam > cuffmerge

comparison cuffmerge_wrapper.py @ 12:1707a530e598 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-devteam/tree/master/tool_collections/cufflinks/cuffmerge commit eb18f691975ef9539b5ebd4f118343c8ad967a1f
author devteam
date Tue, 07 Feb 2017 18:39:42 -0500
parents b6e3849293b1
children cf747d1bd79a
comparison
equal deleted inserted replaced
11:8160c8ea4eb9 12:1707a530e598
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 2
3 import optparse, os, shutil, subprocess, sys, tempfile 3 import optparse
4 import os
5 import shutil
6 import subprocess
7 import sys
8 import tempfile
4 9
5 def stop_err( msg ): 10
6 sys.stderr.write( '%s\n' % msg ) 11 def stop_err(msg):
12 sys.stderr.write('%s\n' % msg)
7 sys.exit() 13 sys.exit()
8 14
15
9 def __main__(): 16 def __main__():
10 #Parse Command Line
11 parser = optparse.OptionParser() 17 parser = optparse.OptionParser()
12 parser.add_option( '-g', dest='ref_annotation', help='An optional "reference" annotation GTF. Each sample is matched against this file, and sample isoforms are tagged as overlapping, matching, or novel where appropriate. See the refmap and tmap output file descriptions below.' ) 18 parser.add_option('-g', dest='ref_annotation', help='An optional "reference" annotation GTF. Each sample is matched against this file, and sample isoforms are tagged as overlapping, matching, or novel where appropriate. See the refmap and tmap output file descriptions below.')
13 parser.add_option( '-s', dest='use_seq_data', action="store_true", help='Causes cuffmerge to look into for fasta files with the underlying genomic sequences (one file per contig) against which your reads were aligned for some optional classification functions. For example, Cufflinks transcripts consisting mostly of lower-case bases are classified as repeats. Note that <seq_dir> must contain one fasta file per reference chromosome, and each file must be named after the chromosome, and have a .fa or .fasta extension.') 19 parser.add_option('-s', dest='use_seq_data', action="store_true", help='Causes cuffmerge to look into for fasta files with the underlying genomic sequences (one file per contig) against which your reads were aligned for some optional classification functions. For example, Cufflinks transcripts consisting mostly of lower-case bases are classified as repeats. Note that <seq_dir> must contain one fasta file per reference chromosome, and each file must be named after the chromosome, and have a .fa or .fasta extension.')
14 parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' ) 20 parser.add_option('-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.')
15 21
16 # Wrapper / Galaxy options. 22 # Wrapper / Galaxy options.
17 parser.add_option( '', '--index', dest='index', help='The path of the reference genome' ) 23 parser.add_option('', '--index', dest='index', help='The path of the reference genome')
18 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) 24 parser.add_option('', '--ref_file', dest='ref_file', help='The reference dataset from the history')
19 25
20 # Outputs. 26 # Outputs.
21 parser.add_option( '', '--merged-transcripts', dest='merged_transcripts' ) 27 parser.add_option('', '--merged-transcripts', dest='merged_transcripts')
22 parser.add_option( '--min-isoform-fraction', dest='min_isoform_fraction' ) 28 parser.add_option('--min-isoform-fraction', dest='min_isoform_fraction')
23 29
24 (options, args) = parser.parse_args() 30 (options, args) = parser.parse_args()
25 31
26 # output version # of tool 32 # output version # of tool
27 try: 33 try:
28 tmp = tempfile.NamedTemporaryFile().name 34 with tempfile.NamedTemporaryFile() as tmp_stdout:
29 tmp_stdout = open( tmp, 'wb' ) 35 returncode = subprocess.call(args='cuffmerge -v 2>&1', stdout=tmp_stdout, shell=True)
30 proc = subprocess.Popen( args='cuffmerge -v 2>&1', shell=True, stdout=tmp_stdout ) 36 stdout = None
31 tmp_stdout.close() 37 with open(tmp_stdout.name) as tmp_stdout2:
32 returncode = proc.wait() 38 for line in tmp_stdout2:
33 stdout = None 39 if line.lower().find('merge_cuff_asms v') >= 0:
34 for line in open( tmp_stdout.name, 'rb' ): 40 stdout = line.strip()
35 if line.lower().find( 'merge_cuff_asms v' ) >= 0: 41 break
36 stdout = line.strip()
37 break
38 if stdout: 42 if stdout:
39 sys.stdout.write( '%s\n' % stdout ) 43 sys.stdout.write('%s\n' % stdout)
40 else: 44 else:
41 raise Exception 45 raise Exception
42 except: 46 except:
43 sys.stdout.write( 'Could not determine Cuffmerge version\n' ) 47 sys.stdout.write('Could not determine Cuffmerge version\n')
44 48
45 # Set/link to sequence file. 49 # Set/link to sequence file.
46 if options.use_seq_data: 50 if options.use_seq_data:
47 if options.ref_file: 51 if options.ref_file:
48 # Sequence data from history. 52 # Sequence data from history.
49 # Create symbolic link to ref_file so that index will be created in working directory. 53 # Create symbolic link to ref_file so that index will be created in working directory.
50 seq_path = "ref.fa" 54 seq_path = "ref.fa"
51 os.symlink( options.ref_file, seq_path ) 55 os.symlink(options.ref_file, seq_path)
52 else: 56 else:
53 if not os.path.exists( options.index ): 57 if not os.path.exists(options.index):
54 stop_err( 'Reference genome %s not present, request it by reporting this error.' % options.index ) 58 stop_err('Reference genome %s not present, request it by reporting this error.' % options.index)
55 seq_path = options.index 59 seq_path = options.index
56 60
57 # Build command. 61 # Build command.
58 62
59 # Base. 63 # Base.
60 cmd = "cuffmerge -o cm_output " 64 cmd = "cuffmerge -o cm_output "
61 65
62 # Add options. 66 # Add options.
63 if options.num_threads: 67 if options.num_threads:
64 cmd += ( " -p %i " % int ( options.num_threads ) ) 68 cmd += (" -p %i " % int(options.num_threads))
65 if options.ref_annotation: 69 if options.ref_annotation:
66 cmd += " -g %s " % options.ref_annotation 70 cmd += " -g %s " % options.ref_annotation
67 if options.use_seq_data: 71 if options.use_seq_data:
68 cmd += " -s %s " % seq_path 72 cmd += " -s %s " % seq_path
69 if options.min_isoform_fraction: 73 if options.min_isoform_fraction:
70 cmd += " --min-isoform-fraction %s " % (options.min_isoform_fraction) 74 cmd += " --min-isoform-fraction %s " % (options.min_isoform_fraction)
71 # Add input files to a file. 75 # Add input files to a file.
72 inputs_file_name = tempfile.NamedTemporaryFile( dir="." ).name 76 with tempfile.NamedTemporaryFile(mode='w', dir=".", delete=False) as inputs_file:
73 inputs_file = open( inputs_file_name, 'w' ) 77 for arg in args:
74 for arg in args: 78 inputs_file.write(arg + "\n")
75 inputs_file.write( arg + "\n" ) 79 cmd += inputs_file.name
76 inputs_file.close()
77 cmd += inputs_file_name
78 80
79 # Debugging.
80 print cmd
81
82 # Run command. 81 # Run command.
83 try: 82 try:
84 tmp_name = tempfile.NamedTemporaryFile( dir="." ).name 83 with tempfile.NamedTemporaryFile(dir=".") as tmp_stderr:
85 tmp_stderr = open( tmp_name, 'wb' ) 84 returncode = subprocess.call(args=cmd, stderr=tmp_stderr, shell=True)
86 proc = subprocess.Popen( args=cmd, shell=True, stderr=tmp_stderr.fileno() ) 85
87 returncode = proc.wait() 86 # Error checking.
88 tmp_stderr.close() 87 if returncode != 0:
89 88 # Get stderr, allowing for case where it's very large.
90 # Get stderr, allowing for case where it's very large. 89 buffsize = 1048576
91 tmp_stderr = open( tmp_name, 'rb' ) 90 stderr = ''
92 stderr = '' 91 with open(tmp_stderr.name, 'r') as tmp_stderr2:
93 buffsize = 1048576 92 try:
94 try: 93 while True:
95 while True: 94 stderr += tmp_stderr2.read(buffsize)
96 stderr += tmp_stderr.read( buffsize ) 95 if not stderr or len(stderr) % buffsize != 0:
97 if not stderr or len( stderr ) % buffsize != 0: 96 break
98 break 97 except OverflowError:
99 except OverflowError: 98 pass
100 pass 99 raise Exception(stderr)
101 tmp_stderr.close() 100
102 101 if len(open("cm_output/merged.gtf", 'r').read().strip()) == 0:
103 # Error checking. 102 raise Exception('The output file is empty, there may be an error with your input file or settings.')
104 if returncode != 0: 103
105 raise Exception, stderr
106
107 if len( open( "cm_output/merged.gtf", 'rb' ).read().strip() ) == 0:
108 raise Exception, 'The output file is empty, there may be an error with your input file or settings.'
109
110 # Copy outputs. 104 # Copy outputs.
111 shutil.copyfile( "cm_output/merged.gtf" , options.merged_transcripts ) 105 shutil.move("cm_output/merged.gtf", options.merged_transcripts)
112 106 except Exception as e:
113 except Exception, e: 107 stop_err('Error running cuffmerge: %s' % e)
114 stop_err( 'Error running cuffmerge. ' + str( e ) ) 108
115 109
116 if __name__=="__main__": __main__() 110 if __name__ == "__main__":
111 __main__()