Mercurial > repos > devteam > cuffcompare

comparison cuffcompare_wrapper.py @ 8:1322b73ffe44 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/cuffcompare commit eb18f691975ef9539b5ebd4f118343c8ad967a1f
author devteam
date Tue, 07 Feb 2017 18:38:41 -0500
parents 8e534225baa9
children e66b9b5b8580
comparison
equal deleted inserted replaced
7:b77178f66fc3 8:1322b73ffe44
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 2
3 # Supports Cuffcompare versions v1.3.0 and newer. 3 # Supports Cuffcompare versions v1.3.0 and newer.
4 4
5 import optparse, os, shutil, subprocess, sys, tempfile 5 import optparse
6 import os
7 import shutil
8 import subprocess
9 import sys
10 import tempfile
6 11
7 def stop_err( msg ): 12
8 sys.stderr.write( '%s\n' % msg ) 13 def stop_err(msg):
14 sys.stderr.write('%s\n' % msg)
9 sys.exit() 15 sys.exit()
10 16
17
11 def __main__(): 18 def __main__():
12 #Parse Command Line
13 parser = optparse.OptionParser() 19 parser = optparse.OptionParser()
14 parser.add_option( '-r', 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.' ) 20 parser.add_option('-r', 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( '-R', action="store_true", dest='ignore_nonoverlap_reference', help='If -r was specified, this option causes cuffcompare to ignore reference transcripts that are not overlapped by any transcript in one of cuff1.gtf,...,cuffN.gtf. Useful for ignoring annotated transcripts that are not present in your RNA-Seq samples and thus adjusting the "sensitivity" calculation in the accuracy report written in the transcripts accuracy file' ) 21 parser.add_option('-R', action="store_true", dest='ignore_nonoverlap_reference', help='If -r was specified, this option causes cuffcompare to ignore reference transcripts that are not overlapped by any transcript in one of cuff1.gtf,...,cuffN.gtf. Useful for ignoring annotated transcripts that are not present in your RNA-Seq samples and thus adjusting the "sensitivity" calculation in the accuracy report written in the transcripts accuracy file')
16 parser.add_option( '-Q', action="store_true", dest='ignore_nonoverlap_transfrag', help='If -r was specified, this option causes cuffcompare to consider only the input transcripts that overlap any of the reference transcripts (Sp correction); Warning: this will discard all "novel" loci!)' ) 22 parser.add_option('-Q', action="store_true", dest='ignore_nonoverlap_transfrag', help='If -r was specified, this option causes cuffcompare to consider only the input transcripts that overlap any of the reference transcripts (Sp correction); Warning: this will discard all "novel" loci!)')
17 23
18 parser.add_option( '-s', dest='use_seq_data', action="store_true", help='Causes cuffcompare 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.') 24 parser.add_option('-s', dest='use_seq_data', action="store_true", help='Causes cuffcompare 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 25
20 parser.add_option( '-M', action="store_true", dest='discard_single_exon_all', help='discard (ignore) single-exon transfrags and reference transcript') 26 parser.add_option('-M', action="store_true", dest='discard_single_exon_all', help='discard (ignore) single-exon transfrags and reference transcript')
21 parser.add_option( '-N', action="store_true", dest='discard_single_exon_ref', help='discard (ignore) single-exon reference transcripts') 27 parser.add_option('-N', action="store_true", dest='discard_single_exon_ref', help='discard (ignore) single-exon reference transcripts')
22 parser.add_option( '-e', dest='max_dist_exon', help='Max. Distance for assessing exon accuracy" help="max. distance (range) allowed from free ends of terminal exons of reference transcripts when assessing exon accuracy. Default: 100') 28 parser.add_option('-e', dest='max_dist_exon', help='Max. Distance for assessing exon accuracy" help="max. distance (range) allowed from free ends of terminal exons of reference transcripts when assessing exon accuracy. Default: 100')
23 parser.add_option( '-d', dest='max_dist_group', help='Max.Distance for transcript grouping" help="max. distance (range) for grouping transcript start sites. Default: 100') 29 parser.add_option('-d', dest='max_dist_group', help='Max.Distance for transcript grouping" help="max. distance (range) for grouping transcript start sites. Default: 100')
24 parser.add_option( '-F', action="store_true", dest='discard_redundant_intron_transfrags', help='Discard intron-redundant transfrags if they share the 5-prime end (if they differ only at the 3-prime end)') 30 parser.add_option('-F', action="store_true", dest='discard_redundant_intron_transfrags', help='Discard intron-redundant transfrags if they share the 5-prime end (if they differ only at the 3-prime end)')
25 31
26 # Wrapper / Galaxy options. 32 # Wrapper / Galaxy options.
27 parser.add_option( '', '--index', dest='index', help='The path of the reference genome' ) 33 parser.add_option('', '--index', dest='index', help='The path of the reference genome')
28 parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' ) 34 parser.add_option('', '--ref_file', dest='ref_file', help='The reference dataset from the history')
29 35
30 # Outputs. 36 # Outputs.
31 parser.add_option( '', '--combined-transcripts', dest='combined_transcripts' ) 37 parser.add_option('', '--combined-transcripts', dest='combined_transcripts')
32 38
33 (options, args) = parser.parse_args() 39 (options, args) = parser.parse_args()
34 40
35 # output version # of tool 41 # output version # of tool
36 try: 42 try:
37 tmp = tempfile.NamedTemporaryFile().name 43 with tempfile.NamedTemporaryFile() as tmp_stdout:
38 tmp_stdout = open( tmp, 'wb' ) 44 returncode = subprocess.call(args='cuffcompare 2>&1', stdout=tmp_stdout, shell=True)
39 proc = subprocess.Popen( args='cuffcompare 2>&1', shell=True, stdout=tmp_stdout ) 45 stdout = None
40 tmp_stdout.close() 46 with open(tmp_stdout.name) as tmp_stdout2:
41 returncode = proc.wait() 47 for line in tmp_stdout2:
42 stdout = None 48 if line.lower().find('cuffcompare v') >= 0:
43 for line in open( tmp_stdout.name, 'rb' ): 49 stdout = line.strip()
44 if line.lower().find( 'cuffcompare v' ) >= 0: 50 break
45 stdout = line.strip()
46 break
47 if stdout: 51 if stdout:
48 sys.stdout.write( '%s\n' % stdout ) 52 sys.stdout.write('%s\n' % stdout)
49 else: 53 else:
50 raise Exception 54 raise Exception
51 except: 55 except:
52 sys.stdout.write( 'Could not determine Cuffcompare version\n' ) 56 sys.stdout.write('Could not determine Cuffcompare version\n')
53 57
54 # Set/link to sequence file. 58 # Set/link to sequence file.
55 if options.use_seq_data: 59 if options.use_seq_data:
56 if options.ref_file: 60 if options.ref_file:
57 # Sequence data from history. 61 # Sequence data from history.
58 # Create symbolic link to ref_file so that index will be created in working directory. 62 # Create symbolic link to ref_file so that index will be created in working directory.
59 seq_path = "ref.fa" 63 seq_path = "ref.fa"
60 os.symlink( options.ref_file, seq_path ) 64 os.symlink(options.ref_file, seq_path)
61 else: 65 else:
62 if not os.path.exists( options.index ): 66 if not os.path.exists(options.index):
63 stop_err( 'Reference genome %s not present, request it by reporting this error.' % options.index ) 67 stop_err('Reference genome %s not present, request it by reporting this error.' % options.index)
64 seq_path = options.index 68 seq_path = options.index
65 69
66 # Build command. 70 # Build command.
67 71
68 # Base. 72 # Base.
69 cmd = "cuffcompare -o cc_output " 73 cmd = "cuffcompare -o cc_output "
70 74
71 # Add options. 75 # Add options.
72 if options.ref_annotation: 76 if options.ref_annotation:
73 cmd += " -r %s " % options.ref_annotation 77 cmd += " -r %s " % options.ref_annotation
74 if options.ignore_nonoverlap_reference: 78 if options.ignore_nonoverlap_reference:
75 cmd += " -R " 79 cmd += " -R "
76 if options.ignore_nonoverlap_transfrag: 80 if options.ignore_nonoverlap_transfrag:
77 cmd += " -Q " 81 cmd += " -Q "
78 if options.use_seq_data: 82 if options.use_seq_data:
79 cmd += " -s %s " % seq_path 83 cmd += " -s %s " % seq_path
80 if options.discard_single_exon_all: 84 if options.discard_single_exon_all:
81 cmd += " -M " 85 cmd += " -M "
82 if options.discard_single_exon_ref: 86 if options.discard_single_exon_ref:
83 cmd += " -N " 87 cmd += " -N "
84 if options.max_dist_exon: 88 if options.max_dist_exon:
85 cmd += " -e %i " % int( options.max_dist_exon ) 89 cmd += " -e %i " % int(options.max_dist_exon)
86 if options.max_dist_group: 90 if options.max_dist_group:
87 cmd += " -d %i " % int( options.max_dist_group ) 91 cmd += " -d %i " % int(options.max_dist_group)
88 if options.discard_redundant_intron_transfrags: 92 if options.discard_redundant_intron_transfrags:
89 cmd += " -F " 93 cmd += " -F "
90 # Add input files. 94 # Add input files.
91 95
92 # Need to symlink inputs so that output files are written to temp directory. 96 # Need to symlink inputs so that output files are written to the current directory.
93 for i, arg in enumerate( args ): 97 for i, arg in enumerate(args):
94 input_file_name = "./input%i" % ( i+1 ) 98 input_file_name = "./input%i" % (i + 1)
95 os.symlink( arg, input_file_name ) 99 os.symlink(arg, input_file_name)
96 cmd += "%s " % input_file_name 100 cmd += "%s " % input_file_name
97 101
98 # Debugging.
99 print cmd
100
101 # Run command. 102 # Run command.
102 try: 103 try:
103 tmp_name = tempfile.NamedTemporaryFile( dir="." ).name 104 with tempfile.NamedTemporaryFile(dir=".") as tmp_stderr:
104 tmp_stderr = open( tmp_name, 'wb' ) 105 returncode = subprocess.call(args=cmd, stderr=tmp_stderr, shell=True)
105 proc = subprocess.Popen( args=cmd, shell=True, stderr=tmp_stderr.fileno() ) 106
106 returncode = proc.wait() 107 # Error checking.
107 tmp_stderr.close() 108 if returncode != 0:
108 109 # Get stderr, allowing for case where it's very large.
109 # Get stderr, allowing for case where it's very large. 110 buffsize = 1048576
110 tmp_stderr = open( tmp_name, 'rb' ) 111 stderr = ''
111 stderr = '' 112 with open(tmp_stderr.name) as tmp_stderr2:
112 buffsize = 1048576 113 try:
113 try: 114 while True:
114 while True: 115 stderr += tmp_stderr2.read(buffsize)
115 stderr += tmp_stderr.read( buffsize ) 116 if not stderr or len(stderr) % buffsize != 0:
116 if not stderr or len( stderr ) % buffsize != 0: 117 break
117 break 118 except OverflowError:
118 except OverflowError: 119 pass
119 pass 120 raise Exception(stderr)
120 tmp_stderr.close() 121
121
122 # Error checking.
123 if returncode != 0:
124 raise Exception, stderr
125
126 # Copy outputs. 122 # Copy outputs.
127 shutil.copyfile( "cc_output.combined.gtf" , options.combined_transcripts ) 123 shutil.move("cc_output.combined.gtf", options.combined_transcripts)
128 124
129 # check that there are results in the output file 125 # check that there are results in the output file
130 cc_output_fname = "cc_output.stats" 126 cc_output_fname = "cc_output.stats"
131 if len( open( cc_output_fname, 'rb' ).read().strip() ) == 0: 127 if len(open(cc_output_fname).read().strip()) == 0:
132 raise Exception, 'The main output file is empty, there may be an error with your input file or settings.' 128 raise Exception('The main output file is empty, there may be an error with your input file or settings.')
133 except Exception, e: 129 except Exception as e:
134 stop_err( 'Error running cuffcompare. ' + str( e ) ) 130 stop_err('Error running cuffcompare: %s' % e)
135 131
136 if __name__=="__main__": __main__() 132
133 if __name__ == "__main__":
134 __main__()