Mercurial > repos > devteam > cuffmerge
annotate cuffmerge_wrapper.py @ 2:5b285b6e4ee3
Update to the new data table specification.
| author | Dave Bouvier <dave@bx.psu.edu> |
|---|---|
| date | Wed, 04 Dec 2013 13:24:59 -0500 |
| parents | dbbd37e013aa |
| children | b6e3849293b1 |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 | |
| 3 # Supports Cuffmerge versions 1.3 and newer. | |
| 4 | |
| 5 import optparse, os, shutil, subprocess, sys, tempfile | |
| 6 | |
| 7 def stop_err( msg ): | |
| 8 sys.stderr.write( '%s\n' % msg ) | |
| 9 sys.exit() | |
| 10 | |
| 11 def __main__(): | |
| 12 #Parse Command Line | |
| 13 parser = optparse.OptionParser() | |
| 14 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.' ) | |
| 15 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.') | |
| 16 parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' ) | |
| 17 | |
| 18 | |
| 19 # Wrapper / Galaxy options. | |
|
2
5b285b6e4ee3
Update to the new data table specification.
Dave Bouvier <dave@bx.psu.edu>
parents:
0
diff
changeset
|
20 parser.add_option( '', '--index', dest='index', help='The path of the reference genome' ) |
| 0 | 21 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) |
| 22 | |
| 23 # Outputs. | |
| 24 parser.add_option( '', '--merged-transcripts', dest='merged_transcripts' ) | |
| 25 | |
| 26 (options, args) = parser.parse_args() | |
| 27 | |
| 28 # output version # of tool | |
| 29 try: | |
| 30 tmp = tempfile.NamedTemporaryFile().name | |
| 31 tmp_stdout = open( tmp, 'wb' ) | |
| 32 proc = subprocess.Popen( args='cuffmerge -v 2>&1', shell=True, stdout=tmp_stdout ) | |
| 33 tmp_stdout.close() | |
| 34 returncode = proc.wait() | |
| 35 stdout = None | |
| 36 for line in open( tmp_stdout.name, 'rb' ): | |
| 37 if line.lower().find( 'merge_cuff_asms v' ) >= 0: | |
| 38 stdout = line.strip() | |
| 39 break | |
| 40 if stdout: | |
| 41 sys.stdout.write( '%s\n' % stdout ) | |
| 42 else: | |
| 43 raise Exception | |
| 44 except: | |
| 45 sys.stdout.write( 'Could not determine Cuffmerge version\n' ) | |
| 46 | |
| 47 # Set/link to sequence file. | |
| 48 if options.use_seq_data: | |
|
2
5b285b6e4ee3
Update to the new data table specification.
Dave Bouvier <dave@bx.psu.edu>
parents:
0
diff
changeset
|
49 if options.ref_file: |
| 0 | 50 # Sequence data from history. |
| 51 # Create symbolic link to ref_file so that index will be created in working directory. | |
| 52 seq_path = "ref.fa" | |
| 53 os.symlink( options.ref_file, seq_path ) | |
| 54 else: | |
|
2
5b285b6e4ee3
Update to the new data table specification.
Dave Bouvier <dave@bx.psu.edu>
parents:
0
diff
changeset
|
55 if not os.path.exists( options.index ): |
|
5b285b6e4ee3
Update to the new data table specification.
Dave Bouvier <dave@bx.psu.edu>
parents:
0
diff
changeset
|
56 stop_err( 'Reference genome %s not present, request it by reporting this error.' % options.index ) |
|
5b285b6e4ee3
Update to the new data table specification.
Dave Bouvier <dave@bx.psu.edu>
parents:
0
diff
changeset
|
57 seq_path = options.index |
| 0 | 58 |
| 59 # Build command. | |
| 60 | |
| 61 # Base. | |
| 62 cmd = "cuffmerge -o cm_output " | |
| 63 | |
| 64 # Add options. | |
| 65 if options.num_threads: | |
| 66 cmd += ( " -p %i " % int ( options.num_threads ) ) | |
| 67 if options.ref_annotation: | |
| 68 cmd += " -g %s " % options.ref_annotation | |
| 69 if options.use_seq_data: | |
| 70 cmd += " -s %s " % seq_path | |
| 71 | |
| 72 # Add input files to a file. | |
| 73 inputs_file_name = tempfile.NamedTemporaryFile( dir="." ).name | |
| 74 inputs_file = open( inputs_file_name, 'w' ) | |
| 75 for arg in args: | |
| 76 inputs_file.write( arg + "\n" ) | |
| 77 inputs_file.close() | |
| 78 cmd += inputs_file_name | |
| 79 | |
| 80 # Debugging. | |
| 81 print cmd | |
| 82 | |
| 83 # Run command. | |
| 84 try: | |
| 85 tmp_name = tempfile.NamedTemporaryFile( dir="." ).name | |
| 86 tmp_stderr = open( tmp_name, 'wb' ) | |
| 87 proc = subprocess.Popen( args=cmd, shell=True, stderr=tmp_stderr.fileno() ) | |
| 88 returncode = proc.wait() | |
| 89 tmp_stderr.close() | |
| 90 | |
| 91 # Get stderr, allowing for case where it's very large. | |
| 92 tmp_stderr = open( tmp_name, 'rb' ) | |
| 93 stderr = '' | |
| 94 buffsize = 1048576 | |
| 95 try: | |
| 96 while True: | |
| 97 stderr += tmp_stderr.read( buffsize ) | |
| 98 if not stderr or len( stderr ) % buffsize != 0: | |
| 99 break | |
| 100 except OverflowError: | |
| 101 pass | |
| 102 tmp_stderr.close() | |
| 103 | |
| 104 # Error checking. | |
| 105 if returncode != 0: | |
| 106 raise Exception, stderr | |
| 107 | |
| 108 if len( open( "cm_output/merged.gtf", 'rb' ).read().strip() ) == 0: | |
| 109 raise Exception, 'The output file is empty, there may be an error with your input file or settings.' | |
| 110 | |
| 111 # Copy outputs. | |
| 112 shutil.copyfile( "cm_output/merged.gtf" , options.merged_transcripts ) | |
| 113 | |
| 114 except Exception, e: | |
| 115 stop_err( 'Error running cuffmerge. ' + str( e ) ) | |
| 116 | |
| 117 if __name__=="__main__": __main__() |