root/galaxy-central/tools/ngs_rna/cuffcompare_wrapper.py

リビジョン 2, 4.9 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
1#!/usr/bin/env python
2
3import optparse, os, shutil, subprocess, sys, tempfile
4
5def stop_err( msg ):
6    sys.stderr.write( "%s\n" % msg )
7    sys.exit()
8
9def __main__():
10    #Parse Command Line
11    parser = optparse.OptionParser()
12    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.' )
13    parser.add_option( '-R', action="store_true", dest='ignore_nonoverlap', 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' )
14   
15    # Wrapper / Galaxy options.
16    parser.add_option( '-1', dest='input1')
17    parser.add_option( '-2', dest='input2')
18    parser.add_option( '-A', '--transcripts-accuracy-output', dest='transcripts_accuracy_output_file', help='' )
19    parser.add_option( '-B', '--transcripts-combined-output', dest='transcripts_combined_output_file', help='' )
20    parser.add_option( '-C', '--transcripts-tracking-output', dest='transcripts_tracking_output_file', help='' )
21    parser.add_option( '', '--input1-tmap-output', dest='input1_tmap_output_file', help='' )
22    parser.add_option( '', '--input1-refmap-output', dest='input1_refmap_output_file', help='' )
23    parser.add_option( '', '--input2-tmap-output', dest='input2_tmap_output_file', help='' )
24    parser.add_option( '', '--input2-refmap-output', dest='input2_refmap_output_file', help='' )
25   
26   
27    (options, args) = parser.parse_args()
28   
29    # Make temp directory for output.
30    tmp_output_dir = tempfile.mkdtemp()
31   
32    # Build command.
33   
34    # Base.
35    cmd = "cuffcompare -o cc_output"
36   
37    # Add options.
38    if options.ref_annotation:
39        cmd += " -r %s" % options.ref_annotation
40    if options.ignore_nonoverlap:
41        cmd += " -R "
42       
43    # Add input files.
44       
45    # Need to symlink inputs so that output files are written to temp directory.
46    input1_file_name = tmp_output_dir + "/input1"
47    os.symlink( options.input1,  input1_file_name )
48    cmd += " %s" % input1_file_name
49    two_inputs = ( options.input2 != None)
50    if two_inputs:
51        input2_file_name = tmp_output_dir + "/input2"
52        os.symlink( options.input2, input2_file_name )
53        cmd += " %s" % input2_file_name
54   
55    # Run command.
56    try:       
57        tmp_name = tempfile.NamedTemporaryFile( dir=tmp_output_dir ).name
58        tmp_stderr = open( tmp_name, 'wb' )
59        proc = subprocess.Popen( args=cmd, shell=True, cwd=tmp_output_dir, stderr=tmp_stderr.fileno() )
60        returncode = proc.wait()
61        tmp_stderr.close()
62       
63        # Get stderr, allowing for case where it's very large.
64        tmp_stderr = open( tmp_name, 'rb' )
65        stderr = ''
66        buffsize = 1048576
67        try:
68            while True:
69                stderr += tmp_stderr.read( buffsize )
70                if not stderr or len( stderr ) % buffsize != 0:
71                    break
72        except OverflowError:
73            pass
74        tmp_stderr.close()
75       
76        # Error checking.
77        if returncode != 0:
78            raise Exception, stderr
79           
80        # check that there are results in the output file
81        if len( open( tmp_output_dir + "/cc_output", 'rb' ).read().strip() ) == 0:
82            raise Exception, 'The main output file is empty, there may be an error with your input file or settings.'
83    except Exception, e:
84        stop_err( 'Error running cuffcompare. ' + str( e ) )
85       
86    # Copy output files from tmp directory to specified files.
87    try:
88        try:
89            shutil.copyfile( tmp_output_dir + "/cc_output", options.transcripts_accuracy_output_file )
90            shutil.copyfile( tmp_output_dir + "/input1.tmap", options.input1_tmap_output_file )
91            shutil.copyfile( tmp_output_dir + "/input1.refmap", options.input1_refmap_output_file )
92            if two_inputs:
93                shutil.copyfile( tmp_output_dir + "/cc_output.combined.gtf", options.transcripts_combined_output_file )
94                shutil.copyfile( tmp_output_dir + "/cc_output.tracking", options.transcripts_tracking_output_file )
95                shutil.copyfile( tmp_output_dir + "/input2.tmap", options.input2_tmap_output_file )
96                shutil.copyfile( tmp_output_dir + "/input2.refmap", options.input2_refmap_output_file )
97        except Exception, e:
98            stop_err( 'Error in cuffcompare:\n' + str( e ) )
99    finally:
100        # Clean up temp dirs
101        if not os.path.exists( tmp_output_dir ):
102            shutil.rmtree( tmp_output_dir )
103
104if __name__=="__main__": __main__()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。