|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 # Supports Cuffcompare versions v1.3.0 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( '-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' )
|
|
|
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!)' )
|
|
|
17
|
|
|
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.')
|
|
|
19
|
|
|
20 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')
|
|
|
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')
|
|
|
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')
|
|
|
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)')
|
|
|
25
|
|
|
26 # Wrapper / Galaxy options.
|
|
|
27 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' )
|
|
|
29
|
|
|
30 # Outputs.
|
|
|
31 parser.add_option( '', '--combined-transcripts', dest='combined_transcripts' )
|
|
|
32
|
|
|
33 (options, args) = parser.parse_args()
|
|
|
34
|
|
|
35 # output version # of tool
|
|
|
36 try:
|
|
|
37 tmp = tempfile.NamedTemporaryFile().name
|
|
|
38 tmp_stdout = open( tmp, 'wb' )
|
|
|
39 proc = subprocess.Popen( args='cuffcompare 2>&1', shell=True, stdout=tmp_stdout )
|
|
|
40 tmp_stdout.close()
|
|
|
41 returncode = proc.wait()
|
|
|
42 stdout = None
|
|
|
43 for line in open( tmp_stdout.name, 'rb' ):
|
|
|
44 if line.lower().find( 'cuffcompare v' ) >= 0:
|
|
|
45 stdout = line.strip()
|
|
|
46 break
|
|
|
47 if stdout:
|
|
|
48 sys.stdout.write( '%s\n' % stdout )
|
|
|
49 else:
|
|
|
50 raise Exception
|
|
|
51 except:
|
|
|
52 sys.stdout.write( 'Could not determine Cuffcompare version\n' )
|
|
|
53
|
|
|
54 # Set/link to sequence file.
|
|
|
55 if options.use_seq_data:
|
|
|
56 if options.ref_file:
|
|
|
57 # Sequence data from history.
|
|
|
58 # Create symbolic link to ref_file so that index will be created in working directory.
|
|
|
59 seq_path = "ref.fa"
|
|
|
60 os.symlink( options.ref_file, seq_path )
|
|
|
61 else:
|
|
|
62 if not os.path.exists( options.index ):
|
|
|
63 stop_err( 'Reference genome %s not present, request it by reporting this error.' % options.index )
|
|
|
64 seq_path = options.index
|
|
|
65
|
|
|
66 # Build command.
|
|
|
67
|
|
|
68 # Base.
|
|
|
69 cmd = "cuffcompare -o cc_output "
|
|
|
70
|
|
|
71 # Add options.
|
|
|
72 if options.ref_annotation:
|
|
|
73 cmd += " -r %s " % options.ref_annotation
|
|
|
74 if options.ignore_nonoverlap_reference:
|
|
|
75 cmd += " -R "
|
|
|
76 if options.ignore_nonoverlap_transfrag:
|
|
|
77 cmd += " -Q "
|
|
|
78 if options.use_seq_data:
|
|
|
79 cmd += " -s %s " % seq_path
|
|
|
80 if options.discard_single_exon_all:
|
|
|
81 cmd += " -M "
|
|
|
82 if options.discard_single_exon_ref:
|
|
|
83 cmd += " -N "
|
|
|
84 if options.max_dist_exon:
|
|
|
85 cmd += " -e %i " % int( options.max_dist_exon )
|
|
|
86 if options.max_dist_group:
|
|
|
87 cmd += " -d %i " % int( options.max_dist_group )
|
|
|
88 if options.discard_redundant_intron_transfrags:
|
|
|
89 cmd += " -F "
|
|
|
90 # Add input files.
|
|
|
91
|
|
|
92 # Need to symlink inputs so that output files are written to temp directory.
|
|
|
93 for i, arg in enumerate( args ):
|
|
|
94 input_file_name = "./input%i" % ( i+1 )
|
|
|
95 os.symlink( arg, input_file_name )
|
|
|
96 cmd += "%s " % input_file_name
|
|
|
97
|
|
|
98 # Debugging.
|
|
|
99 print cmd
|
|
|
100
|
|
|
101 # Run command.
|
|
|
102 try:
|
|
|
103 tmp_name = tempfile.NamedTemporaryFile( dir="." ).name
|
|
|
104 tmp_stderr = open( tmp_name, 'wb' )
|
|
|
105 proc = subprocess.Popen( args=cmd, shell=True, stderr=tmp_stderr.fileno() )
|
|
|
106 returncode = proc.wait()
|
|
|
107 tmp_stderr.close()
|
|
|
108
|
|
|
109 # Get stderr, allowing for case where it's very large.
|
|
|
110 tmp_stderr = open( tmp_name, 'rb' )
|
|
|
111 stderr = ''
|
|
|
112 buffsize = 1048576
|
|
|
113 try:
|
|
|
114 while True:
|
|
|
115 stderr += tmp_stderr.read( buffsize )
|
|
|
116 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
117 break
|
|
|
118 except OverflowError:
|
|
|
119 pass
|
|
|
120 tmp_stderr.close()
|
|
|
121
|
|
|
122 # Error checking.
|
|
|
123 if returncode != 0:
|
|
|
124 raise Exception, stderr
|
|
|
125
|
|
|
126 # Copy outputs.
|
|
|
127 shutil.copyfile( "cc_output.combined.gtf" , options.combined_transcripts )
|
|
|
128
|
|
|
129 # check that there are results in the output file
|
|
|
130 cc_output_fname = "cc_output.stats"
|
|
|
131 if len( open( cc_output_fname, 'rb' ).read().strip() ) == 0:
|
|
|
132 raise Exception, 'The main output file is empty, there may be an error with your input file or settings.'
|
|
|
133 except Exception, e:
|
|
|
134 stop_err( 'Error running cuffcompare. ' + str( e ) )
|
|
|
135
|
|
|
136 if __name__=="__main__": __main__()
|
