PageRenderTime 29ms CodeModel.GetById 13ms app.highlight 11ms RepoModel.GetById 1ms app.codeStats 0ms

/tools/ngs_rna/cuffcompare_wrapper.py

https://bitbucket.org/cistrome/cistrome-harvard/
Python | 131 lines | 122 code | 5 blank | 4 comment | 5 complexity | 3e0873f4bc3e9c9cc9deb4a2f297c6fb MD5 | raw file
  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
  9# Copied from sam_to_bam.py:
 10def check_seq_file( dbkey, cached_seqs_pointer_file ):
 11    seq_path = ''
 12    for line in open( cached_seqs_pointer_file ):
 13        line = line.rstrip( '\r\n' )
 14        if line and not line.startswith( '#' ) and line.startswith( 'index' ):
 15            fields = line.split( '\t' )
 16            if len( fields ) < 3:
 17                continue
 18            if fields[1] == dbkey:
 19                seq_path = fields[2].strip()
 20                break
 21    return seq_path
 22
 23def __main__():
 24    #Parse Command Line
 25    parser = optparse.OptionParser()
 26    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.' )
 27    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' )
 28    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.')
 29    
 30    # Wrapper / Galaxy options.
 31    parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' )
 32    parser.add_option( '', '--index_dir', dest='index_dir', help='GALAXY_DATA_INDEX_DIR' )
 33    parser.add_option( '', '--ref_file', dest='ref_file', help='The reference dataset from the history' )
 34    
 35    # Outputs.
 36    parser.add_option( '', '--combined-transcripts', dest='combined_transcripts' )
 37    
 38    (options, args) = parser.parse_args()
 39    
 40    # output version # of tool
 41    try:
 42        tmp = tempfile.NamedTemporaryFile().name
 43        tmp_stdout = open( tmp, 'wb' )
 44        proc = subprocess.Popen( args='cuffcompare 2>&1', shell=True, stdout=tmp_stdout )
 45        tmp_stdout.close()
 46        returncode = proc.wait()
 47        stdout = None
 48        for line in open( tmp_stdout.name, 'rb' ):
 49            if line.lower().find( 'cuffcompare v' ) >= 0:
 50                stdout = line.strip()
 51                break
 52        if stdout:
 53            sys.stdout.write( '%s\n' % stdout )
 54        else:
 55            raise Exception
 56    except:
 57        sys.stdout.write( 'Could not determine Cuffcompare version\n' )
 58        
 59    # Set/link to sequence file.
 60    if options.use_seq_data:
 61        cached_seqs_pointer_file = os.path.join( options.index_dir, 'sam_fa_indices.loc' )
 62        if not os.path.exists( cached_seqs_pointer_file ):
 63            stop_err( 'The required file (%s) does not exist.' % cached_seqs_pointer_file )
 64        # If found for the dbkey, seq_path will look something like /galaxy/data/equCab2/sam_index/equCab2.fa,
 65        # and the equCab2.fa file will contain fasta sequences.
 66        seq_path = check_seq_file( options.dbkey, cached_seqs_pointer_file )
 67        if options.ref_file != 'None':
 68            # Create symbolic link to ref_file so that index will be created in working directory.
 69            seq_path = "ref.fa"
 70            os.symlink( options.ref_file, seq_path  )
 71    
 72    # Build command.
 73    
 74    # Base.
 75    cmd = "cuffcompare -o cc_output "
 76    
 77    # Add options.
 78    if options.ref_annotation:
 79        cmd += " -r %s " % options.ref_annotation
 80    if options.ignore_nonoverlap:
 81        cmd += " -R "
 82    if options.use_seq_data:
 83        cmd += " -s %s " % seq_path
 84        
 85    # Add input files.
 86        
 87    # Need to symlink inputs so that output files are written to temp directory.
 88    for i, arg in enumerate( args ):
 89        input_file_name = "./input%i" % ( i+1 )
 90        os.symlink( arg, input_file_name )
 91        cmd += "%s " % input_file_name
 92
 93    # Debugging.
 94    print cmd
 95    
 96    # Run command.
 97    try:        
 98        tmp_name = tempfile.NamedTemporaryFile( dir="." ).name
 99        tmp_stderr = open( tmp_name, 'wb' )
100        proc = subprocess.Popen( args=cmd, shell=True, stderr=tmp_stderr.fileno() )
101        returncode = proc.wait()
102        tmp_stderr.close()
103        
104        # Get stderr, allowing for case where it's very large.
105        tmp_stderr = open( tmp_name, 'rb' )
106        stderr = ''
107        buffsize = 1048576
108        try:
109            while True:
110                stderr += tmp_stderr.read( buffsize )
111                if not stderr or len( stderr ) % buffsize != 0:
112                    break
113        except OverflowError:
114            pass
115        tmp_stderr.close()
116        
117        # Error checking.
118        if returncode != 0:
119            raise Exception, stderr
120            
121        # Copy outputs.
122        shutil.copyfile( "cc_output.combined.gtf" , options.combined_transcripts )    
123            
124        # check that there are results in the output file
125        cc_output_fname = "cc_output.stats"
126        if len( open( cc_output_fname, 'rb' ).read().strip() ) == 0:
127            raise Exception, 'The main output file is empty, there may be an error with your input file or settings.'
128    except Exception, e:
129        stop_err( 'Error running cuffcompare. ' + str( e ) )
130        
131if __name__=="__main__": __main__()